AliPhysics  d497547 (d497547)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskDmesonJets.cxx
Go to the documentation of this file.
1 /*************************************************************************
2 * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15 
16 // Root
17 #include <TClonesArray.h>
18 #include <TDatabasePDG.h>
19 #include <TParticlePDG.h>
20 #include <TVector3.h>
21 #include <THnSparse.h>
22 #include <TParticle.h>
23 #include <TMath.h>
24 #include <THashList.h>
25 #include <TFile.h>
26 #include <TRandom3.h>
27 
28 // Aliroot general
29 #include "AliLog.h"
30 #include "AliEMCALGeometry.h"
31 #include "AliAnalysisManager.h"
32 #include "AliVEventHandler.h"
33 
34 // Aliroot HF
35 #include "AliAODRecoCascadeHF.h"
37 #include "AliRDHFCutsD0toKpi.h"
40 #include "AliHFTrackContainer.h"
41 #include "AliAnalysisVertexingHF.h"
42 
43 // Aliroot EMCal jet framework
44 #include "AliEmcalJetTask.h"
45 #include "AliEmcalJet.h"
46 #include "AliJetContainer.h"
47 #include "AliParticleContainer.h"
48 #include "AliEmcalParticle.h"
49 #include "AliFJWrapper.h"
50 #include "AliRhoParameter.h"
51 
53 
54 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfo
55 
59 
67 {
68  dphi = TVector2::Phi_mpi_pi(fMomentum.Phi() - jet.Phi());;
69  deta = fMomentum.Eta() - jet.Eta();
70  return TMath::Sqrt(dphi*dphi + deta*deta);
71 }
72 
78 {
79  Double_t deta = 0;
80  Double_t dphi = 0;
81  return GetDistance(jet, deta, dphi);
82 }
83 
84 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonJetInfo
85 
89 
92  fDmesonParticle(0),
93  fD(),
94  fSoftPionPt(0),
95  fInvMass2Prong(0),
96  fJets(),
97  fMCLabel(-1),
98  fReconstructed(kFALSE),
99  fParton(0),
100  fPartonType(0),
101  fAncestor(0),
102  fSelectionType(0)
103 {
104 }
105 
110  fDmesonParticle(source.fDmesonParticle),
111  fD(source.fD),
112  fSoftPionPt(source.fSoftPionPt),
113  fInvMass2Prong(source.fInvMass2Prong),
114  fJets(source.fJets),
115  fMCLabel(source.fMCLabel),
116  fReconstructed(source.fReconstructed),
117  fParton(source.fParton),
118  fPartonType(source.fPartonType),
119  fAncestor(source.fAncestor),
120  fSelectionType(source.fSelectionType)
121 {
122 }
123 
128 {
129  new (this) AliDmesonJetInfo(source);
130  return *this;
131 }
132 
135 {
136  fD.SetPtEtaPhiE(0,0,0,0);
137  fSoftPionPt = 0;
138  fInvMass2Prong = 0;
139  fDmesonParticle = 0;
140  fMCLabel = -1;
141  fReconstructed = kFALSE;
142  fParton = 0;
143  fPartonType = 0;
144  fAncestor = 0;
145  for (auto &jet : fJets) {
146  jet.second.fMomentum.SetPtEtaPhiE(0,0,0,0);
147  jet.second.fNConstituents = 0;
148  jet.second.fNEF = 0;
149  jet.second.fMaxChargedPt = 0;
150  jet.second.fMaxNeutralPt = 0;
151  }
152 }
153 
156 {
157  Printf("Printing D Meson Jet object.");
158  Printf("D Meson: pT = %.3f, eta = %.3f, phi = %.3f, inv. mass = %.3f", fD.Pt(), fD.Eta(), fD.Phi_0_2pi(), fD.M());
159  Printf("Soft pion pT: %.3f. 2-Prong Invariant mass = %.3f", fSoftPionPt, fInvMass2Prong);
160  for (auto &jet : fJets) {
161  Printf("Jet %s: pT = %.3f, eta = %.3f, phi = %.3f", jet.first.c_str(), jet.second.Pt(), jet.second.Eta(), jet.second.Phi_0_2pi());
162  Printf("Jet N Consituents = %d", jet.second.fNConstituents);
163  }
164 }
165 
170 {
171  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
172  if (it == fJets.end()) return 0;
173 
174  Double_t z = 0;
175 
176  if ((*it).second.Pt() > 0) {
177  TVector3 dvect = fD.Vect();
178  TVector3 jvect = (*it).second.fMomentum.Vect();
179 
180  Double_t jetMom = jvect * jvect;
181 
182  if (jetMom < 1e-6) {
183  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetZ", "Zero jet momentum!");
184  z = 0.999;
185  }
186  else {
187  z = (dvect * jvect) / jetMom;
188  }
189 
190  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin 0.9-1 rather than 1-1.1
191  }
192 
193  return z;
194 }
195 
200 {
201  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
202  if (it == fJets.end()) return 0;
203 
204  Double_t z = 0;
205 
206  if ((*it).second.Pt() > 0) {
207  TVector3 dvect = fD.Vect();
208  TVector3 jvect = (*it).second.fMomentum.Vect();
209  // If the corr pt is < 0, assign 0.
210  Double_t corrpt = (*it).second.fCorrPt > 0 ? (*it).second.fCorrPt : 0.;
211  jvect.SetPerp(corrpt);
212 
213  Double_t jetMom = jvect * jvect;
214 
215  if (jetMom < 1e-6) {
216  z = 1.0;
217  }
218  else {
219  z = (dvect * jvect) / jetMom;
220  }
221  }
222 
223  return z;
224 }
225 
233 {
234  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
235  if (it == fJets.end()) return 0;
236 
237  return GetDistance((*it).second, deta, dphi);
238 }
239 
245 {
246  Double_t deta = 0;
247  Double_t dphi = 0;
248  return GetDistance(n, deta, dphi);
249 }
250 
258 {
259  dphi = TVector2::Phi_mpi_pi(fD.Phi() - jet.Phi());;
260  deta = fD.Eta() - jet.Eta();
261  return TMath::Sqrt(dphi*dphi + deta*deta);
262 }
263 
269 {
270  Double_t deta = 0;
271  Double_t dphi = 0;
272  return GetDistance(jet, deta, dphi);
273 }
274 
280 {
281  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
282  if (it == fJets.end()) {
283  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
284  n.c_str());
285  return 0;
286  }
287  return &((*it).second);
288 }
289 
295 {
296  std::map<std::string, AliJetInfo>::iterator it = fJets.find(n);
297  if (it == fJets.end()) {
298  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
299  n.c_str());
300  return 0;
301  }
302  return &((*it).second);
303 }
304 
305 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoSummary
306 
310 
316  fPt(0),
317  fEta(0),
318  fPhi(0),
319  fR(0),
320  fZ(0),
321  fN(0)
322 {
323  Set(source, n);
324 }
325 
328 {
329  fPt = 0;
330  fEta = 0;
331  fPhi = 0;
332  fR = 0;
333  fZ = 0;
334  fN = 0;
335 }
336 
342 {
343  std::map<std::string, AliJetInfo>::const_iterator it = source.fJets.find(n);
344  if (it == source.fJets.end()) return;
345 
346  Set((*it).second);
347 
348  fR = source.GetDistance(n);
349  fZ = source.GetZ(n);
350 }
351 
356 {
357  fPt = source.Pt();
358  fEta = source.Eta();
359  fPhi = source.Phi_0_2pi();
360  fN = source.GetNConstituents();
361  fR = 0;
362  fZ = 0;
363 }
364 
365 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoPbPbSummary
366 
370 
377  fCorrPt(0),
378  fCorrZ(0),
379  fArea(0)
380 {
381  Set(source, n);
382 }
383 
386 {
388  fCorrPt = 0;
389  fCorrZ = 0;
390  fArea = 0;
391 }
392 
397 {
398  AliJetInfoSummary::Set(source);
399  fArea = source.fArea;
400  fCorrPt = source.fCorrPt;
401 }
402 
408 {
409  AliJetInfoSummary::Set(source, n);
410  fCorrZ = source.GetCorrZ(n);
411 }
412 
413 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonInfoSummary
414 
418 
423  fPt(0),
424  fEta(0),
425  fPhi(0)
426 {
427  Set(source);
428 }
429 
434 {
435  fPt = source.fD.Pt();
436  fEta = source.fD.Eta();
437  fPhi = source.fD.Phi_0_2pi();
438 }
439 
442 {
443  fPt = 0;
444  fEta = 0;
445  fPhi = 0;
446 }
447 
448 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary
449 
453 
458  AliDmesonInfoSummary(source),
459  fPartonType(0),
460  fPartonPt(0),
461  fAncestorPDG(0)
462 {
463  Set(source);
464 }
465 
470 {
472 
473  fPartonType = source.fPartonType;
474 
475  if (source.fParton) {
476  fPartonPt = source.fParton->Pt();
477  }
478  else {
479  fPartonPt = 0.;
480  }
481 
482  if (source.fAncestor) {
483  fAncestorPDG = (UShort_t)((UInt_t)(TMath::Abs(source.fAncestor->GetPdgCode())));
484  }
485 }
486 
489 {
491  fPartonType = 0,
492  fPartonPt = 0.;
493  fAncestorPDG = 0;
494 }
495 
496 // Definitions of class AliAnalysisTaskDmesonJets::AliD0InfoSummary
497 
501 
506  AliDmesonInfoSummary(source),
507  fInvMass(source.fD.M()),
508  fSelectionType(0)
509 {
510 }
511 
516 {
517  fInvMass = source.fD.M();
518  fSelectionType = source.fSelectionType;
520 }
521 
524 {
526  fSelectionType = 0;
527  fInvMass = 0;
528 }
529 
530 // Definitions of class AliAnalysisTaskDmesonJets::AliDStarInfoSummary
531 
535 
540  AliDmesonInfoSummary(source),
541  f2ProngInvMass(source.fInvMass2Prong),
542  fDeltaInvMass(source.fD.M() - source.fInvMass2Prong)
543 {
544 }
545 
550 {
551  f2ProngInvMass = source.fInvMass2Prong;
552  fDeltaInvMass = source.fD.M() - source.fInvMass2Prong;
554 }
555 
558 {
560 
561  f2ProngInvMass = 0;
562  fDeltaInvMass = 0;
563 }
564 
565 // Definitions of class AliAnalysisTaskDmesonJets::AliJetDefinition
566 
570 
573  TObject(),
574  fJetType(AliJetContainer::kChargedJet),
575  fRadius(0),
576  fJetAlgo(AliJetContainer::antikt_algorithm),
577  fRecoScheme(AliJetContainer::pt_scheme),
578  fMinJetPt(0.),
579  fMaxJetPt(500.),
580  fMinJetPhi(0.),
581  fMaxJetPhi(0.),
582  fMinJetEta(-1.),
583  fMaxJetEta(1.),
584  fMinChargedPt(0.),
585  fMaxChargedPt(0.),
586  fMinNeutralPt(0.),
587  fMaxNeutralPt(0.),
588  fRhoName(),
589  fRho(0)
590 {
591 }
592 
600  TObject(),
601  fJetType(type),
602  fRadius(r),
603  fJetAlgo(algo),
604  fRecoScheme(reco),
605  fMinJetPt(0.),
606  fMaxJetPt(500.),
607  fMinJetPhi(0.),
608  fMaxJetPhi(0.),
609  fMinJetEta(-1.),
610  fMaxJetEta(1.),
611  fMinChargedPt(0.),
612  fMaxChargedPt(0.),
613  fMinNeutralPt(0.),
614  fMaxNeutralPt(0.),
615  fRhoName(),
616  fRho(0)
617 {
618 }
619 
627  TObject(),
628  fJetType(type),
629  fRadius(r),
630  fJetAlgo(algo),
631  fRecoScheme(reco),
632  fMinJetPt(0.),
633  fMaxJetPt(500.),
634  fMinJetPhi(0.),
635  fMaxJetPhi(0.),
636  fMinJetEta(-1.),
637  fMaxJetEta(1.),
638  fMinChargedPt(0.),
639  fMaxChargedPt(0.),
640  fMinNeutralPt(0.),
641  fMaxNeutralPt(0.),
642  fRhoName(rhoName),
643  fRho(0)
644 {
645 }
646 
651  TObject(),
652  fJetType(source.fJetType),
653  fRadius(source.fRadius),
654  fJetAlgo(source.fJetAlgo),
655  fRecoScheme(source.fRecoScheme),
656  fMinJetPt(source.fMinJetPt),
657  fMaxJetPt(source.fMaxJetPt),
658  fMinJetPhi(source.fMinJetPhi),
659  fMaxJetPhi(source.fMaxJetPhi),
660  fMinJetEta(source.fMinJetEta),
661  fMaxJetEta(source.fMaxJetEta),
662  fMinChargedPt(source.fMinChargedPt),
663  fMaxChargedPt(source.fMaxChargedPt),
664  fMinNeutralPt(source.fMinNeutralPt),
665  fMaxNeutralPt(source.fMaxNeutralPt),
666  fRhoName(source.fRhoName),
667  fRho(0)
668 {
669 }
670 
675 {
676  new (this) AliHFJetDefinition(source);
677  return *this;
678 }
679 
682 {
683  static TString name;
684 
685  name = AliJetContainer::GenerateJetName(fJetType, fJetAlgo, fRecoScheme, fRadius, 0, 0, "Jet");
686 
687  return name.Data();
688 }
689 
695 {
696  if (fMinJetEta < fMaxJetEta && (jet.Eta() < fMinJetEta || jet.Eta() > fMaxJetEta)) return kFALSE;
697  if (fMinJetPhi < fMaxJetPhi && (jet.Phi() < fMinJetPhi || jet.Phi() > fMaxJetPhi)) return kFALSE;
698  if (jet.Pt() > fMaxJetPt || jet.Pt() < fMinJetPt) return kFALSE;
699  if (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt) return kFALSE;
700  if (jet.fMaxNeutralPt < fMinNeutralPt || jet.fMaxNeutralPt > fMaxNeutralPt) return kFALSE;
701 
702  return kTRUE;
703 }
704 
710 {
711  const AliJetInfo* jet = dMesonJet.GetJet(n);
712  if (!jet) return kFALSE;
713  return IsJetInAcceptance((*jet));
714 }
715 
722 {
723  if (lhs.fJetType > rhs.fJetType) return false;
724  else if (lhs.fJetType < rhs.fJetType) return true;
725  else {
726  if (lhs.fRadius > rhs.fRadius) return false;
727  else if (lhs.fRadius < rhs.fRadius) return true;
728  else {
729  if (lhs.fJetAlgo > rhs.fJetAlgo) return false;
730  else if (lhs.fJetAlgo < rhs.fJetAlgo) return true;
731  else {
732  if (lhs.fRecoScheme < rhs.fRecoScheme) return true;
733  else return false;
734  }
735  }
736  }
737 }
738 
745 {
746  if (lhs.fJetType != rhs.fJetType) return false;
747  if (lhs.fRadius != rhs.fRadius) return false;
748  if (lhs.fJetAlgo != rhs.fJetAlgo) return false;
749  if (lhs.fRecoScheme != rhs.fRecoScheme) return false;
750  return true;
751 }
752 
753 // Definitions of class AliAnalysisTaskDmesonJets::AnalysisEngine
754 
758 
761  TObject(),
762  fPartons(),
763  fCandidateType(kD0toKpi),
764  fCandidateName(),
765  fCandidatePDG(0),
766  fNDaughters(0),
767  fPDGdaughters(),
768  fBranchName(),
769  fMCMode(kNoMC),
770  fNMassBins(0),
771  fMinMass(0),
772  fMaxMass(0),
773  fRDHFCuts(0),
774  fRejectedOrigin(0),
775  fAcceptedDecay(0),
776  fInhibit(kFALSE),
777  fJetDefinitions(),
778  fPtBinWidth(0.5),
779  fMaxPt(100),
780  fRandomGen(0),
781  fTrackEfficiency(0),
782  fRejectISR(kFALSE),
783  fDataSlotNumber(-1),
784  fTree(0),
785  fCurrentDmesonJetInfo(0),
786  fCurrentJetInfo(0),
787  fCandidateArray(0),
788  fMCContainer(),
789  fTrackContainers(),
790  fClusterContainers(),
791  fAodEvent(0),
792  fFastJetWrapper(0),
793  fHistManager(0),
794  fCent(-1)
795 {
796 }
797 
806  TObject(),
807  fPartons(),
808  fCandidateType(type),
809  fCandidateName(),
810  fCandidatePDG(0),
811  fNDaughters(0),
812  fPDGdaughters(),
813  fBranchName(),
814  fMCMode(MCmode),
815  fNMassBins(nMassBins),
816  fMinMass(0),
817  fMaxMass(0),
818  fRDHFCuts(cuts),
819  fRejectedOrigin(0),
820  fAcceptedDecay(kAnyDecay),
821  fInhibit(kFALSE),
822  fJetDefinitions(),
823  fPtBinWidth(0.5),
824  fMaxPt(100),
825  fRandomGen(0),
826  fTrackEfficiency(0),
827  fDataSlotNumber(-1),
828  fTree(0),
829  fCurrentDmesonJetInfo(0),
830  fCurrentJetInfo(0),
831  fCandidateArray(0),
832  fMCContainer(),
833  fTrackContainers(),
834  fClusterContainers(),
835  fAodEvent(0),
836  fFastJetWrapper(0),
837  fHistManager(0),
838  fCent(-1)
839 {
840  SetCandidateProperties(range);
841 }
842 
847  TObject(source),
848  fPartons(source.fPartons),
849  fCandidateType(source.fCandidateType),
850  fCandidateName(source.fCandidateName),
851  fCandidatePDG(source.fCandidatePDG),
852  fNDaughters(source.fNDaughters),
853  fPDGdaughters(source.fPDGdaughters),
854  fBranchName(source.fBranchName),
855  fMCMode(source.fMCMode),
856  fNMassBins(source.fNMassBins),
857  fMinMass(source.fMinMass),
858  fMaxMass(source.fMaxMass),
859  fRDHFCuts(),
860  fRejectedOrigin(source.fRejectedOrigin),
861  fAcceptedDecay(source.fAcceptedDecay),
862  fInhibit(source.fInhibit),
863  fJetDefinitions(source.fJetDefinitions),
864  fPtBinWidth(source.fPtBinWidth),
865  fMaxPt(source.fMaxPt),
866  fRandomGen(source.fRandomGen),
868  fDataSlotNumber(-1),
869  fTree(0),
870  fCurrentDmesonJetInfo(0),
871  fCurrentJetInfo(0),
872  fCandidateArray(source.fCandidateArray),
873  fMCContainer(source.fMCContainer),
874  fTrackContainers(source.fTrackContainers),
875  fClusterContainers(source.fClusterContainers),
876  fAodEvent(source.fAodEvent),
878  fHistManager(source.fHistManager),
879  fCent(-1)
880 {
881  SetRDHFCuts(source.fRDHFCuts);
882 }
883 
884 // Destructor
886 {
887  delete fRDHFCuts;
888 }
889 
894 {
895  new (this) AnalysisEngine(source);
896  return *this;
897 }
898 
903 {
904  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
905  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
906  }
907 
908  return kFALSE;
909 }
910 
912 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const /*geom*/, Int_t /*runNumber*/)
913 {
914 }
915 
920 {
921  switch (fCandidateType) {
922  case kD0toKpi:
923  fCandidatePDG = 421;
924  fCandidateName = "D0";
925  fNDaughters = 2;
926  fPDGdaughters.Set(fNDaughters);
927  fPDGdaughters.Reset();
928  fPDGdaughters[0] = 211; // pi
929  fPDGdaughters[1] = 321; // K
930  fBranchName = "D0toKpi";
931  fAcceptedDecay = kDecayD0toKpi;
932  break;
933  case kD0toKpiLikeSign:
934  fCandidatePDG = 421;
935  fCandidateName = "2ProngLikeSign";
936  fNDaughters = 2;
937  fPDGdaughters.Set(fNDaughters);
938  fPDGdaughters.Reset();
939  fPDGdaughters[0] = 211; // pi
940  fPDGdaughters[1] = 321; // K
941  fBranchName = "LikeSign2Prong";
942  break;
943  case kDstartoKpipi:
944  fCandidatePDG = 413;
945  fCandidateName = "DStar";
946  fNDaughters = 3;
947  fPDGdaughters.Set(fNDaughters);
948  fPDGdaughters.Reset();
949  fPDGdaughters[0] = 211; // pi soft
950  fPDGdaughters[1] = 211; // pi fromD0
951  fPDGdaughters[2] = 321; // K from D0
952  fBranchName = "Dstar";
953  fAcceptedDecay = kDecayDStartoKpipi;
954  break;
955  default:
956  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
957  }
958 
959  CalculateMassLimits(range, fCandidatePDG, fNMassBins, fMinMass, fMaxMass);
960 }
961 
966 {
967  if (fRDHFCuts) delete fRDHFCuts;
968  fRDHFCuts = cuts;
969 }
970 
975 {
976  if (!cuts) return;
977  if (fRDHFCuts) delete fRDHFCuts;
978  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
979 }
980 
985 {
986  static TString name;
987 
988  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
989 
990  return name.Data();
991 }
992 
997 {
998  static TString name;
999 
1000  name = fCandidateName;
1001  switch (fMCMode) {
1002  case kBackgroundOnly:
1003  name += "_kBackgroundOnly";
1004  break;
1005  case kSignalOnly:
1006  name += "_kSignalOnly";
1007  break;
1008  case kMCTruth:
1009  name += "_MCTruth";
1010  break;
1011  case kWrongPID:
1012  name += "_WrongPID";
1013  break;
1014  default:
1015  break;
1016  }
1017 
1018  return name.Data();
1019 }
1020 
1028 {
1029  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
1030 
1031  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
1032  fJetDefinitions.push_back(def);
1033  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
1034  "Total number of jet definitions is now %lu.",
1035  def.GetName(), GetName(), fJetDefinitions.size());
1036  // For detector level set maximum track pt to 100 GeV/c
1037  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
1038  }
1039  else {
1040  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
1041  }
1042 
1043  return &(*it);
1044 }
1045 
1057 {
1058  AliHFJetDefinition def(type, r, algo, reco);
1059 
1060  return AddJetDefinition(def);
1061 }
1062 
1068 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
1069 {
1070  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
1071  while (it != fJetDefinitions.end() && (*it) != def) it++;
1072  return it;
1073 }
1074 
1081 {
1082  if (lhs.fCandidateType > rhs.fCandidateType) return false;
1083  else if (lhs.fCandidateType < rhs.fCandidateType) return true;
1084  else {
1085  if (lhs.fMCMode < rhs.fMCMode) return true;
1086  else return false;
1087  }
1088 }
1089 
1096 {
1097  if (lhs.fCandidateType != rhs.fCandidateType) return false;
1098  if (lhs.fMCMode != rhs.fMCMode) return false;
1099  return true;
1100 }
1101 
1110 {
1111  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) { // D0 candidate
1112  return ExtractD0Attributes(Dcand, DmesonJet, i);
1113  }
1114  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
1115  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
1116  }
1117  else {
1118  return kFALSE;
1119  }
1120 }
1121 
1130 {
1131  AliDebug(10,"Checking if D0 meson is selected");
1132  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
1133  if (isSelected == 0) return kFALSE;
1134 
1135  Int_t MCtruthPdgCode = 0;
1136 
1137  Double_t invMassD = 0;
1138 
1139  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
1140  // Checks also the origin, and if it matches the rejected origin mask, return false
1141  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly || fMCMode == kWrongPID) {
1142  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
1143  DmesonJet.fMCLabel = mcLab;
1144 
1145  // Retrieve the generated particle (if exists) and its PDG code
1146  if (mcLab >= 0) {
1147  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1148 
1149  if (aodMcPart) {
1150  // Check origin and return false if it matches the rejected origin mask
1151  if (fRejectedOrigin) {
1152  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
1153  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1154  }
1155  MCtruthPdgCode = aodMcPart->PdgCode();
1156  }
1157  }
1158  }
1159 
1160  if (isSelected == 1) { // selected as a D0
1161  if (i != 0) return kFALSE; // only one mass hypothesis thanks to PID
1162 
1163  if (fMCMode == kNoMC ||
1164  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1165  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly) ||
1166  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kWrongPID)) {
1167  // both background and signal are requested OR (it is a true D0 AND signal is requested) OR (it is NOT a D0 and background is requested)
1168  AliDebug(10,"Selected as D0");
1169  invMassD = Dcand->InvMassD0();
1170  }
1171  else { // conditions above not passed, so return FALSE
1172  return kFALSE;
1173  }
1174  }
1175  else if (isSelected == 2) { // selected as a D0bar
1176  if (i != 1) return kFALSE; // only one mass hypothesis thanks to PID
1177 
1178  if (fMCMode == kNoMC ||
1179  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1180  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly) ||
1181  (MCtruthPdgCode == fCandidatePDG && fMCMode == kWrongPID)) {
1182  // both background and signal are requested OR (it is a true D0bar AND signal is requested) OR (it is NOT a D0bar and background is requested)
1183  AliDebug(10,"Selected as D0bar");
1184  invMassD = Dcand->InvMassD0bar();
1185  }
1186  else { // conditions above not passed, so return FALSE
1187  return kFALSE;
1188  }
1189  }
1190  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
1191  AliDebug(10,"Selected as either D0 or D0bar");
1192 
1193  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
1194  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1195  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kWrongPID))) {
1196  if (i != 0) return kFALSE;
1197  AliDebug(10, "MC truth is D0");
1198  invMassD = Dcand->InvMassD0();
1199  }
1200  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1201  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kWrongPID))) {
1202  if (i != 1) return kFALSE;
1203  AliDebug(10, "MC truth is D0bar");
1204  invMassD = Dcand->InvMassD0bar();
1205  }
1206  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
1207 
1208  // Only accept it if background-only OR background-and-signal was requested
1209  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
1210  // Select D0 or D0bar depending on the i-parameter
1211  if (i == 0) {
1212  AliDebug(10, "Returning invariant mass with D0 hypothesis");
1213  invMassD = Dcand->InvMassD0();
1214  }
1215  else if (i == 1) {
1216  AliDebug(10, "Returning invariant mass with D0bar hypothesis");
1217  invMassD = Dcand->InvMassD0bar();
1218  }
1219  else { // i > 1
1220  return kFALSE;
1221  }
1222  }
1223  else { // signal-only was requested but this is not a true D0
1224  return kFALSE;
1225  }
1226  }
1227  }
1228  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
1229  return kTRUE;
1230 }
1231 
1240 {
1241  AliDebug(10,"Checking if D* meson is selected");
1242  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoCascadeHF*>(DstarCand), AliRDHFCuts::kAll, fAodEvent);
1243  if (isSelected == 0) return kFALSE;
1244 
1245  if ((i == 1 && DstarCand->Charge()>0) || (i == 0 && DstarCand->Charge()<0) || i > 1) return kFALSE; // only one mass hypothesis for the D*
1246 
1247  Int_t MCtruthPdgCode = 0;
1248 
1249  Double_t invMassD = 0;
1250 
1251  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1252  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1253  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1254 
1255  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1256  AliDebug(10, Form("MC label is %d", mcLab));
1257  DmesonJet.fMCLabel = mcLab;
1258  if (mcLab >= 0) {
1259  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1260 
1261  if (aodMcPart) {
1262  if (fRejectedOrigin) {
1263  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
1264  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1265  }
1266 
1267  MCtruthPdgCode = aodMcPart->PdgCode();
1268  AliDebug(10, Form("MC truth pdg code is %d",MCtruthPdgCode));
1269  }
1270  }
1271  }
1272 
1273  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1274  if (fMCMode == kNoMC ||
1275  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1276  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1277  // both background and signal are requested OR (it is a true D*/D*bar AND signal is requested) OR (it is NOT a D*/D*bar and background is requested)
1278  invMassD = DstarCand->InvMassDstarKpipi();
1279  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1280  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1281  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1282  return kTRUE;
1283  }
1284  else { // conditions above not passed, so return FALSE
1285  return kFALSE;
1286  }
1287 }
1288 
1296 {
1297  if (!part) return kUnknownDecay;
1298  if (!mcArray) return kUnknownDecay;
1299 
1301 
1302  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1303 
1304  if (part->GetNDaughters() == 2) {
1305 
1306  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1307  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1308 
1309  if (!d1 || !d2) {
1310  return decay;
1311  }
1312 
1313  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1314  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1315 
1316  if (absPdgPart == 421) { // D0 -> K pi
1317  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1318  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1319  decay = kDecayD0toKpi;
1320  }
1321  }
1322 
1323  if (absPdgPart == 413) { // D* -> D0 pi
1324  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1325  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1326  if (D0decay == kDecayD0toKpi) {
1327  decay = kDecayDStartoKpipi;
1328  }
1329  }
1330  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1331  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1332  if (D0decay == kDecayD0toKpi) {
1333  decay = kDecayDStartoKpipi;
1334  }
1335  }
1336  }
1337  }
1338 
1339  return decay;
1340 }
1341 
1348 std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> AliAnalysisTaskDmesonJets::AnalysisEngine::IsPromptCharm(const AliAODMCParticle* part, TClonesArray* mcArray)
1349 {
1350  std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> result(kUnknownQuark, 0);
1351 
1352  if (!part) return result;
1353  if (!mcArray) return result;
1354 
1355  static std::set<UInt_t> partons = { 4, 5 };
1356 
1357  AliAODMCParticle* parton = FindParticleOrigin(part, mcArray, kFindLast, partons);
1358  if (parton) {
1359  result.second = parton;
1360  UInt_t absPdgParton = TMath::Abs(parton->GetPdgCode());
1361  if (absPdgParton == 4) result.first = kFromCharm;
1362  else if (absPdgParton == 5) result.first = kFromBottom;
1363  }
1364 
1365  return result;
1366 }
1367 
1375 
1376 AliAODMCParticle* AliAnalysisTaskDmesonJets::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode)
1377 {
1378  static std::set<UInt_t> pdgSet;
1379 
1380  return FindParticleOrigin(part, mcArray, mode, pdgSet);
1381 }
1382 
1391 
1392 AliAODMCParticle* AliAnalysisTaskDmesonJets::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode, const std::set<UInt_t>& pdgSet)
1393 {
1394  AliAODMCParticle* result = nullptr;
1395 
1396  Int_t mother = part->GetMother();
1397  while (mother >= 0) {
1398  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1399  if (mcGranma) {
1400  UInt_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
1401 
1402  // If the current particle is one of the particle types that is being searched assign it to the result pointer
1403  if (pdgSet.empty() || pdgSet.count(abspdgGranma) > 0) {
1404  result = mcGranma;
1405  // If the last particle in the fragmentation tree (first when going reverse) was requested then stop the loop
1406  if (mode == kFindLast) break;
1407  }
1408  mother = mcGranma->GetMother();
1409  }
1410  else {
1411  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::FindParticleOrigin", "Could not retrieve mother particle %d!", mother);
1412  break;
1413  }
1414  }
1415 
1416  return result;
1417 }
1418 
1421 {
1422  for (auto& jetDef : fJetDefinitions) {
1423  jetDef.fJets.clear();
1424  }
1425 
1426  if (fMCMode == kMCTruth) {
1427  RunParticleLevelAnalysis();
1428  }
1429  else {
1430  RunDetectorLevelAnalysis();
1431  }
1432 }
1433 
1436 {
1437  // Fill the vertex info of the candidates
1438  // Needed for reduced delta AOD, where the vertex info has been deleted
1439  // to reduce the delta AOD file size
1441 
1442  const Int_t nD = fCandidateArray->GetEntriesFast();
1443 
1444  AliDmesonJetInfo DmesonJet;
1445 
1446  std::map<AliHFJetDefinition*,Double_t> maxJetPt;
1447  for (auto& def : fJetDefinitions) maxJetPt[&def] = 0;
1448  Double_t maxDPt = 0;
1449 
1450  Int_t nAccCharm[3] = {0};
1451  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1452  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1453  if (!charmCand) continue;
1454  if(!(vHF.FillRecoCand(fAodEvent,charmCand))) continue;
1455 
1456  // region of interest + cuts
1457  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1458  Int_t nMassHypo = 0; // number of mass hypothesis accepted for this D meson
1459  if (charmCand->Pt() > maxDPt) maxDPt = charmCand->Pt();
1460  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1461  DmesonJet.Reset();
1462  DmesonJet.fDmesonParticle = charmCand;
1463  DmesonJet.fSelectionType = im + 1;
1464  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1465  for (auto& def : fJetDefinitions) {
1466  if (FindJet(charmCand, DmesonJet, def)) {
1467  Double_t jetPt = DmesonJet.fJets[def.GetName()].fMomentum.Pt();
1468  if (jetPt > maxJetPt[&def]) maxJetPt[&def] = jetPt;
1469  }
1470  else {
1471  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1472  def.GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1473  }
1474  }
1475  fDmesonJets[(icharm+1)*(1-(im*2))] = DmesonJet;
1476  nMassHypo++;
1477  nAccCharm[im]++;
1478  }
1479  }
1480  if (nMassHypo == 2) {
1481  nAccCharm[0]--;
1482  nAccCharm[1]--;
1483  nAccCharm[2] += 2;
1484  }
1485  if (nMassHypo == 2) { // both mass hypothesis accepted
1486  fDmesonJets[(icharm+1)].fSelectionType = 3;
1487  fDmesonJets[-(icharm+1)].fSelectionType = 3;
1488  }
1489  } // end of D cand loop
1490 
1491  TString hname;
1492 
1493  Int_t ntracks = 0;
1494 
1495  for (auto track_cont : fTrackContainers) {
1496  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
1497  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(nullptr);
1498  ntracks += track_cont->GetNAcceptEntries();
1499  }
1500 
1501  for (auto& def : fJetDefinitions) {
1502  if (!def.fRho) continue;
1503  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
1504  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
1505 
1506  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
1507  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
1508 
1509  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
1510  fHistManager->FillTH2(hname, fCent, def.fRho->GetVal());
1511 
1512  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
1513  fHistManager->FillTH2(hname, fCent, maxJetPt[&def]);
1514 
1515  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
1516  fHistManager->FillTH2(hname, fCent, maxDPt);
1517 
1518  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
1519  fHistManager->FillTH2(hname, ntracks, def.fRho->GetVal());
1520 
1521  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
1522  fHistManager->FillTH2(hname, ntracks, maxJetPt[&def]);
1523 
1524  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
1525  fHistManager->FillTH2(hname, ntracks, maxDPt);
1526  }
1527 
1528  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1529  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1530  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1531  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1532 
1533  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1534  fHistManager->FillTH2(hname, ntracks, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1535 
1536  hname = TString::Format("%s/fHistNDmesons", GetName());
1537  fHistManager->FillTH1(hname, nD);
1538 }
1539 
1550 {
1551  TString hname;
1552 
1553  Double_t rho = 0;
1554  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
1555 
1557  fFastJetWrapper->SetR(jetDef.fRadius);
1560 
1561  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1562 
1563  if (jetDef.fJetType != AliJetContainer::kNeutralJet) {
1564  for (auto track_cont : fTrackContainers) {
1565  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
1566  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(Dcand);
1567  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1568  AddInputVectors(track_cont, 100, static_cast<TH2*>(fHistManager->FindObject(hname)), fTrackEfficiency);
1569 
1570  if (hftrack_cont) {
1571  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1572  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1573  const TObjArray& daughters = hftrack_cont->GetDaughterList();
1574  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1575  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1576  if (!hftrack_cont->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1577  }
1578  }
1579  }
1580  }
1581 
1582  if (jetDef.fJetType != AliJetContainer::kChargedJet) {
1583  for (auto clus_cont : fClusterContainers) {
1584  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1585  AddInputVectors(clus_cont, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1586  }
1587  }
1588 
1589  // run jet finder
1590  fFastJetWrapper->Run();
1591 
1592  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1593 
1594  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1595  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1596 
1597  Bool_t isDmesonJet = kFALSE;
1598 
1599  Double_t maxChPt = 0;
1600  Double_t maxNePt = 0;
1601  Double_t totalNeutralPt = 0;
1602  Int_t nConst = 1;
1603 
1604  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1605  if (constituents[ic].user_index() == 0) {
1606  isDmesonJet = kTRUE;
1607  }
1608  else if (constituents[ic].user_index() >= 100) {
1609  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1610  nConst++;
1611  }
1612  else if (constituents[ic].user_index() <= -100) {
1613  totalNeutralPt += constituents[ic].pt();
1614  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1615  nConst++;
1616  }
1617  }
1618 
1619  if (isDmesonJet) {
1620  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1621  DmesonJet.fJets[jetDef.GetName()].fNConstituents = nConst;
1622  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1623  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1624  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1625  DmesonJet.fJets[jetDef.GetName()].fArea = jets_incl[ijet].area();
1626  DmesonJet.fJets[jetDef.GetName()].fCorrPt = DmesonJet.fJets[jetDef.GetName()].fMomentum.Pt() - jets_incl[ijet].area() * rho;
1627 
1628  return kTRUE;
1629  }
1630  }
1631 
1632  return kFALSE;
1633 }
1634 
1638 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist, Double_t eff)
1639 {
1640  auto itcont = cont->all_momentum();
1641  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1642  UInt_t rejectionReason = 0;
1643  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1644  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1645  continue;
1646  }
1647  if (fRandomGen && eff > 0 && eff < 1) {
1648  Double_t rnd = fRandomGen->Rndm();
1649  if (eff < rnd) {
1650  if (rejectHist) rejectHist->Fill(6, it->first.Pt());
1651  continue;
1652  }
1653  }
1654  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1655  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1656  }
1657 }
1658 
1661 {
1662  TString hname;
1663 
1664  fMCContainer->SetSpecialPDG(fCandidatePDG);
1665  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1666  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1668 
1669  if (!fMCContainer->IsSpecialPDGFound()) return;
1670 
1671  Int_t nAccCharm[3] = {0};
1672 
1673  std::map<AliHFJetDefinition*, Double_t> maxJetPt;
1674  Double_t maxDPt = 0;
1675 
1676  for (auto &jetDef : fJetDefinitions) {
1677  maxJetPt[&jetDef] = 0;
1678  Double_t rho = 0;
1679  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
1680  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", GetName(), jetDef.GetName());
1681  TH1* histNDmesonsVsNconstituents = static_cast<TH1*>(fHistManager->FindObject(hname));
1682 
1683  switch (jetDef.fJetType) {
1685  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
1686  break;
1688  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
1689  break;
1691  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
1692  break;
1693  }
1694 
1696  fFastJetWrapper->SetR(jetDef.fRadius);
1699 
1700  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1701  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1702 
1703  fFastJetWrapper->Run();
1704 
1705  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1706 
1707  for (auto jet : jets_incl) {
1708  Int_t nDmesonsInJet = 0;
1709 
1710  for (auto constituent : jet.constituents()) {
1711  Int_t iPart = constituent.user_index() - 100;
1712  if (constituent.perp() < 1e-6) continue; // reject ghost particles
1713  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
1714  if (!part) {
1715  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1716  continue;
1717  }
1718  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1719  nDmesonsInJet++;
1720  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1721  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1722  if (part->Pt() > maxDPt) maxDPt = part->Pt();
1723  std::pair<int, AliDmesonJetInfo> element;
1724  element.first = iPart;
1725  dMesonJetIt = fDmesonJets.insert(element).first;
1726  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1727  (*dMesonJetIt).second.fDmesonParticle = part;
1728  (*dMesonJetIt).second.fSelectionType = part->PdgCode() > 0 ? 1 : 2;
1729 
1730  UShort_t p = 0;
1731  UInt_t rs = 0;
1732 
1733  auto origin = IsPromptCharm(part, fMCContainer->GetArray());
1734  p = 0;
1735  rs = origin.first;
1736  while (rs >>= 1) { p++; }
1737  (*dMesonJetIt).second.fPartonType = p;
1738  (*dMesonJetIt).second.fParton = origin.second;
1739 
1740  (*dMesonJetIt).second.fAncestor = FindParticleOrigin(part, fMCContainer->GetArray(), kFindFirst);
1741 
1742  if (part->PdgCode() > 0) { // D0
1743  nAccCharm[0]++;
1744  }
1745  else { // D0bar
1746  nAccCharm[1]++;
1747  }
1748  }
1749 
1750  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
1751  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = jet.constituents().size();
1752  (*dMesonJetIt).second.fJets[jetDef.GetName()].fArea = jet.area();
1753  (*dMesonJetIt).second.fJets[jetDef.GetName()].fCorrPt = (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.Pt() - jet.area() * rho;
1754  if (jet.perp() > maxJetPt[&jetDef]) maxJetPt[&jetDef] = jet.perp();
1755  } // if constituent is a D meson
1756  } // for each constituent
1757  if (nDmesonsInJet > 0) histNDmesonsVsNconstituents->Fill(jet.constituents().size(), nDmesonsInJet);
1758  } // for each jet
1759  } // for each jet definition
1760 
1762 
1763  for (auto& def : fJetDefinitions) {
1764  if (!def.fRho) continue;
1765  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
1766  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
1767 
1768  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
1769  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
1770 
1771  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
1772  fHistManager->FillTH2(hname, fCent, def.fRho->GetVal());
1773 
1774  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
1775  fHistManager->FillTH2(hname, fCent, maxJetPt[&def]);
1776 
1777  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
1778  fHistManager->FillTH2(hname, fCent, maxDPt);
1779 
1780  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
1781  fHistManager->FillTH2(hname, npart, def.fRho->GetVal());
1782 
1783  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
1784  fHistManager->FillTH2(hname, npart, maxJetPt[&def]);
1785 
1786  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
1787  fHistManager->FillTH2(hname, npart, maxDPt);
1788  }
1789 
1790  if (fDmesonJets.size() != nAccCharm[0]+nAccCharm[1]) AliError(Form("I found %lu mesons (%d)?", fDmesonJets.size(), nAccCharm[0]+nAccCharm[1]));
1791  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1792  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1793  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1794  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1795 
1796  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1797  fHistManager->FillTH2(hname, npart, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1798 
1799  hname = TString::Format("%s/fHistNDmesons", GetName());
1800  fHistManager->FillTH1(hname, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]); // same as the number of accepted D mesons, since no selection is performed
1801 }
1802 
1807 {
1808  TString classname;
1809  if (fMCMode == kMCTruth) {
1810  classname = "AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary";
1811  fCurrentDmesonJetInfo = new AliDmesonMCInfoSummary();
1812  }
1813  else {
1814  switch (fCandidateType) {
1815  case kD0toKpi:
1816  case kD0toKpiLikeSign:
1817  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1818  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1819  break;
1820  case kDstartoKpipi:
1821  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1822  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1823  break;
1824  }
1825  }
1826  TString treeName = TString::Format("%s_%s", taskName, GetName());
1827  fTree = new TTree(treeName, treeName);
1828  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
1829  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1830  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1831  if (fJetDefinitions[i].fRhoName.IsNull()) {
1832  fCurrentJetInfo[i] = new AliJetInfoSummary();
1833  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1834  }
1835  else {
1836  fCurrentJetInfo[i] = new AliJetInfoPbPbSummary();
1837  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoPbPbSummary", &fCurrentJetInfo[i]);
1838  }
1839  }
1840 
1841  return fTree;
1842 }
1843 
1848 {
1849  TString hname;
1850 
1851  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
1852 
1853  for (auto &jetDef : fJetDefinitions) {
1854 
1855  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
1856 
1857  Double_t radius = jetDef.fRadius;
1858 
1859  TString title[30] = {""};
1860  Int_t nbins[30] = {0 };
1861  Double_t min [30] = {0.};
1862  Double_t max [30] = {0.};
1863  Int_t dim = 0 ;
1864 
1865  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1866  nbins[dim] = nPtBins;
1867  min[dim] = 0;
1868  max[dim] = fMaxPt;
1869  dim++;
1870 
1871  if ((enabledAxis & kPositionD) != 0) {
1872  title[dim] = "#eta_{D}";
1873  nbins[dim] = 50;
1874  min[dim] = -1;
1875  max[dim] = 1;
1876  dim++;
1877 
1878  title[dim] = "#phi_{D} (rad)";
1879  nbins[dim] = 150;
1880  min[dim] = 0;
1881  max[dim] = TMath::TwoPi();
1882  dim++;
1883  }
1884 
1885  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1886  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1887  nbins[dim] = fNMassBins;
1888  min[dim] = fMinMass;
1889  max[dim] = fMaxMass;
1890  dim++;
1891  }
1892 
1893  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1894  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1895  nbins[dim] = fNMassBins;
1896  min[dim] = fMinMass;
1897  max[dim] = fMaxMass;
1898  dim++;
1899  }
1900 
1901  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1902  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1903  nbins[dim] = fNMassBins;
1904  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1905  dim++;
1906  }
1907 
1908  if (fCandidateType == kDstartoKpipi) {
1909  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1910  nbins[dim] = fNMassBins*6;
1911  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1912 
1913  // subtract mass of D0
1914  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1915  min[dim] -= D0mass;
1916  max[dim] -= D0mass;
1917 
1918  dim++;
1919  }
1920 
1921  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1922  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1923  nbins[dim] = 100;
1924  min[dim] = 0;
1925  max[dim] = 25;
1926  dim++;
1927  }
1928 
1929  title[dim] = "#it{z}_{D}";
1930  nbins[dim] = 110;
1931  min[dim] = 0;
1932  max[dim] = 1.10;
1933  dim++;
1934 
1935  if ((enabledAxis & kDeltaR) != 0) {
1936  title[dim] = "#Delta R_{D-jet}";
1937  nbins[dim] = 100;
1938  min[dim] = 0;
1939  max[dim] = radius * 1.5;
1940  dim++;
1941  }
1942 
1943  if ((enabledAxis & kDeltaEta) != 0) {
1944  title[dim] = "#eta_{D} - #eta_{jet}";
1945  nbins[dim] = 100;
1946  min[dim] = -radius * 1.2;
1947  max[dim] = radius * 1.2;
1948  dim++;
1949  }
1950 
1951  if ((enabledAxis & kDeltaPhi) != 0) {
1952  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1953  nbins[dim] = 100;
1954  min[dim] = -radius * 1.2;
1955  max[dim] = radius * 1.2;
1956  dim++;
1957  }
1958 
1959  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1960  nbins[dim] = nPtBins;
1961  min[dim] = 0;
1962  max[dim] = fMaxPt;
1963  dim++;
1964 
1965  if ((enabledAxis & kPositionJet) != 0) {
1966  title[dim] = "#eta_{jet}";
1967  nbins[dim] = 50;
1968  min[dim] = -1;
1969  max[dim] = 1;
1970  dim++;
1971 
1972  title[dim] = "#phi_{jet} (rad)";
1973  nbins[dim] = 150;
1974  min[dim] = 0;
1975  max[dim] = TMath::TwoPi();
1976  dim++;
1977  }
1978 
1979  if ((enabledAxis & kJetConstituents) != 0) {
1980  title[dim] = "No. of constituents";
1981  nbins[dim] = 50;
1982  min[dim] = -0.5;
1983  max[dim] = 49.5;
1984  dim++;
1985  }
1986 
1987  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1988  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1989  for (Int_t j = 0; j < dim; j++) {
1990  h->GetAxis(j)->SetTitle(title[j]);
1991  }
1992  }
1993 }
1994 
1999 {
2000  TString hname;
2001  fPartons.clear();
2002 
2003  TH1* histAncestor = nullptr;
2004  TH1* histPrompt = nullptr;
2005 
2006  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2007  hname = TString::Format("%s/fHistPrompt", GetName());
2008  histPrompt = static_cast<TH1*>(fHistManager->FindObject(hname));
2009 
2010  hname = TString::Format("%s/fHistAncestor", GetName());
2011  histAncestor = static_cast<TH1*>(fHistManager->FindObject(hname));
2012  }
2013 
2014  for (auto& dmeson_pair : fDmesonJets) {
2015  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
2016  Int_t accJets = 0;
2017  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
2018  fCurrentJetInfo[ij]->Reset();
2019  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
2020  if (!jet) continue;
2021  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
2022  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
2023  fHistManager->FillTH1(hname, jet->Pt());
2024  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
2025  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2026  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
2027  fHistManager->FillTH1(hname, jet->Eta());
2028  continue;
2029  }
2030  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
2031  accJets++;
2032  }
2033  if (accJets > 0) {
2034  if (histPrompt) {
2035  if (dmeson_pair.second.fParton) {
2036  fPartons[dmeson_pair.second.fParton] = dmeson_pair.second.fPartonType;
2037  UInt_t absPdgParton = TMath::Abs(dmeson_pair.second.fParton->GetPdgCode());
2038  if (absPdgParton == 4) {
2039  histPrompt->Fill("Prompt", 1);
2040  }
2041  else if (absPdgParton == 5) {
2042  histPrompt->Fill("Non-Prompt", 1);
2043  }
2044  else {
2045  histPrompt->Fill("Unknown", 1);
2046  }
2047  }
2048  else {
2049  histPrompt->Fill("Unknown", 1);
2050  }
2051  }
2052 
2053  if (histAncestor) {
2054  if (dmeson_pair.second.fAncestor) {
2055  UInt_t absPdgAncestor = TMath::Abs(dmeson_pair.second.fAncestor->GetPdgCode());
2056  if (absPdgAncestor == 4) {
2057  histAncestor->Fill("Charm", 1);
2058  }
2059  else if (absPdgAncestor == 5) {
2060  histAncestor->Fill("Bottom", 1);
2061  }
2062  else if (absPdgAncestor == 2212) {
2063  histAncestor->Fill("Proton", 1);
2064  }
2065  else {
2066  histAncestor->Fill("Unknown", 1);
2067  }
2068  }
2069  else {
2070  histAncestor->Fill("Unknown", 1);
2071  }
2072  }
2073 
2074  fTree->Fill();
2075  }
2076  else {
2077  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2078  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2079  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2080  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2081  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2082  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2083  if (fMCMode != kMCTruth) {
2084  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
2085  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
2086  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
2087  }
2088  else if (fCandidateType == kDstartoKpipi) {
2089  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
2090  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
2091 
2092  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
2093  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
2094  }
2095  }
2096  }
2097  }
2098 
2099  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2100  hname = TString::Format("%s/fHistPartonPt", GetName());
2101  TH1* histPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
2102  hname = TString::Format("%s/fHistPartonEta", GetName());
2103  TH1* histPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
2104  hname = TString::Format("%s/fHistPartonPhi", GetName());
2105  TH1* histPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
2106  hname = TString::Format("%s/fHistPartonType", GetName());
2107  TH1* histPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
2108 
2109  for (auto parton : fPartons) {
2110  if (!parton.first) continue;
2111  histPartonPt->Fill(parton.first->Pt());
2112  histPartonEta->Fill(parton.first->Eta());
2113  histPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
2114  histPartonType->Fill(parton.second);
2115  }
2116  }
2117 
2118  return kTRUE;
2119 }
2120 
2126 {
2127  TString hname;
2128 
2129  TH1* histAncestor = nullptr;
2130  TH1* histPrompt = nullptr;
2131 
2132  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2133  hname = TString::Format("%s/fHistPrompt", GetName());
2134  histPrompt = static_cast<TH1*>(fHistManager->FindObject(hname));
2135 
2136  hname = TString::Format("%s/fHistAncestor", GetName());
2137  histAncestor = static_cast<TH1*>(fHistManager->FindObject(hname));
2138  }
2139 
2140  fPartons.clear();
2141  for (auto& dmeson_pair : fDmesonJets) {
2142  Int_t accJets = 0;
2143  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
2144  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
2145  if (!jet) continue;
2146  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
2147  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
2148  fHistManager->FillTH1(hname, jet->Pt());
2149  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
2150  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2151  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
2152  fHistManager->FillTH1(hname, jet->Eta());
2153  continue;
2154  }
2155  accJets++;
2156  }
2157  if (accJets > 0) {
2158  if (histPrompt) {
2159  if (dmeson_pair.second.fParton) {
2160  fPartons[dmeson_pair.second.fParton] = dmeson_pair.second.fPartonType;
2161  UInt_t absPdgParton = TMath::Abs(dmeson_pair.second.fParton->GetPdgCode());
2162  if (absPdgParton == 4) {
2163  histPrompt->Fill("Prompt", 1);
2164  }
2165  else if (absPdgParton == 5) {
2166  histPrompt->Fill("Non-Prompt", 1);
2167  }
2168  else {
2169  histPrompt->Fill("Unknown", 1);
2170  }
2171  }
2172  else {
2173  histPrompt->Fill("Unknown", 1);
2174  }
2175  }
2176 
2177  if (histAncestor) {
2178  if (dmeson_pair.second.fAncestor) {
2179  UInt_t absPdgAncestor = TMath::Abs(dmeson_pair.second.fAncestor->GetPdgCode());
2180  if (absPdgAncestor == 4) {
2181  histAncestor->Fill("Charm", 1);
2182  }
2183  else if (absPdgAncestor == 5) {
2184  histAncestor->Fill("Bottom", 1);
2185  }
2186  else if (absPdgAncestor == 2212) {
2187  histAncestor->Fill("Proton", 1);
2188  }
2189  else {
2190  histAncestor->Fill("Unknown", 1);
2191  }
2192  }
2193  else {
2194  histAncestor->Fill("Unknown", 1);
2195  }
2196  }
2197  }
2198  else {
2199  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2200  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2201  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2202  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2203  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2204  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2205  if (fMCMode != kMCTruth) {
2206  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
2207  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
2208  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
2209  }
2210  else if (fCandidateType == kDstartoKpipi) {
2211  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
2212  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
2213 
2214  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
2215  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
2216  }
2217  }
2218  }
2219  }
2220 
2221  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2222  hname = TString::Format("%s/fHistPartonPt", GetName());
2223  TH1* histPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
2224  hname = TString::Format("%s/fHistPartonEta", GetName());
2225  TH1* histPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
2226  hname = TString::Format("%s/fHistPartonPhi", GetName());
2227  TH1* histPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
2228  hname = TString::Format("%s/fHistPartonType", GetName());
2229  TH1* histPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
2230 
2231  for (auto parton : fPartons) {
2232  if (!parton.first) continue;
2233  histPartonPt->Fill(parton.first->Pt());
2234  histPartonEta->Fill(parton.first->Eta());
2235  histPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
2236  histPartonType->Fill(parton.second);
2237  }
2238  }
2239 
2240  return kTRUE;
2241 }
2242 
2247 {
2248  TString hname;
2249 
2250  for (auto& dmeson_pair : fDmesonJets) {
2251  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
2252  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2253  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2254  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2255  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2256  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2257  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2258  }
2259  }
2260 
2261  for (auto &jetDef : fJetDefinitions) {
2262 
2263  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
2264  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
2265 
2266  for (auto& dmeson_pair : fDmesonJets) {
2267  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
2268  if (!jet) continue;
2269  if (!jetDef.IsJetInAcceptance(*jet)) {
2270  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
2271  fHistManager->FillTH1(hname, jet->Pt());
2272  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
2273  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2274  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
2275  fHistManager->FillTH1(hname, jet->Eta());
2276  continue;
2277  }
2278  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
2279  }
2280  }
2281 
2282  return kTRUE;
2283 }
2284 
2291 {
2292  // Fill the THnSparse histogram.
2293 
2294  Double_t contents[30] = {0.};
2295 
2296  Double_t z = DmesonJet.GetZ(n);
2297  Double_t deltaPhi = 0;
2298  Double_t deltaEta = 0;
2299  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
2300 
2301  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
2302  if (it == DmesonJet.fJets.end()) return kFALSE;
2303 
2304  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
2305  TString title(h->GetAxis(i)->GetTitle());
2306  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
2307  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
2308  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
2309  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
2310  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
2311  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
2312  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
2313  else if (title=="#it{z}_{D}") contents[i] = z;
2314  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
2315  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
2316  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
2317  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
2318  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
2319  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
2320  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
2321  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
2322  }
2323 
2324  h->Fill(contents);
2325 
2326  return kTRUE;
2327 }
2328 
2329 // Definitions of class AliAnalysisTaskDmesonJets
2330 
2332 ClassImp(AliAnalysisTaskDmesonJets);
2334 
2338  fAnalysisEngines(),
2339  fEnabledAxis(0),
2341  fHistManager(),
2342  fApplyKinematicCuts(kTRUE),
2343  fNOutputTrees(0),
2344  fTrackEfficiency(0),
2345  fRejectISR(kFALSE),
2346  fJetAreaType(fastjet::active_area),
2347  fJetGhostArea(0.005),
2348  fMCContainer(0),
2349  fAodEvent(0),
2350  fFastJetWrapper(0)
2351 {
2352  SetMakeGeneralHistograms(kTRUE);
2353 }
2354 
2359  AliAnalysisTaskEmcalLight(name, kTRUE),
2360  fAnalysisEngines(),
2361  fEnabledAxis(k2ProngInvMass),
2362  fOutputType(kTreeOutput),
2363  fHistManager(name),
2364  fApplyKinematicCuts(kTRUE),
2365  fNOutputTrees(nOutputTrees),
2366  fTrackEfficiency(0),
2367  fRejectISR(kFALSE),
2368  fJetAreaType(fastjet::active_area),
2369  fJetGhostArea(0.005),
2370  fMCContainer(0),
2371  fAodEvent(0),
2372  fFastJetWrapper(0)
2373 {
2374  SetMakeGeneralHistograms(kTRUE);
2375  for (Int_t i = 0; i < nOutputTrees; i++){
2376  DefineOutput(2+i, TTree::Class());
2377  }
2378 }
2379 
2382 {
2383  if (fFastJetWrapper) delete fFastJetWrapper;
2384 }
2385 
2393 {
2394  AliRDHFCuts* analysiscuts = 0;
2395  TFile* filecuts = TFile::Open(cutfname);
2396  if (!filecuts || filecuts->IsZombie()) {
2397  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Input file not found: will use std cuts.");
2398  filecuts = 0;
2399  }
2400 
2401  if (filecuts) analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
2402 
2403  if (!analysiscuts) {
2404  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Could not find analysis cuts '%s' in '%s'.", cutsname.Data(), cutfname.Data());
2405  if (filecuts) {
2406  filecuts->ls();
2407  }
2408  }
2409  else {
2410  ::Info("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Cuts '%s' loaded from file '%s'", cutsname.Data(), cutfname.Data());
2411  }
2412 
2413  return analysiscuts;
2414 }
2415 
2425 {
2427  return AddAnalysisEngine(type, cutfname, MCmode, jetDef, rhoName);
2428 }
2429 
2439 {
2440  AliRDHFCuts* cuts = 0;
2441 
2442  if (!cutfname.IsNull()) {
2443  TString cutsname;
2444 
2445  switch (type) {
2446  case kD0toKpi:
2447  case kD0toKpiLikeSign:
2448  cutsname = "D0toKpiCuts";
2449  break;
2450  case kDstartoKpipi:
2451  cutsname = "DStartoKpipiCuts";
2452  break;
2453  default:
2454  return 0;
2455  }
2456 
2457  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
2458  if (cuts) cuts->PrintAll();
2459  }
2460 
2461  AnalysisEngine eng(type, MCmode, cuts);
2462 
2463  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
2464 
2465  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
2466  eng.AddJetDefinition(jetDef);
2467  it = fAnalysisEngines.insert(it, eng);
2468  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
2469  }
2470  else {
2471  AnalysisEngine* found_eng = &(*it);
2472  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "An analysis engine '%s' with %lu jet definitions has been found. The total number of analysis engines is %lu. A new jet definition '%s' is being added.", found_eng->GetName(), found_eng->fJetDefinitions.size(), fAnalysisEngines.size(), jetDef.GetName());
2473  found_eng->AddJetDefinition(jetDef);
2474 
2475  if (cuts) {
2476  if (found_eng->fRDHFCuts != 0) ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
2477  found_eng->SetRDHFCuts(cuts);
2478  }
2479  }
2480 
2481  return &(*it);
2482 }
2483 
2484 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
2485 {
2486  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
2487  while (it != fAnalysisEngines.end() && (*it) < eng) it++;
2488  return it;
2489 }
2490 
2493 {
2494  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
2495 
2497 
2498  // Define histograms
2499  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
2500 
2501  TString hname;
2502  TString htitle;
2503  TH1* h = 0;
2504  Int_t treeSlot = 0;
2505 
2506  Int_t maxTracks = 6000;
2507  Double_t maxRho = 500;
2508  if (fForceBeamType == kpp) {
2509  maxRho = 50;
2510  maxTracks = 200;
2511  }
2512  else if (fForceBeamType == kpA) {
2513  maxRho = 200;
2514  maxTracks = 500;
2515  }
2516 
2517  hname = "fHistCharmPt";
2518  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2519  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2520 
2521  hname = "fHistCharmEta";
2522  htitle = hname + ";#eta_{charm};counts";
2523  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2524 
2525  hname = "fHistCharmPhi";
2526  htitle = hname + ";#phi_{charm};counts";
2527  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2528 
2529  hname = "fHistCharmPt_Eta05";
2530  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2531  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2532 
2533  hname = "fHistBottomPt";
2534  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2535  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2536 
2537  hname = "fHistBottomEta";
2538  htitle = hname + ";#eta_{bottom};counts";
2539  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2540 
2541  hname = "fHistBottomPhi";
2542  htitle = hname + ";#phi_{bottom};counts";
2543  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2544 
2545  hname = "fHistBottomPt_Eta05";
2546  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2547  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2548 
2549  hname = "fHistHighestPartonPt";
2550  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2551  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2552 
2553  hname = "fHistHighestPartonType";
2554  htitle = hname + ";type;counts";
2555  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2556 
2557  hname = "fHistNHeavyQuarks";
2558  htitle = hname + ";number of heavy-quarks;counts";
2559  fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
2560 
2561  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
2562  for (auto &param : fAnalysisEngines) {
2563  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
2564 
2565  param.fHistManager = &fHistManager;
2566 
2567  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", param.GetName());
2568  htitle = hname + ";#it{N}_{tracks};#it{N}_{D};events";
2569  h = fHistManager.CreateTH2(hname, htitle, 251, -0.5, 250.5, 21, -0.5, 20.5);
2570 
2571  hname = TString::Format("%s/fHistNTotAcceptedDmesons", param.GetName());
2572  htitle = hname + ";;#it{N}_{D}";
2573  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
2574 
2575  hname = TString::Format("%s/fHistNDmesons", param.GetName());
2576  htitle = hname + ";#it{N}_{D};events";
2577  h = fHistManager.CreateTH1(hname, htitle, 501, -0.5, 500.5);
2578 
2579  hname = TString::Format("%s/fHistNEvents", param.GetName());
2580  htitle = hname + ";Event status;counts";
2581  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
2582  h->GetXaxis()->SetBinLabel(1, "Accepted");
2583  h->GetXaxis()->SetBinLabel(2, "Rejected");
2584 
2585  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
2586  htitle = hname + ";Rejection reason;counts";
2587  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
2588 
2589  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
2590  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
2591  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2592 
2593  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
2594  htitle = hname + ";#it{#eta}_{D};counts";
2595  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2596 
2597  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
2598  htitle = hname + ";#it{#phi}_{D};counts";
2599  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2600 
2601  if (param.fMCMode != kMCTruth) {
2602  if (param.fCandidateType == kD0toKpi || param.fCandidateType == kD0toKpiLikeSign) {
2603  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
2604  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2605  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
2606  }
2607  else if (param.fCandidateType == kDstartoKpipi) {
2608  Double_t min = 0;
2609  Double_t max = 0;
2610 
2611  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
2612  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2613  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
2614  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
2615 
2616  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2617  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
2618  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2619  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
2620  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
2621  }
2622  }
2623 
2624  if (param.fMCMode == kMCTruth) {
2625  hname = TString::Format("%s/fHistPartonPt", param.GetName());
2626  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
2627  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2628 
2629  hname = TString::Format("%s/fHistPartonEta", param.GetName());
2630  htitle = hname + ";#eta_{parton};counts";
2631  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2632 
2633  hname = TString::Format("%s/fHistPartonPhi", param.GetName());
2634  htitle = hname + ";#phi_{parton};counts";
2635  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2636 
2637  hname = TString::Format("%s/fHistPartonType", param.GetName());
2638  htitle = hname + ";type;counts";
2639  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2640 
2641  hname = TString::Format("%s/fHistPrompt", param.GetName());
2642  htitle = hname + ";Type;counts";
2643  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
2644  h->GetXaxis()->SetBinLabel(1, "Unknown");
2645  h->GetXaxis()->SetBinLabel(2, "Prompt");
2646  h->GetXaxis()->SetBinLabel(3, "Non-Prompt");
2647 
2648  hname = TString::Format("%s/fHistAncestor", param.GetName());
2649  htitle = hname + ";Ancestor;counts";
2650  h = fHistManager.CreateTH1(hname, htitle, 4, 0, 4);
2651  h->GetXaxis()->SetBinLabel(1, "Unknown");
2652  h->GetXaxis()->SetBinLabel(2, "Charm");
2653  h->GetXaxis()->SetBinLabel(3, "Bottom");
2654  h->GetXaxis()->SetBinLabel(4, "Proton");
2655  }
2656 
2657  for (auto& jetDef : param.fJetDefinitions) {
2658  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef.GetName());
2659 
2660  if (param.fMCMode == kMCTruth) {
2661  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", param.GetName(), jetDef.GetName());
2662  htitle = hname + ";#it{N}_{constituents};#it{N}_{D};counts";
2663  h = fHistManager.CreateTH2(hname, htitle, 51, -0.5, 50.5, 10, 0.5, 10.5);
2664  }
2665 
2666  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef.GetName());
2667  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2668  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2669  SetRejectionReasonLabels(h->GetXaxis());
2670 
2671  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef.GetName());
2672  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2673  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2674  SetRejectionReasonLabels(h->GetXaxis());
2675 
2676  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef.GetName());
2677  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
2678  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2679  SetRejectionReasonLabels(h->GetXaxis());
2680 
2681  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef.GetName());
2682  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
2683  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
2684 
2685  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef.GetName());
2686  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
2687  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2688 
2689  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef.GetName());
2690  htitle = hname + ";#it{#eta}_{jet};counts";
2691  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2692 
2693  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef.GetName());
2694  htitle = hname + ";#it{#phi}_{jet};counts";
2695  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2696 
2697  if (!jetDef.fRhoName.IsNull()) {
2698  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", param.GetName(), jetDef.GetName());
2699  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});#rho (GeV/#it{c} #times rad^{-1});counts";
2700  fHistManager.CreateTH2(hname, htitle, 300, 0, 150, 1000, 0, maxRho);
2701 
2702  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", param.GetName(), jetDef.GetName());
2703  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});#rho (GeV/#it{c} #times rad^{-1});counts";
2704  fHistManager.CreateTH2(hname, htitle, 300, 0, 150, 1000, 0, maxRho);
2705 
2706  hname = TString::Format("%s/%s/fHistRhoVsCent", param.GetName(), jetDef.GetName());
2707  htitle = hname + ";Centrality (%);#rho (GeV/#it{c} #times rad^{-1});counts";
2708  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 1000, 0, maxRho);
2709 
2710  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", param.GetName(), jetDef.GetName());
2711  htitle = hname + ";Centrality (%);#it{p}_{T,jet} (GeV/#it{c});counts";
2712  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
2713 
2714  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", param.GetName(), jetDef.GetName());
2715  htitle = hname + ";Centrality (%);#it{p}_{T,D} (GeV/#it{c});counts";
2716  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
2717 
2718  hname = TString::Format("%s/%s/fHistRhoVsNTracks", param.GetName(), jetDef.GetName());
2719  htitle = hname + ";no. of tracks;#rho (GeV/#it{c} #times rad^{-1});counts";
2720  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 1000, 0, maxRho);
2721 
2722  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", param.GetName(), jetDef.GetName());
2723  htitle = hname + ";no. of tracks;#it{p}_{T,jet} (GeV/#it{c});counts";
2724  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
2725 
2726  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", param.GetName(), jetDef.GetName());
2727  htitle = hname + ";no. of tracks;#it{p}_{T,D} (GeV/#it{c});counts";
2728  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
2729  }
2730  }
2731  switch (fOutputType) {
2732  case kTreeOutput:
2733  param.BuildTree(GetName());
2734  if (treeSlot < fNOutputTrees) {
2735  param.AssignDataSlot(treeSlot+2);
2736  treeSlot++;
2738  }
2739  else {
2740  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
2741  }
2742  break;
2743  case kTHnOutput:
2744  param.BuildHnSparse(fEnabledAxis);
2745  break;
2746  case kNoOutput:
2747  break;
2748  }
2749  }
2750 
2752 
2753  PostData(1, fOutput);
2754 }
2755 
2759 {
2761 
2762  // Load the event
2763  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
2764 
2765  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
2766 
2767  fFastJetWrapper->SetAreaType((fastjet::AreaType)fJetAreaType);
2769 
2770  if (!fAodEvent) {
2771  AliError(Form("This task need an AOD event (Task '%s'). Expect troubles...", GetName()));
2772  //return;
2773  }
2774 
2775  TRandom* rnd = 0;
2776  if (fTrackEfficiency > 0 && fTrackEfficiency < 1) rnd = new TRandom3(0);
2777 
2778  for (auto cont_it : fParticleCollArray) {
2779  AliHFAODMCParticleContainer* part_cont = dynamic_cast<AliHFAODMCParticleContainer*>(cont_it.second);
2780  if (part_cont) fMCContainer = part_cont;
2781  }
2782 
2783  for (auto &params : fAnalysisEngines) {
2784 
2785  params.fAodEvent = fAodEvent;
2786  params.fFastJetWrapper = fFastJetWrapper;
2787  params.fTrackEfficiency = fTrackEfficiency;
2788  params.fRejectISR = fRejectISR;
2789  params.fRandomGen = rnd;
2790 
2791  for (auto &jetdef: params.fJetDefinitions) {
2792  if (!jetdef.fRhoName.IsNull()) {
2793  jetdef.fRho = dynamic_cast<AliRhoParameter*>(fInputEvent->FindListObject(jetdef.fRhoName));
2794  if (!jetdef.fRho) {
2795  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2796  "%s: Could not find rho object '%s' for engine '%s'",
2797  GetName(), jetdef.fRhoName.Data(), params.GetName());
2798  }
2799  }
2800  }
2801 
2802  if (!params.fRDHFCuts) {
2803  if (params.fMCMode == kMCTruth) {
2804  ::Warning("AliAnalysisTaskDmesonJets::ExecOnce",
2805  "%s: RDHF cuts not provided for engine '%s'.",
2806  GetName(), params.GetName());
2807  }
2808  else {
2809  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2810  "%s: RDHF cuts not provided. Engine '%s' disabled.",
2811  GetName(), params.GetName());
2812  params.fInhibit = kTRUE;
2813  }
2814  }
2815 
2816  params.fMCContainer = fMCContainer;
2817 
2818  for (auto cont_it : fParticleCollArray) {
2819  AliTrackContainer* track_cont = dynamic_cast<AliTrackContainer*>(cont_it.second);
2820  if (track_cont) params.fTrackContainers.push_back(track_cont);
2821  }
2822 
2823  for (auto cont_it : fClusterCollArray) params.fClusterContainers.push_back(cont_it.second);
2824 
2825  if (fAodEvent) params.Init(fGeom, fAodEvent->GetRunNumber());
2826 
2827  if (params.fMCMode != kMCTruth && fAodEvent) {
2828  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
2829 
2830  if (params.fCandidateArray) {
2831  TString className;
2832  if (params.fCandidateType == kD0toKpi || params.fCandidateType == kD0toKpiLikeSign) {
2833  className = "AliAODRecoDecayHF2Prong";
2834  }
2835  else if (params.fCandidateType == kDstartoKpipi) {
2836  className = "AliAODRecoCascadeHF";
2837  }
2838  if (!params.fCandidateArray->GetClass()->InheritsFrom(className)) {
2839  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2840  "%s: Objects of type %s in %s are not inherited from %s! Task will be disabled!",
2841  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName(), className.Data());
2842  params.fCandidateArray = 0;
2843  params.fInhibit = kTRUE;
2844  }
2845  }
2846  else {
2847  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2848  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
2849  params.fBranchName.Data(), params.GetName());
2850  params.fInhibit = kTRUE;
2851  }
2852  }
2853 
2854  if (params.fMCMode != kNoMC) {
2855  if (!params.fMCContainer) {
2856  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2857  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
2858  params.GetName());
2859  params.fInhibit = kTRUE;
2860  }
2861  }
2862 
2863  if (params.fMCMode != kMCTruth) {
2864  if (params.fTrackContainers.size() == 0 && params.fClusterContainers.size() == 0) {
2865  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2866  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
2867  params.GetName());
2868  params.fInhibit = kTRUE;
2869  }
2870  }
2871  }
2872 }
2873 
2878 {
2879  TString hname;
2880 
2881  // Fix for temporary bug in ESDfilter
2882  // The AODs with null vertex pointer didn't pass the PhysSel
2883  // Now adding an entry in the histogram so as to check that this is actually cutting anything out
2884  if (fAodEvent && (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001)) {
2885  for (auto &eng : fAnalysisEngines) {
2886  if (eng.fInhibit) continue;
2887  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2888  fHistManager.FillTH1(hname, "ESDfilterBug");
2889  }
2890  return kFALSE;
2891  }
2892 
2893  if (fAodEvent) {
2894  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
2895  if (matchingAODdeltaAODlevel <= 0) {
2896  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
2897  for (auto &eng : fAnalysisEngines) {
2898  if (eng.fInhibit) continue;
2899  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2900  fHistManager.FillTH1(hname, "MismatchDeltaAOD");
2901  }
2902  return kFALSE;
2903  }
2904  }
2905 
2906  for (auto &eng : fAnalysisEngines) {
2907  eng.fDmesonJets.clear();
2908  if (eng.fInhibit) continue;
2909 
2910  eng.fCent = fCent;
2911 
2912  //Event selection
2913  hname = TString::Format("%s/fHistNEvents", eng.GetName());
2914  if (fAodEvent) {
2915  Bool_t iseventselected = kTRUE;
2916  if (eng.fRDHFCuts) iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
2917  if (!iseventselected) {
2918  fHistManager.FillTH1(hname, "Rejected");
2919  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2920  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
2921  TString label;
2922  do {
2923  label = GetHFEventRejectionReasonLabel(bitmap);
2924  if (label.IsNull()) break;
2925  fHistManager.FillTH1(hname, label);
2926  } while (true);
2927 
2928  continue;
2929  }
2930  }
2931 
2932  fHistManager.FillTH1(hname, "Accepted");
2933 
2934  AliDebug(2, "Event selected");
2935 
2936  eng.RunAnalysis();
2937  }
2938  return kTRUE;
2939 }
2940 
2945 {
2946  for (auto &param : fAnalysisEngines) {
2947  if (param.fInhibit) continue;
2948 
2949  if (fOutputType == kTreeOutput) {
2950  param.FillTree(fApplyKinematicCuts);
2952  }
2953  else if (fOutputType == kTHnOutput) {
2954  param.FillHnSparse(fApplyKinematicCuts);
2955  }
2956  }
2958  return kTRUE;
2959 }
2960 
2963 {
2964  auto itcont = fMCContainer->all_momentum();
2965  Int_t nHQ = 0;
2966  Double_t highestPartonPt = 0;
2967  Int_t absPdgHighParton = 0;
2968  for (auto part : itcont) {
2969  Int_t absPdgCode = TMath::Abs(part.second->GetPdgCode());
2970 
2971  // Skip all particles that are not either quarks or gluons
2972  if (absPdgCode > 9 && absPdgCode != 21) continue;
2973 
2974  // Look for highest momentum parton
2975  if (highestPartonPt < part.first.Pt()) {
2976  highestPartonPt = part.first.Pt();
2977  absPdgHighParton = absPdgCode;
2978  }
2979  /*
2980  // Look for the mother PDG code
2981  Int_t motherIndex = part.second->GetMother();
2982  AliAODMCParticle *mother = 0;
2983  Int_t motherPdg = 0;
2984  Double_t motherPt = 0;
2985  if (motherIndex >= 0) {
2986  mother = fMCContainer->GetMCParticle(motherIndex);
2987  if (motherIndex) {
2988  motherPdg = TMath::Abs(mother->GetPdgCode());
2989  motherPt = mother->Pt();
2990  }
2991  }
2992  */
2993  if (absPdgCode != 4 && absPdgCode != 5) continue;
2994  Bool_t notLastInPartonShower = kFALSE;
2995  for (Int_t idaugh = 0; idaugh < 2; idaugh++){
2996  Int_t daughterIndex = part.second->GetDaughter(idaugh);
2997  if (daughterIndex < 0) {
2998  AliDebug(10, Form("Could not find daughter of heavy quark (pdg=%d, pt=%.3f)!", absPdgCode, part.first.Pt()));
2999  continue;
3000  }
3001  AliAODMCParticle *daughter = fMCContainer->GetMCParticle(daughterIndex);
3002  if (!daughter) {
3003  AliDebug(10, Form("Could not find daughter %d of heavy quark (pdg=%d, pt=%.3f)!", daughterIndex, absPdgCode, part.first.Pt()));
3004  continue;
3005  }
3006  Int_t daughterAbsPdgCode = TMath::Abs(daughter->GetPdgCode());
3007  if (daughterAbsPdgCode <= 9 || daughterAbsPdgCode == 21) notLastInPartonShower = kTRUE; // this parton is not the last parton in the shower
3008  AliDebug(10, Form("Found daughter with PDG=%d, pt=%.3f", daughterAbsPdgCode, daughter->Pt()));
3009  }
3010  if (notLastInPartonShower) continue;
3011 
3012  if (absPdgCode == 4) {
3013  fHistManager.FillTH1("fHistCharmPt", part.first.Pt());
3014  fHistManager.FillTH1("fHistCharmEta", part.first.Eta());
3015  fHistManager.FillTH1("fHistCharmPhi", part.first.Phi_0_2pi());
3016  if (TMath::Abs(part.first.Eta()) < 0.5) fHistManager.FillTH1("fHistCharmPt_Eta05", part.first.Pt());
3017  }
3018  else if (absPdgCode == 5) {
3019  fHistManager.FillTH1("fHistBottomPt", part.first.Pt());
3020  fHistManager.FillTH1("fHistBottomEta", part.first.Eta());
3021  fHistManager.FillTH1("fHistBottomPhi", part.first.Phi_0_2pi());
3022  if (TMath::Abs(part.first.Eta()) < 0.5) fHistManager.FillTH1("fHistBottomPt_Eta05", part.first.Pt());
3023  }
3024  nHQ++;
3025  }
3026  fHistManager.FillTH1("fHistNHeavyQuarks", nHQ);
3027  fHistManager.FillTH1("fHistHighestPartonPt",highestPartonPt);
3028  Int_t partonType = 0;
3029  if (absPdgHighParton == 9 || absPdgHighParton == 21) {
3030  partonType = 7;
3031  }
3032  else {
3033  partonType = absPdgHighParton;
3034  }
3035  fHistManager.FillTH1("fHistHighestPartonType",partonType);
3036 }
3037 
3046 {
3047  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
3048 
3049  Double_t mass = part->Mass();
3050 
3051  // To make sure that the PDG mass value is not at the edge of a bin
3052  if (nbins % 2 == 0) {
3053  minMass = mass - range / 2 - range / nbins / 2;
3054  maxMass = mass + range / 2 - range / nbins / 2;
3055  }
3056  else {
3057  minMass = mass - range / 2;
3058  maxMass = mass + range / 2;
3059  }
3060 }
3061 
3068 {
3069  static TString label;
3070  label = "";
3071 
3072  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
3073  label = "NotSelTrigger";
3074  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
3075  return label.Data();
3076  }
3077  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
3078  label = "NoVertex";
3079  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
3080  return label.Data();
3081  }
3082  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
3083  label = "TooFewVtxContrib";
3084  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
3085  return label.Data();
3086  }
3087  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
3088  label = "ZVtxOutFid";
3089  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
3090  return label.Data();
3091  }
3092  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
3093  label = "Pileup";
3094  bitmap &= ~BIT(AliRDHFCuts::kPileup);
3095  return label.Data();
3096  }
3097  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
3098  label = "OutsideCentrality";
3099  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
3100  return label.Data();
3101  }
3102  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
3103  label = "PhysicsSelection";
3104  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
3105  return label.Data();
3106  }
3107  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
3108  label = "BadSPDVertex";
3109  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
3110  return label.Data();
3111  }
3112  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
3113  label = "ZVtxSPDOutFid";
3114  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
3115  return label.Data();
3116  }
3117  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
3118  label = "CentralityFlattening";
3119  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
3120  return label.Data();
3121  }
3122  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
3123  label = "BadTrackV0Correl";
3124  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
3125  return label.Data();
3126  }
3127 
3128  return label.Data();
3129 }
3130 
3137 {
3138  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
3139  PostData(eng.GetDataSlotNumber(), eng.GetTree());
3140  return eng.GetDataSlotNumber();
3141  }
3142  else {
3143  return -1;
3144  }
3145 }
3146 
3156 {
3157  // Get the pointer to the existing analysis manager via the static access method
3158  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
3159  if (!mgr) {
3160  ::Error("AddTaskDmesonJets", "No analysis manager to connect to.");
3161  return NULL;
3162  }
3163 
3164  // Check the analysis type using the event handlers connected to the analysis manager
3165  AliVEventHandler* handler = mgr->GetInputEventHandler();
3166  if (!handler) {
3167  ::Error("AddTaskEmcalJetSpectraQA", "This task requires an input event handler");
3168  return NULL;
3169  }
3170 
3171  EDataType_t dataType = kUnknownDataType;
3172 
3173  if (handler->InheritsFrom("AliESDInputHandler")) {
3174  dataType = kESD;
3175  }
3176  else if (handler->InheritsFrom("AliAODInputHandler")) {
3177  dataType = kAOD;
3178  }
3179 
3180  // Init the task and do settings
3181  if (ntracks == "usedefault") {
3182  if (dataType == kESD) {
3183  ntracks = "Tracks";
3184  }
3185  else if (dataType == kAOD) {
3186  ntracks = "tracks";
3187  }
3188  else {
3189  ntracks = "";
3190  }
3191  }
3192 
3193  if (nclusters == "usedefault") {
3194  if (dataType == kESD) {
3195  nclusters = "CaloClusters";
3196  }
3197  else if (dataType == kAOD) {
3198  nclusters = "caloClusters";
3199  }
3200  else {
3201  nclusters = "";
3202  }
3203  }
3204 
3205  if (nMCpart == "usedefault") {
3206  nMCpart = "mcparticles"; // Always needs AliAODMCParticle objects
3207  }
3208 
3209  TString name("AliAnalysisTaskDmesonJets");
3210  if (strcmp(suffix, "") != 0) {
3211  name += TString::Format("_%s", suffix.Data());
3212  }
3213 
3214  AliAnalysisTaskDmesonJets* jetTask = new AliAnalysisTaskDmesonJets(name, nMaxTrees);
3215 
3216  if (!ntracks.IsNull()) {
3217  AliHFTrackContainer* trackCont = new AliHFTrackContainer(ntracks);
3218  jetTask->AdoptParticleContainer(trackCont);
3219  }
3220 
3221  if (!nMCpart.IsNull()) {
3222  AliMCParticleContainer* partCont = new AliHFAODMCParticleContainer(nMCpart);
3223  partCont->SetEtaLimits(-1.5, 1.5);
3224  partCont->SetPtLimits(0, 1000);
3225  jetTask->AdoptParticleContainer(partCont);
3226  }
3227 
3228  jetTask->AddClusterContainer(nclusters.Data());
3229 
3230  // Final settings, pass to manager and set the containers
3231  mgr->AddTask(jetTask);
3232 
3233  // Create containers for input/output
3234  AliAnalysisDataContainer* cinput1 = mgr->GetCommonInputContainer();
3235  TString contname1(name);
3236  contname1 += "_histos";
3237  AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(contname1.Data(),
3238  TList::Class(), AliAnalysisManager::kOutputContainer,
3239  Form("%s", AliAnalysisManager::GetCommonFileName()));
3240 
3241  mgr->ConnectInput(jetTask, 0, cinput1);
3242  mgr->ConnectOutput(jetTask, 1, coutput1);
3243 
3244  for (Int_t i = 0; i < nMaxTrees; i++) {
3245  TString contname = TString::Format("%s_tree_%d", name.Data(), i);
3246  AliAnalysisDataContainer *coutput = mgr->CreateContainer(contname.Data(),
3247  TTree::Class(),AliAnalysisManager::kOutputContainer,
3248  Form("%s", AliAnalysisManager::GetCommonFileName()));
3249  mgr->ConnectOutput(jetTask, 2+i, coutput);
3250  }
3251  return jetTask;
3252 }
3253 
Int_t pdg
AliClusterContainer * AddClusterContainer(std::string branchName, std::string contName="")
void Print() const
Prints the content of this object in the standard output.
AliAODMCParticle * fParton
! pointer to the parton in the shower tree of the D meson (only for particle level D mesons) ...
std::list< AnalysisEngine > fAnalysisEngines
Array of analysis parameters.
AliEMCALGeometry * fGeom
!emcal geometry
void SetDMesonCandidate(AliAODRecoDecay *c)
Analysis task for D meson jets.
virtual void Set(const AliDmesonJetInfo &source)
Invariant mass of the D0 meson candidate in GeV/c2.
virtual Int_t Run()
double Double_t
Definition: External.C:58
Bool_t IsJetInAcceptance(const AliJetInfo &jet) const
Double_t fSoftPionPt
! Transverse momentum of the soft pion of the D* candidate
EDataType_t
Switch for the data type.
virtual AliAODMCParticle * GetMCParticle(Int_t i=-1) const
const char * title
Definition: MakeQAPdf.C:27
AliRDHFCuts * fRDHFCuts
D meson candidates cuts.
virtual void Set(const AliDmesonJetInfo &source)
Lightweight class that encapsulates D meson jets.
Double_t fTrackEfficiency
Artificial tracking inefficiency (0...1)
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
AliHFAODMCParticleContainer * fMCContainer
! MC particle container
ERecoScheme_t fRecoScheme
Jet recombination scheme (pt scheme, E scheme, ...)
const TObjArray & GetDaughterList() const
Double_t fMaxNeutralPt
Transverse momentum of the leading neutral particle (or cluster)
void RunDetectorLevelAnalysis()
Run a detector level analysis.
Double_t mass
Container with name, TClonesArray and cuts for particles.
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
virtual void UserCreateOutputObjects()
Creates the output containers.
void FillPartonLevelHistograms()
Fill histograms with parton-level information.
static Int_t CheckMatchingAODdeltaAODevents()
Lightweight class that encapsulates D meson jets.
Double_t fInvMass2Prong
! 2-prong mass of the D* candidate (w/o the soft pion)
static void CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t &minMass, Double_t &maxMass)
AliAODEvent * fAodEvent
! AOD event
Bool_t FindJet(AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, AliHFJetDefinition &jetDef)
UInt_t fEnabledAxis
Use bit defined in EAxis_t to enable axis in the THnSparse.
AliRhoParameter * fRho
Object that holds the average background value.
EMCMode_t fMCMode
MC mode: No MC (data and MC detector level), background-only (MC), signal-only (MC), MC truth (particle level)
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
Double_t fJetGhostArea
Area of the ghost particles.
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:122
AnalysisEngine()
This is the default constructor, used for ROOT I/O purposes.
void Init(const AliEMCALGeometry *const geom, Int_t runNumber)
Initialize the analysis engine.
static FJRecoScheme ConvertToFJRecoScheme(ERecoScheme_t reco)
AliVParticle * fDmesonParticle
! pointer to the particle object
static TString GenerateJetName(EJetType_t jetType, EJetAlgo_t jetAlgo, ERecoScheme_t recoScheme, Double_t radius, AliParticleContainer *partCont, AliClusterContainer *clusCont, TString tag)
static std::pair< AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle * > IsPromptCharm(const AliAODMCParticle *part, TClonesArray *mcArray)
void SetR(Double_t r)
Definition: AliFJWrapper.h:127
Select tracks based on specific prescriptions of HF analysis.
ECandidateType_t fCandidateType
Candidate type.
Double_t GetDistance(std::string n, Double_t &deta, Double_t &dphi) const
Int_t fMCLabel
! MC label, i.e. index of the generator level D meson (only for detector level D meson candidates) ...
EBeamType_t fForceBeamType
forced beam type
const Int_t nPtBins
virtual void Set(const AliDmesonJetInfo &source)
void AddInputVectors(AliEmcalContainer *cont, Int_t offset, TH2 *rejectHist=0, Double_t eff=0.)
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:34
EJetType_t fJetType
Jet type (charged, full, neutral)
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
Create a new TH2 within the container.
virtual void Set(const AliDmesonJetInfo &source)
AliTLorentzVector fD
! 4-momentum of the D meson candidate
Double_t GetDistance(const AliJetInfo &jet, Double_t &deta, Double_t &dphi) const
TObject * FindObject(const char *name) const
Find an object inside the container.
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
AliAnalysisTaskDmesonJets()
This is the default constructor, used for ROOT I/O purposes.
virtual void Reset()
Reset the current object.
AnalysisEngine & operator=(const AnalysisEngine &source)
Lightweight class that encapsulates D meson jets for PbPb analysis.
Struct that encapsulates analysis parameters.
std::vector< AliHFJetDefinition >::iterator FindJetDefinition(const AliHFJetDefinition &eng)
Base task in the EMCAL framework (lighter version of AliAnalysisTaskEmcal)
Double_t Phi_0_2pi() const
static AliAODMCParticle * FindParticleOrigin(const AliAODMCParticle *part, TClonesArray *mcArray, EFindParticleOriginMode_t mode, const std::set< UInt_t > &pdgSet)
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:121
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
AliAODTrack * GetBachelor() const
AliRDHFCuts * LoadDMesonCutsFromFile(TString cutfname, TString cutsname)
bool operator<(const AliAnalysisTaskDmesonJets::AliHFJetDefinition &lhs, const AliAnalysisTaskDmesonJets::AliHFJetDefinition &rhs)
std::list< AnalysisEngine >::iterator FindAnalysisEngine(const AnalysisEngine &eng)
std::vector< AliHFJetDefinition > fJetDefinitions
Jet definitions.
EJetAlgo_t fJetAlgo
Jet algorithm (kt, anti-kt,...)
static const char * GetHFEventRejectionReasonLabel(UInt_t &bitmap)
Lightweight class that encapsulates D*.
virtual void Clear(const Option_t *="")
Int_t mode
Definition: anaM.C:41
Bool_t ExtractDstarAttributes(const AliAODRecoCascadeHF *DstarCand, AliDmesonJetInfo &DmesonJet, UInt_t i)
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
const AliJetInfo * GetJet(std::string n) const
Short_t fPartonType
! type of the parton in the shower tree (only for particle level D mesons)
const AliMCParticleIterableMomentumContainer all_momentum() const
virtual ~AliAnalysisTaskDmesonJets()
This is the standard destructor.
Byte_t fSelectionType
! for D0: 0=not selected, 1=D0, 2=D0bar, 3=both
Bool_t ExtractD0Attributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
Select MC particles based on specific prescriptions of HF analysis.
static EMesonDecayChannel_t CheckDecayChannel(const AliAODMCParticle *part, TClonesArray *mcArray)
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:125
AliDmesonJetInfo & operator=(const AliDmesonJetInfo &source)
Double_t fMaxChargedPt
Transverse momentum of the leading charged particle (or track)
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
static FJJetAlgo ConvertToFJAlgo(EJetAlgo_t algo)
decay
Definition: HFPtSpectrum.C:41
std::map< std::string, AliJetInfo > fJets
! list of jets
AliAODMCParticle * fAncestor
! pointer to the ancestor particle in the shower tree of the D meson (only for particle level D meson...
virtual void PrintAll() const
Definition: External.C:220
Double_t minMass
void RunParticleLevelAnalysis()
Run a particle level analysis.
std::map< std::string, AliParticleContainer * > fParticleCollArray
particle/track collection array
Bool_t fApplyKinematicCuts
Apply jet kinematic cuts.
AnalysisEngine * AddAnalysisEngine(ECandidateType_t type, TString cutfname, EMCMode_t bkgMode, EJetType_t jettype, Double_t jetradius, TString rhoName="")
const char * GetName() const
Generate a name for this jet definition.
virtual void Set(const AliDmesonJetInfo &source, std::string n)
unsigned short UShort_t
Definition: External.C:28
Bool_t ExtractRecoDecayAttributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
Double_t InvMassD0() const
Bool_t fRejectISR
Reject initial state radiation.
EOutputType_t fOutputType
Output type: none, TTree or THnSparse.
virtual void AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
bool operator==(const AliAnalysisTaskDmesonJets::AliHFJetDefinition &lhs, const AliAnalysisTaskDmesonJets::AliHFJetDefinition &rhs)
Double_t InvMassDstarKpipi() const
const Int_t nbins
Double_t maxMass
bool Bool_t
Definition: External.C:53
Int_t fNOutputTrees
Maximum number of output trees.
AliHFJetDefinition()
This is the default constructor, used for ROOT I/O purposes.
AliHFJetDefinition & operator=(const AliHFJetDefinition &source)
Int_t GetNAcceptedParticles() const
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:123
AliTLorentzVector fMomentum
4-momentum of the jet
std::map< std::string, AliClusterContainer * > fClusterCollArray
cluster collection array
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
Create a new THnSparse within the container.
Double_t fCorrPt
Transverse momentum of the jet after subtracting the average background.
AliHFJetDefinition * AddJetDefinition(EJetType_t type, Double_t r, EJetAlgo_t algo, ERecoScheme_t reco)
Container for MC-true particles within the EMCAL framework.
Lightweight class that encapsulates D0.
Bool_t IsAnyJetInAcceptance(const AliDmesonJetInfo &dMesonJet) const
void AdoptParticleContainer(AliParticleContainer *cont)
Double_t fCent
!event centrality
void SetCharge(EChargeCut_t c)
Container for jet within the EMCAL jet framework.
Definition: External.C:196
Int_t PostDataFromAnalysisEngine(const AnalysisEngine &eng)
AliFJWrapper * fFastJetWrapper
! Fastjet wrapper
THistManager fHistManager
Histogram manager.
void Reset()
Reset all fields to their default values.
static AliAnalysisTaskDmesonJets * AddTaskDmesonJets(TString ntracks="usedefault", TString nclusters="usedefault", TString nMCpart="", Int_t nMaxTrees=2, TString suffix="")
virtual void Set(const AliDmesonJetInfo &source, std::string n)