AliPhysics  068200c (068200c)
 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 //C++
17 #include <sstream>
18 
19 // Root
20 #include <TClonesArray.h>
21 #include <TDatabasePDG.h>
22 #include <TParticlePDG.h>
23 #include <TVector3.h>
24 #include <THnSparse.h>
25 #include <TParticle.h>
26 #include <TMath.h>
27 #include <THashList.h>
28 #include <TFile.h>
29 #include <TRandom3.h>
30 
31 // Aliroot general
32 #include "AliLog.h"
33 #include "AliEMCALGeometry.h"
34 #include "AliAnalysisManager.h"
35 #include "AliVEventHandler.h"
36 
37 // Aliroot HF
38 #include "AliAODRecoCascadeHF.h"
40 #include "AliRDHFCutsD0toKpi.h"
43 #include "AliHFTrackContainer.h"
44 #include "AliAnalysisVertexingHF.h"
45 
46 // Aliroot EMCal jet framework
47 #include "AliEmcalJetTask.h"
48 #include "AliEmcalJet.h"
49 #include "AliJetContainer.h"
50 #include "AliParticleContainer.h"
51 #include "AliEmcalParticle.h"
52 #include "AliFJWrapper.h"
53 #include "AliRhoParameter.h"
54 
56 
57 AliAnalysisTaskDmesonJets::AliEventNotFound::AliEventNotFound(const std::string& class_name, const std::string& method_name) :
58  std::exception(),
59  fClassName(class_name),
60  fAccessMethodName(method_name)
61 {
62  std::stringstream what_str;
63  what_str << "ALICE event not found in class '" << fClassName << "' using method '" << method_name << "'.";
64  fWhat = what_str.str();
65 }
66 
67 #if !(defined(__CINT__) || defined(__MAKECINT__))
69 {
70  return fWhat.c_str();
71 }
72 #endif
73 
74 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfo
75 
79 
87 {
88  dphi = TVector2::Phi_mpi_pi(fMomentum.Phi() - jet.Phi());;
89  deta = fMomentum.Eta() - jet.Eta();
90  return TMath::Sqrt(dphi*dphi + deta*deta);
91 }
92 
98 {
99  Double_t deta = 0;
100  Double_t dphi = 0;
101  return GetDistance(jet, deta, dphi);
102 }
103 
104 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonJetInfo
105 
109 
112  fDmesonParticle(0),
113  fD(),
114  fSoftPionPt(0),
115  fInvMass2Prong(0),
116  fJets(),
117  fMCLabel(-1),
118  fReconstructed(kFALSE),
119  fParton(0),
120  fPartonType(0),
121  fAncestor(0),
122  fD0D0bar(kFALSE),
123  fSelectionType(0),
124  fEvent(nullptr)
125 {
126 }
127 
132  fDmesonParticle(source.fDmesonParticle),
133  fD(source.fD),
134  fSoftPionPt(source.fSoftPionPt),
135  fInvMass2Prong(source.fInvMass2Prong),
136  fJets(source.fJets),
137  fMCLabel(source.fMCLabel),
138  fReconstructed(source.fReconstructed),
139  fParton(source.fParton),
140  fPartonType(source.fPartonType),
141  fAncestor(source.fAncestor),
142  fD0D0bar(source.fD0D0bar),
143  fSelectionType(source.fSelectionType),
144  fEvent(source.fEvent)
145 {
146 }
147 
152 {
153  new (this) AliDmesonJetInfo(source);
154  return *this;
155 }
156 
159 {
160  fD.SetPtEtaPhiE(0,0,0,0);
161  fSoftPionPt = 0;
162  fInvMass2Prong = 0;
163  fDmesonParticle = 0;
164  fMCLabel = -1;
165  fReconstructed = kFALSE;
166  fParton = 0;
167  fPartonType = 0;
168  fAncestor = 0;
169  fD0D0bar = kFALSE;
170  for (auto &jet : fJets) {
171  jet.second.fMomentum.SetPtEtaPhiE(0,0,0,0);
172  jet.second.fNConstituents = 0;
173  jet.second.fNEF = 0;
174  jet.second.fMaxChargedPt = 0;
175  jet.second.fMaxNeutralPt = 0;
176  }
177 }
178 
181 {
182  Printf("Printing D Meson Jet object.");
183  Printf("D Meson: pT = %.3f, eta = %.3f, phi = %.3f, inv. mass = %.3f", fD.Pt(), fD.Eta(), fD.Phi_0_2pi(), fD.M());
184  Printf("Soft pion pT: %.3f. 2-Prong Invariant mass = %.3f", fSoftPionPt, fInvMass2Prong);
185  for (auto &jet : fJets) {
186  Printf("Jet %s: pT = %.3f, eta = %.3f, phi = %.3f", jet.first.c_str(), jet.second.Pt(), jet.second.Eta(), jet.second.Phi_0_2pi());
187  Printf("Jet N Consituents = %d", jet.second.fNConstituents);
188  }
189 }
190 
195 {
196  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
197  if (it == fJets.end()) return 0;
198 
199  Double_t z = 0;
200 
201  if ((*it).second.Pt() > 0) {
202  TVector3 dvect = fD.Vect();
203  TVector3 jvect = (*it).second.fMomentum.Vect();
204 
205  Double_t jetMom = jvect * jvect;
206 
207  if (jetMom < 1e-6) {
208  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetZ", "Zero jet momentum!");
209  z = 0.999;
210  }
211  else {
212  z = (dvect * jvect) / jetMom;
213  }
214 
215  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
216  }
217 
218  return z;
219 }
220 
225 {
226  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
227  if (it == fJets.end()) return 0;
228 
229  Double_t z = 0;
230 
231  if ((*it).second.Pt() > 0) {
232  TVector3 dvect = fD.Vect();
233  TVector3 jvect = (*it).second.fMomentum.Vect();
234  // If the corr pt is < 0, assign 0.
235  Double_t corrpt = (*it).second.fCorrPt > 0 ? (*it).second.fCorrPt : 0.;
236  jvect.SetPerp(corrpt);
237 
238  Double_t jetMom = jvect * jvect;
239 
240  if (jetMom < 1e-6) {
241  z = 1.0;
242  }
243  else {
244  z = (dvect * jvect) / jetMom;
245  }
246  }
247 
248  return z;
249 }
250 
258 {
259  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
260  if (it == fJets.end()) return 0;
261 
262  return GetDistance((*it).second, deta, dphi);
263 }
264 
270 {
271  Double_t deta = 0;
272  Double_t dphi = 0;
273  return GetDistance(n, deta, dphi);
274 }
275 
283 {
284  dphi = TVector2::Phi_mpi_pi(fD.Phi() - jet.Phi());;
285  deta = fD.Eta() - jet.Eta();
286  return TMath::Sqrt(dphi*dphi + deta*deta);
287 }
288 
294 {
295  Double_t deta = 0;
296  Double_t dphi = 0;
297  return GetDistance(jet, deta, dphi);
298 }
299 
305 {
306  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
307  if (it == fJets.end()) {
308  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
309  n.c_str());
310  return 0;
311  }
312  return &((*it).second);
313 }
314 
320 {
321  std::map<std::string, AliJetInfo>::iterator it = fJets.find(n);
322  if (it == fJets.end()) {
323  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
324  n.c_str());
325  return 0;
326  }
327  return &((*it).second);
328 }
329 
330 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoSummary
331 
335 
341  fPt(0),
342  fEta(0),
343  fPhi(0),
344  fR(0),
345  fZ(0),
346  fN(0)
347 {
348  Set(source, n);
349 }
350 
353 {
354  fPt = 0;
355  fEta = 0;
356  fPhi = 0;
357  fR = 0;
358  fZ = 0;
359  fN = 0;
360 }
361 
367 {
368  std::map<std::string, AliJetInfo>::const_iterator it = source.fJets.find(n);
369  if (it == source.fJets.end()) return;
370 
371  Set((*it).second);
372 
373  fR = source.GetDistance(n);
374  fZ = source.GetZ(n);
375 }
376 
381 {
382  fPt = source.Pt();
383  fEta = source.Eta();
384  fPhi = source.Phi_0_2pi();
385  fN = source.GetNConstituents();
386  fR = 0;
387  fZ = 0;
388 }
389 
390 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoPbPbSummary
391 
395 
402  fCorrPt(0),
403  fCorrZ(0),
404  fArea(0)
405 {
406  Set(source, n);
407 }
408 
411 {
413  fCorrPt = 0;
414  fCorrZ = 0;
415  fArea = 0;
416 }
417 
422 {
423  AliJetInfoSummary::Set(source);
424  fArea = source.fArea;
425  fCorrPt = source.fCorrPt;
426 }
427 
433 {
434  AliJetInfoSummary::Set(source, n);
435  fCorrZ = source.GetCorrZ(n);
436 }
437 
438 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonInfoSummary
439 
443 
448  fPt(0),
449  fEta(0),
450  fPhi(0)
451 {
452  Set(source);
453 }
454 
459 {
460  fPt = source.fD.Pt();
461  fEta = source.fD.Eta();
462  fPhi = source.fD.Phi_0_2pi();
463 }
464 
467 {
468  fPt = 0;
469  fEta = 0;
470  fPhi = 0;
471 }
472 
473 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary
474 
478 
483  AliDmesonInfoSummary(source),
484  fPartonType(0),
485  fPartonPt(0),
486  fAncestorPDG(0)
487 {
488  Set(source);
489 }
490 
495 {
497 
498  fPartonType = source.fPartonType;
499 
500  if (source.fParton) {
501  fPartonPt = source.fParton->Pt();
502  }
503  else {
504  fPartonPt = 0.;
505  }
506 
507  if (source.fAncestor) {
508  fAncestorPDG = (UShort_t)((UInt_t)(TMath::Abs(source.fAncestor->GetPdgCode())));
509  }
510 }
511 
514 {
516  fPartonType = 0,
517  fPartonPt = 0.;
518  fAncestorPDG = 0;
519 }
520 
521 // Definitions of class AliAnalysisTaskDmesonJets::AliD0InfoSummary
522 
526 
531  AliDmesonInfoSummary(source),
532  fInvMass(source.fD.M()),
533  fSelectionType(0)
534 {
535 }
536 
541 {
542  fInvMass = source.fD.M();
543  fSelectionType = source.GetSelectionTypeSummary();
545 }
546 
549 {
551  fSelectionType = 0;
552  fInvMass = 0;
553 }
554 
555 // Definitions of class AliAnalysisTaskDmesonJets::AliD0ExtendedInfoSummary
556 
560 
565  AliD0InfoSummary(source),
566  fDCA(0),
567  fCosThetaStar(0),
568  fd0K(0),
569  fd0Pi(0),
570  fd0d0(0),
571  fCosPointing(0),
572  fMaxNormd0(0)
573 {
574  Set(source);
575 }
576 
581 {
582  AliD0InfoSummary::Set(source);
583 
584  AliAODRecoDecayHF2Prong* recoDecay = dynamic_cast<AliAODRecoDecayHF2Prong*>(source.fDmesonParticle);
585  if (recoDecay) {
586  fDCA = recoDecay->GetDCA();
587  if (source.fSelectionType == 1) { // D0
588  fCosThetaStar = recoDecay->CosThetaStarD0();
589  fPtK = recoDecay->PtProng(0);
590  fPtPi = recoDecay->PtProng(1);
591  fd0K = recoDecay->Getd0Prong(0);
592  fd0Pi = recoDecay->Getd0Prong(1);
593  }
594  else { //D0bar
595  fCosThetaStar = recoDecay->CosThetaStarD0bar();
596  fPtK = recoDecay->PtProng(1);
597  fPtPi = recoDecay->PtProng(0);
598  fd0K = recoDecay->Getd0Prong(1);
599  fd0Pi = recoDecay->Getd0Prong(0);
600  }
601 
602  fMaxNormd0 = 0.;
603  // Based on Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent* aod)
604  // Line 480 and following
605  if (source.fEvent) {
606  for (Int_t ipr=0; ipr < 2; ipr++) {
607  Double_t diffIP = 0., errdiffIP = 0.;
608  recoDecay->Getd0MeasMinusExpProng(ipr, source.fEvent->GetMagneticField(), diffIP, errdiffIP);
609  Double_t normdd0 = 0.;
610  if (errdiffIP > 0.) {
611  normdd0 = diffIP / errdiffIP;
612  }
613  else {
614  if (diffIP == 0) {
615  normdd0 = 0;
616  }
617  else {
618  normdd0 = diffIP > 0 ? 9999. : -9999.;
619  }
620  }
621  if (TMath::Abs(normdd0) > TMath::Abs(fMaxNormd0)) {
622  fMaxNormd0 = normdd0;
623  }
624  }
625  }
626  else {
627  throw AliAnalysisTaskDmesonJets::AliEventNotFound("AliAnalysisTaskDmesonJets::AliDmesonJetInfo", "fEvent");
628  }
629 
630  fd0d0 = recoDecay->Prodd0d0();
631  fCosPointing = recoDecay->CosPointingAngle();
632  }
633 }
634 
637 {
639  fDCA = 0;
640  fCosThetaStar = 0;
641  fd0K = 0;
642  fd0Pi = 0;
643  fd0d0 = 0;
644  fCosPointing = 0;
645  fMaxNormd0 = 0;
646 }
647 
648 // Definitions of class AliAnalysisTaskDmesonJets::AliDStarInfoSummary
649 
653 
658  AliDmesonInfoSummary(source),
659  f2ProngInvMass(source.fInvMass2Prong),
660  fDeltaInvMass(source.fD.M() - source.fInvMass2Prong)
661 {
662 }
663 
668 {
669  f2ProngInvMass = source.fInvMass2Prong;
670  fDeltaInvMass = source.fD.M() - source.fInvMass2Prong;
672 }
673 
676 {
678 
679  f2ProngInvMass = 0;
680  fDeltaInvMass = 0;
681 }
682 
683 // Definitions of class AliAnalysisTaskDmesonJets::AliJetDefinition
684 
688 
691  TObject(),
692  fJetType(AliJetContainer::kChargedJet),
693  fRadius(0),
694  fJetAlgo(AliJetContainer::antikt_algorithm),
695  fRecoScheme(AliJetContainer::pt_scheme),
696  fMinJetPt(0.),
697  fMaxJetPt(500.),
698  fMinJetPhi(0.),
699  fMaxJetPhi(0.),
700  fMinJetEta(-1.),
701  fMaxJetEta(1.),
702  fMinChargedPt(0.),
703  fMaxChargedPt(0.),
704  fMinNeutralPt(0.),
705  fMaxNeutralPt(0.),
706  fRhoName(),
707  fRho(0)
708 {
709 }
710 
718  TObject(),
719  fJetType(type),
720  fRadius(r),
721  fJetAlgo(algo),
722  fRecoScheme(reco),
723  fMinJetPt(0.),
724  fMaxJetPt(500.),
725  fMinJetPhi(0.),
726  fMaxJetPhi(0.),
727  fMinJetEta(-1.),
728  fMaxJetEta(1.),
729  fMinChargedPt(0.),
730  fMaxChargedPt(0.),
731  fMinNeutralPt(0.),
732  fMaxNeutralPt(0.),
733  fRhoName(),
734  fRho(0)
735 {
736 }
737 
745  TObject(),
746  fJetType(type),
747  fRadius(r),
748  fJetAlgo(algo),
749  fRecoScheme(reco),
750  fMinJetPt(0.),
751  fMaxJetPt(0.),
752  fMinJetPhi(0.),
753  fMaxJetPhi(0.),
754  fMinJetEta(0.),
755  fMaxJetEta(0.),
756  fMinChargedPt(0.),
757  fMaxChargedPt(0.),
758  fMinNeutralPt(0.),
759  fMaxNeutralPt(0.),
760  fRhoName(rhoName),
761  fRho(0)
762 {
763 }
764 
769  TObject(),
770  fJetType(source.fJetType),
771  fRadius(source.fRadius),
772  fJetAlgo(source.fJetAlgo),
773  fRecoScheme(source.fRecoScheme),
774  fMinJetPt(source.fMinJetPt),
775  fMaxJetPt(source.fMaxJetPt),
776  fMinJetPhi(source.fMinJetPhi),
777  fMaxJetPhi(source.fMaxJetPhi),
778  fMinJetEta(source.fMinJetEta),
779  fMaxJetEta(source.fMaxJetEta),
780  fMinChargedPt(source.fMinChargedPt),
781  fMaxChargedPt(source.fMaxChargedPt),
782  fMinNeutralPt(source.fMinNeutralPt),
783  fMaxNeutralPt(source.fMaxNeutralPt),
784  fRhoName(source.fRhoName),
785  fRho(0)
786 {
787 }
788 
793 {
794  new (this) AliHFJetDefinition(source);
795  return *this;
796 }
797 
800 {
801  static TString name;
802 
803  name = AliJetContainer::GenerateJetName(fJetType, fJetAlgo, fRecoScheme, fRadius, 0, 0, "Jet");
804 
805  return name.Data();
806 }
807 
813 {
814  if (fMinJetEta < fMaxJetEta && (jet.Eta() < fMinJetEta || jet.Eta() > fMaxJetEta)) return kFALSE;
815  if (fMinJetPhi < fMaxJetPhi && (jet.Phi() < fMinJetPhi || jet.Phi() > fMaxJetPhi)) return kFALSE;
816  if (fMinJetPt < fMaxJetPt && (jet.Pt() > fMaxJetPt || jet.Pt() < fMinJetPt)) return kFALSE;
817  if (fMinChargedPt < fMaxChargedPt && (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt)) return kFALSE;
818  if (fMinNeutralPt < fMaxNeutralPt && (jet.fMaxNeutralPt < fMinNeutralPt || jet.fMaxNeutralPt > fMaxNeutralPt)) return kFALSE;
819 
820  return kTRUE;
821 }
822 
828 {
829  const AliJetInfo* jet = dMesonJet.GetJet(n);
830  if (!jet) return kFALSE;
831  return IsJetInAcceptance((*jet));
832 }
833 
840 {
841  if (lhs.fJetType > rhs.fJetType) return false;
842  else if (lhs.fJetType < rhs.fJetType) return true;
843  else {
844  if (lhs.fRadius > rhs.fRadius) return false;
845  else if (lhs.fRadius < rhs.fRadius) return true;
846  else {
847  if (lhs.fJetAlgo > rhs.fJetAlgo) return false;
848  else if (lhs.fJetAlgo < rhs.fJetAlgo) return true;
849  else {
850  if (lhs.fRecoScheme < rhs.fRecoScheme) return true;
851  else return false;
852  }
853  }
854  }
855 }
856 
863 {
864  if (lhs.fJetType != rhs.fJetType) return false;
865  if (lhs.fRadius != rhs.fRadius) return false;
866  if (lhs.fJetAlgo != rhs.fJetAlgo) return false;
867  if (lhs.fRecoScheme != rhs.fRecoScheme) return false;
868  return true;
869 }
870 
871 // Definitions of class AliAnalysisTaskDmesonJets::AnalysisEngine
872 
876 
879  TObject(),
880  fPartons(),
881  fCandidateType(kD0toKpi),
882  fCandidateName(),
883  fCandidatePDG(0),
884  fNDaughters(0),
885  fPDGdaughters(),
886  fBranchName(),
887  fMCMode(kNoMC),
888  fNMassBins(0),
889  fMinMass(0),
890  fMaxMass(0),
891  fRDHFCuts(0),
892  fRejectedOrigin(0),
893  fAcceptedDecay(0),
894  fInhibit(kFALSE),
895  fJetDefinitions(),
896  fPtBinWidth(0.5),
897  fMaxPt(100),
898  fD0Extended(kFALSE),
899  fRandomGen(0),
900  fTrackEfficiency(0),
901  fRejectISR(kFALSE),
902  fDataSlotNumber(-1),
903  fTree(0),
904  fCurrentDmesonJetInfo(0),
905  fCurrentJetInfo(0),
906  fCandidateArray(0),
907  fMCContainer(),
908  fTrackContainers(),
909  fClusterContainers(),
910  fAodEvent(0),
911  fFastJetWrapper(0),
912  fHistManager(0),
913  fCent(-1)
914 {
915 }
916 
925  TObject(),
926  fPartons(),
927  fCandidateType(type),
928  fCandidateName(),
929  fCandidatePDG(0),
930  fNDaughters(0),
931  fPDGdaughters(),
932  fBranchName(),
933  fMCMode(MCmode),
934  fNMassBins(nMassBins),
935  fMinMass(0),
936  fMaxMass(0),
937  fRDHFCuts(cuts),
938  fRejectedOrigin(0),
939  fAcceptedDecay(kAnyDecay),
940  fInhibit(kFALSE),
941  fJetDefinitions(),
942  fPtBinWidth(0.5),
943  fMaxPt(100),
944  fD0Extended(kFALSE),
945  fRandomGen(0),
946  fTrackEfficiency(0),
947  fDataSlotNumber(-1),
948  fTree(0),
949  fCurrentDmesonJetInfo(0),
950  fCurrentJetInfo(0),
951  fCandidateArray(0),
952  fMCContainer(),
953  fTrackContainers(),
954  fClusterContainers(),
955  fAodEvent(0),
956  fFastJetWrapper(0),
957  fHistManager(0),
958  fCent(-1)
959 {
960  SetCandidateProperties(range);
961 }
962 
967  TObject(source),
968  fPartons(source.fPartons),
969  fCandidateType(source.fCandidateType),
970  fCandidateName(source.fCandidateName),
971  fCandidatePDG(source.fCandidatePDG),
972  fNDaughters(source.fNDaughters),
973  fPDGdaughters(source.fPDGdaughters),
974  fBranchName(source.fBranchName),
975  fMCMode(source.fMCMode),
976  fNMassBins(source.fNMassBins),
977  fMinMass(source.fMinMass),
978  fMaxMass(source.fMaxMass),
979  fRDHFCuts(),
980  fRejectedOrigin(source.fRejectedOrigin),
981  fAcceptedDecay(source.fAcceptedDecay),
982  fInhibit(source.fInhibit),
983  fJetDefinitions(source.fJetDefinitions),
984  fPtBinWidth(source.fPtBinWidth),
985  fMaxPt(source.fMaxPt),
986  fRandomGen(source.fRandomGen),
988  fDataSlotNumber(-1),
989  fTree(0),
990  fCurrentDmesonJetInfo(0),
991  fCurrentJetInfo(0),
992  fCandidateArray(source.fCandidateArray),
993  fMCContainer(source.fMCContainer),
994  fTrackContainers(source.fTrackContainers),
995  fClusterContainers(source.fClusterContainers),
996  fAodEvent(source.fAodEvent),
998  fHistManager(source.fHistManager),
999  fCent(-1)
1000 {
1001  SetRDHFCuts(source.fRDHFCuts);
1002 }
1003 
1004 // Destructor
1006 {
1007  delete fRDHFCuts;
1008 }
1009 
1014 {
1015  new (this) AnalysisEngine(source);
1016  return *this;
1017 }
1018 
1023 {
1024  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
1025  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
1026  }
1027 
1028  return kFALSE;
1029 }
1030 
1032 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const /*geom*/, Int_t /*runNumber*/)
1033 {
1034 }
1035 
1040 {
1041  switch (fCandidateType) {
1042  case kD0toKpi:
1043  fCandidatePDG = 421;
1044  fCandidateName = "D0";
1045  fNDaughters = 2;
1046  fPDGdaughters.Set(fNDaughters);
1047  fPDGdaughters.Reset();
1048  fPDGdaughters[0] = 211; // pi
1049  fPDGdaughters[1] = 321; // K
1050  fBranchName = "D0toKpi";
1051  fAcceptedDecay = kDecayD0toKpi;
1052  break;
1053  case kD0toKpiLikeSign:
1054  fCandidatePDG = 421;
1055  fCandidateName = "2ProngLikeSign";
1056  fNDaughters = 2;
1057  fPDGdaughters.Set(fNDaughters);
1058  fPDGdaughters.Reset();
1059  fPDGdaughters[0] = 211; // pi
1060  fPDGdaughters[1] = 321; // K
1061  fBranchName = "LikeSign2Prong";
1062  break;
1063  case kDstartoKpipi:
1064  fCandidatePDG = 413;
1065  fCandidateName = "DStar";
1066  fNDaughters = 3;
1067  fPDGdaughters.Set(fNDaughters);
1068  fPDGdaughters.Reset();
1069  fPDGdaughters[0] = 211; // pi soft
1070  fPDGdaughters[1] = 211; // pi fromD0
1071  fPDGdaughters[2] = 321; // K from D0
1072  fBranchName = "Dstar";
1073  fAcceptedDecay = kDecayDStartoKpipi;
1074  break;
1075  default:
1076  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
1077  }
1078 
1079  CalculateMassLimits(range, fCandidatePDG, fNMassBins, fMinMass, fMaxMass);
1080 }
1081 
1086 {
1087  if (fRDHFCuts) delete fRDHFCuts;
1088  fRDHFCuts = cuts;
1089 }
1090 
1095 {
1096  if (!cuts) return;
1097  if (fRDHFCuts) delete fRDHFCuts;
1098  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
1099 }
1100 
1105 {
1106  static TString name;
1107 
1108  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
1109 
1110  return name.Data();
1111 }
1112 
1117 {
1118  fName = fCandidateName;
1119  switch (fMCMode) {
1120  case kBackgroundOnly:
1121  fName += "_kBackgroundOnly";
1122  break;
1123  case kSignalOnly:
1124  fName += "_kSignalOnly";
1125  break;
1126  case kMCTruth:
1127  fName += "_MCTruth";
1128  break;
1129  case kWrongPID:
1130  fName += "_WrongPID";
1131  break;
1132  default:
1133  break;
1134  }
1135 
1136  if (fRDHFCuts) fName += TString::Format("_%s", fRDHFCuts->GetName());
1137 
1138  return fName.Data();
1139 }
1140 
1148 {
1149  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
1150 
1151  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
1152  fJetDefinitions.push_back(def);
1153  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
1154  "Total number of jet definitions is now %lu.",
1155  def.GetName(), GetName(), fJetDefinitions.size());
1156  // For detector level set maximum track pt to 100 GeV/c
1157  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
1158  }
1159  else {
1160  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
1161  }
1162 
1163  return &(*it);
1164 }
1165 
1177 {
1178  AliHFJetDefinition def(type, r, algo, reco);
1179 
1180  return AddJetDefinition(def);
1181 }
1182 
1188 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
1189 {
1190  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
1191  while (it != fJetDefinitions.end() && (*it) != def) it++;
1192  return it;
1193 }
1194 
1199 {
1200  for (auto &jetdef : fJetDefinitions) jetdef.SetJetPhiRange(min, max);
1201 }
1202 
1207 {
1208  for (auto &jetdef : fJetDefinitions) jetdef.SetJetEtaRange(min, max);
1209 }
1210 
1215 {
1216  for (auto &jetdef : fJetDefinitions) jetdef.SetJetPtRange(min, max);
1217 }
1218 
1223 {
1224  for (auto &jetdef : fJetDefinitions) jetdef.SetChargedPtRange(min, max);
1225 }
1226 
1231 {
1232  for (auto &jetdef : fJetDefinitions) jetdef.SetNeutralPtRange(min, max);
1233 }
1234 
1241 {
1242  if (lhs.fCandidateType < rhs.fCandidateType) {
1243  return true;
1244  }
1245  else if (lhs.fCandidateType > rhs.fCandidateType) {
1246  return false;
1247  }
1248  else if (lhs.fMCMode < rhs.fMCMode) {
1249  return true;
1250  }
1251  else if (lhs.fMCMode > rhs.fMCMode) {
1252  return false;
1253  }
1254  else if (lhs.fRDHFCuts && !rhs.fRDHFCuts) {
1255  return true;
1256  }
1257  else if (lhs.fRDHFCuts && rhs.fRDHFCuts && strcmp(lhs.fRDHFCuts->GetName(), rhs.fRDHFCuts->GetName()) < 0) {
1258  return true;
1259  }
1260  else {
1261  return false;
1262  }
1263 }
1264 
1271 {
1272  if (lhs.fCandidateType != rhs.fCandidateType) return false;
1273  if (lhs.fMCMode != rhs.fMCMode) return false;
1274  if (lhs.fRDHFCuts == nullptr && rhs.fRDHFCuts != nullptr) return false;
1275  if (lhs.fRDHFCuts != nullptr && rhs.fRDHFCuts == nullptr) return false;
1276  if (lhs.fRDHFCuts && rhs.fRDHFCuts && strcmp(lhs.fRDHFCuts->GetName(), rhs.fRDHFCuts->GetName()) != 0) return false;
1277  return true;
1278 }
1279 
1288 {
1289  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) { // D0 candidate
1290  return ExtractD0Attributes(Dcand, DmesonJet, i);
1291  }
1292  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
1293  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
1294  }
1295  else {
1296  return kFALSE;
1297  }
1298 }
1299 
1308 {
1309  AliDebug(10,"Checking if D0 meson is selected");
1310  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
1311  if (isSelected == 0) return kFALSE;
1312 
1313  Int_t MCtruthPdgCode = 0;
1314 
1315  Double_t invMassD = 0;
1316 
1317  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
1318  // Checks also the origin, and if it matches the rejected origin mask, return false
1319  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly || fMCMode == kWrongPID) {
1320  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
1321  DmesonJet.fMCLabel = mcLab;
1322 
1323  // Retrieve the generated particle (if exists) and its PDG code
1324  if (mcLab >= 0) {
1325  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1326 
1327  if (aodMcPart) {
1328  // Check origin and return false if it matches the rejected origin mask
1329  if (fRejectedOrigin) {
1330  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
1331  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1332  }
1333  MCtruthPdgCode = aodMcPart->PdgCode();
1334  }
1335  }
1336  }
1337 
1338  if (isSelected == 1) { // selected as a D0
1339  if (i != 0) return kFALSE; // only one mass hypothesis thanks to PID
1340 
1341  if (fMCMode == kNoMC ||
1342  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1343  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly) ||
1344  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kWrongPID)) {
1345  // 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)
1346  AliDebug(10,"Selected as D0");
1347  invMassD = Dcand->InvMassD0();
1348  }
1349  else { // conditions above not passed, so return FALSE
1350  return kFALSE;
1351  }
1352  }
1353  else if (isSelected == 2) { // selected as a D0bar
1354  if (i != 1) return kFALSE; // only one mass hypothesis thanks to PID
1355 
1356  if (fMCMode == kNoMC ||
1357  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1358  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly) ||
1359  (MCtruthPdgCode == fCandidatePDG && fMCMode == kWrongPID)) {
1360  // 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)
1361  AliDebug(10,"Selected as D0bar");
1362  invMassD = Dcand->InvMassD0bar();
1363  }
1364  else { // conditions above not passed, so return FALSE
1365  return kFALSE;
1366  }
1367  }
1368  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
1369  AliDebug(10,"Selected as either D0 or D0bar");
1370 
1371  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
1372  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1373  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kWrongPID))) {
1374  if (i != 0) return kFALSE;
1375  AliDebug(10, "MC truth is D0");
1376  invMassD = Dcand->InvMassD0();
1377  }
1378  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1379  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kWrongPID))) {
1380  if (i != 1) return kFALSE;
1381  AliDebug(10, "MC truth is D0bar");
1382  invMassD = Dcand->InvMassD0bar();
1383  }
1384  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
1385 
1386  // Only accept it if background-only OR background-and-signal was requested
1387  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
1388  // Select D0 or D0bar depending on the i-parameter
1389  if (i == 0) {
1390  AliDebug(10, "Returning invariant mass with D0 hypothesis");
1391  invMassD = Dcand->InvMassD0();
1392  }
1393  else if (i == 1) {
1394  AliDebug(10, "Returning invariant mass with D0bar hypothesis");
1395  invMassD = Dcand->InvMassD0bar();
1396  }
1397  else { // i > 1
1398  return kFALSE;
1399  }
1400  }
1401  else { // signal-only was requested but this is not a true D0
1402  return kFALSE;
1403  }
1404  }
1405  }
1406  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
1407  return kTRUE;
1408 }
1409 
1418 {
1419  AliDebug(10,"Checking if D* meson is selected");
1420  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoCascadeHF*>(DstarCand), AliRDHFCuts::kAll, fAodEvent);
1421  if (isSelected == 0) return kFALSE;
1422 
1423  if ((i == 1 && DstarCand->Charge()>0) || (i == 0 && DstarCand->Charge()<0) || i > 1) return kFALSE; // only one mass hypothesis for the D*
1424 
1425  Int_t MCtruthPdgCode = 0;
1426 
1427  Double_t invMassD = 0;
1428 
1429  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1430  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1431  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1432 
1433  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1434  AliDebug(10, Form("MC label is %d", mcLab));
1435  DmesonJet.fMCLabel = mcLab;
1436  if (mcLab >= 0) {
1437  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1438 
1439  if (aodMcPart) {
1440  if (fRejectedOrigin) {
1441  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
1442  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1443  }
1444 
1445  MCtruthPdgCode = aodMcPart->PdgCode();
1446  AliDebug(10, Form("MC truth pdg code is %d",MCtruthPdgCode));
1447  }
1448  }
1449  }
1450 
1451  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1452  if (fMCMode == kNoMC ||
1453  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1454  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1455  // 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)
1456  invMassD = DstarCand->InvMassDstarKpipi();
1457  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1458  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1459  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1460  return kTRUE;
1461  }
1462  else { // conditions above not passed, so return FALSE
1463  return kFALSE;
1464  }
1465 }
1466 
1474 {
1475  if (!part) return kUnknownDecay;
1476  if (!mcArray) return kUnknownDecay;
1477 
1479 
1480  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1481 
1482  if (part->GetNDaughters() == 2) {
1483 
1484  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1485  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1486 
1487  if (!d1 || !d2) {
1488  return decay;
1489  }
1490 
1491  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1492  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1493 
1494  if (absPdgPart == 421) { // D0 -> K pi
1495  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1496  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1497  decay = kDecayD0toKpi;
1498  }
1499  }
1500 
1501  if (absPdgPart == 413) { // D* -> D0 pi
1502  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1503  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1504  if (D0decay == kDecayD0toKpi) {
1505  decay = kDecayDStartoKpipi;
1506  }
1507  }
1508  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1509  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1510  if (D0decay == kDecayD0toKpi) {
1511  decay = kDecayDStartoKpipi;
1512  }
1513  }
1514  }
1515  }
1516 
1517  return decay;
1518 }
1519 
1526 std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> AliAnalysisTaskDmesonJets::AnalysisEngine::IsPromptCharm(const AliAODMCParticle* part, TClonesArray* mcArray)
1527 {
1528  std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> result(kUnknownQuark, 0);
1529 
1530  if (!part) return result;
1531  if (!mcArray) return result;
1532 
1533  static std::set<UInt_t> partons = { 4, 5 };
1534 
1535  AliAODMCParticle* parton = FindParticleOrigin(part, mcArray, kFindLast, partons);
1536  if (parton) {
1537  result.second = parton;
1538  UInt_t absPdgParton = TMath::Abs(parton->GetPdgCode());
1539  if (absPdgParton == 4) result.first = kFromCharm;
1540  else if (absPdgParton == 5) result.first = kFromBottom;
1541  }
1542 
1543  return result;
1544 }
1545 
1553 
1554 AliAODMCParticle* AliAnalysisTaskDmesonJets::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode)
1555 {
1556  static std::set<UInt_t> pdgSet;
1557 
1558  return FindParticleOrigin(part, mcArray, mode, pdgSet);
1559 }
1560 
1569 
1570 AliAODMCParticle* AliAnalysisTaskDmesonJets::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode, const std::set<UInt_t>& pdgSet)
1571 {
1572  AliAODMCParticle* result = nullptr;
1573 
1574  Int_t mother = part->GetMother();
1575  while (mother >= 0) {
1576  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1577  if (mcGranma) {
1578  UInt_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
1579 
1580  // If the current particle is one of the particle types that is being searched assign it to the result pointer
1581  if (pdgSet.empty() || pdgSet.count(abspdgGranma) > 0) {
1582  result = mcGranma;
1583  // If the last particle in the fragmentation tree (first when going reverse) was requested then stop the loop
1584  if (mode == kFindLast) break;
1585  }
1586  mother = mcGranma->GetMother();
1587  }
1588  else {
1589  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::FindParticleOrigin", "Could not retrieve mother particle %d!", mother);
1590  break;
1591  }
1592  }
1593 
1594  return result;
1595 }
1596 
1599 {
1600  for (auto& jetDef : fJetDefinitions) {
1601  jetDef.fJets.clear();
1602  }
1603 
1604  if (fMCMode == kMCTruth) {
1605  RunParticleLevelAnalysis();
1606  }
1607  else {
1608  RunDetectorLevelAnalysis();
1609  }
1610 }
1611 
1614 {
1615  // Fill the vertex info of the candidates
1616  // Needed for reduced delta AOD, where the vertex info has been deleted
1617  // to reduce the delta AOD file size
1619 
1620  const Int_t nD = fCandidateArray->GetEntriesFast();
1621 
1622  AliDmesonJetInfo DmesonJet;
1623  DmesonJet.fEvent = this->fAodEvent;
1624 
1625  std::map<AliHFJetDefinition*,Double_t> maxJetPt;
1626  for (auto& def : fJetDefinitions) maxJetPt[&def] = 0;
1627  Double_t maxDPt = 0;
1628 
1629  Int_t nAccCharm[3] = {0};
1630  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1631  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1632  if (!charmCand) continue;
1633  if(!(vHF.FillRecoCand(fAodEvent,charmCand))) continue;
1634 
1635  // region of interest + cuts
1636  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1637  Int_t nMassHypo = 0; // number of mass hypothesis accepted for this D meson
1638  if (charmCand->Pt() > maxDPt) maxDPt = charmCand->Pt();
1639  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1640  DmesonJet.Reset();
1641  DmesonJet.fDmesonParticle = charmCand;
1642  DmesonJet.fSelectionType = im + 1;
1643  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1644  for (auto& def : fJetDefinitions) {
1645  if (FindJet(charmCand, DmesonJet, def)) {
1646  Double_t jetPt = DmesonJet.fJets[def.GetName()].fMomentum.Pt();
1647  if (jetPt > maxJetPt[&def]) maxJetPt[&def] = jetPt;
1648  }
1649  else {
1650  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1651  def.GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1652  }
1653  }
1654  fDmesonJets[(icharm+1)*(1-(im*2))] = DmesonJet;
1655  nMassHypo++;
1656  nAccCharm[im]++;
1657  }
1658  }
1659  if (nMassHypo == 2) {
1660  nAccCharm[0]--;
1661  nAccCharm[1]--;
1662  nAccCharm[2] += 2;
1663  }
1664  if (nMassHypo == 2) { // both mass hypothesis accepted
1665  fDmesonJets[(icharm+1)].fD0D0bar = kTRUE;
1666  fDmesonJets[-(icharm+1)].fD0D0bar = kTRUE;
1667  }
1668  } // end of D cand loop
1669 
1670  TString hname;
1671 
1672  Int_t ntracks = 0;
1673 
1674  for (auto track_cont : fTrackContainers) {
1675  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
1676  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(nullptr);
1677  ntracks += track_cont->GetNAcceptEntries();
1678  }
1679 
1680  for (auto& def : fJetDefinitions) {
1681  if (!def.fRho) continue;
1682  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
1683  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
1684 
1685  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
1686  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
1687 
1688  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
1689  fHistManager->FillTH2(hname, fCent, def.fRho->GetVal());
1690 
1691  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
1692  fHistManager->FillTH2(hname, fCent, maxJetPt[&def]);
1693 
1694  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
1695  fHistManager->FillTH2(hname, fCent, maxDPt);
1696 
1697  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
1698  fHistManager->FillTH2(hname, ntracks, def.fRho->GetVal());
1699 
1700  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
1701  fHistManager->FillTH2(hname, ntracks, maxJetPt[&def]);
1702 
1703  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
1704  fHistManager->FillTH2(hname, ntracks, maxDPt);
1705  }
1706 
1707  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1708  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1709  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1710  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1711 
1712  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1713  fHistManager->FillTH2(hname, ntracks, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1714 
1715  hname = TString::Format("%s/fHistNDmesons", GetName());
1716  fHistManager->FillTH1(hname, nD);
1717 }
1718 
1729 {
1730  TString hname;
1731 
1732  Double_t rho = 0;
1733  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
1734 
1736  fFastJetWrapper->SetR(jetDef.fRadius);
1739 
1740  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1741 
1742  if (jetDef.fJetType != AliJetContainer::kNeutralJet) {
1743  for (auto track_cont : fTrackContainers) {
1744  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
1745  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(Dcand);
1746  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1747  AddInputVectors(track_cont, 100, static_cast<TH2*>(fHistManager->FindObject(hname)), fTrackEfficiency);
1748 
1749  if (hftrack_cont) {
1750  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1751  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1752  const TObjArray& daughters = hftrack_cont->GetDaughterList();
1753  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1754  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1755  if (!hftrack_cont->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1756  }
1757  }
1758  }
1759  }
1760 
1761  if (jetDef.fJetType != AliJetContainer::kChargedJet) {
1762  for (auto clus_cont : fClusterContainers) {
1763  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1764  AddInputVectors(clus_cont, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1765  }
1766  }
1767 
1768  // run jet finder
1769  fFastJetWrapper->Run();
1770 
1771  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1772 
1773  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1774  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1775 
1776  Bool_t isDmesonJet = kFALSE;
1777 
1778  Double_t maxChPt = 0;
1779  Double_t maxNePt = 0;
1780  Double_t totalNeutralPt = 0;
1781  Int_t nConst = 1;
1782 
1783  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1784  if (constituents[ic].user_index() == 0) {
1785  isDmesonJet = kTRUE;
1786  }
1787  else if (constituents[ic].user_index() >= 100) {
1788  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1789  nConst++;
1790  }
1791  else if (constituents[ic].user_index() <= -100) {
1792  totalNeutralPt += constituents[ic].pt();
1793  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1794  nConst++;
1795  }
1796  }
1797 
1798  if (isDmesonJet) {
1799  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1800  DmesonJet.fJets[jetDef.GetName()].fNConstituents = nConst;
1801  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1802  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1803  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1804  DmesonJet.fJets[jetDef.GetName()].fArea = jets_incl[ijet].area();
1805  DmesonJet.fJets[jetDef.GetName()].fCorrPt = DmesonJet.fJets[jetDef.GetName()].fMomentum.Pt() - jets_incl[ijet].area() * rho;
1806 
1807  return kTRUE;
1808  }
1809  }
1810 
1811  return kFALSE;
1812 }
1813 
1817 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist, Double_t eff)
1818 {
1819  auto itcont = cont->all_momentum();
1820  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1821  UInt_t rejectionReason = 0;
1822  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1823  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1824  continue;
1825  }
1826  if (fRandomGen && eff > 0 && eff < 1) {
1827  Double_t rnd = fRandomGen->Rndm();
1828  if (eff < rnd) {
1829  if (rejectHist) rejectHist->Fill(6, it->first.Pt());
1830  continue;
1831  }
1832  }
1833  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1834  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1835  }
1836 }
1837 
1840 {
1841  TString hname;
1842 
1843  fMCContainer->SetSpecialPDG(fCandidatePDG);
1844  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1845  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1847 
1848  if (!fMCContainer->IsSpecialPDGFound()) return;
1849 
1850  Int_t nAccCharm[3] = {0};
1851 
1852  std::map<AliHFJetDefinition*, Double_t> maxJetPt;
1853  Double_t maxDPt = 0;
1854 
1855  for (auto &jetDef : fJetDefinitions) {
1856  maxJetPt[&jetDef] = 0;
1857  Double_t rho = 0;
1858  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
1859  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", GetName(), jetDef.GetName());
1860  TH1* histNDmesonsVsNconstituents = static_cast<TH1*>(fHistManager->FindObject(hname));
1861 
1862  switch (jetDef.fJetType) {
1864  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
1865  break;
1867  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
1868  break;
1870  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
1871  break;
1872  }
1873 
1875  fFastJetWrapper->SetR(jetDef.fRadius);
1878 
1879  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1880  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1881 
1882  fFastJetWrapper->Run();
1883 
1884  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1885 
1886  for (auto jet : jets_incl) {
1887  Int_t nDmesonsInJet = 0;
1888 
1889  for (auto constituent : jet.constituents()) {
1890  Int_t iPart = constituent.user_index() - 100;
1891  if (constituent.perp() < 1e-6) continue; // reject ghost particles
1892  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
1893  if (!part) {
1894  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1895  continue;
1896  }
1897  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1898  nDmesonsInJet++;
1899  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1900  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1901  if (part->Pt() > maxDPt) maxDPt = part->Pt();
1902  std::pair<int, AliDmesonJetInfo> element;
1903  element.first = iPart;
1904  dMesonJetIt = fDmesonJets.insert(element).first;
1905  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1906  (*dMesonJetIt).second.fDmesonParticle = part;
1907  (*dMesonJetIt).second.fSelectionType = part->PdgCode() > 0 ? 1 : 2;
1908 
1909  UShort_t p = 0;
1910  UInt_t rs = 0;
1911 
1912  auto origin = IsPromptCharm(part, fMCContainer->GetArray());
1913  p = 0;
1914  rs = origin.first;
1915  while (rs >>= 1) { p++; }
1916  (*dMesonJetIt).second.fPartonType = p;
1917  (*dMesonJetIt).second.fParton = origin.second;
1918 
1919  (*dMesonJetIt).second.fAncestor = FindParticleOrigin(part, fMCContainer->GetArray(), kFindFirst);
1920 
1921  if (part->PdgCode() > 0) { // D0
1922  nAccCharm[0]++;
1923  }
1924  else { // D0bar
1925  nAccCharm[1]++;
1926  }
1927  }
1928 
1929  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
1930  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = jet.constituents().size();
1931  (*dMesonJetIt).second.fJets[jetDef.GetName()].fArea = jet.area();
1932  (*dMesonJetIt).second.fJets[jetDef.GetName()].fCorrPt = (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.Pt() - jet.area() * rho;
1933  if (jet.perp() > maxJetPt[&jetDef]) maxJetPt[&jetDef] = jet.perp();
1934  } // if constituent is a D meson
1935  } // for each constituent
1936  if (nDmesonsInJet > 0) histNDmesonsVsNconstituents->Fill(jet.constituents().size(), nDmesonsInJet);
1937  } // for each jet
1938  } // for each jet definition
1939 
1941 
1942  for (auto& def : fJetDefinitions) {
1943  if (!def.fRho) continue;
1944  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
1945  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
1946 
1947  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
1948  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
1949 
1950  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
1951  fHistManager->FillTH2(hname, fCent, def.fRho->GetVal());
1952 
1953  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
1954  fHistManager->FillTH2(hname, fCent, maxJetPt[&def]);
1955 
1956  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
1957  fHistManager->FillTH2(hname, fCent, maxDPt);
1958 
1959  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
1960  fHistManager->FillTH2(hname, npart, def.fRho->GetVal());
1961 
1962  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
1963  fHistManager->FillTH2(hname, npart, maxJetPt[&def]);
1964 
1965  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
1966  fHistManager->FillTH2(hname, npart, maxDPt);
1967  }
1968 
1969  if (fDmesonJets.size() != nAccCharm[0]+nAccCharm[1]) AliError(Form("I found %lu mesons (%d)?", fDmesonJets.size(), nAccCharm[0]+nAccCharm[1]));
1970  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1971  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1972  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1973  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1974 
1975  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1976  fHistManager->FillTH2(hname, npart, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1977 
1978  hname = TString::Format("%s/fHistNDmesons", GetName());
1979  fHistManager->FillTH1(hname, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]); // same as the number of accepted D mesons, since no selection is performed
1980 }
1981 
1986 {
1987  TString classname;
1988  if (fMCMode == kMCTruth) {
1989  classname = "AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary";
1990  fCurrentDmesonJetInfo = new AliDmesonMCInfoSummary();
1991  }
1992  else {
1993  switch (fCandidateType) {
1994  case kD0toKpi:
1995  case kD0toKpiLikeSign:
1996  if (fD0Extended) {
1997  classname = "AliAnalysisTaskDmesonJets::AliD0ExtendedInfoSummary";
1998  fCurrentDmesonJetInfo = new AliD0ExtendedInfoSummary();
1999  }
2000  else {
2001  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
2002  fCurrentDmesonJetInfo = new AliD0InfoSummary();
2003  }
2004  break;
2005  case kDstartoKpipi:
2006  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
2007  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
2008  break;
2009  }
2010  }
2011  TString treeName = TString::Format("%s_%s", taskName, GetName());
2012  fTree = new TTree(treeName, treeName);
2013  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
2014  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
2015  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
2016  if (fJetDefinitions[i].fRhoName.IsNull()) {
2017  fCurrentJetInfo[i] = new AliJetInfoSummary();
2018  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
2019  }
2020  else {
2021  fCurrentJetInfo[i] = new AliJetInfoPbPbSummary();
2022  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoPbPbSummary", &fCurrentJetInfo[i]);
2023  }
2024  }
2025 
2026  return fTree;
2027 }
2028 
2033 {
2034  TString hname;
2035 
2036  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
2037 
2038  for (auto &jetDef : fJetDefinitions) {
2039 
2040  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
2041 
2042  Double_t radius = jetDef.fRadius;
2043 
2044  TString title[30] = {""};
2045  Int_t nbins[30] = {0 };
2046  Double_t min [30] = {0.};
2047  Double_t max [30] = {0.};
2048  Int_t dim = 0 ;
2049 
2050  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
2051  nbins[dim] = nPtBins;
2052  min[dim] = 0;
2053  max[dim] = fMaxPt;
2054  dim++;
2055 
2056  if ((enabledAxis & kPositionD) != 0) {
2057  title[dim] = "#eta_{D}";
2058  nbins[dim] = 50;
2059  min[dim] = -1;
2060  max[dim] = 1;
2061  dim++;
2062 
2063  title[dim] = "#phi_{D} (rad)";
2064  nbins[dim] = 150;
2065  min[dim] = 0;
2066  max[dim] = TMath::TwoPi();
2067  dim++;
2068  }
2069 
2070  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
2071  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
2072  nbins[dim] = fNMassBins;
2073  min[dim] = fMinMass;
2074  max[dim] = fMaxMass;
2075  dim++;
2076  }
2077 
2078  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
2079  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
2080  nbins[dim] = fNMassBins;
2081  min[dim] = fMinMass;
2082  max[dim] = fMaxMass;
2083  dim++;
2084  }
2085 
2086  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
2087  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
2088  nbins[dim] = fNMassBins;
2089  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
2090  dim++;
2091  }
2092 
2093  if (fCandidateType == kDstartoKpipi) {
2094  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
2095  nbins[dim] = fNMassBins*6;
2096  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
2097 
2098  // subtract mass of D0
2099  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2100  min[dim] -= D0mass;
2101  max[dim] -= D0mass;
2102 
2103  dim++;
2104  }
2105 
2106  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
2107  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
2108  nbins[dim] = 100;
2109  min[dim] = 0;
2110  max[dim] = 25;
2111  dim++;
2112  }
2113 
2114  title[dim] = "#it{z}_{D}";
2115  nbins[dim] = 110;
2116  min[dim] = 0;
2117  max[dim] = 1.10;
2118  dim++;
2119 
2120  if ((enabledAxis & kDeltaR) != 0) {
2121  title[dim] = "#Delta R_{D-jet}";
2122  nbins[dim] = 100;
2123  min[dim] = 0;
2124  max[dim] = radius * 1.5;
2125  dim++;
2126  }
2127 
2128  if ((enabledAxis & kDeltaEta) != 0) {
2129  title[dim] = "#eta_{D} - #eta_{jet}";
2130  nbins[dim] = 100;
2131  min[dim] = -radius * 1.2;
2132  max[dim] = radius * 1.2;
2133  dim++;
2134  }
2135 
2136  if ((enabledAxis & kDeltaPhi) != 0) {
2137  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
2138  nbins[dim] = 100;
2139  min[dim] = -radius * 1.2;
2140  max[dim] = radius * 1.2;
2141  dim++;
2142  }
2143 
2144  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
2145  nbins[dim] = nPtBins;
2146  min[dim] = 0;
2147  max[dim] = fMaxPt;
2148  dim++;
2149 
2150  if ((enabledAxis & kPositionJet) != 0) {
2151  title[dim] = "#eta_{jet}";
2152  nbins[dim] = 50;
2153  min[dim] = -1;
2154  max[dim] = 1;
2155  dim++;
2156 
2157  title[dim] = "#phi_{jet} (rad)";
2158  nbins[dim] = 150;
2159  min[dim] = 0;
2160  max[dim] = TMath::TwoPi();
2161  dim++;
2162  }
2163 
2164  if ((enabledAxis & kJetConstituents) != 0) {
2165  title[dim] = "No. of constituents";
2166  nbins[dim] = 50;
2167  min[dim] = -0.5;
2168  max[dim] = 49.5;
2169  dim++;
2170  }
2171 
2172  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
2173  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
2174  for (Int_t j = 0; j < dim; j++) {
2175  h->GetAxis(j)->SetTitle(title[j]);
2176  }
2177  }
2178 }
2179 
2184 {
2185  TString hname;
2186  fPartons.clear();
2187 
2188  TH1* histAncestor = nullptr;
2189  TH1* histPrompt = nullptr;
2190 
2191  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2192  hname = TString::Format("%s/fHistPrompt", GetName());
2193  histPrompt = static_cast<TH1*>(fHistManager->FindObject(hname));
2194 
2195  hname = TString::Format("%s/fHistAncestor", GetName());
2196  histAncestor = static_cast<TH1*>(fHistManager->FindObject(hname));
2197  }
2198 
2199  for (auto& dmeson_pair : fDmesonJets) {
2200  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
2201  Int_t accJets = 0;
2202  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
2203  fCurrentJetInfo[ij]->Reset();
2204  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
2205  if (!jet) continue;
2206  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
2207  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
2208  fHistManager->FillTH1(hname, jet->Pt());
2209  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
2210  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2211  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
2212  fHistManager->FillTH1(hname, jet->Eta());
2213  continue;
2214  }
2215  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
2216  accJets++;
2217  }
2218  if (accJets > 0) {
2219  if (histPrompt) {
2220  if (dmeson_pair.second.fParton) {
2221  fPartons[dmeson_pair.second.fParton] = dmeson_pair.second.fPartonType;
2222  UInt_t absPdgParton = TMath::Abs(dmeson_pair.second.fParton->GetPdgCode());
2223  if (absPdgParton == 4) {
2224  histPrompt->Fill("Prompt", 1);
2225  }
2226  else if (absPdgParton == 5) {
2227  histPrompt->Fill("Non-Prompt", 1);
2228  }
2229  else {
2230  histPrompt->Fill("Unknown", 1);
2231  }
2232  }
2233  else {
2234  histPrompt->Fill("Unknown", 1);
2235  }
2236  }
2237 
2238  if (histAncestor) {
2239  if (dmeson_pair.second.fAncestor) {
2240  UInt_t absPdgAncestor = TMath::Abs(dmeson_pair.second.fAncestor->GetPdgCode());
2241  if (absPdgAncestor == 4) {
2242  histAncestor->Fill("Charm", 1);
2243  }
2244  else if (absPdgAncestor == 5) {
2245  histAncestor->Fill("Bottom", 1);
2246  }
2247  else if (absPdgAncestor == 2212) {
2248  histAncestor->Fill("Proton", 1);
2249  }
2250  else {
2251  histAncestor->Fill("Unknown", 1);
2252  }
2253  }
2254  else {
2255  histAncestor->Fill("Unknown", 1);
2256  }
2257  }
2258 
2259  fTree->Fill();
2260  }
2261  else {
2262  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2263  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2264  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2265  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2266  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2267  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2268  if (fMCMode != kMCTruth) {
2269  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
2270  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
2271  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
2272  }
2273  else if (fCandidateType == kDstartoKpipi) {
2274  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
2275  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
2276 
2277  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
2278  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
2279  }
2280  }
2281  }
2282  }
2283 
2284  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2285  hname = TString::Format("%s/fHistPartonPt", GetName());
2286  TH1* histPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
2287  hname = TString::Format("%s/fHistPartonEta", GetName());
2288  TH1* histPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
2289  hname = TString::Format("%s/fHistPartonPhi", GetName());
2290  TH1* histPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
2291  hname = TString::Format("%s/fHistPartonType", GetName());
2292  TH1* histPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
2293 
2294  for (auto parton : fPartons) {
2295  if (!parton.first) continue;
2296  histPartonPt->Fill(parton.first->Pt());
2297  histPartonEta->Fill(parton.first->Eta());
2298  histPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
2299  histPartonType->Fill(parton.second);
2300  }
2301  }
2302 
2303  return kTRUE;
2304 }
2305 
2311 {
2312  TString hname;
2313 
2314  TH1* histAncestor = nullptr;
2315  TH1* histPrompt = nullptr;
2316 
2317  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2318  hname = TString::Format("%s/fHistPrompt", GetName());
2319  histPrompt = static_cast<TH1*>(fHistManager->FindObject(hname));
2320 
2321  hname = TString::Format("%s/fHistAncestor", GetName());
2322  histAncestor = static_cast<TH1*>(fHistManager->FindObject(hname));
2323  }
2324 
2325  fPartons.clear();
2326  for (auto& dmeson_pair : fDmesonJets) {
2327  Int_t accJets = 0;
2328  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
2329  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
2330  if (!jet) continue;
2331  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
2332  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
2333  fHistManager->FillTH1(hname, jet->Pt());
2334  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
2335  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2336  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
2337  fHistManager->FillTH1(hname, jet->Eta());
2338  continue;
2339  }
2340  accJets++;
2341  }
2342  if (accJets > 0) {
2343  if (histPrompt) {
2344  if (dmeson_pair.second.fParton) {
2345  fPartons[dmeson_pair.second.fParton] = dmeson_pair.second.fPartonType;
2346  UInt_t absPdgParton = TMath::Abs(dmeson_pair.second.fParton->GetPdgCode());
2347  if (absPdgParton == 4) {
2348  histPrompt->Fill("Prompt", 1);
2349  }
2350  else if (absPdgParton == 5) {
2351  histPrompt->Fill("Non-Prompt", 1);
2352  }
2353  else {
2354  histPrompt->Fill("Unknown", 1);
2355  }
2356  }
2357  else {
2358  histPrompt->Fill("Unknown", 1);
2359  }
2360  }
2361 
2362  if (histAncestor) {
2363  if (dmeson_pair.second.fAncestor) {
2364  UInt_t absPdgAncestor = TMath::Abs(dmeson_pair.second.fAncestor->GetPdgCode());
2365  if (absPdgAncestor == 4) {
2366  histAncestor->Fill("Charm", 1);
2367  }
2368  else if (absPdgAncestor == 5) {
2369  histAncestor->Fill("Bottom", 1);
2370  }
2371  else if (absPdgAncestor == 2212) {
2372  histAncestor->Fill("Proton", 1);
2373  }
2374  else {
2375  histAncestor->Fill("Unknown", 1);
2376  }
2377  }
2378  else {
2379  histAncestor->Fill("Unknown", 1);
2380  }
2381  }
2382  }
2383  else {
2384  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2385  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2386  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2387  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2388  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2389  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2390  if (fMCMode != kMCTruth) {
2391  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
2392  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
2393  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
2394  }
2395  else if (fCandidateType == kDstartoKpipi) {
2396  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
2397  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
2398 
2399  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
2400  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
2401  }
2402  }
2403  }
2404  }
2405 
2406  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2407  hname = TString::Format("%s/fHistPartonPt", GetName());
2408  TH1* histPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
2409  hname = TString::Format("%s/fHistPartonEta", GetName());
2410  TH1* histPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
2411  hname = TString::Format("%s/fHistPartonPhi", GetName());
2412  TH1* histPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
2413  hname = TString::Format("%s/fHistPartonType", GetName());
2414  TH1* histPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
2415 
2416  for (auto parton : fPartons) {
2417  if (!parton.first) continue;
2418  histPartonPt->Fill(parton.first->Pt());
2419  histPartonEta->Fill(parton.first->Eta());
2420  histPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
2421  histPartonType->Fill(parton.second);
2422  }
2423  }
2424 
2425  return kTRUE;
2426 }
2427 
2432 {
2433  TString hname;
2434 
2435  for (auto& dmeson_pair : fDmesonJets) {
2436  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
2437  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2438  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2439  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2440  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2441  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2442  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2443  }
2444  }
2445 
2446  for (auto &jetDef : fJetDefinitions) {
2447 
2448  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
2449  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
2450 
2451  for (auto& dmeson_pair : fDmesonJets) {
2452  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
2453  if (!jet) continue;
2454  if (!jetDef.IsJetInAcceptance(*jet)) {
2455  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
2456  fHistManager->FillTH1(hname, jet->Pt());
2457  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
2458  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2459  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
2460  fHistManager->FillTH1(hname, jet->Eta());
2461  continue;
2462  }
2463  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
2464  }
2465  }
2466 
2467  return kTRUE;
2468 }
2469 
2476 {
2477  // Fill the THnSparse histogram.
2478 
2479  Double_t contents[30] = {0.};
2480 
2481  Double_t z = DmesonJet.GetZ(n);
2482  Double_t deltaPhi = 0;
2483  Double_t deltaEta = 0;
2484  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
2485 
2486  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
2487  if (it == DmesonJet.fJets.end()) return kFALSE;
2488 
2489  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
2490  TString title(h->GetAxis(i)->GetTitle());
2491  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
2492  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
2493  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
2494  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
2495  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
2496  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
2497  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
2498  else if (title=="#it{z}_{D}") contents[i] = z;
2499  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
2500  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
2501  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
2502  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
2503  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
2504  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
2505  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
2506  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
2507  }
2508 
2509  h->Fill(contents);
2510 
2511  return kTRUE;
2512 }
2513 
2514 // Definitions of class AliAnalysisTaskDmesonJets
2515 
2517 ClassImp(AliAnalysisTaskDmesonJets);
2519 
2523  fAnalysisEngines(),
2524  fEnabledAxis(0),
2526  fHistManager(),
2527  fApplyKinematicCuts(kTRUE),
2528  fNOutputTrees(0),
2529  fTrackEfficiency(0),
2530  fRejectISR(kFALSE),
2531  fJetAreaType(fastjet::active_area),
2532  fJetGhostArea(0.005),
2533  fMCContainer(0),
2534  fAodEvent(0),
2535  fFastJetWrapper(0)
2536 {
2537  SetMakeGeneralHistograms(kTRUE);
2538 }
2539 
2544  AliAnalysisTaskEmcalLight(name, kTRUE),
2545  fAnalysisEngines(),
2546  fEnabledAxis(k2ProngInvMass),
2547  fOutputType(kTreeOutput),
2548  fHistManager(name),
2549  fApplyKinematicCuts(kTRUE),
2550  fNOutputTrees(nOutputTrees),
2551  fTrackEfficiency(0),
2552  fRejectISR(kFALSE),
2553  fJetAreaType(fastjet::active_area),
2554  fJetGhostArea(0.005),
2555  fMCContainer(0),
2556  fAodEvent(0),
2557  fFastJetWrapper(0)
2558 {
2559  SetMakeGeneralHistograms(kTRUE);
2560  for (Int_t i = 0; i < nOutputTrees; i++){
2561  DefineOutput(2+i, TTree::Class());
2562  }
2563 }
2564 
2567 {
2568  if (fFastJetWrapper) delete fFastJetWrapper;
2569 }
2570 
2578 {
2579  AliRDHFCuts* analysiscuts = 0;
2580  TFile* filecuts = TFile::Open(cutfname);
2581  if (!filecuts || filecuts->IsZombie()) {
2582  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Input file not found: will use std cuts.");
2583  filecuts = 0;
2584  }
2585 
2586  if (filecuts) analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
2587 
2588  if (!analysiscuts) {
2589  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Could not find analysis cuts '%s' in '%s'.", cutsname.Data(), cutfname.Data());
2590  if (filecuts) {
2591  filecuts->ls();
2592  }
2593  }
2594  else {
2595  ::Info("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Cuts '%s' loaded from file '%s'", cutsname.Data(), cutfname.Data());
2596  }
2597 
2598  return analysiscuts;
2599 }
2600 
2613 {
2615  return AddAnalysisEngine(type, cutfname, cuttype, MCmode, jetDef, rhoName);
2616 }
2617 
2629 {
2630  AliRDHFCuts* cuts = 0;
2631 
2632  if (!cutfname.IsNull()) {
2633  TString cutsname;
2634 
2635  switch (type) {
2636  case kD0toKpi:
2637  case kD0toKpiLikeSign:
2638  cutsname = "D0toKpiCuts";
2639  break;
2640  case kDstartoKpipi:
2641  cutsname = "DStartoKpipiCuts";
2642  break;
2643  default:
2644  return 0;
2645  }
2646 
2647  if (!cuttype.IsNull()) {
2648  cutsname += TString::Format("_%s", cuttype.Data());
2649  }
2650 
2651  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
2652  if (cuts) cuts->PrintAll();
2653  }
2654 
2655  AnalysisEngine eng(type, MCmode, cuts);
2656 
2657  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
2658 
2659  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
2660  eng.AddJetDefinition(jetDef);
2661  it = fAnalysisEngines.insert(it, eng);
2662  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
2663  }
2664  else {
2665  AnalysisEngine* found_eng = &(*it);
2666  ::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());
2667  found_eng->AddJetDefinition(jetDef);
2668 
2669  if (cuts) {
2670  if (found_eng->fRDHFCuts != 0) ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
2671  found_eng->SetRDHFCuts(cuts);
2672  }
2673  }
2674 
2675  return &(*it);
2676 }
2677 
2678 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
2679 {
2680  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
2681  while (it != fAnalysisEngines.end() && (*it) != eng) it++;
2682  return it;
2683 }
2684 
2687 {
2688  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
2689 
2691 
2692  // Define histograms
2693  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
2694 
2695  TString hname;
2696  TString htitle;
2697  TH1* h = 0;
2698  Int_t treeSlot = 0;
2699 
2700  Int_t maxTracks = 6000;
2701  Double_t maxRho = 500;
2702  if (fForceBeamType == kpp) {
2703  maxRho = 50;
2704  maxTracks = 200;
2705  }
2706  else if (fForceBeamType == kpA) {
2707  maxRho = 200;
2708  maxTracks = 500;
2709  }
2710 
2711  hname = "fHistCharmPt";
2712  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2713  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2714 
2715  hname = "fHistCharmEta";
2716  htitle = hname + ";#eta_{charm};counts";
2717  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2718 
2719  hname = "fHistCharmPhi";
2720  htitle = hname + ";#phi_{charm};counts";
2721  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2722 
2723  hname = "fHistCharmPt_Eta05";
2724  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2725  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2726 
2727  hname = "fHistBottomPt";
2728  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2729  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2730 
2731  hname = "fHistBottomEta";
2732  htitle = hname + ";#eta_{bottom};counts";
2733  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2734 
2735  hname = "fHistBottomPhi";
2736  htitle = hname + ";#phi_{bottom};counts";
2737  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2738 
2739  hname = "fHistBottomPt_Eta05";
2740  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2741  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2742 
2743  hname = "fHistHighestPartonPt";
2744  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2745  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2746 
2747  hname = "fHistHighestPartonType";
2748  htitle = hname + ";type;counts";
2749  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2750 
2751  hname = "fHistNHeavyQuarks";
2752  htitle = hname + ";number of heavy-quarks;counts";
2753  fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
2754 
2755  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
2756  for (auto &param : fAnalysisEngines) {
2757  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
2758 
2759  param.fHistManager = &fHistManager;
2760 
2761  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", param.GetName());
2762  htitle = hname + ";#it{N}_{tracks};#it{N}_{D};events";
2763  h = fHistManager.CreateTH2(hname, htitle, 251, -0.5, 250.5, 21, -0.5, 20.5);
2764 
2765  hname = TString::Format("%s/fHistNTotAcceptedDmesons", param.GetName());
2766  htitle = hname + ";;#it{N}_{D}";
2767  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
2768 
2769  hname = TString::Format("%s/fHistNDmesons", param.GetName());
2770  htitle = hname + ";#it{N}_{D};events";
2771  h = fHistManager.CreateTH1(hname, htitle, 501, -0.5, 500.5);
2772 
2773  hname = TString::Format("%s/fHistNEvents", param.GetName());
2774  htitle = hname + ";Event status;counts";
2775  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
2776  h->GetXaxis()->SetBinLabel(1, "Accepted");
2777  h->GetXaxis()->SetBinLabel(2, "Rejected");
2778 
2779  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
2780  htitle = hname + ";Rejection reason;counts";
2781  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
2782 
2783  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
2784  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
2785  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2786 
2787  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
2788  htitle = hname + ";#it{#eta}_{D};counts";
2789  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2790 
2791  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
2792  htitle = hname + ";#it{#phi}_{D};counts";
2793  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2794 
2795  if (param.fMCMode != kMCTruth) {
2796  if (param.fCandidateType == kD0toKpi || param.fCandidateType == kD0toKpiLikeSign) {
2797  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
2798  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2799  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
2800  }
2801  else if (param.fCandidateType == kDstartoKpipi) {
2802  Double_t min = 0;
2803  Double_t max = 0;
2804 
2805  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
2806  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2807  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
2808  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
2809 
2810  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2811  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
2812  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2813  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
2814  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
2815  }
2816  }
2817 
2818  if (param.fMCMode == kMCTruth) {
2819  hname = TString::Format("%s/fHistPartonPt", param.GetName());
2820  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
2821  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2822 
2823  hname = TString::Format("%s/fHistPartonEta", param.GetName());
2824  htitle = hname + ";#eta_{parton};counts";
2825  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2826 
2827  hname = TString::Format("%s/fHistPartonPhi", param.GetName());
2828  htitle = hname + ";#phi_{parton};counts";
2829  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2830 
2831  hname = TString::Format("%s/fHistPartonType", param.GetName());
2832  htitle = hname + ";type;counts";
2833  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2834 
2835  hname = TString::Format("%s/fHistPrompt", param.GetName());
2836  htitle = hname + ";Type;counts";
2837  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
2838  h->GetXaxis()->SetBinLabel(1, "Unknown");
2839  h->GetXaxis()->SetBinLabel(2, "Prompt");
2840  h->GetXaxis()->SetBinLabel(3, "Non-Prompt");
2841 
2842  hname = TString::Format("%s/fHistAncestor", param.GetName());
2843  htitle = hname + ";Ancestor;counts";
2844  h = fHistManager.CreateTH1(hname, htitle, 4, 0, 4);
2845  h->GetXaxis()->SetBinLabel(1, "Unknown");
2846  h->GetXaxis()->SetBinLabel(2, "Charm");
2847  h->GetXaxis()->SetBinLabel(3, "Bottom");
2848  h->GetXaxis()->SetBinLabel(4, "Proton");
2849  }
2850 
2851  for (auto& jetDef : param.fJetDefinitions) {
2852  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef.GetName());
2853 
2854  if (param.fMCMode == kMCTruth) {
2855  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", param.GetName(), jetDef.GetName());
2856  htitle = hname + ";#it{N}_{constituents};#it{N}_{D};counts";
2857  h = fHistManager.CreateTH2(hname, htitle, 51, -0.5, 50.5, 10, 0.5, 10.5);
2858  }
2859 
2860  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef.GetName());
2861  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2862  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2863  SetRejectionReasonLabels(h->GetXaxis());
2864 
2865  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef.GetName());
2866  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2867  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2868  SetRejectionReasonLabels(h->GetXaxis());
2869 
2870  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef.GetName());
2871  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
2872  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2873  SetRejectionReasonLabels(h->GetXaxis());
2874 
2875  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef.GetName());
2876  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
2877  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
2878 
2879  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef.GetName());
2880  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
2881  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2882 
2883  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef.GetName());
2884  htitle = hname + ";#it{#eta}_{jet};counts";
2885  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2886 
2887  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef.GetName());
2888  htitle = hname + ";#it{#phi}_{jet};counts";
2889  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2890 
2891  if (!jetDef.fRhoName.IsNull()) {
2892  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", param.GetName(), jetDef.GetName());
2893  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});#rho (GeV/#it{c} #times rad^{-1});counts";
2894  fHistManager.CreateTH2(hname, htitle, 300, 0, 150, 1000, 0, maxRho);
2895 
2896  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", param.GetName(), jetDef.GetName());
2897  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});#rho (GeV/#it{c} #times rad^{-1});counts";
2898  fHistManager.CreateTH2(hname, htitle, 300, 0, 150, 1000, 0, maxRho);
2899 
2900  hname = TString::Format("%s/%s/fHistRhoVsCent", param.GetName(), jetDef.GetName());
2901  htitle = hname + ";Centrality (%);#rho (GeV/#it{c} #times rad^{-1});counts";
2902  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 1000, 0, maxRho);
2903 
2904  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", param.GetName(), jetDef.GetName());
2905  htitle = hname + ";Centrality (%);#it{p}_{T,jet} (GeV/#it{c});counts";
2906  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
2907 
2908  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", param.GetName(), jetDef.GetName());
2909  htitle = hname + ";Centrality (%);#it{p}_{T,D} (GeV/#it{c});counts";
2910  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
2911 
2912  hname = TString::Format("%s/%s/fHistRhoVsNTracks", param.GetName(), jetDef.GetName());
2913  htitle = hname + ";no. of tracks;#rho (GeV/#it{c} #times rad^{-1});counts";
2914  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 1000, 0, maxRho);
2915 
2916  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", param.GetName(), jetDef.GetName());
2917  htitle = hname + ";no. of tracks;#it{p}_{T,jet} (GeV/#it{c});counts";
2918  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
2919 
2920  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", param.GetName(), jetDef.GetName());
2921  htitle = hname + ";no. of tracks;#it{p}_{T,D} (GeV/#it{c});counts";
2922  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
2923  }
2924  }
2925  switch (fOutputType) {
2926  case kTreeOutput:
2927  param.BuildTree(GetName());
2928  if (treeSlot < fNOutputTrees) {
2929  param.AssignDataSlot(treeSlot+2);
2930  treeSlot++;
2932  }
2933  else {
2934  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
2935  }
2936  break;
2937  case kTHnOutput:
2938  param.BuildHnSparse(fEnabledAxis);
2939  break;
2940  case kNoOutput:
2941  break;
2942  }
2943  }
2944 
2946 
2947  PostData(1, fOutput);
2948 }
2949 
2953 {
2955 
2956  // Load the event
2957  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
2958 
2959  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
2960 
2961  fFastJetWrapper->SetAreaType((fastjet::AreaType)fJetAreaType);
2963 
2964  if (!fAodEvent) {
2965  AliError(Form("This task need an AOD event (Task '%s'). Expect troubles...", GetName()));
2966  //return;
2967  }
2968 
2969  TRandom* rnd = 0;
2970  if (fTrackEfficiency > 0 && fTrackEfficiency < 1) rnd = new TRandom3(0);
2971 
2972  for (auto cont_it : fParticleCollArray) {
2973  AliHFAODMCParticleContainer* part_cont = dynamic_cast<AliHFAODMCParticleContainer*>(cont_it.second);
2974  if (part_cont) fMCContainer = part_cont;
2975  }
2976 
2977  for (auto &params : fAnalysisEngines) {
2978 
2979  params.fAodEvent = fAodEvent;
2980  params.fFastJetWrapper = fFastJetWrapper;
2981  params.fTrackEfficiency = fTrackEfficiency;
2982  params.fRejectISR = fRejectISR;
2983  params.fRandomGen = rnd;
2984 
2985  for (auto &jetdef: params.fJetDefinitions) {
2986  if (!jetdef.fRhoName.IsNull()) {
2987  jetdef.fRho = dynamic_cast<AliRhoParameter*>(fInputEvent->FindListObject(jetdef.fRhoName));
2988  if (!jetdef.fRho) {
2989  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2990  "%s: Could not find rho object '%s' for engine '%s'",
2991  GetName(), jetdef.fRhoName.Data(), params.GetName());
2992  }
2993  }
2994  }
2995 
2996  if (!params.fRDHFCuts) {
2997  if (params.fMCMode == kMCTruth) {
2998  ::Warning("AliAnalysisTaskDmesonJets::ExecOnce",
2999  "%s: RDHF cuts not provided for engine '%s'.",
3000  GetName(), params.GetName());
3001  }
3002  else {
3003  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3004  "%s: RDHF cuts not provided. Engine '%s' disabled.",
3005  GetName(), params.GetName());
3006  params.fInhibit = kTRUE;
3007  }
3008  }
3009 
3010  params.fMCContainer = fMCContainer;
3011 
3012  for (auto cont_it : fParticleCollArray) {
3013  AliTrackContainer* track_cont = dynamic_cast<AliTrackContainer*>(cont_it.second);
3014  if (track_cont) params.fTrackContainers.push_back(track_cont);
3015  }
3016 
3017  for (auto cont_it : fClusterCollArray) params.fClusterContainers.push_back(cont_it.second);
3018 
3019  if (fAodEvent) params.Init(fGeom, fAodEvent->GetRunNumber());
3020 
3021  if (params.fMCMode != kMCTruth && fAodEvent) {
3022  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
3023 
3024  if (params.fCandidateArray) {
3025  TString className;
3026  if (params.fCandidateType == kD0toKpi || params.fCandidateType == kD0toKpiLikeSign) {
3027  className = "AliAODRecoDecayHF2Prong";
3028  }
3029  else if (params.fCandidateType == kDstartoKpipi) {
3030  className = "AliAODRecoCascadeHF";
3031  }
3032  if (!params.fCandidateArray->GetClass()->InheritsFrom(className)) {
3033  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3034  "%s: Objects of type %s in %s are not inherited from %s! Task will be disabled!",
3035  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName(), className.Data());
3036  params.fCandidateArray = 0;
3037  params.fInhibit = kTRUE;
3038  }
3039  }
3040  else {
3041  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3042  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
3043  params.fBranchName.Data(), params.GetName());
3044  params.fInhibit = kTRUE;
3045  }
3046  }
3047 
3048  if (params.fMCMode != kNoMC) {
3049  if (!params.fMCContainer) {
3050  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3051  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
3052  params.GetName());
3053  params.fInhibit = kTRUE;
3054  }
3055  }
3056 
3057  if (params.fMCMode != kMCTruth) {
3058  if (params.fTrackContainers.size() == 0 && params.fClusterContainers.size() == 0) {
3059  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3060  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
3061  params.GetName());
3062  params.fInhibit = kTRUE;
3063  }
3064  }
3065  }
3066 }
3067 
3072 {
3073  TString hname;
3074 
3075  // Fix for temporary bug in ESDfilter
3076  // The AODs with null vertex pointer didn't pass the PhysSel
3077  // Now adding an entry in the histogram so as to check that this is actually cutting anything out
3078  if (fAodEvent && (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001)) {
3079  for (auto &eng : fAnalysisEngines) {
3080  if (eng.fInhibit) continue;
3081  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
3082  fHistManager.FillTH1(hname, "ESDfilterBug");
3083  }
3084  return kFALSE;
3085  }
3086 
3087  if (fAodEvent) {
3088  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
3089  if (matchingAODdeltaAODlevel <= 0) {
3090  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
3091  for (auto &eng : fAnalysisEngines) {
3092  if (eng.fInhibit) continue;
3093  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
3094  fHistManager.FillTH1(hname, "MismatchDeltaAOD");
3095  }
3096  return kFALSE;
3097  }
3098  }
3099 
3100  for (auto &eng : fAnalysisEngines) {
3101  eng.fDmesonJets.clear();
3102  if (eng.fInhibit) continue;
3103 
3104  eng.fCent = fCent;
3105 
3106  //Event selection
3107  hname = TString::Format("%s/fHistNEvents", eng.GetName());
3108  if (fAodEvent) {
3109  Bool_t iseventselected = kTRUE;
3110  if (eng.fRDHFCuts) iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
3111  if (!iseventselected) {
3112  fHistManager.FillTH1(hname, "Rejected");
3113  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
3114  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
3115  TString label;
3116  do {
3117  label = GetHFEventRejectionReasonLabel(bitmap);
3118  if (label.IsNull()) break;
3119  fHistManager.FillTH1(hname, label);
3120  } while (true);
3121 
3122  continue;
3123  }
3124  }
3125 
3126  fHistManager.FillTH1(hname, "Accepted");
3127 
3128  AliDebug(2, "Event selected");
3129 
3130  eng.RunAnalysis();
3131  }
3132  return kTRUE;
3133 }
3134 
3139 {
3140  for (auto &param : fAnalysisEngines) {
3141  if (param.fInhibit) continue;
3142 
3143  if (fOutputType == kTreeOutput) {
3144  param.FillTree(fApplyKinematicCuts);
3146  }
3147  else if (fOutputType == kTHnOutput) {
3148  param.FillHnSparse(fApplyKinematicCuts);
3149  }
3150  }
3152  return kTRUE;
3153 }
3154 
3157 {
3158  auto itcont = fMCContainer->all_momentum();
3159  Int_t nHQ = 0;
3160  Double_t highestPartonPt = 0;
3161  Int_t absPdgHighParton = 0;
3162  for (auto part : itcont) {
3163  Int_t absPdgCode = TMath::Abs(part.second->GetPdgCode());
3164 
3165  // Skip all particles that are not either quarks or gluons
3166  if (absPdgCode > 9 && absPdgCode != 21) continue;
3167 
3168  // Look for highest momentum parton
3169  if (highestPartonPt < part.first.Pt()) {
3170  highestPartonPt = part.first.Pt();
3171  absPdgHighParton = absPdgCode;
3172  }
3173  /*
3174  // Look for the mother PDG code
3175  Int_t motherIndex = part.second->GetMother();
3176  AliAODMCParticle *mother = 0;
3177  Int_t motherPdg = 0;
3178  Double_t motherPt = 0;
3179  if (motherIndex >= 0) {
3180  mother = fMCContainer->GetMCParticle(motherIndex);
3181  if (motherIndex) {
3182  motherPdg = TMath::Abs(mother->GetPdgCode());
3183  motherPt = mother->Pt();
3184  }
3185  }
3186  */
3187  if (absPdgCode != 4 && absPdgCode != 5) continue;
3188  Bool_t notLastInPartonShower = kFALSE;
3189  for (Int_t idaugh = 0; idaugh < 2; idaugh++){
3190  Int_t daughterIndex = part.second->GetDaughter(idaugh);
3191  if (daughterIndex < 0) {
3192  AliDebug(10, Form("Could not find daughter of heavy quark (pdg=%d, pt=%.3f)!", absPdgCode, part.first.Pt()));
3193  continue;
3194  }
3195  AliAODMCParticle *daughter = fMCContainer->GetMCParticle(daughterIndex);
3196  if (!daughter) {
3197  AliDebug(10, Form("Could not find daughter %d of heavy quark (pdg=%d, pt=%.3f)!", daughterIndex, absPdgCode, part.first.Pt()));
3198  continue;
3199  }
3200  Int_t daughterAbsPdgCode = TMath::Abs(daughter->GetPdgCode());
3201  if (daughterAbsPdgCode <= 9 || daughterAbsPdgCode == 21) notLastInPartonShower = kTRUE; // this parton is not the last parton in the shower
3202  AliDebug(10, Form("Found daughter with PDG=%d, pt=%.3f", daughterAbsPdgCode, daughter->Pt()));
3203  }
3204  if (notLastInPartonShower) continue;
3205 
3206  if (absPdgCode == 4) {
3207  fHistManager.FillTH1("fHistCharmPt", part.first.Pt());
3208  fHistManager.FillTH1("fHistCharmEta", part.first.Eta());
3209  fHistManager.FillTH1("fHistCharmPhi", part.first.Phi_0_2pi());
3210  if (TMath::Abs(part.first.Eta()) < 0.5) fHistManager.FillTH1("fHistCharmPt_Eta05", part.first.Pt());
3211  }
3212  else if (absPdgCode == 5) {
3213  fHistManager.FillTH1("fHistBottomPt", part.first.Pt());
3214  fHistManager.FillTH1("fHistBottomEta", part.first.Eta());
3215  fHistManager.FillTH1("fHistBottomPhi", part.first.Phi_0_2pi());
3216  if (TMath::Abs(part.first.Eta()) < 0.5) fHistManager.FillTH1("fHistBottomPt_Eta05", part.first.Pt());
3217  }
3218  nHQ++;
3219  }
3220  fHistManager.FillTH1("fHistNHeavyQuarks", nHQ);
3221  fHistManager.FillTH1("fHistHighestPartonPt",highestPartonPt);
3222  Int_t partonType = 0;
3223  if (absPdgHighParton == 9 || absPdgHighParton == 21) {
3224  partonType = 7;
3225  }
3226  else {
3227  partonType = absPdgHighParton;
3228  }
3229  fHistManager.FillTH1("fHistHighestPartonType",partonType);
3230 }
3231 
3240 {
3241  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
3242 
3243  Double_t mass = part->Mass();
3244 
3245  // To make sure that the PDG mass value is not at the edge of a bin
3246  if (nbins % 2 == 0) {
3247  minMass = mass - range / 2 - range / nbins / 2;
3248  maxMass = mass + range / 2 - range / nbins / 2;
3249  }
3250  else {
3251  minMass = mass - range / 2;
3252  maxMass = mass + range / 2;
3253  }
3254 }
3255 
3262 {
3263  static TString label;
3264  label = "";
3265 
3266  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
3267  label = "NotSelTrigger";
3268  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
3269  return label.Data();
3270  }
3271  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
3272  label = "NoVertex";
3273  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
3274  return label.Data();
3275  }
3276  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
3277  label = "TooFewVtxContrib";
3278  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
3279  return label.Data();
3280  }
3281  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
3282  label = "ZVtxOutFid";
3283  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
3284  return label.Data();
3285  }
3286  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
3287  label = "Pileup";
3288  bitmap &= ~BIT(AliRDHFCuts::kPileup);
3289  return label.Data();
3290  }
3291  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
3292  label = "OutsideCentrality";
3293  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
3294  return label.Data();
3295  }
3296  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
3297  label = "PhysicsSelection";
3298  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
3299  return label.Data();
3300  }
3301  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
3302  label = "BadSPDVertex";
3303  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
3304  return label.Data();
3305  }
3306  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
3307  label = "ZVtxSPDOutFid";
3308  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
3309  return label.Data();
3310  }
3311  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
3312  label = "CentralityFlattening";
3313  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
3314  return label.Data();
3315  }
3316  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
3317  label = "BadTrackV0Correl";
3318  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
3319  return label.Data();
3320  }
3321 
3322  return label.Data();
3323 }
3324 
3331 {
3332  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
3333  PostData(eng.GetDataSlotNumber(), eng.GetTree());
3334  return eng.GetDataSlotNumber();
3335  }
3336  else {
3337  return -1;
3338  }
3339 }
3340 
3350 {
3351  // Get the pointer to the existing analysis manager via the static access method
3352  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
3353  if (!mgr) {
3354  ::Error("AddTaskDmesonJets", "No analysis manager to connect to.");
3355  return NULL;
3356  }
3357 
3358  // Check the analysis type using the event handlers connected to the analysis manager
3359  AliVEventHandler* handler = mgr->GetInputEventHandler();
3360  if (!handler) {
3361  ::Error("AddTaskEmcalJetSpectraQA", "This task requires an input event handler");
3362  return NULL;
3363  }
3364 
3365  EDataType_t dataType = kUnknownDataType;
3366 
3367  if (handler->InheritsFrom("AliESDInputHandler")) {
3368  dataType = kESD;
3369  }
3370  else if (handler->InheritsFrom("AliAODInputHandler")) {
3371  dataType = kAOD;
3372  }
3373 
3374  // Init the task and do settings
3375  if (ntracks == "usedefault") {
3376  if (dataType == kESD) {
3377  ntracks = "Tracks";
3378  }
3379  else if (dataType == kAOD) {
3380  ntracks = "tracks";
3381  }
3382  else {
3383  ntracks = "";
3384  }
3385  }
3386 
3387  if (nclusters == "usedefault") {
3388  if (dataType == kESD) {
3389  nclusters = "CaloClusters";
3390  }
3391  else if (dataType == kAOD) {
3392  nclusters = "caloClusters";
3393  }
3394  else {
3395  nclusters = "";
3396  }
3397  }
3398 
3399  if (nMCpart == "usedefault") {
3400  nMCpart = "mcparticles"; // Always needs AliAODMCParticle objects
3401  }
3402 
3403  TString name("AliAnalysisTaskDmesonJets");
3404  if (strcmp(suffix, "") != 0) {
3405  name += TString::Format("_%s", suffix.Data());
3406  }
3407 
3408  AliAnalysisTaskDmesonJets* jetTask = new AliAnalysisTaskDmesonJets(name, nMaxTrees);
3409 
3410  if (!ntracks.IsNull()) {
3411  AliHFTrackContainer* trackCont = new AliHFTrackContainer(ntracks);
3412  jetTask->AdoptParticleContainer(trackCont);
3413  }
3414 
3415  if (!nMCpart.IsNull()) {
3416  AliMCParticleContainer* partCont = new AliHFAODMCParticleContainer(nMCpart);
3417  partCont->SetEtaLimits(-1.5, 1.5);
3418  partCont->SetPtLimits(0, 1000);
3419  jetTask->AdoptParticleContainer(partCont);
3420  }
3421 
3422  jetTask->AddClusterContainer(nclusters.Data());
3423 
3424  // Final settings, pass to manager and set the containers
3425  mgr->AddTask(jetTask);
3426 
3427  // Create containers for input/output
3428  AliAnalysisDataContainer* cinput1 = mgr->GetCommonInputContainer();
3429  TString contname1(name);
3430  contname1 += "_histos";
3431  AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(contname1.Data(),
3432  TList::Class(), AliAnalysisManager::kOutputContainer,
3433  Form("%s", AliAnalysisManager::GetCommonFileName()));
3434 
3435  mgr->ConnectInput(jetTask, 0, cinput1);
3436  mgr->ConnectOutput(jetTask, 1, coutput1);
3437 
3438  for (Int_t i = 0; i < nMaxTrees; i++) {
3439  TString contname = TString::Format("%s_tree_%d", name.Data(), i);
3440  AliAnalysisDataContainer *coutput = mgr->CreateContainer(contname.Data(),
3441  TTree::Class(),AliAnalysisManager::kOutputContainer,
3442  Form("%s", AliAnalysisManager::GetCommonFileName()));
3443  mgr->ConnectOutput(jetTask, 2+i, coutput);
3444  }
3445  return jetTask;
3446 }
Lightweight class that encapsulates D meson jets.
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.
std::string fClassName
Class name where the event was not found.
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
void Getd0MeasMinusExpProng(Int_t ip, Double_t magf, Double_t &d0diff, Double_t &errd0diff) 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)
virtual void Set(const AliDmesonJetInfo &source)
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
Double_t CosThetaStarD0bar() const
angle of K
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
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.
AliEventNotFound(const std::string &class_name, const std::string &method_name)
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 AliVEvent * fEvent
! pointer to the ESD/AOD event
const Int_t nbins
Double_t maxMass
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
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
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)