AliPhysics  6eb37c2 (6eb37c2)
 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(0.),
634  fMinJetPhi(0.),
635  fMaxJetPhi(0.),
636  fMinJetEta(0.),
637  fMaxJetEta(0.),
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 (fMinJetPt < fMaxJetPt && (jet.Pt() > fMaxJetPt || jet.Pt() < fMinJetPt)) return kFALSE;
699  if (fMinChargedPt < fMaxChargedPt && (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt)) return kFALSE;
700  if (fMinNeutralPt < fMaxNeutralPt && (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  fName = fCandidateName;
999  switch (fMCMode) {
1000  case kBackgroundOnly:
1001  fName += "_kBackgroundOnly";
1002  break;
1003  case kSignalOnly:
1004  fName += "_kSignalOnly";
1005  break;
1006  case kMCTruth:
1007  fName += "_MCTruth";
1008  break;
1009  case kWrongPID:
1010  fName += "_WrongPID";
1011  break;
1012  default:
1013  break;
1014  }
1015 
1016  if (fRDHFCuts) fName += TString::Format("_%s", fRDHFCuts->GetName());
1017 
1018  return fName.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 
1079 {
1080  for (auto &jetdef : fJetDefinitions) jetdef.SetJetPhiRange(min, max);
1081 }
1082 
1087 {
1088  for (auto &jetdef : fJetDefinitions) jetdef.SetJetEtaRange(min, max);
1089 }
1090 
1095 {
1096  for (auto &jetdef : fJetDefinitions) jetdef.SetJetPtRange(min, max);
1097 }
1098 
1103 {
1104  for (auto &jetdef : fJetDefinitions) jetdef.SetChargedPtRange(min, max);
1105 }
1106 
1111 {
1112  for (auto &jetdef : fJetDefinitions) jetdef.SetNeutralPtRange(min, max);
1113 }
1114 
1121 {
1122  if (lhs.fCandidateType < rhs.fCandidateType) {
1123  return true;
1124  }
1125  else if (lhs.fCandidateType > rhs.fCandidateType) {
1126  return false;
1127  }
1128  else if (lhs.fMCMode < rhs.fMCMode) {
1129  return true;
1130  }
1131  else if (lhs.fMCMode > rhs.fMCMode) {
1132  return false;
1133  }
1134  else if (lhs.fRDHFCuts && !rhs.fRDHFCuts) {
1135  return true;
1136  }
1137  else if (lhs.fRDHFCuts && rhs.fRDHFCuts && strcmp(lhs.fRDHFCuts->GetName(), rhs.fRDHFCuts->GetName()) < 0) {
1138  return true;
1139  }
1140  else {
1141  return false;
1142  }
1143 }
1144 
1151 {
1152  if (lhs.fCandidateType != rhs.fCandidateType) return false;
1153  if (lhs.fMCMode != rhs.fMCMode) return false;
1154  if (lhs.fRDHFCuts == nullptr && rhs.fRDHFCuts != nullptr) return false;
1155  if (lhs.fRDHFCuts != nullptr && rhs.fRDHFCuts == nullptr) return false;
1156  if (lhs.fRDHFCuts && rhs.fRDHFCuts && strcmp(lhs.fRDHFCuts->GetName(), rhs.fRDHFCuts->GetName()) != 0) return false;
1157  return true;
1158 }
1159 
1168 {
1169  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) { // D0 candidate
1170  return ExtractD0Attributes(Dcand, DmesonJet, i);
1171  }
1172  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
1173  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
1174  }
1175  else {
1176  return kFALSE;
1177  }
1178 }
1179 
1188 {
1189  AliDebug(10,"Checking if D0 meson is selected");
1190  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
1191  if (isSelected == 0) return kFALSE;
1192 
1193  Int_t MCtruthPdgCode = 0;
1194 
1195  Double_t invMassD = 0;
1196 
1197  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
1198  // Checks also the origin, and if it matches the rejected origin mask, return false
1199  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly || fMCMode == kWrongPID) {
1200  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
1201  DmesonJet.fMCLabel = mcLab;
1202 
1203  // Retrieve the generated particle (if exists) and its PDG code
1204  if (mcLab >= 0) {
1205  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1206 
1207  if (aodMcPart) {
1208  // Check origin and return false if it matches the rejected origin mask
1209  if (fRejectedOrigin) {
1210  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
1211  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1212  }
1213  MCtruthPdgCode = aodMcPart->PdgCode();
1214  }
1215  }
1216  }
1217 
1218  if (isSelected == 1) { // selected as a D0
1219  if (i != 0) return kFALSE; // only one mass hypothesis thanks to PID
1220 
1221  if (fMCMode == kNoMC ||
1222  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1223  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly) ||
1224  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kWrongPID)) {
1225  // 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)
1226  AliDebug(10,"Selected as D0");
1227  invMassD = Dcand->InvMassD0();
1228  }
1229  else { // conditions above not passed, so return FALSE
1230  return kFALSE;
1231  }
1232  }
1233  else if (isSelected == 2) { // selected as a D0bar
1234  if (i != 1) return kFALSE; // only one mass hypothesis thanks to PID
1235 
1236  if (fMCMode == kNoMC ||
1237  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1238  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly) ||
1239  (MCtruthPdgCode == fCandidatePDG && fMCMode == kWrongPID)) {
1240  // 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)
1241  AliDebug(10,"Selected as D0bar");
1242  invMassD = Dcand->InvMassD0bar();
1243  }
1244  else { // conditions above not passed, so return FALSE
1245  return kFALSE;
1246  }
1247  }
1248  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
1249  AliDebug(10,"Selected as either D0 or D0bar");
1250 
1251  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
1252  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1253  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kWrongPID))) {
1254  if (i != 0) return kFALSE;
1255  AliDebug(10, "MC truth is D0");
1256  invMassD = Dcand->InvMassD0();
1257  }
1258  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1259  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kWrongPID))) {
1260  if (i != 1) return kFALSE;
1261  AliDebug(10, "MC truth is D0bar");
1262  invMassD = Dcand->InvMassD0bar();
1263  }
1264  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
1265 
1266  // Only accept it if background-only OR background-and-signal was requested
1267  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
1268  // Select D0 or D0bar depending on the i-parameter
1269  if (i == 0) {
1270  AliDebug(10, "Returning invariant mass with D0 hypothesis");
1271  invMassD = Dcand->InvMassD0();
1272  }
1273  else if (i == 1) {
1274  AliDebug(10, "Returning invariant mass with D0bar hypothesis");
1275  invMassD = Dcand->InvMassD0bar();
1276  }
1277  else { // i > 1
1278  return kFALSE;
1279  }
1280  }
1281  else { // signal-only was requested but this is not a true D0
1282  return kFALSE;
1283  }
1284  }
1285  }
1286  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
1287  return kTRUE;
1288 }
1289 
1298 {
1299  AliDebug(10,"Checking if D* meson is selected");
1300  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoCascadeHF*>(DstarCand), AliRDHFCuts::kAll, fAodEvent);
1301  if (isSelected == 0) return kFALSE;
1302 
1303  if ((i == 1 && DstarCand->Charge()>0) || (i == 0 && DstarCand->Charge()<0) || i > 1) return kFALSE; // only one mass hypothesis for the D*
1304 
1305  Int_t MCtruthPdgCode = 0;
1306 
1307  Double_t invMassD = 0;
1308 
1309  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1310  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1311  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1312 
1313  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1314  AliDebug(10, Form("MC label is %d", mcLab));
1315  DmesonJet.fMCLabel = mcLab;
1316  if (mcLab >= 0) {
1317  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1318 
1319  if (aodMcPart) {
1320  if (fRejectedOrigin) {
1321  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
1322  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1323  }
1324 
1325  MCtruthPdgCode = aodMcPart->PdgCode();
1326  AliDebug(10, Form("MC truth pdg code is %d",MCtruthPdgCode));
1327  }
1328  }
1329  }
1330 
1331  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1332  if (fMCMode == kNoMC ||
1333  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1334  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1335  // 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)
1336  invMassD = DstarCand->InvMassDstarKpipi();
1337  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1338  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1339  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1340  return kTRUE;
1341  }
1342  else { // conditions above not passed, so return FALSE
1343  return kFALSE;
1344  }
1345 }
1346 
1354 {
1355  if (!part) return kUnknownDecay;
1356  if (!mcArray) return kUnknownDecay;
1357 
1359 
1360  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1361 
1362  if (part->GetNDaughters() == 2) {
1363 
1364  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1365  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1366 
1367  if (!d1 || !d2) {
1368  return decay;
1369  }
1370 
1371  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1372  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1373 
1374  if (absPdgPart == 421) { // D0 -> K pi
1375  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1376  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1377  decay = kDecayD0toKpi;
1378  }
1379  }
1380 
1381  if (absPdgPart == 413) { // D* -> D0 pi
1382  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1383  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1384  if (D0decay == kDecayD0toKpi) {
1385  decay = kDecayDStartoKpipi;
1386  }
1387  }
1388  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1389  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1390  if (D0decay == kDecayD0toKpi) {
1391  decay = kDecayDStartoKpipi;
1392  }
1393  }
1394  }
1395  }
1396 
1397  return decay;
1398 }
1399 
1406 std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> AliAnalysisTaskDmesonJets::AnalysisEngine::IsPromptCharm(const AliAODMCParticle* part, TClonesArray* mcArray)
1407 {
1408  std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> result(kUnknownQuark, 0);
1409 
1410  if (!part) return result;
1411  if (!mcArray) return result;
1412 
1413  static std::set<UInt_t> partons = { 4, 5 };
1414 
1415  AliAODMCParticle* parton = FindParticleOrigin(part, mcArray, kFindLast, partons);
1416  if (parton) {
1417  result.second = parton;
1418  UInt_t absPdgParton = TMath::Abs(parton->GetPdgCode());
1419  if (absPdgParton == 4) result.first = kFromCharm;
1420  else if (absPdgParton == 5) result.first = kFromBottom;
1421  }
1422 
1423  return result;
1424 }
1425 
1433 
1434 AliAODMCParticle* AliAnalysisTaskDmesonJets::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode)
1435 {
1436  static std::set<UInt_t> pdgSet;
1437 
1438  return FindParticleOrigin(part, mcArray, mode, pdgSet);
1439 }
1440 
1449 
1450 AliAODMCParticle* AliAnalysisTaskDmesonJets::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode, const std::set<UInt_t>& pdgSet)
1451 {
1452  AliAODMCParticle* result = nullptr;
1453 
1454  Int_t mother = part->GetMother();
1455  while (mother >= 0) {
1456  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1457  if (mcGranma) {
1458  UInt_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
1459 
1460  // If the current particle is one of the particle types that is being searched assign it to the result pointer
1461  if (pdgSet.empty() || pdgSet.count(abspdgGranma) > 0) {
1462  result = mcGranma;
1463  // If the last particle in the fragmentation tree (first when going reverse) was requested then stop the loop
1464  if (mode == kFindLast) break;
1465  }
1466  mother = mcGranma->GetMother();
1467  }
1468  else {
1469  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::FindParticleOrigin", "Could not retrieve mother particle %d!", mother);
1470  break;
1471  }
1472  }
1473 
1474  return result;
1475 }
1476 
1479 {
1480  for (auto& jetDef : fJetDefinitions) {
1481  jetDef.fJets.clear();
1482  }
1483 
1484  if (fMCMode == kMCTruth) {
1485  RunParticleLevelAnalysis();
1486  }
1487  else {
1488  RunDetectorLevelAnalysis();
1489  }
1490 }
1491 
1494 {
1495  // Fill the vertex info of the candidates
1496  // Needed for reduced delta AOD, where the vertex info has been deleted
1497  // to reduce the delta AOD file size
1499 
1500  const Int_t nD = fCandidateArray->GetEntriesFast();
1501 
1502  AliDmesonJetInfo DmesonJet;
1503 
1504  std::map<AliHFJetDefinition*,Double_t> maxJetPt;
1505  for (auto& def : fJetDefinitions) maxJetPt[&def] = 0;
1506  Double_t maxDPt = 0;
1507 
1508  Int_t nAccCharm[3] = {0};
1509  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1510  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1511  if (!charmCand) continue;
1512  if(!(vHF.FillRecoCand(fAodEvent,charmCand))) continue;
1513 
1514  // region of interest + cuts
1515  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1516  Int_t nMassHypo = 0; // number of mass hypothesis accepted for this D meson
1517  if (charmCand->Pt() > maxDPt) maxDPt = charmCand->Pt();
1518  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1519  DmesonJet.Reset();
1520  DmesonJet.fDmesonParticle = charmCand;
1521  DmesonJet.fSelectionType = im + 1;
1522  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1523  for (auto& def : fJetDefinitions) {
1524  if (FindJet(charmCand, DmesonJet, def)) {
1525  Double_t jetPt = DmesonJet.fJets[def.GetName()].fMomentum.Pt();
1526  if (jetPt > maxJetPt[&def]) maxJetPt[&def] = jetPt;
1527  }
1528  else {
1529  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1530  def.GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1531  }
1532  }
1533  fDmesonJets[(icharm+1)*(1-(im*2))] = DmesonJet;
1534  nMassHypo++;
1535  nAccCharm[im]++;
1536  }
1537  }
1538  if (nMassHypo == 2) {
1539  nAccCharm[0]--;
1540  nAccCharm[1]--;
1541  nAccCharm[2] += 2;
1542  }
1543  if (nMassHypo == 2) { // both mass hypothesis accepted
1544  fDmesonJets[(icharm+1)].fSelectionType = 3;
1545  fDmesonJets[-(icharm+1)].fSelectionType = 3;
1546  }
1547  } // end of D cand loop
1548 
1549  TString hname;
1550 
1551  Int_t ntracks = 0;
1552 
1553  for (auto track_cont : fTrackContainers) {
1554  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
1555  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(nullptr);
1556  ntracks += track_cont->GetNAcceptEntries();
1557  }
1558 
1559  for (auto& def : fJetDefinitions) {
1560  if (!def.fRho) continue;
1561  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
1562  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
1563 
1564  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
1565  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
1566 
1567  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
1568  fHistManager->FillTH2(hname, fCent, def.fRho->GetVal());
1569 
1570  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
1571  fHistManager->FillTH2(hname, fCent, maxJetPt[&def]);
1572 
1573  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
1574  fHistManager->FillTH2(hname, fCent, maxDPt);
1575 
1576  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
1577  fHistManager->FillTH2(hname, ntracks, def.fRho->GetVal());
1578 
1579  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
1580  fHistManager->FillTH2(hname, ntracks, maxJetPt[&def]);
1581 
1582  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
1583  fHistManager->FillTH2(hname, ntracks, maxDPt);
1584  }
1585 
1586  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1587  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1588  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1589  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1590 
1591  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1592  fHistManager->FillTH2(hname, ntracks, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1593 
1594  hname = TString::Format("%s/fHistNDmesons", GetName());
1595  fHistManager->FillTH1(hname, nD);
1596 }
1597 
1608 {
1609  TString hname;
1610 
1611  Double_t rho = 0;
1612  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
1613 
1615  fFastJetWrapper->SetR(jetDef.fRadius);
1618 
1619  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1620 
1621  if (jetDef.fJetType != AliJetContainer::kNeutralJet) {
1622  for (auto track_cont : fTrackContainers) {
1623  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
1624  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(Dcand);
1625  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1626  AddInputVectors(track_cont, 100, static_cast<TH2*>(fHistManager->FindObject(hname)), fTrackEfficiency);
1627 
1628  if (hftrack_cont) {
1629  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1630  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1631  const TObjArray& daughters = hftrack_cont->GetDaughterList();
1632  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1633  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1634  if (!hftrack_cont->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1635  }
1636  }
1637  }
1638  }
1639 
1640  if (jetDef.fJetType != AliJetContainer::kChargedJet) {
1641  for (auto clus_cont : fClusterContainers) {
1642  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1643  AddInputVectors(clus_cont, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1644  }
1645  }
1646 
1647  // run jet finder
1648  fFastJetWrapper->Run();
1649 
1650  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1651 
1652  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1653  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1654 
1655  Bool_t isDmesonJet = kFALSE;
1656 
1657  Double_t maxChPt = 0;
1658  Double_t maxNePt = 0;
1659  Double_t totalNeutralPt = 0;
1660  Int_t nConst = 1;
1661 
1662  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1663  if (constituents[ic].user_index() == 0) {
1664  isDmesonJet = kTRUE;
1665  }
1666  else if (constituents[ic].user_index() >= 100) {
1667  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1668  nConst++;
1669  }
1670  else if (constituents[ic].user_index() <= -100) {
1671  totalNeutralPt += constituents[ic].pt();
1672  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1673  nConst++;
1674  }
1675  }
1676 
1677  if (isDmesonJet) {
1678  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1679  DmesonJet.fJets[jetDef.GetName()].fNConstituents = nConst;
1680  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1681  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1682  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1683  DmesonJet.fJets[jetDef.GetName()].fArea = jets_incl[ijet].area();
1684  DmesonJet.fJets[jetDef.GetName()].fCorrPt = DmesonJet.fJets[jetDef.GetName()].fMomentum.Pt() - jets_incl[ijet].area() * rho;
1685 
1686  return kTRUE;
1687  }
1688  }
1689 
1690  return kFALSE;
1691 }
1692 
1696 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist, Double_t eff)
1697 {
1698  auto itcont = cont->all_momentum();
1699  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1700  UInt_t rejectionReason = 0;
1701  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1702  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1703  continue;
1704  }
1705  if (fRandomGen && eff > 0 && eff < 1) {
1706  Double_t rnd = fRandomGen->Rndm();
1707  if (eff < rnd) {
1708  if (rejectHist) rejectHist->Fill(6, it->first.Pt());
1709  continue;
1710  }
1711  }
1712  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1713  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1714  }
1715 }
1716 
1719 {
1720  TString hname;
1721 
1722  fMCContainer->SetSpecialPDG(fCandidatePDG);
1723  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1724  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1726 
1727  if (!fMCContainer->IsSpecialPDGFound()) return;
1728 
1729  Int_t nAccCharm[3] = {0};
1730 
1731  std::map<AliHFJetDefinition*, Double_t> maxJetPt;
1732  Double_t maxDPt = 0;
1733 
1734  for (auto &jetDef : fJetDefinitions) {
1735  maxJetPt[&jetDef] = 0;
1736  Double_t rho = 0;
1737  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
1738  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", GetName(), jetDef.GetName());
1739  TH1* histNDmesonsVsNconstituents = static_cast<TH1*>(fHistManager->FindObject(hname));
1740 
1741  switch (jetDef.fJetType) {
1743  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
1744  break;
1746  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
1747  break;
1749  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
1750  break;
1751  }
1752 
1754  fFastJetWrapper->SetR(jetDef.fRadius);
1757 
1758  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1759  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1760 
1761  fFastJetWrapper->Run();
1762 
1763  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1764 
1765  for (auto jet : jets_incl) {
1766  Int_t nDmesonsInJet = 0;
1767 
1768  for (auto constituent : jet.constituents()) {
1769  Int_t iPart = constituent.user_index() - 100;
1770  if (constituent.perp() < 1e-6) continue; // reject ghost particles
1771  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
1772  if (!part) {
1773  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1774  continue;
1775  }
1776  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1777  nDmesonsInJet++;
1778  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1779  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1780  if (part->Pt() > maxDPt) maxDPt = part->Pt();
1781  std::pair<int, AliDmesonJetInfo> element;
1782  element.first = iPart;
1783  dMesonJetIt = fDmesonJets.insert(element).first;
1784  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1785  (*dMesonJetIt).second.fDmesonParticle = part;
1786  (*dMesonJetIt).second.fSelectionType = part->PdgCode() > 0 ? 1 : 2;
1787 
1788  UShort_t p = 0;
1789  UInt_t rs = 0;
1790 
1791  auto origin = IsPromptCharm(part, fMCContainer->GetArray());
1792  p = 0;
1793  rs = origin.first;
1794  while (rs >>= 1) { p++; }
1795  (*dMesonJetIt).second.fPartonType = p;
1796  (*dMesonJetIt).second.fParton = origin.second;
1797 
1798  (*dMesonJetIt).second.fAncestor = FindParticleOrigin(part, fMCContainer->GetArray(), kFindFirst);
1799 
1800  if (part->PdgCode() > 0) { // D0
1801  nAccCharm[0]++;
1802  }
1803  else { // D0bar
1804  nAccCharm[1]++;
1805  }
1806  }
1807 
1808  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
1809  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = jet.constituents().size();
1810  (*dMesonJetIt).second.fJets[jetDef.GetName()].fArea = jet.area();
1811  (*dMesonJetIt).second.fJets[jetDef.GetName()].fCorrPt = (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.Pt() - jet.area() * rho;
1812  if (jet.perp() > maxJetPt[&jetDef]) maxJetPt[&jetDef] = jet.perp();
1813  } // if constituent is a D meson
1814  } // for each constituent
1815  if (nDmesonsInJet > 0) histNDmesonsVsNconstituents->Fill(jet.constituents().size(), nDmesonsInJet);
1816  } // for each jet
1817  } // for each jet definition
1818 
1820 
1821  for (auto& def : fJetDefinitions) {
1822  if (!def.fRho) continue;
1823  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
1824  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
1825 
1826  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
1827  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
1828 
1829  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
1830  fHistManager->FillTH2(hname, fCent, def.fRho->GetVal());
1831 
1832  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
1833  fHistManager->FillTH2(hname, fCent, maxJetPt[&def]);
1834 
1835  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
1836  fHistManager->FillTH2(hname, fCent, maxDPt);
1837 
1838  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
1839  fHistManager->FillTH2(hname, npart, def.fRho->GetVal());
1840 
1841  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
1842  fHistManager->FillTH2(hname, npart, maxJetPt[&def]);
1843 
1844  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
1845  fHistManager->FillTH2(hname, npart, maxDPt);
1846  }
1847 
1848  if (fDmesonJets.size() != nAccCharm[0]+nAccCharm[1]) AliError(Form("I found %lu mesons (%d)?", fDmesonJets.size(), nAccCharm[0]+nAccCharm[1]));
1849  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1850  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1851  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1852  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1853 
1854  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1855  fHistManager->FillTH2(hname, npart, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1856 
1857  hname = TString::Format("%s/fHistNDmesons", GetName());
1858  fHistManager->FillTH1(hname, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]); // same as the number of accepted D mesons, since no selection is performed
1859 }
1860 
1865 {
1866  TString classname;
1867  if (fMCMode == kMCTruth) {
1868  classname = "AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary";
1869  fCurrentDmesonJetInfo = new AliDmesonMCInfoSummary();
1870  }
1871  else {
1872  switch (fCandidateType) {
1873  case kD0toKpi:
1874  case kD0toKpiLikeSign:
1875  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1876  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1877  break;
1878  case kDstartoKpipi:
1879  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1880  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1881  break;
1882  }
1883  }
1884  TString treeName = TString::Format("%s_%s", taskName, GetName());
1885  fTree = new TTree(treeName, treeName);
1886  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
1887  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1888  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1889  if (fJetDefinitions[i].fRhoName.IsNull()) {
1890  fCurrentJetInfo[i] = new AliJetInfoSummary();
1891  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1892  }
1893  else {
1894  fCurrentJetInfo[i] = new AliJetInfoPbPbSummary();
1895  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoPbPbSummary", &fCurrentJetInfo[i]);
1896  }
1897  }
1898 
1899  return fTree;
1900 }
1901 
1906 {
1907  TString hname;
1908 
1909  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
1910 
1911  for (auto &jetDef : fJetDefinitions) {
1912 
1913  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
1914 
1915  Double_t radius = jetDef.fRadius;
1916 
1917  TString title[30] = {""};
1918  Int_t nbins[30] = {0 };
1919  Double_t min [30] = {0.};
1920  Double_t max [30] = {0.};
1921  Int_t dim = 0 ;
1922 
1923  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1924  nbins[dim] = nPtBins;
1925  min[dim] = 0;
1926  max[dim] = fMaxPt;
1927  dim++;
1928 
1929  if ((enabledAxis & kPositionD) != 0) {
1930  title[dim] = "#eta_{D}";
1931  nbins[dim] = 50;
1932  min[dim] = -1;
1933  max[dim] = 1;
1934  dim++;
1935 
1936  title[dim] = "#phi_{D} (rad)";
1937  nbins[dim] = 150;
1938  min[dim] = 0;
1939  max[dim] = TMath::TwoPi();
1940  dim++;
1941  }
1942 
1943  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1944  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1945  nbins[dim] = fNMassBins;
1946  min[dim] = fMinMass;
1947  max[dim] = fMaxMass;
1948  dim++;
1949  }
1950 
1951  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1952  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1953  nbins[dim] = fNMassBins;
1954  min[dim] = fMinMass;
1955  max[dim] = fMaxMass;
1956  dim++;
1957  }
1958 
1959  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1960  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1961  nbins[dim] = fNMassBins;
1962  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1963  dim++;
1964  }
1965 
1966  if (fCandidateType == kDstartoKpipi) {
1967  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1968  nbins[dim] = fNMassBins*6;
1969  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1970 
1971  // subtract mass of D0
1972  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1973  min[dim] -= D0mass;
1974  max[dim] -= D0mass;
1975 
1976  dim++;
1977  }
1978 
1979  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1980  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1981  nbins[dim] = 100;
1982  min[dim] = 0;
1983  max[dim] = 25;
1984  dim++;
1985  }
1986 
1987  title[dim] = "#it{z}_{D}";
1988  nbins[dim] = 110;
1989  min[dim] = 0;
1990  max[dim] = 1.10;
1991  dim++;
1992 
1993  if ((enabledAxis & kDeltaR) != 0) {
1994  title[dim] = "#Delta R_{D-jet}";
1995  nbins[dim] = 100;
1996  min[dim] = 0;
1997  max[dim] = radius * 1.5;
1998  dim++;
1999  }
2000 
2001  if ((enabledAxis & kDeltaEta) != 0) {
2002  title[dim] = "#eta_{D} - #eta_{jet}";
2003  nbins[dim] = 100;
2004  min[dim] = -radius * 1.2;
2005  max[dim] = radius * 1.2;
2006  dim++;
2007  }
2008 
2009  if ((enabledAxis & kDeltaPhi) != 0) {
2010  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
2011  nbins[dim] = 100;
2012  min[dim] = -radius * 1.2;
2013  max[dim] = radius * 1.2;
2014  dim++;
2015  }
2016 
2017  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
2018  nbins[dim] = nPtBins;
2019  min[dim] = 0;
2020  max[dim] = fMaxPt;
2021  dim++;
2022 
2023  if ((enabledAxis & kPositionJet) != 0) {
2024  title[dim] = "#eta_{jet}";
2025  nbins[dim] = 50;
2026  min[dim] = -1;
2027  max[dim] = 1;
2028  dim++;
2029 
2030  title[dim] = "#phi_{jet} (rad)";
2031  nbins[dim] = 150;
2032  min[dim] = 0;
2033  max[dim] = TMath::TwoPi();
2034  dim++;
2035  }
2036 
2037  if ((enabledAxis & kJetConstituents) != 0) {
2038  title[dim] = "No. of constituents";
2039  nbins[dim] = 50;
2040  min[dim] = -0.5;
2041  max[dim] = 49.5;
2042  dim++;
2043  }
2044 
2045  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
2046  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
2047  for (Int_t j = 0; j < dim; j++) {
2048  h->GetAxis(j)->SetTitle(title[j]);
2049  }
2050  }
2051 }
2052 
2057 {
2058  TString hname;
2059  fPartons.clear();
2060 
2061  TH1* histAncestor = nullptr;
2062  TH1* histPrompt = nullptr;
2063 
2064  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2065  hname = TString::Format("%s/fHistPrompt", GetName());
2066  histPrompt = static_cast<TH1*>(fHistManager->FindObject(hname));
2067 
2068  hname = TString::Format("%s/fHistAncestor", GetName());
2069  histAncestor = static_cast<TH1*>(fHistManager->FindObject(hname));
2070  }
2071 
2072  for (auto& dmeson_pair : fDmesonJets) {
2073  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
2074  Int_t accJets = 0;
2075  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
2076  fCurrentJetInfo[ij]->Reset();
2077  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
2078  if (!jet) continue;
2079  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
2080  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
2081  fHistManager->FillTH1(hname, jet->Pt());
2082  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
2083  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2084  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
2085  fHistManager->FillTH1(hname, jet->Eta());
2086  continue;
2087  }
2088  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
2089  accJets++;
2090  }
2091  if (accJets > 0) {
2092  if (histPrompt) {
2093  if (dmeson_pair.second.fParton) {
2094  fPartons[dmeson_pair.second.fParton] = dmeson_pair.second.fPartonType;
2095  UInt_t absPdgParton = TMath::Abs(dmeson_pair.second.fParton->GetPdgCode());
2096  if (absPdgParton == 4) {
2097  histPrompt->Fill("Prompt", 1);
2098  }
2099  else if (absPdgParton == 5) {
2100  histPrompt->Fill("Non-Prompt", 1);
2101  }
2102  else {
2103  histPrompt->Fill("Unknown", 1);
2104  }
2105  }
2106  else {
2107  histPrompt->Fill("Unknown", 1);
2108  }
2109  }
2110 
2111  if (histAncestor) {
2112  if (dmeson_pair.second.fAncestor) {
2113  UInt_t absPdgAncestor = TMath::Abs(dmeson_pair.second.fAncestor->GetPdgCode());
2114  if (absPdgAncestor == 4) {
2115  histAncestor->Fill("Charm", 1);
2116  }
2117  else if (absPdgAncestor == 5) {
2118  histAncestor->Fill("Bottom", 1);
2119  }
2120  else if (absPdgAncestor == 2212) {
2121  histAncestor->Fill("Proton", 1);
2122  }
2123  else {
2124  histAncestor->Fill("Unknown", 1);
2125  }
2126  }
2127  else {
2128  histAncestor->Fill("Unknown", 1);
2129  }
2130  }
2131 
2132  fTree->Fill();
2133  }
2134  else {
2135  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2136  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2137  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2138  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2139  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2140  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2141  if (fMCMode != kMCTruth) {
2142  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
2143  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
2144  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
2145  }
2146  else if (fCandidateType == kDstartoKpipi) {
2147  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
2148  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
2149 
2150  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
2151  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
2152  }
2153  }
2154  }
2155  }
2156 
2157  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2158  hname = TString::Format("%s/fHistPartonPt", GetName());
2159  TH1* histPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
2160  hname = TString::Format("%s/fHistPartonEta", GetName());
2161  TH1* histPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
2162  hname = TString::Format("%s/fHistPartonPhi", GetName());
2163  TH1* histPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
2164  hname = TString::Format("%s/fHistPartonType", GetName());
2165  TH1* histPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
2166 
2167  for (auto parton : fPartons) {
2168  if (!parton.first) continue;
2169  histPartonPt->Fill(parton.first->Pt());
2170  histPartonEta->Fill(parton.first->Eta());
2171  histPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
2172  histPartonType->Fill(parton.second);
2173  }
2174  }
2175 
2176  return kTRUE;
2177 }
2178 
2184 {
2185  TString hname;
2186 
2187  TH1* histAncestor = nullptr;
2188  TH1* histPrompt = nullptr;
2189 
2190  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2191  hname = TString::Format("%s/fHistPrompt", GetName());
2192  histPrompt = static_cast<TH1*>(fHistManager->FindObject(hname));
2193 
2194  hname = TString::Format("%s/fHistAncestor", GetName());
2195  histAncestor = static_cast<TH1*>(fHistManager->FindObject(hname));
2196  }
2197 
2198  fPartons.clear();
2199  for (auto& dmeson_pair : fDmesonJets) {
2200  Int_t accJets = 0;
2201  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
2202  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
2203  if (!jet) continue;
2204  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
2205  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
2206  fHistManager->FillTH1(hname, jet->Pt());
2207  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
2208  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2209  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
2210  fHistManager->FillTH1(hname, jet->Eta());
2211  continue;
2212  }
2213  accJets++;
2214  }
2215  if (accJets > 0) {
2216  if (histPrompt) {
2217  if (dmeson_pair.second.fParton) {
2218  fPartons[dmeson_pair.second.fParton] = dmeson_pair.second.fPartonType;
2219  UInt_t absPdgParton = TMath::Abs(dmeson_pair.second.fParton->GetPdgCode());
2220  if (absPdgParton == 4) {
2221  histPrompt->Fill("Prompt", 1);
2222  }
2223  else if (absPdgParton == 5) {
2224  histPrompt->Fill("Non-Prompt", 1);
2225  }
2226  else {
2227  histPrompt->Fill("Unknown", 1);
2228  }
2229  }
2230  else {
2231  histPrompt->Fill("Unknown", 1);
2232  }
2233  }
2234 
2235  if (histAncestor) {
2236  if (dmeson_pair.second.fAncestor) {
2237  UInt_t absPdgAncestor = TMath::Abs(dmeson_pair.second.fAncestor->GetPdgCode());
2238  if (absPdgAncestor == 4) {
2239  histAncestor->Fill("Charm", 1);
2240  }
2241  else if (absPdgAncestor == 5) {
2242  histAncestor->Fill("Bottom", 1);
2243  }
2244  else if (absPdgAncestor == 2212) {
2245  histAncestor->Fill("Proton", 1);
2246  }
2247  else {
2248  histAncestor->Fill("Unknown", 1);
2249  }
2250  }
2251  else {
2252  histAncestor->Fill("Unknown", 1);
2253  }
2254  }
2255  }
2256  else {
2257  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2258  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2259  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2260  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2261  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2262  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2263  if (fMCMode != kMCTruth) {
2264  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
2265  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
2266  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
2267  }
2268  else if (fCandidateType == kDstartoKpipi) {
2269  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
2270  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
2271 
2272  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
2273  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
2274  }
2275  }
2276  }
2277  }
2278 
2279  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2280  hname = TString::Format("%s/fHistPartonPt", GetName());
2281  TH1* histPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
2282  hname = TString::Format("%s/fHistPartonEta", GetName());
2283  TH1* histPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
2284  hname = TString::Format("%s/fHistPartonPhi", GetName());
2285  TH1* histPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
2286  hname = TString::Format("%s/fHistPartonType", GetName());
2287  TH1* histPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
2288 
2289  for (auto parton : fPartons) {
2290  if (!parton.first) continue;
2291  histPartonPt->Fill(parton.first->Pt());
2292  histPartonEta->Fill(parton.first->Eta());
2293  histPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
2294  histPartonType->Fill(parton.second);
2295  }
2296  }
2297 
2298  return kTRUE;
2299 }
2300 
2305 {
2306  TString hname;
2307 
2308  for (auto& dmeson_pair : fDmesonJets) {
2309  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
2310  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2311  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2312  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2313  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2314  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2315  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2316  }
2317  }
2318 
2319  for (auto &jetDef : fJetDefinitions) {
2320 
2321  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
2322  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
2323 
2324  for (auto& dmeson_pair : fDmesonJets) {
2325  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
2326  if (!jet) continue;
2327  if (!jetDef.IsJetInAcceptance(*jet)) {
2328  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
2329  fHistManager->FillTH1(hname, jet->Pt());
2330  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
2331  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2332  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
2333  fHistManager->FillTH1(hname, jet->Eta());
2334  continue;
2335  }
2336  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
2337  }
2338  }
2339 
2340  return kTRUE;
2341 }
2342 
2349 {
2350  // Fill the THnSparse histogram.
2351 
2352  Double_t contents[30] = {0.};
2353 
2354  Double_t z = DmesonJet.GetZ(n);
2355  Double_t deltaPhi = 0;
2356  Double_t deltaEta = 0;
2357  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
2358 
2359  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
2360  if (it == DmesonJet.fJets.end()) return kFALSE;
2361 
2362  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
2363  TString title(h->GetAxis(i)->GetTitle());
2364  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
2365  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
2366  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
2367  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
2368  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
2369  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
2370  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
2371  else if (title=="#it{z}_{D}") contents[i] = z;
2372  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
2373  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
2374  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
2375  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
2376  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
2377  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
2378  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
2379  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
2380  }
2381 
2382  h->Fill(contents);
2383 
2384  return kTRUE;
2385 }
2386 
2387 // Definitions of class AliAnalysisTaskDmesonJets
2388 
2390 ClassImp(AliAnalysisTaskDmesonJets);
2392 
2396  fAnalysisEngines(),
2397  fEnabledAxis(0),
2399  fHistManager(),
2400  fApplyKinematicCuts(kTRUE),
2401  fNOutputTrees(0),
2402  fTrackEfficiency(0),
2403  fRejectISR(kFALSE),
2404  fJetAreaType(fastjet::active_area),
2405  fJetGhostArea(0.005),
2406  fMCContainer(0),
2407  fAodEvent(0),
2408  fFastJetWrapper(0)
2409 {
2410  SetMakeGeneralHistograms(kTRUE);
2411 }
2412 
2417  AliAnalysisTaskEmcalLight(name, kTRUE),
2418  fAnalysisEngines(),
2419  fEnabledAxis(k2ProngInvMass),
2420  fOutputType(kTreeOutput),
2421  fHistManager(name),
2422  fApplyKinematicCuts(kTRUE),
2423  fNOutputTrees(nOutputTrees),
2424  fTrackEfficiency(0),
2425  fRejectISR(kFALSE),
2426  fJetAreaType(fastjet::active_area),
2427  fJetGhostArea(0.005),
2428  fMCContainer(0),
2429  fAodEvent(0),
2430  fFastJetWrapper(0)
2431 {
2432  SetMakeGeneralHistograms(kTRUE);
2433  for (Int_t i = 0; i < nOutputTrees; i++){
2434  DefineOutput(2+i, TTree::Class());
2435  }
2436 }
2437 
2440 {
2441  if (fFastJetWrapper) delete fFastJetWrapper;
2442 }
2443 
2451 {
2452  AliRDHFCuts* analysiscuts = 0;
2453  TFile* filecuts = TFile::Open(cutfname);
2454  if (!filecuts || filecuts->IsZombie()) {
2455  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Input file not found: will use std cuts.");
2456  filecuts = 0;
2457  }
2458 
2459  if (filecuts) analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
2460 
2461  if (!analysiscuts) {
2462  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Could not find analysis cuts '%s' in '%s'.", cutsname.Data(), cutfname.Data());
2463  if (filecuts) {
2464  filecuts->ls();
2465  }
2466  }
2467  else {
2468  ::Info("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Cuts '%s' loaded from file '%s'", cutsname.Data(), cutfname.Data());
2469  }
2470 
2471  return analysiscuts;
2472 }
2473 
2486 {
2488  return AddAnalysisEngine(type, cutfname, cuttype, MCmode, jetDef, rhoName);
2489 }
2490 
2502 {
2503  AliRDHFCuts* cuts = 0;
2504 
2505  if (!cutfname.IsNull()) {
2506  TString cutsname;
2507 
2508  switch (type) {
2509  case kD0toKpi:
2510  case kD0toKpiLikeSign:
2511  cutsname = "D0toKpiCuts";
2512  break;
2513  case kDstartoKpipi:
2514  cutsname = "DStartoKpipiCuts";
2515  break;
2516  default:
2517  return 0;
2518  }
2519 
2520  if (!cuttype.IsNull()) {
2521  cutsname += TString::Format("_%s", cuttype.Data());
2522  }
2523 
2524  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
2525  if (cuts) cuts->PrintAll();
2526  }
2527 
2528  AnalysisEngine eng(type, MCmode, cuts);
2529 
2530  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
2531 
2532  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
2533  eng.AddJetDefinition(jetDef);
2534  it = fAnalysisEngines.insert(it, eng);
2535  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
2536  }
2537  else {
2538  AnalysisEngine* found_eng = &(*it);
2539  ::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());
2540  found_eng->AddJetDefinition(jetDef);
2541 
2542  if (cuts) {
2543  if (found_eng->fRDHFCuts != 0) ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
2544  found_eng->SetRDHFCuts(cuts);
2545  }
2546  }
2547 
2548  return &(*it);
2549 }
2550 
2551 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
2552 {
2553  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
2554  while (it != fAnalysisEngines.end() && (*it) != eng) it++;
2555  return it;
2556 }
2557 
2560 {
2561  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
2562 
2564 
2565  // Define histograms
2566  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
2567 
2568  TString hname;
2569  TString htitle;
2570  TH1* h = 0;
2571  Int_t treeSlot = 0;
2572 
2573  Int_t maxTracks = 6000;
2574  Double_t maxRho = 500;
2575  if (fForceBeamType == kpp) {
2576  maxRho = 50;
2577  maxTracks = 200;
2578  }
2579  else if (fForceBeamType == kpA) {
2580  maxRho = 200;
2581  maxTracks = 500;
2582  }
2583 
2584  hname = "fHistCharmPt";
2585  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2586  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2587 
2588  hname = "fHistCharmEta";
2589  htitle = hname + ";#eta_{charm};counts";
2590  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2591 
2592  hname = "fHistCharmPhi";
2593  htitle = hname + ";#phi_{charm};counts";
2594  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2595 
2596  hname = "fHistCharmPt_Eta05";
2597  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2598  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2599 
2600  hname = "fHistBottomPt";
2601  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2602  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2603 
2604  hname = "fHistBottomEta";
2605  htitle = hname + ";#eta_{bottom};counts";
2606  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2607 
2608  hname = "fHistBottomPhi";
2609  htitle = hname + ";#phi_{bottom};counts";
2610  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2611 
2612  hname = "fHistBottomPt_Eta05";
2613  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2614  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2615 
2616  hname = "fHistHighestPartonPt";
2617  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2618  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2619 
2620  hname = "fHistHighestPartonType";
2621  htitle = hname + ";type;counts";
2622  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2623 
2624  hname = "fHistNHeavyQuarks";
2625  htitle = hname + ";number of heavy-quarks;counts";
2626  fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
2627 
2628  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
2629  for (auto &param : fAnalysisEngines) {
2630  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
2631 
2632  param.fHistManager = &fHistManager;
2633 
2634  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", param.GetName());
2635  htitle = hname + ";#it{N}_{tracks};#it{N}_{D};events";
2636  h = fHistManager.CreateTH2(hname, htitle, 251, -0.5, 250.5, 21, -0.5, 20.5);
2637 
2638  hname = TString::Format("%s/fHistNTotAcceptedDmesons", param.GetName());
2639  htitle = hname + ";;#it{N}_{D}";
2640  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
2641 
2642  hname = TString::Format("%s/fHistNDmesons", param.GetName());
2643  htitle = hname + ";#it{N}_{D};events";
2644  h = fHistManager.CreateTH1(hname, htitle, 501, -0.5, 500.5);
2645 
2646  hname = TString::Format("%s/fHistNEvents", param.GetName());
2647  htitle = hname + ";Event status;counts";
2648  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
2649  h->GetXaxis()->SetBinLabel(1, "Accepted");
2650  h->GetXaxis()->SetBinLabel(2, "Rejected");
2651 
2652  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
2653  htitle = hname + ";Rejection reason;counts";
2654  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
2655 
2656  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
2657  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
2658  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2659 
2660  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
2661  htitle = hname + ";#it{#eta}_{D};counts";
2662  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2663 
2664  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
2665  htitle = hname + ";#it{#phi}_{D};counts";
2666  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2667 
2668  if (param.fMCMode != kMCTruth) {
2669  if (param.fCandidateType == kD0toKpi || param.fCandidateType == kD0toKpiLikeSign) {
2670  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
2671  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2672  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
2673  }
2674  else if (param.fCandidateType == kDstartoKpipi) {
2675  Double_t min = 0;
2676  Double_t max = 0;
2677 
2678  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
2679  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2680  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
2681  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
2682 
2683  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2684  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
2685  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2686  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
2687  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
2688  }
2689  }
2690 
2691  if (param.fMCMode == kMCTruth) {
2692  hname = TString::Format("%s/fHistPartonPt", param.GetName());
2693  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
2694  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2695 
2696  hname = TString::Format("%s/fHistPartonEta", param.GetName());
2697  htitle = hname + ";#eta_{parton};counts";
2698  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2699 
2700  hname = TString::Format("%s/fHistPartonPhi", param.GetName());
2701  htitle = hname + ";#phi_{parton};counts";
2702  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2703 
2704  hname = TString::Format("%s/fHistPartonType", param.GetName());
2705  htitle = hname + ";type;counts";
2706  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2707 
2708  hname = TString::Format("%s/fHistPrompt", param.GetName());
2709  htitle = hname + ";Type;counts";
2710  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
2711  h->GetXaxis()->SetBinLabel(1, "Unknown");
2712  h->GetXaxis()->SetBinLabel(2, "Prompt");
2713  h->GetXaxis()->SetBinLabel(3, "Non-Prompt");
2714 
2715  hname = TString::Format("%s/fHistAncestor", param.GetName());
2716  htitle = hname + ";Ancestor;counts";
2717  h = fHistManager.CreateTH1(hname, htitle, 4, 0, 4);
2718  h->GetXaxis()->SetBinLabel(1, "Unknown");
2719  h->GetXaxis()->SetBinLabel(2, "Charm");
2720  h->GetXaxis()->SetBinLabel(3, "Bottom");
2721  h->GetXaxis()->SetBinLabel(4, "Proton");
2722  }
2723 
2724  for (auto& jetDef : param.fJetDefinitions) {
2725  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef.GetName());
2726 
2727  if (param.fMCMode == kMCTruth) {
2728  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", param.GetName(), jetDef.GetName());
2729  htitle = hname + ";#it{N}_{constituents};#it{N}_{D};counts";
2730  h = fHistManager.CreateTH2(hname, htitle, 51, -0.5, 50.5, 10, 0.5, 10.5);
2731  }
2732 
2733  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef.GetName());
2734  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2735  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2736  SetRejectionReasonLabels(h->GetXaxis());
2737 
2738  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef.GetName());
2739  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2740  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2741  SetRejectionReasonLabels(h->GetXaxis());
2742 
2743  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef.GetName());
2744  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
2745  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2746  SetRejectionReasonLabels(h->GetXaxis());
2747 
2748  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef.GetName());
2749  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
2750  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
2751 
2752  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef.GetName());
2753  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
2754  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2755 
2756  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef.GetName());
2757  htitle = hname + ";#it{#eta}_{jet};counts";
2758  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2759 
2760  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef.GetName());
2761  htitle = hname + ";#it{#phi}_{jet};counts";
2762  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2763 
2764  if (!jetDef.fRhoName.IsNull()) {
2765  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", param.GetName(), jetDef.GetName());
2766  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});#rho (GeV/#it{c} #times rad^{-1});counts";
2767  fHistManager.CreateTH2(hname, htitle, 300, 0, 150, 1000, 0, maxRho);
2768 
2769  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", param.GetName(), jetDef.GetName());
2770  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});#rho (GeV/#it{c} #times rad^{-1});counts";
2771  fHistManager.CreateTH2(hname, htitle, 300, 0, 150, 1000, 0, maxRho);
2772 
2773  hname = TString::Format("%s/%s/fHistRhoVsCent", param.GetName(), jetDef.GetName());
2774  htitle = hname + ";Centrality (%);#rho (GeV/#it{c} #times rad^{-1});counts";
2775  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 1000, 0, maxRho);
2776 
2777  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", param.GetName(), jetDef.GetName());
2778  htitle = hname + ";Centrality (%);#it{p}_{T,jet} (GeV/#it{c});counts";
2779  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
2780 
2781  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", param.GetName(), jetDef.GetName());
2782  htitle = hname + ";Centrality (%);#it{p}_{T,D} (GeV/#it{c});counts";
2783  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
2784 
2785  hname = TString::Format("%s/%s/fHistRhoVsNTracks", param.GetName(), jetDef.GetName());
2786  htitle = hname + ";no. of tracks;#rho (GeV/#it{c} #times rad^{-1});counts";
2787  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 1000, 0, maxRho);
2788 
2789  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", param.GetName(), jetDef.GetName());
2790  htitle = hname + ";no. of tracks;#it{p}_{T,jet} (GeV/#it{c});counts";
2791  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
2792 
2793  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", param.GetName(), jetDef.GetName());
2794  htitle = hname + ";no. of tracks;#it{p}_{T,D} (GeV/#it{c});counts";
2795  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
2796  }
2797  }
2798  switch (fOutputType) {
2799  case kTreeOutput:
2800  param.BuildTree(GetName());
2801  if (treeSlot < fNOutputTrees) {
2802  param.AssignDataSlot(treeSlot+2);
2803  treeSlot++;
2805  }
2806  else {
2807  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
2808  }
2809  break;
2810  case kTHnOutput:
2811  param.BuildHnSparse(fEnabledAxis);
2812  break;
2813  case kNoOutput:
2814  break;
2815  }
2816  }
2817 
2819 
2820  PostData(1, fOutput);
2821 }
2822 
2826 {
2828 
2829  // Load the event
2830  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
2831 
2832  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
2833 
2834  fFastJetWrapper->SetAreaType((fastjet::AreaType)fJetAreaType);
2836 
2837  if (!fAodEvent) {
2838  AliError(Form("This task need an AOD event (Task '%s'). Expect troubles...", GetName()));
2839  //return;
2840  }
2841 
2842  TRandom* rnd = 0;
2843  if (fTrackEfficiency > 0 && fTrackEfficiency < 1) rnd = new TRandom3(0);
2844 
2845  for (auto cont_it : fParticleCollArray) {
2846  AliHFAODMCParticleContainer* part_cont = dynamic_cast<AliHFAODMCParticleContainer*>(cont_it.second);
2847  if (part_cont) fMCContainer = part_cont;
2848  }
2849 
2850  for (auto &params : fAnalysisEngines) {
2851 
2852  params.fAodEvent = fAodEvent;
2853  params.fFastJetWrapper = fFastJetWrapper;
2854  params.fTrackEfficiency = fTrackEfficiency;
2855  params.fRejectISR = fRejectISR;
2856  params.fRandomGen = rnd;
2857 
2858  for (auto &jetdef: params.fJetDefinitions) {
2859  if (!jetdef.fRhoName.IsNull()) {
2860  jetdef.fRho = dynamic_cast<AliRhoParameter*>(fInputEvent->FindListObject(jetdef.fRhoName));
2861  if (!jetdef.fRho) {
2862  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2863  "%s: Could not find rho object '%s' for engine '%s'",
2864  GetName(), jetdef.fRhoName.Data(), params.GetName());
2865  }
2866  }
2867  }
2868 
2869  if (!params.fRDHFCuts) {
2870  if (params.fMCMode == kMCTruth) {
2871  ::Warning("AliAnalysisTaskDmesonJets::ExecOnce",
2872  "%s: RDHF cuts not provided for engine '%s'.",
2873  GetName(), params.GetName());
2874  }
2875  else {
2876  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2877  "%s: RDHF cuts not provided. Engine '%s' disabled.",
2878  GetName(), params.GetName());
2879  params.fInhibit = kTRUE;
2880  }
2881  }
2882 
2883  params.fMCContainer = fMCContainer;
2884 
2885  for (auto cont_it : fParticleCollArray) {
2886  AliTrackContainer* track_cont = dynamic_cast<AliTrackContainer*>(cont_it.second);
2887  if (track_cont) params.fTrackContainers.push_back(track_cont);
2888  }
2889 
2890  for (auto cont_it : fClusterCollArray) params.fClusterContainers.push_back(cont_it.second);
2891 
2892  if (fAodEvent) params.Init(fGeom, fAodEvent->GetRunNumber());
2893 
2894  if (params.fMCMode != kMCTruth && fAodEvent) {
2895  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
2896 
2897  if (params.fCandidateArray) {
2898  TString className;
2899  if (params.fCandidateType == kD0toKpi || params.fCandidateType == kD0toKpiLikeSign) {
2900  className = "AliAODRecoDecayHF2Prong";
2901  }
2902  else if (params.fCandidateType == kDstartoKpipi) {
2903  className = "AliAODRecoCascadeHF";
2904  }
2905  if (!params.fCandidateArray->GetClass()->InheritsFrom(className)) {
2906  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2907  "%s: Objects of type %s in %s are not inherited from %s! Task will be disabled!",
2908  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName(), className.Data());
2909  params.fCandidateArray = 0;
2910  params.fInhibit = kTRUE;
2911  }
2912  }
2913  else {
2914  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2915  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
2916  params.fBranchName.Data(), params.GetName());
2917  params.fInhibit = kTRUE;
2918  }
2919  }
2920 
2921  if (params.fMCMode != kNoMC) {
2922  if (!params.fMCContainer) {
2923  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2924  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
2925  params.GetName());
2926  params.fInhibit = kTRUE;
2927  }
2928  }
2929 
2930  if (params.fMCMode != kMCTruth) {
2931  if (params.fTrackContainers.size() == 0 && params.fClusterContainers.size() == 0) {
2932  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2933  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
2934  params.GetName());
2935  params.fInhibit = kTRUE;
2936  }
2937  }
2938  }
2939 }
2940 
2945 {
2946  TString hname;
2947 
2948  // Fix for temporary bug in ESDfilter
2949  // The AODs with null vertex pointer didn't pass the PhysSel
2950  // Now adding an entry in the histogram so as to check that this is actually cutting anything out
2951  if (fAodEvent && (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001)) {
2952  for (auto &eng : fAnalysisEngines) {
2953  if (eng.fInhibit) continue;
2954  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2955  fHistManager.FillTH1(hname, "ESDfilterBug");
2956  }
2957  return kFALSE;
2958  }
2959 
2960  if (fAodEvent) {
2961  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
2962  if (matchingAODdeltaAODlevel <= 0) {
2963  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
2964  for (auto &eng : fAnalysisEngines) {
2965  if (eng.fInhibit) continue;
2966  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2967  fHistManager.FillTH1(hname, "MismatchDeltaAOD");
2968  }
2969  return kFALSE;
2970  }
2971  }
2972 
2973  for (auto &eng : fAnalysisEngines) {
2974  eng.fDmesonJets.clear();
2975  if (eng.fInhibit) continue;
2976 
2977  eng.fCent = fCent;
2978 
2979  //Event selection
2980  hname = TString::Format("%s/fHistNEvents", eng.GetName());
2981  if (fAodEvent) {
2982  Bool_t iseventselected = kTRUE;
2983  if (eng.fRDHFCuts) iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
2984  if (!iseventselected) {
2985  fHistManager.FillTH1(hname, "Rejected");
2986  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2987  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
2988  TString label;
2989  do {
2990  label = GetHFEventRejectionReasonLabel(bitmap);
2991  if (label.IsNull()) break;
2992  fHistManager.FillTH1(hname, label);
2993  } while (true);
2994 
2995  continue;
2996  }
2997  }
2998 
2999  fHistManager.FillTH1(hname, "Accepted");
3000 
3001  AliDebug(2, "Event selected");
3002 
3003  eng.RunAnalysis();
3004  }
3005  return kTRUE;
3006 }
3007 
3012 {
3013  for (auto &param : fAnalysisEngines) {
3014  if (param.fInhibit) continue;
3015 
3016  if (fOutputType == kTreeOutput) {
3017  param.FillTree(fApplyKinematicCuts);
3019  }
3020  else if (fOutputType == kTHnOutput) {
3021  param.FillHnSparse(fApplyKinematicCuts);
3022  }
3023  }
3025  return kTRUE;
3026 }
3027 
3030 {
3031  auto itcont = fMCContainer->all_momentum();
3032  Int_t nHQ = 0;
3033  Double_t highestPartonPt = 0;
3034  Int_t absPdgHighParton = 0;
3035  for (auto part : itcont) {
3036  Int_t absPdgCode = TMath::Abs(part.second->GetPdgCode());
3037 
3038  // Skip all particles that are not either quarks or gluons
3039  if (absPdgCode > 9 && absPdgCode != 21) continue;
3040 
3041  // Look for highest momentum parton
3042  if (highestPartonPt < part.first.Pt()) {
3043  highestPartonPt = part.first.Pt();
3044  absPdgHighParton = absPdgCode;
3045  }
3046  /*
3047  // Look for the mother PDG code
3048  Int_t motherIndex = part.second->GetMother();
3049  AliAODMCParticle *mother = 0;
3050  Int_t motherPdg = 0;
3051  Double_t motherPt = 0;
3052  if (motherIndex >= 0) {
3053  mother = fMCContainer->GetMCParticle(motherIndex);
3054  if (motherIndex) {
3055  motherPdg = TMath::Abs(mother->GetPdgCode());
3056  motherPt = mother->Pt();
3057  }
3058  }
3059  */
3060  if (absPdgCode != 4 && absPdgCode != 5) continue;
3061  Bool_t notLastInPartonShower = kFALSE;
3062  for (Int_t idaugh = 0; idaugh < 2; idaugh++){
3063  Int_t daughterIndex = part.second->GetDaughter(idaugh);
3064  if (daughterIndex < 0) {
3065  AliDebug(10, Form("Could not find daughter of heavy quark (pdg=%d, pt=%.3f)!", absPdgCode, part.first.Pt()));
3066  continue;
3067  }
3068  AliAODMCParticle *daughter = fMCContainer->GetMCParticle(daughterIndex);
3069  if (!daughter) {
3070  AliDebug(10, Form("Could not find daughter %d of heavy quark (pdg=%d, pt=%.3f)!", daughterIndex, absPdgCode, part.first.Pt()));
3071  continue;
3072  }
3073  Int_t daughterAbsPdgCode = TMath::Abs(daughter->GetPdgCode());
3074  if (daughterAbsPdgCode <= 9 || daughterAbsPdgCode == 21) notLastInPartonShower = kTRUE; // this parton is not the last parton in the shower
3075  AliDebug(10, Form("Found daughter with PDG=%d, pt=%.3f", daughterAbsPdgCode, daughter->Pt()));
3076  }
3077  if (notLastInPartonShower) continue;
3078 
3079  if (absPdgCode == 4) {
3080  fHistManager.FillTH1("fHistCharmPt", part.first.Pt());
3081  fHistManager.FillTH1("fHistCharmEta", part.first.Eta());
3082  fHistManager.FillTH1("fHistCharmPhi", part.first.Phi_0_2pi());
3083  if (TMath::Abs(part.first.Eta()) < 0.5) fHistManager.FillTH1("fHistCharmPt_Eta05", part.first.Pt());
3084  }
3085  else if (absPdgCode == 5) {
3086  fHistManager.FillTH1("fHistBottomPt", part.first.Pt());
3087  fHistManager.FillTH1("fHistBottomEta", part.first.Eta());
3088  fHistManager.FillTH1("fHistBottomPhi", part.first.Phi_0_2pi());
3089  if (TMath::Abs(part.first.Eta()) < 0.5) fHistManager.FillTH1("fHistBottomPt_Eta05", part.first.Pt());
3090  }
3091  nHQ++;
3092  }
3093  fHistManager.FillTH1("fHistNHeavyQuarks", nHQ);
3094  fHistManager.FillTH1("fHistHighestPartonPt",highestPartonPt);
3095  Int_t partonType = 0;
3096  if (absPdgHighParton == 9 || absPdgHighParton == 21) {
3097  partonType = 7;
3098  }
3099  else {
3100  partonType = absPdgHighParton;
3101  }
3102  fHistManager.FillTH1("fHistHighestPartonType",partonType);
3103 }
3104 
3113 {
3114  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
3115 
3116  Double_t mass = part->Mass();
3117 
3118  // To make sure that the PDG mass value is not at the edge of a bin
3119  if (nbins % 2 == 0) {
3120  minMass = mass - range / 2 - range / nbins / 2;
3121  maxMass = mass + range / 2 - range / nbins / 2;
3122  }
3123  else {
3124  minMass = mass - range / 2;
3125  maxMass = mass + range / 2;
3126  }
3127 }
3128 
3135 {
3136  static TString label;
3137  label = "";
3138 
3139  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
3140  label = "NotSelTrigger";
3141  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
3142  return label.Data();
3143  }
3144  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
3145  label = "NoVertex";
3146  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
3147  return label.Data();
3148  }
3149  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
3150  label = "TooFewVtxContrib";
3151  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
3152  return label.Data();
3153  }
3154  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
3155  label = "ZVtxOutFid";
3156  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
3157  return label.Data();
3158  }
3159  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
3160  label = "Pileup";
3161  bitmap &= ~BIT(AliRDHFCuts::kPileup);
3162  return label.Data();
3163  }
3164  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
3165  label = "OutsideCentrality";
3166  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
3167  return label.Data();
3168  }
3169  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
3170  label = "PhysicsSelection";
3171  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
3172  return label.Data();
3173  }
3174  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
3175  label = "BadSPDVertex";
3176  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
3177  return label.Data();
3178  }
3179  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
3180  label = "ZVtxSPDOutFid";
3181  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
3182  return label.Data();
3183  }
3184  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
3185  label = "CentralityFlattening";
3186  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
3187  return label.Data();
3188  }
3189  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
3190  label = "BadTrackV0Correl";
3191  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
3192  return label.Data();
3193  }
3194 
3195  return label.Data();
3196 }
3197 
3204 {
3205  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
3206  PostData(eng.GetDataSlotNumber(), eng.GetTree());
3207  return eng.GetDataSlotNumber();
3208  }
3209  else {
3210  return -1;
3211  }
3212 }
3213 
3223 {
3224  // Get the pointer to the existing analysis manager via the static access method
3225  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
3226  if (!mgr) {
3227  ::Error("AddTaskDmesonJets", "No analysis manager to connect to.");
3228  return NULL;
3229  }
3230 
3231  // Check the analysis type using the event handlers connected to the analysis manager
3232  AliVEventHandler* handler = mgr->GetInputEventHandler();
3233  if (!handler) {
3234  ::Error("AddTaskEmcalJetSpectraQA", "This task requires an input event handler");
3235  return NULL;
3236  }
3237 
3238  EDataType_t dataType = kUnknownDataType;
3239 
3240  if (handler->InheritsFrom("AliESDInputHandler")) {
3241  dataType = kESD;
3242  }
3243  else if (handler->InheritsFrom("AliAODInputHandler")) {
3244  dataType = kAOD;
3245  }
3246 
3247  // Init the task and do settings
3248  if (ntracks == "usedefault") {
3249  if (dataType == kESD) {
3250  ntracks = "Tracks";
3251  }
3252  else if (dataType == kAOD) {
3253  ntracks = "tracks";
3254  }
3255  else {
3256  ntracks = "";
3257  }
3258  }
3259 
3260  if (nclusters == "usedefault") {
3261  if (dataType == kESD) {
3262  nclusters = "CaloClusters";
3263  }
3264  else if (dataType == kAOD) {
3265  nclusters = "caloClusters";
3266  }
3267  else {
3268  nclusters = "";
3269  }
3270  }
3271 
3272  if (nMCpart == "usedefault") {
3273  nMCpart = "mcparticles"; // Always needs AliAODMCParticle objects
3274  }
3275 
3276  TString name("AliAnalysisTaskDmesonJets");
3277  if (strcmp(suffix, "") != 0) {
3278  name += TString::Format("_%s", suffix.Data());
3279  }
3280 
3281  AliAnalysisTaskDmesonJets* jetTask = new AliAnalysisTaskDmesonJets(name, nMaxTrees);
3282 
3283  if (!ntracks.IsNull()) {
3284  AliHFTrackContainer* trackCont = new AliHFTrackContainer(ntracks);
3285  jetTask->AdoptParticleContainer(trackCont);
3286  }
3287 
3288  if (!nMCpart.IsNull()) {
3289  AliMCParticleContainer* partCont = new AliHFAODMCParticleContainer(nMCpart);
3290  partCont->SetEtaLimits(-1.5, 1.5);
3291  partCont->SetPtLimits(0, 1000);
3292  jetTask->AdoptParticleContainer(partCont);
3293  }
3294 
3295  jetTask->AddClusterContainer(nclusters.Data());
3296 
3297  // Final settings, pass to manager and set the containers
3298  mgr->AddTask(jetTask);
3299 
3300  // Create containers for input/output
3301  AliAnalysisDataContainer* cinput1 = mgr->GetCommonInputContainer();
3302  TString contname1(name);
3303  contname1 += "_histos";
3304  AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(contname1.Data(),
3305  TList::Class(), AliAnalysisManager::kOutputContainer,
3306  Form("%s", AliAnalysisManager::GetCommonFileName()));
3307 
3308  mgr->ConnectInput(jetTask, 0, cinput1);
3309  mgr->ConnectOutput(jetTask, 1, coutput1);
3310 
3311  for (Int_t i = 0; i < nMaxTrees; i++) {
3312  TString contname = TString::Format("%s_tree_%d", name.Data(), i);
3313  AliAnalysisDataContainer *coutput = mgr->CreateContainer(contname.Data(),
3314  TTree::Class(),AliAnalysisManager::kOutputContainer,
3315  Form("%s", AliAnalysisManager::GetCommonFileName()));
3316  mgr->ConnectOutput(jetTask, 2+i, coutput);
3317  }
3318  return jetTask;
3319 }
3320 
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.
AnalysisEngine * AddAnalysisEngine(ECandidateType_t type, TString cutfname, TString cuttype, EMCMode_t bkgMode, EJetType_t jettype, Double_t jetradius, TString rhoName="")
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.
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)