AliPhysics  d565ceb (d565ceb)
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 
133  fD(source.fD),
134  fSoftPionPt(source.fSoftPionPt),
136  fJets(source.fJets),
137  fMCLabel(source.fMCLabel),
139  fParton(source.fParton),
140  fPartonType(source.fPartonType),
141  fAncestor(source.fAncestor),
142  fD0D0bar(source.fD0D0bar),
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();
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 {
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),
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 
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),
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),
950  fCurrentJetInfo(0),
951  fCandidateArray(0),
952  fMCContainer(),
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),
972  fNDaughters(source.fNDaughters),
974  fBranchName(source.fBranchName),
975  fMCMode(source.fMCMode),
976  fNMassBins(source.fNMassBins),
977  fMinMass(source.fMinMass),
978  fMaxMass(source.fMaxMass),
979  fRDHFCuts(),
982  fInhibit(source.fInhibit),
984  fPtBinWidth(source.fPtBinWidth),
985  fMaxPt(source.fMaxPt),
986  fD0Extended(source.fD0Extended),
987  fRandomGen(source.fRandomGen),
989  fDataSlotNumber(-1),
990  fTree(0),
992  fCurrentJetInfo(0),
994  fMCContainer(source.fMCContainer),
997  fAodEvent(source.fAodEvent),
999  fHistManager(source.fHistManager),
1000  fCent(-1)
1001 {
1002  SetRDHFCuts(source.fRDHFCuts);
1003 }
1004 
1005 // Destructor
1007 {
1008  delete fRDHFCuts;
1009 }
1010 
1015 {
1016  new (this) AnalysisEngine(source);
1017  return *this;
1018 }
1019 
1024 {
1025  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
1026  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
1027  }
1028 
1029  return kFALSE;
1030 }
1031 
1033 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const /*geom*/, Int_t /*runNumber*/)
1034 {
1035 }
1036 
1041 {
1042  switch (fCandidateType) {
1043  case kD0toKpi:
1044  fCandidatePDG = 421;
1045  fCandidateName = "D0";
1046  fNDaughters = 2;
1048  fPDGdaughters.Reset();
1049  fPDGdaughters[0] = 211; // pi
1050  fPDGdaughters[1] = 321; // K
1051  fBranchName = "D0toKpi";
1053  break;
1054  case kD0toKpiLikeSign:
1055  fCandidatePDG = 421;
1056  fCandidateName = "2ProngLikeSign";
1057  fNDaughters = 2;
1059  fPDGdaughters.Reset();
1060  fPDGdaughters[0] = 211; // pi
1061  fPDGdaughters[1] = 321; // K
1062  fBranchName = "LikeSign2Prong";
1063  break;
1064  case kDstartoKpipi:
1065  fCandidatePDG = 413;
1066  fCandidateName = "DStar";
1067  fNDaughters = 3;
1069  fPDGdaughters.Reset();
1070  fPDGdaughters[0] = 211; // pi soft
1071  fPDGdaughters[1] = 211; // pi fromD0
1072  fPDGdaughters[2] = 321; // K from D0
1073  fBranchName = "Dstar";
1075  break;
1076  default:
1077  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
1078  }
1079 
1081 }
1082 
1087 {
1088  if (fRDHFCuts) delete fRDHFCuts;
1089  fRDHFCuts = cuts;
1090 }
1091 
1096 {
1097  if (!cuts) return;
1098  if (fRDHFCuts) delete fRDHFCuts;
1099  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
1100 }
1101 
1106 {
1107  static TString name;
1108 
1109  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
1110 
1111  return name.Data();
1112 }
1113 
1118 {
1120  switch (fMCMode) {
1121  case kBackgroundOnly:
1122  fName += "_BackgroundOnly";
1123  break;
1124  case kSignalOnly:
1125  fName += "_SignalOnly";
1126  break;
1127  case kMCTruth:
1128  fName += "_MCTruth";
1129  break;
1130  case kD0Reflection:
1131  fName += "_D0Reflection";
1132  break;
1133  case kOnlyWrongPIDAccepted:
1134  fName += "_OnlyWrongPIDAccepted";
1135  break;
1136  default:
1137  break;
1138  }
1139 
1140  if (fRDHFCuts) fName += TString::Format("_%s", fRDHFCuts->GetName());
1141 
1142  return fName.Data();
1143 }
1144 
1152 {
1153  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
1154 
1155  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
1156  fJetDefinitions.push_back(def);
1157  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
1158  "Total number of jet definitions is now %lu.",
1159  def.GetName(), GetName(), fJetDefinitions.size());
1160  // For detector level set maximum track pt to 100 GeV/c
1161  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
1162  }
1163  else {
1164  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
1165  }
1166 
1167  return &(*it);
1168 }
1169 
1181 {
1182  AliHFJetDefinition def(type, r, algo, reco);
1183 
1184  return AddJetDefinition(def);
1185 }
1186 
1192 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
1193 {
1194  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
1195  while (it != fJetDefinitions.end() && (*it) != def) it++;
1196  return it;
1197 }
1198 
1203 {
1204  for (auto &jetdef : fJetDefinitions) jetdef.SetJetPhiRange(min, max);
1205 }
1206 
1211 {
1212  for (auto &jetdef : fJetDefinitions) jetdef.SetJetEtaRange(min, max);
1213 }
1214 
1219 {
1220  for (auto &jetdef : fJetDefinitions) jetdef.SetJetPtRange(min, max);
1221 }
1222 
1227 {
1228  for (auto &jetdef : fJetDefinitions) jetdef.SetChargedPtRange(min, max);
1229 }
1230 
1235 {
1236  for (auto &jetdef : fJetDefinitions) jetdef.SetNeutralPtRange(min, max);
1237 }
1238 
1245 {
1246  if (lhs.fCandidateType < rhs.fCandidateType) {
1247  return true;
1248  }
1249  else if (lhs.fCandidateType > rhs.fCandidateType) {
1250  return false;
1251  }
1252  else if (lhs.fMCMode < rhs.fMCMode) {
1253  return true;
1254  }
1255  else if (lhs.fMCMode > rhs.fMCMode) {
1256  return false;
1257  }
1258  else if (lhs.fRDHFCuts && !rhs.fRDHFCuts) {
1259  return true;
1260  }
1261  else if (lhs.fRDHFCuts && rhs.fRDHFCuts && strcmp(lhs.fRDHFCuts->GetName(), rhs.fRDHFCuts->GetName()) < 0) {
1262  return true;
1263  }
1264  else {
1265  return false;
1266  }
1267 }
1268 
1275 {
1276  if (lhs.fCandidateType != rhs.fCandidateType) return false;
1277  if (lhs.fMCMode != rhs.fMCMode) return false;
1278  if (lhs.fRDHFCuts == nullptr && rhs.fRDHFCuts != nullptr) return false;
1279  if (lhs.fRDHFCuts != nullptr && rhs.fRDHFCuts == nullptr) return false;
1280  if (lhs.fRDHFCuts && rhs.fRDHFCuts && strcmp(lhs.fRDHFCuts->GetName(), rhs.fRDHFCuts->GetName()) != 0) return false;
1281  return true;
1282 }
1283 
1292 {
1293  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) { // D0 candidate
1294  return ExtractD0Attributes(Dcand, DmesonJet, i);
1295  }
1296  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
1297  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
1298  }
1299  else {
1300  return kFALSE;
1301  }
1302 }
1303 
1312 {
1313  AliDebug(10,"Checking if D0 meson is selected");
1314  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
1315  if (isSelected == 0) return kFALSE;
1316 
1317  Int_t MCtruthPdgCode = 0;
1318 
1319  Double_t invMassD = 0;
1320 
1321  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
1322  // Checks also the origin, and if it matches the rejected origin mask, return false
1323  if (fMCMode != kNoMC) {
1324  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
1325  DmesonJet.fMCLabel = mcLab;
1326 
1327  // Retrieve the generated particle (if exists) and its PDG code
1328  if (mcLab >= 0) {
1329  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1330 
1331  if (aodMcPart) {
1332  // Check origin and return false if it matches the rejected origin mask
1333  if (fRejectedOrigin) {
1334  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
1335  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1336  }
1337  MCtruthPdgCode = aodMcPart->PdgCode();
1338  }
1339  }
1340  }
1341 
1342  if (isSelected == 1) { // selected as a D0
1343  if (i != 0) return kFALSE; // only one mass hypothesis thanks to PID
1344 
1345  if (fMCMode == kNoMC ||
1346  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1347  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly) ||
1348  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kD0Reflection || fMCMode == kOnlyWrongPIDAccepted))) {
1349  // 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)
1350  AliDebug(10,"Selected as D0");
1351  invMassD = Dcand->InvMassD0();
1352  }
1353  else { // conditions above not passed, so return FALSE
1354  return kFALSE;
1355  }
1356  }
1357  else if (isSelected == 2) { // selected as a D0bar
1358  if (i != 1) return kFALSE; // only one mass hypothesis thanks to PID
1359 
1360  if (fMCMode == kNoMC ||
1361  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1362  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly) ||
1363  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kD0Reflection || fMCMode == kOnlyWrongPIDAccepted))) {
1364  // 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)
1365  AliDebug(10,"Selected as D0bar");
1366  invMassD = Dcand->InvMassD0bar();
1367  }
1368  else { // conditions above not passed, so return FALSE
1369  return kFALSE;
1370  }
1371  }
1372  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
1373  AliDebug(10,"Selected as either D0 or D0bar");
1374 
1375  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
1376  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1377  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kD0Reflection))) {
1378  if (i != 0) return kFALSE;
1379  AliDebug(10, "MC truth is D0");
1380  invMassD = Dcand->InvMassD0();
1381  }
1382  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1383  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kD0Reflection))) {
1384  if (i != 1) return kFALSE;
1385  AliDebug(10, "MC truth is D0bar");
1386  invMassD = Dcand->InvMassD0bar();
1387  }
1388  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
1389 
1390  // Only accept it if background-only OR background-and-signal was requested
1391  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
1392  // Select D0 or D0bar depending on the i-parameter
1393  if (i == 0) {
1394  AliDebug(10, "Returning invariant mass with D0 hypothesis");
1395  invMassD = Dcand->InvMassD0();
1396  }
1397  else if (i == 1) {
1398  AliDebug(10, "Returning invariant mass with D0bar hypothesis");
1399  invMassD = Dcand->InvMassD0bar();
1400  }
1401  else { // i > 1
1402  return kFALSE;
1403  }
1404  }
1405  else { // signal-only was requested but this is not a true D0
1406  return kFALSE;
1407  }
1408  }
1409  }
1410  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
1411  return kTRUE;
1412 }
1413 
1422 {
1423  AliDebug(10,"Checking if D* meson is selected");
1424  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoCascadeHF*>(DstarCand), AliRDHFCuts::kAll, fAodEvent);
1425  if (isSelected == 0) return kFALSE;
1426 
1427  if ((i == 1 && DstarCand->Charge()>0) || (i == 0 && DstarCand->Charge()<0) || i > 1) return kFALSE; // only one mass hypothesis for the D*
1428 
1429  Int_t MCtruthPdgCode = 0;
1430 
1431  Double_t invMassD = 0;
1432 
1433  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1434  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1435  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1436 
1437  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1438  AliDebug(10, Form("MC label is %d", mcLab));
1439  DmesonJet.fMCLabel = mcLab;
1440  if (mcLab >= 0) {
1441  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1442 
1443  if (aodMcPart) {
1444  if (fRejectedOrigin) {
1445  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
1446  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1447  }
1448 
1449  MCtruthPdgCode = aodMcPart->PdgCode();
1450  AliDebug(10, Form("MC truth pdg code is %d",MCtruthPdgCode));
1451  }
1452  }
1453  }
1454 
1455  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1456  if (fMCMode == kNoMC ||
1457  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1458  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1459  // 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)
1460  invMassD = DstarCand->InvMassDstarKpipi();
1461  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1462  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1463  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1464  return kTRUE;
1465  }
1466  else { // conditions above not passed, so return FALSE
1467  return kFALSE;
1468  }
1469 }
1470 
1478 {
1479  if (!part) return kUnknownDecay;
1480  if (!mcArray) return kUnknownDecay;
1481 
1483 
1484  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1485 
1486  if (part->GetNDaughters() == 2) {
1487 
1488  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1489  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1490 
1491  if (!d1 || !d2) {
1492  return decay;
1493  }
1494 
1495  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1496  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1497 
1498  if (absPdgPart == 421) { // D0 -> K pi
1499  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1500  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1501  decay = kDecayD0toKpi;
1502  }
1503  }
1504 
1505  if (absPdgPart == 413) { // D* -> D0 pi
1506  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1507  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1508  if (D0decay == kDecayD0toKpi) {
1509  decay = kDecayDStartoKpipi;
1510  }
1511  }
1512  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1513  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1514  if (D0decay == kDecayD0toKpi) {
1515  decay = kDecayDStartoKpipi;
1516  }
1517  }
1518  }
1519  }
1520 
1521  return decay;
1522 }
1523 
1530 std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> AliAnalysisTaskDmesonJets::AnalysisEngine::IsPromptCharm(const AliAODMCParticle* part, TClonesArray* mcArray)
1531 {
1532  std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> result(kUnknownQuark, 0);
1533 
1534  if (!part) return result;
1535  if (!mcArray) return result;
1536 
1537  static std::set<UInt_t> partons = { 4, 5 };
1538 
1539  AliAODMCParticle* parton = FindParticleOrigin(part, mcArray, kFindLast, partons);
1540  if (parton) {
1541  result.second = parton;
1542  UInt_t absPdgParton = TMath::Abs(parton->GetPdgCode());
1543  if (absPdgParton == 4) result.first = kFromCharm;
1544  else if (absPdgParton == 5) result.first = kFromBottom;
1545  }
1546 
1547  return result;
1548 }
1549 
1557 
1558 AliAODMCParticle* AliAnalysisTaskDmesonJets::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode)
1559 {
1560  static std::set<UInt_t> pdgSet;
1561 
1562  return FindParticleOrigin(part, mcArray, mode, pdgSet);
1563 }
1564 
1573 
1574 AliAODMCParticle* AliAnalysisTaskDmesonJets::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode, const std::set<UInt_t>& pdgSet)
1575 {
1576  AliAODMCParticle* result = nullptr;
1577 
1578  Int_t mother = part->GetMother();
1579  while (mother >= 0) {
1580  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1581  if (mcGranma) {
1582  UInt_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
1583 
1584  // If the current particle is one of the particle types that is being searched assign it to the result pointer
1585  if (pdgSet.empty() || pdgSet.count(abspdgGranma) > 0) {
1586  result = mcGranma;
1587  // If the last particle in the fragmentation tree (first when going reverse) was requested then stop the loop
1588  if (mode == kFindLast) break;
1589  }
1590  if (mother == mcGranma->GetMother()) { // avoid infinite loop!
1591  AliWarningClassStream() << "Particle " << mother << " (PDG=" << mcGranma->PdgCode() << ") is the mother of itself!?" << std::endl;
1592  break;
1593  }
1594  mother = mcGranma->GetMother();
1595  }
1596  else {
1597  AliErrorClassStream() << "Could not retrieve mother particle " << mother << "!" << std::endl;
1598  break;
1599  }
1600  }
1601 
1602  return result;
1603 }
1604 
1607 {
1608  for (auto& jetDef : fJetDefinitions) {
1609  jetDef.fJets.clear();
1610  }
1611 
1612  if (fMCMode == kMCTruth) {
1614  }
1615  else {
1617  }
1618 }
1619 
1622 {
1623  // Fill the vertex info of the candidates
1624  // Needed for reduced delta AOD, where the vertex info has been deleted
1625  // to reduce the delta AOD file size
1627 
1628  const Int_t nD = fCandidateArray->GetEntriesFast();
1629 
1630  AliDmesonJetInfo DmesonJet;
1631  DmesonJet.fEvent = this->fAodEvent;
1632 
1633  std::map<AliHFJetDefinition*,Double_t> maxJetPt;
1634  for (auto& def : fJetDefinitions) maxJetPt[&def] = 0;
1635  Double_t maxDPt = 0;
1636 
1637  Int_t nAccCharm[3] = {0};
1638  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1639  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1640  if (!charmCand) continue;
1641  if(!(vHF.FillRecoCand(fAodEvent,charmCand))) continue;
1642 
1643  // region of interest + cuts
1644  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1645  Int_t nMassHypo = 0; // number of mass hypothesis accepted for this D meson
1646  if (charmCand->Pt() > maxDPt) maxDPt = charmCand->Pt();
1647  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1648  DmesonJet.Reset();
1649  DmesonJet.fDmesonParticle = charmCand;
1650  DmesonJet.fSelectionType = im + 1;
1651  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1652  for (auto& def : fJetDefinitions) {
1653  if (FindJet(charmCand, DmesonJet, def)) {
1654  Double_t jetPt = DmesonJet.fJets[def.GetName()].fMomentum.Pt();
1655  if (jetPt > maxJetPt[&def]) maxJetPt[&def] = jetPt;
1656  }
1657  else {
1658  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1659  def.GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1660  }
1661  }
1662  fDmesonJets[(icharm+1)*(1-(im*2))] = DmesonJet;
1663  nMassHypo++;
1664  nAccCharm[im]++;
1665  }
1666  }
1667  if (nMassHypo == 2) {
1668  nAccCharm[0]--;
1669  nAccCharm[1]--;
1670  nAccCharm[2] += 2;
1671  }
1672  if (nMassHypo == 2) { // both mass hypothesis accepted
1673  fDmesonJets[(icharm+1)].fD0D0bar = kTRUE;
1674  fDmesonJets[-(icharm+1)].fD0D0bar = kTRUE;
1675  }
1676  } // end of D cand loop
1677 
1678  TString hname;
1679 
1680  Int_t ntracks = 0;
1681 
1682  for (auto track_cont : fTrackContainers) {
1683  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
1684  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(nullptr);
1685  ntracks += track_cont->GetNAcceptEntries();
1686  }
1687 
1688  for (auto& def : fJetDefinitions) {
1689  if (!def.fRho) continue;
1690  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
1691  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
1692 
1693  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
1694  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
1695 
1696  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
1697  fHistManager->FillTH2(hname, fCent, def.fRho->GetVal());
1698 
1699  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
1700  fHistManager->FillTH2(hname, fCent, maxJetPt[&def]);
1701 
1702  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
1703  fHistManager->FillTH2(hname, fCent, maxDPt);
1704 
1705  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
1706  fHistManager->FillTH2(hname, ntracks, def.fRho->GetVal());
1707 
1708  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
1709  fHistManager->FillTH2(hname, ntracks, maxJetPt[&def]);
1710 
1711  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
1712  fHistManager->FillTH2(hname, ntracks, maxDPt);
1713  }
1714 
1715  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1716  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1717  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1718  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1719 
1720  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1721  fHistManager->FillTH2(hname, ntracks, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1722 
1723  hname = TString::Format("%s/fHistNDmesons", GetName());
1724  fHistManager->FillTH1(hname, nD);
1725 }
1726 
1737 {
1738  TString hname;
1739 
1740  Double_t rho = 0;
1741  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
1742 
1744  fFastJetWrapper->SetR(jetDef.fRadius);
1747 
1748  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1749 
1750  if (jetDef.fJetType != AliJetContainer::kNeutralJet) {
1751  for (auto track_cont : fTrackContainers) {
1752  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
1753  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(Dcand);
1754  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1755  AddInputVectors(track_cont, 100, static_cast<TH2*>(fHistManager->FindObject(hname)), fTrackEfficiency);
1756 
1757  if (hftrack_cont) {
1758  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1759  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1760  const TObjArray& daughters = hftrack_cont->GetDaughterList();
1761  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1762  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1763  if (!hftrack_cont->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1764  }
1765  }
1766  }
1767  }
1768 
1769  if (jetDef.fJetType != AliJetContainer::kChargedJet) {
1770  for (auto clus_cont : fClusterContainers) {
1771  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1772  AddInputVectors(clus_cont, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1773  }
1774  }
1775 
1776  // run jet finder
1777  fFastJetWrapper->Run();
1778 
1779  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1780 
1781  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1782  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1783 
1784  Bool_t isDmesonJet = kFALSE;
1785 
1786  Double_t maxChPt = 0;
1787  Double_t maxNePt = 0;
1788  Double_t totalNeutralPt = 0;
1789  Int_t nConst = 1;
1790 
1791  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1792  if (constituents[ic].user_index() == 0) {
1793  isDmesonJet = kTRUE;
1794  }
1795  else if (constituents[ic].user_index() >= 100) {
1796  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1797  nConst++;
1798  }
1799  else if (constituents[ic].user_index() <= -100) {
1800  totalNeutralPt += constituents[ic].pt();
1801  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1802  nConst++;
1803  }
1804  }
1805 
1806  if (isDmesonJet) {
1807  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1808  DmesonJet.fJets[jetDef.GetName()].fNConstituents = nConst;
1809  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1810  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1811  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1812  DmesonJet.fJets[jetDef.GetName()].fArea = jets_incl[ijet].area();
1813  DmesonJet.fJets[jetDef.GetName()].fCorrPt = DmesonJet.fJets[jetDef.GetName()].fMomentum.Pt() - jets_incl[ijet].area() * rho;
1814 
1815  return kTRUE;
1816  }
1817  }
1818 
1819  return kFALSE;
1820 }
1821 
1825 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist, Double_t eff)
1826 {
1827  auto itcont = cont->all_momentum();
1828  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1829  UInt_t rejectionReason = 0;
1830  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1831  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1832  continue;
1833  }
1834  if (fRandomGen && eff > 0 && eff < 1) {
1835  Double_t rnd = fRandomGen->Rndm();
1836  if (eff < rnd) {
1837  if (rejectHist) rejectHist->Fill(6, it->first.Pt());
1838  continue;
1839  }
1840  }
1841  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1842  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1843  }
1844 }
1845 
1848 {
1849  TString hname;
1850 
1855 
1856  if (!fMCContainer->IsSpecialPDGFound()) return;
1857 
1858  Int_t nAccCharm[3] = {0};
1859 
1860  std::map<AliHFJetDefinition*, Double_t> maxJetPt;
1861  Double_t maxDPt = 0;
1862 
1863  for (auto &jetDef : fJetDefinitions) {
1864  maxJetPt[&jetDef] = 0;
1865  Double_t rho = 0;
1866  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
1867  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", GetName(), jetDef.GetName());
1868  TH1* histNDmesonsVsNconstituents = static_cast<TH1*>(fHistManager->FindObject(hname));
1869 
1870  switch (jetDef.fJetType) {
1872  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
1873  break;
1875  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
1876  break;
1878  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
1879  break;
1880  }
1881 
1883  fFastJetWrapper->SetR(jetDef.fRadius);
1886 
1887  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1888  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1889 
1890  fFastJetWrapper->Run();
1891 
1892  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1893 
1894  for (auto jet : jets_incl) {
1895  Int_t nDmesonsInJet = 0;
1896 
1897  for (auto constituent : jet.constituents()) {
1898  Int_t iPart = constituent.user_index() - 100;
1899  if (constituent.perp() < 1e-6) continue; // reject ghost particles
1900  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
1901  if (!part) {
1902  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1903  continue;
1904  }
1905  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1906  nDmesonsInJet++;
1907  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1908  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1909  if (part->Pt() > maxDPt) maxDPt = part->Pt();
1910  std::pair<int, AliDmesonJetInfo> element;
1911  element.first = iPart;
1912  dMesonJetIt = fDmesonJets.insert(element).first;
1913  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1914  (*dMesonJetIt).second.fDmesonParticle = part;
1915  (*dMesonJetIt).second.fSelectionType = part->PdgCode() > 0 ? 1 : 2;
1916 
1917  UShort_t p = 0;
1918  UInt_t rs = 0;
1919 
1920  auto origin = IsPromptCharm(part, fMCContainer->GetArray());
1921  p = 0;
1922  rs = origin.first;
1923  while (rs >>= 1) { p++; }
1924  (*dMesonJetIt).second.fPartonType = p;
1925  (*dMesonJetIt).second.fParton = origin.second;
1926 
1927  (*dMesonJetIt).second.fAncestor = FindParticleOrigin(part, fMCContainer->GetArray(), kFindFirst);
1928 
1929  if (part->PdgCode() > 0) { // D0
1930  nAccCharm[0]++;
1931  }
1932  else { // D0bar
1933  nAccCharm[1]++;
1934  }
1935  }
1936 
1937  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
1938  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = jet.constituents().size();
1939  (*dMesonJetIt).second.fJets[jetDef.GetName()].fArea = jet.area();
1940  (*dMesonJetIt).second.fJets[jetDef.GetName()].fCorrPt = (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.Pt() - jet.area() * rho;
1941  if (jet.perp() > maxJetPt[&jetDef]) maxJetPt[&jetDef] = jet.perp();
1942  } // if constituent is a D meson
1943  } // for each constituent
1944  if (nDmesonsInJet > 0) histNDmesonsVsNconstituents->Fill(jet.constituents().size(), nDmesonsInJet);
1945  } // for each jet
1946  } // for each jet definition
1947 
1949 
1950  for (auto& def : fJetDefinitions) {
1951  if (!def.fRho) continue;
1952  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
1953  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
1954 
1955  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
1956  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
1957 
1958  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
1959  fHistManager->FillTH2(hname, fCent, def.fRho->GetVal());
1960 
1961  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
1962  fHistManager->FillTH2(hname, fCent, maxJetPt[&def]);
1963 
1964  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
1965  fHistManager->FillTH2(hname, fCent, maxDPt);
1966 
1967  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
1968  fHistManager->FillTH2(hname, npart, def.fRho->GetVal());
1969 
1970  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
1971  fHistManager->FillTH2(hname, npart, maxJetPt[&def]);
1972 
1973  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
1974  fHistManager->FillTH2(hname, npart, maxDPt);
1975  }
1976 
1977  if (fDmesonJets.size() != nAccCharm[0]+nAccCharm[1]) AliError(Form("I found %lu mesons (%d)?", fDmesonJets.size(), nAccCharm[0]+nAccCharm[1]));
1978  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1979  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1980  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1981  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1982 
1983  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1984  fHistManager->FillTH2(hname, npart, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1985 
1986  hname = TString::Format("%s/fHistNDmesons", GetName());
1987  fHistManager->FillTH1(hname, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]); // same as the number of accepted D mesons, since no selection is performed
1988 }
1989 
1994 {
1995  TString classname;
1996  if (fMCMode == kMCTruth) {
1997  classname = "AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary";
1999  }
2000  else {
2001  switch (fCandidateType) {
2002  case kD0toKpi:
2003  case kD0toKpiLikeSign:
2004  if (fD0Extended) {
2005  classname = "AliAnalysisTaskDmesonJets::AliD0ExtendedInfoSummary";
2007  }
2008  else {
2009  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
2011  }
2012  break;
2013  case kDstartoKpipi:
2014  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
2016  break;
2017  }
2018  }
2019  TString treeName = TString::Format("%s_%s", taskName, GetName());
2020  fTree = new TTree(treeName, treeName);
2021  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
2023  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
2024  if (fJetDefinitions[i].fRhoName.IsNull()) {
2026  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
2027  }
2028  else {
2030  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoPbPbSummary", &fCurrentJetInfo[i]);
2031  }
2032  }
2033 
2034  return fTree;
2035 }
2036 
2041 {
2042  TString hname;
2043 
2044  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
2045 
2046  for (auto &jetDef : fJetDefinitions) {
2047 
2048  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
2049 
2050  Double_t radius = jetDef.fRadius;
2051 
2052  TString title[30] = {""};
2053  Int_t nbins[30] = {0 };
2054  Double_t min [30] = {0.};
2055  Double_t max [30] = {0.};
2056  Int_t dim = 0 ;
2057 
2058  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
2059  nbins[dim] = nPtBins;
2060  min[dim] = 0;
2061  max[dim] = fMaxPt;
2062  dim++;
2063 
2064  if ((enabledAxis & kPositionD) != 0) {
2065  title[dim] = "#eta_{D}";
2066  nbins[dim] = 50;
2067  min[dim] = -1;
2068  max[dim] = 1;
2069  dim++;
2070 
2071  title[dim] = "#phi_{D} (rad)";
2072  nbins[dim] = 150;
2073  min[dim] = 0;
2074  max[dim] = TMath::TwoPi();
2075  dim++;
2076  }
2077 
2078  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
2079  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
2080  nbins[dim] = fNMassBins;
2081  min[dim] = fMinMass;
2082  max[dim] = fMaxMass;
2083  dim++;
2084  }
2085 
2087  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
2088  nbins[dim] = fNMassBins;
2089  min[dim] = fMinMass;
2090  max[dim] = fMaxMass;
2091  dim++;
2092  }
2093 
2094  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
2095  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
2096  nbins[dim] = fNMassBins;
2097  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
2098  dim++;
2099  }
2100 
2101  if (fCandidateType == kDstartoKpipi) {
2102  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
2103  nbins[dim] = fNMassBins*6;
2104  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
2105 
2106  // subtract mass of D0
2107  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2108  min[dim] -= D0mass;
2109  max[dim] -= D0mass;
2110 
2111  dim++;
2112  }
2113 
2114  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
2115  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
2116  nbins[dim] = 100;
2117  min[dim] = 0;
2118  max[dim] = 25;
2119  dim++;
2120  }
2121 
2122  title[dim] = "#it{z}_{D}";
2123  nbins[dim] = 110;
2124  min[dim] = 0;
2125  max[dim] = 1.10;
2126  dim++;
2127 
2128  if ((enabledAxis & kDeltaR) != 0) {
2129  title[dim] = "#Delta R_{D-jet}";
2130  nbins[dim] = 100;
2131  min[dim] = 0;
2132  max[dim] = radius * 1.5;
2133  dim++;
2134  }
2135 
2136  if ((enabledAxis & kDeltaEta) != 0) {
2137  title[dim] = "#eta_{D} - #eta_{jet}";
2138  nbins[dim] = 100;
2139  min[dim] = -radius * 1.2;
2140  max[dim] = radius * 1.2;
2141  dim++;
2142  }
2143 
2144  if ((enabledAxis & kDeltaPhi) != 0) {
2145  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
2146  nbins[dim] = 100;
2147  min[dim] = -radius * 1.2;
2148  max[dim] = radius * 1.2;
2149  dim++;
2150  }
2151 
2152  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
2153  nbins[dim] = nPtBins;
2154  min[dim] = 0;
2155  max[dim] = fMaxPt;
2156  dim++;
2157 
2158  if ((enabledAxis & kPositionJet) != 0) {
2159  title[dim] = "#eta_{jet}";
2160  nbins[dim] = 50;
2161  min[dim] = -1;
2162  max[dim] = 1;
2163  dim++;
2164 
2165  title[dim] = "#phi_{jet} (rad)";
2166  nbins[dim] = 150;
2167  min[dim] = 0;
2168  max[dim] = TMath::TwoPi();
2169  dim++;
2170  }
2171 
2172  if ((enabledAxis & kJetConstituents) != 0) {
2173  title[dim] = "No. of constituents";
2174  nbins[dim] = 50;
2175  min[dim] = -0.5;
2176  max[dim] = 49.5;
2177  dim++;
2178  }
2179 
2180  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
2181  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
2182  for (Int_t j = 0; j < dim; j++) {
2183  h->GetAxis(j)->SetTitle(title[j]);
2184  }
2185  }
2186 }
2187 
2192 {
2193  TString hname;
2194  fPartons.clear();
2195 
2196  TH1* histAncestor = nullptr;
2197  TH1* histPrompt = nullptr;
2198 
2199  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2200  hname = TString::Format("%s/fHistPrompt", GetName());
2201  histPrompt = static_cast<TH1*>(fHistManager->FindObject(hname));
2202 
2203  hname = TString::Format("%s/fHistAncestor", GetName());
2204  histAncestor = static_cast<TH1*>(fHistManager->FindObject(hname));
2205  }
2206 
2207  for (auto& dmeson_pair : fDmesonJets) {
2208  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
2209  Int_t accJets = 0;
2210  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
2211  fCurrentJetInfo[ij]->Reset();
2212  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
2213  if (!jet) continue;
2214  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
2215  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
2216  fHistManager->FillTH1(hname, jet->Pt());
2217  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
2218  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2219  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
2220  fHistManager->FillTH1(hname, jet->Eta());
2221  continue;
2222  }
2223  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
2224  accJets++;
2225  }
2226  if (accJets > 0) {
2227  if (histPrompt) {
2228  if (dmeson_pair.second.fParton) {
2229  fPartons[dmeson_pair.second.fParton] = dmeson_pair.second.fPartonType;
2230  UInt_t absPdgParton = TMath::Abs(dmeson_pair.second.fParton->GetPdgCode());
2231  if (absPdgParton == 4) {
2232  histPrompt->Fill("Prompt", 1);
2233  }
2234  else if (absPdgParton == 5) {
2235  histPrompt->Fill("Non-Prompt", 1);
2236  }
2237  else {
2238  histPrompt->Fill("Unknown", 1);
2239  }
2240  }
2241  else {
2242  histPrompt->Fill("Unknown", 1);
2243  }
2244  }
2245 
2246  if (histAncestor) {
2247  if (dmeson_pair.second.fAncestor) {
2248  UInt_t absPdgAncestor = TMath::Abs(dmeson_pair.second.fAncestor->GetPdgCode());
2249  if (absPdgAncestor == 4) {
2250  histAncestor->Fill("Charm", 1);
2251  }
2252  else if (absPdgAncestor == 5) {
2253  histAncestor->Fill("Bottom", 1);
2254  }
2255  else if (absPdgAncestor == 2212) {
2256  histAncestor->Fill("Proton", 1);
2257  }
2258  else {
2259  histAncestor->Fill("Unknown", 1);
2260  }
2261  }
2262  else {
2263  histAncestor->Fill("Unknown", 1);
2264  }
2265  }
2266 
2267  fTree->Fill();
2268  }
2269  else {
2270  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2271  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2272  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2273  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2274  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2275  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2276  if (fMCMode != kMCTruth) {
2278  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
2279  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
2280  }
2281  else if (fCandidateType == kDstartoKpipi) {
2282  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
2283  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
2284 
2285  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
2286  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
2287  }
2288  }
2289  }
2290  }
2291 
2292  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2293  hname = TString::Format("%s/fHistPartonPt", GetName());
2294  TH1* histPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
2295  hname = TString::Format("%s/fHistPartonEta", GetName());
2296  TH1* histPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
2297  hname = TString::Format("%s/fHistPartonPhi", GetName());
2298  TH1* histPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
2299  hname = TString::Format("%s/fHistPartonType", GetName());
2300  TH1* histPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
2301 
2302  for (auto parton : fPartons) {
2303  if (!parton.first) continue;
2304  histPartonPt->Fill(parton.first->Pt());
2305  histPartonEta->Fill(parton.first->Eta());
2306  histPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
2307  histPartonType->Fill(parton.second);
2308  }
2309  }
2310 
2311  return kTRUE;
2312 }
2313 
2319 {
2320  TString hname;
2321 
2322  TH1* histAncestor = nullptr;
2323  TH1* histPrompt = nullptr;
2324 
2325  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2326  hname = TString::Format("%s/fHistPrompt", GetName());
2327  histPrompt = static_cast<TH1*>(fHistManager->FindObject(hname));
2328 
2329  hname = TString::Format("%s/fHistAncestor", GetName());
2330  histAncestor = static_cast<TH1*>(fHistManager->FindObject(hname));
2331  }
2332 
2333  fPartons.clear();
2334  for (auto& dmeson_pair : fDmesonJets) {
2335  Int_t accJets = 0;
2336  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
2337  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
2338  if (!jet) continue;
2339  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
2340  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
2341  fHistManager->FillTH1(hname, jet->Pt());
2342  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
2343  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2344  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
2345  fHistManager->FillTH1(hname, jet->Eta());
2346  continue;
2347  }
2348  accJets++;
2349  }
2350  if (accJets > 0) {
2351  if (histPrompt) {
2352  if (dmeson_pair.second.fParton) {
2353  fPartons[dmeson_pair.second.fParton] = dmeson_pair.second.fPartonType;
2354  UInt_t absPdgParton = TMath::Abs(dmeson_pair.second.fParton->GetPdgCode());
2355  if (absPdgParton == 4) {
2356  histPrompt->Fill("Prompt", 1);
2357  }
2358  else if (absPdgParton == 5) {
2359  histPrompt->Fill("Non-Prompt", 1);
2360  }
2361  else {
2362  histPrompt->Fill("Unknown", 1);
2363  }
2364  }
2365  else {
2366  histPrompt->Fill("Unknown", 1);
2367  }
2368  }
2369 
2370  if (histAncestor) {
2371  if (dmeson_pair.second.fAncestor) {
2372  UInt_t absPdgAncestor = TMath::Abs(dmeson_pair.second.fAncestor->GetPdgCode());
2373  if (absPdgAncestor == 4) {
2374  histAncestor->Fill("Charm", 1);
2375  }
2376  else if (absPdgAncestor == 5) {
2377  histAncestor->Fill("Bottom", 1);
2378  }
2379  else if (absPdgAncestor == 2212) {
2380  histAncestor->Fill("Proton", 1);
2381  }
2382  else {
2383  histAncestor->Fill("Unknown", 1);
2384  }
2385  }
2386  else {
2387  histAncestor->Fill("Unknown", 1);
2388  }
2389  }
2390  }
2391  else {
2392  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2393  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2394  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2395  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2396  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2397  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2398  if (fMCMode != kMCTruth) {
2400  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
2401  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
2402  }
2403  else if (fCandidateType == kDstartoKpipi) {
2404  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
2405  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
2406 
2407  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
2408  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
2409  }
2410  }
2411  }
2412  }
2413 
2414  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2415  hname = TString::Format("%s/fHistPartonPt", GetName());
2416  TH1* histPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
2417  hname = TString::Format("%s/fHistPartonEta", GetName());
2418  TH1* histPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
2419  hname = TString::Format("%s/fHistPartonPhi", GetName());
2420  TH1* histPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
2421  hname = TString::Format("%s/fHistPartonType", GetName());
2422  TH1* histPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
2423 
2424  for (auto parton : fPartons) {
2425  if (!parton.first) continue;
2426  histPartonPt->Fill(parton.first->Pt());
2427  histPartonEta->Fill(parton.first->Eta());
2428  histPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
2429  histPartonType->Fill(parton.second);
2430  }
2431  }
2432 
2433  return kTRUE;
2434 }
2435 
2440 {
2441  TString hname;
2442 
2443  for (auto& dmeson_pair : fDmesonJets) {
2444  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
2445  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2446  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2447  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2448  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2449  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2450  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2451  }
2452  }
2453 
2454  for (auto &jetDef : fJetDefinitions) {
2455 
2456  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
2457  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
2458 
2459  for (auto& dmeson_pair : fDmesonJets) {
2460  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
2461  if (!jet) continue;
2462  if (!jetDef.IsJetInAcceptance(*jet)) {
2463  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
2464  fHistManager->FillTH1(hname, jet->Pt());
2465  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
2466  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2467  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
2468  fHistManager->FillTH1(hname, jet->Eta());
2469  continue;
2470  }
2471  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
2472  }
2473  }
2474 
2475  return kTRUE;
2476 }
2477 
2484 {
2485  // Fill the THnSparse histogram.
2486 
2487  Double_t contents[30] = {0.};
2488 
2489  Double_t z = DmesonJet.GetZ(n);
2490  Double_t deltaPhi = 0;
2491  Double_t deltaEta = 0;
2492  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
2493 
2494  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
2495  if (it == DmesonJet.fJets.end()) return kFALSE;
2496 
2497  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
2498  TString title(h->GetAxis(i)->GetTitle());
2499  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
2500  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
2501  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
2502  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
2503  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
2504  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
2505  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
2506  else if (title=="#it{z}_{D}") contents[i] = z;
2507  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
2508  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
2509  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
2510  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
2511  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
2512  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
2513  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
2514  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
2515  }
2516 
2517  h->Fill(contents);
2518 
2519  return kTRUE;
2520 }
2521 
2522 // Definitions of class AliAnalysisTaskDmesonJets
2523 
2525 ClassImp(AliAnalysisTaskDmesonJets);
2527 
2531  fAnalysisEngines(),
2532  fEnabledAxis(0),
2534  fHistManager(),
2535  fApplyKinematicCuts(kTRUE),
2536  fNOutputTrees(0),
2537  fTrackEfficiency(0),
2538  fRejectISR(kFALSE),
2539  fJetAreaType(fastjet::active_area),
2540  fJetGhostArea(0.005),
2541  fMCContainer(0),
2542  fAodEvent(0),
2543  fFastJetWrapper(0)
2544 {
2545  SetMakeGeneralHistograms(kTRUE);
2546 }
2547 
2552  AliAnalysisTaskEmcalLight(name, kTRUE),
2553  fAnalysisEngines(),
2556  fHistManager(name),
2557  fApplyKinematicCuts(kTRUE),
2558  fNOutputTrees(nOutputTrees),
2559  fTrackEfficiency(0),
2560  fRejectISR(kFALSE),
2561  fJetAreaType(fastjet::active_area),
2562  fJetGhostArea(0.005),
2563  fMCContainer(0),
2564  fAodEvent(0),
2565  fFastJetWrapper(0)
2566 {
2567  SetMakeGeneralHistograms(kTRUE);
2568  for (Int_t i = 0; i < nOutputTrees; i++){
2569  DefineOutput(2+i, TTree::Class());
2570  }
2571 }
2572 
2575 {
2576  if (fFastJetWrapper) delete fFastJetWrapper;
2577 }
2578 
2586 {
2587  AliRDHFCuts* analysiscuts = 0;
2588  TFile* filecuts = TFile::Open(cutfname);
2589  if (!filecuts || filecuts->IsZombie()) {
2590  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Input file not found: will use std cuts.");
2591  filecuts = 0;
2592  }
2593 
2594  if (filecuts) analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
2595 
2596  if (!analysiscuts) {
2597  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Could not find analysis cuts '%s' in '%s'.", cutsname.Data(), cutfname.Data());
2598  if (filecuts) {
2599  filecuts->ls();
2600  }
2601  }
2602  else {
2603  ::Info("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Cuts '%s' loaded from file '%s'", cutsname.Data(), cutfname.Data());
2604  }
2605 
2606  return analysiscuts;
2607 }
2608 
2621 {
2623  return AddAnalysisEngine(type, cutfname, cuttype, MCmode, jetDef, rhoName);
2624 }
2625 
2637 {
2638  AliRDHFCuts* cuts = 0;
2639 
2640  if (!cutfname.IsNull()) {
2641  TString cutsname;
2642 
2643  switch (type) {
2644  case kD0toKpi:
2645  case kD0toKpiLikeSign:
2646  cutsname = "D0toKpiCuts";
2647  break;
2648  case kDstartoKpipi:
2649  cutsname = "DStartoKpipiCuts";
2650  break;
2651  default:
2652  return 0;
2653  }
2654 
2655  if (!cuttype.IsNull()) {
2656  cutsname += TString::Format("_%s", cuttype.Data());
2657  }
2658 
2659  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
2660  if (cuts) cuts->PrintAll();
2661  }
2662 
2663  AnalysisEngine eng(type, MCmode, cuts);
2664 
2665  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
2666 
2667  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
2668  eng.AddJetDefinition(jetDef);
2669  it = fAnalysisEngines.insert(it, eng);
2670  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
2671  }
2672  else {
2673  AnalysisEngine* found_eng = &(*it);
2674  ::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());
2675  found_eng->AddJetDefinition(jetDef);
2676 
2677  if (cuts) {
2678  if (found_eng->fRDHFCuts != 0) ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
2679  found_eng->SetRDHFCuts(cuts);
2680  }
2681  }
2682 
2683  return &(*it);
2684 }
2685 
2686 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
2687 {
2688  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
2689  while (it != fAnalysisEngines.end() && (*it) != eng) it++;
2690  return it;
2691 }
2692 
2695 {
2696  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
2697 
2699 
2700  // Define histograms
2701  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
2702 
2703  TString hname;
2704  TString htitle;
2705  TH1* h = 0;
2706  Int_t treeSlot = 0;
2707 
2708  Int_t maxTracks = 6000;
2709  Double_t maxRho = 500;
2710  if (fForceBeamType == kpp) {
2711  maxRho = 50;
2712  maxTracks = 200;
2713  }
2714  else if (fForceBeamType == kpA) {
2715  maxRho = 200;
2716  maxTracks = 500;
2717  }
2718 
2719  hname = "fHistCharmPt";
2720  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2721  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2722 
2723  hname = "fHistCharmEta";
2724  htitle = hname + ";#eta_{charm};counts";
2725  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2726 
2727  hname = "fHistCharmPhi";
2728  htitle = hname + ";#phi_{charm};counts";
2729  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2730 
2731  hname = "fHistBottomPt";
2732  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2733  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2734 
2735  hname = "fHistBottomEta";
2736  htitle = hname + ";#eta_{bottom};counts";
2737  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2738 
2739  hname = "fHistBottomPhi";
2740  htitle = hname + ";#phi_{bottom};counts";
2741  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
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()); // @suppress("Ambiguous problem")
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 PdgHighParton = 0;
3162  for (auto it = itcont.begin(); it != itcont.end(); it++) {
3163  auto part = *it;
3164  if (part.first.Pt() == 0) continue;
3165 
3166  Int_t PdgCode = part.second->GetPdgCode();
3167 
3168  // Skip all particles that are not either quarks or gluons
3169  if ((PdgCode < -9 || PdgCode > 9) && PdgCode != 21 && PdgCode != -21) continue;
3170 
3171  AliDebugStream(5) << "Parton " << it.current_index() <<
3172  " with pdg=" << PdgCode <<
3173  ", px=" << part.first.Px() <<
3174  ", py=" << part.first.Py() <<
3175  ", pz=" << part.first.Pz() <<
3176  ", n daughters = " << part.second->GetNDaughters() <<
3177  std::endl;
3178 
3179  // Skip partons that do not have any children
3180  // Unclear how this can happen, it would seem that this parton were not fragmented by the generator
3181  if (part.second->GetNDaughters() == 0) continue;
3182 
3183  // Look for highest momentum parton
3184  if (highestPartonPt < part.first.Pt()) {
3185  highestPartonPt = part.first.Pt();
3186  PdgHighParton = PdgCode;
3187  }
3188 
3189  // Skip partons that are not HF
3190  if (PdgCode != 4 && PdgCode != 5 && PdgCode != -4 && PdgCode != -5) continue;
3191 
3192  Bool_t lastInPartonShower = kTRUE;
3193  Bool_t hadronDaughter = kFALSE;
3194  for (Int_t daughterIndex = part.second->GetFirstDaughter(); daughterIndex <= part.second->GetLastDaughter(); daughterIndex++){
3195  if (daughterIndex < 0) {
3196  AliDebugStream(5) << "Could not find daughter index!" << std::endl;
3197  continue;
3198  }
3199  AliAODMCParticle *daughter = fMCContainer->GetMCParticle(daughterIndex);
3200  if (!daughter) {
3201  AliDebugStream(5) << "Could not find particle with index " << daughterIndex << "!" << std::endl;
3202  continue;
3203  }
3204  Int_t daughterPdgCode = daughter->GetPdgCode();
3205  if (daughter->GetMother() != it.current_index()) {
3206  AliDebugStream(5) << "Particle " << daughterIndex << " with pdg=" << daughterPdgCode <<
3207  ", px=" << daughter->Px() <<
3208  ", py=" << daughter->Py() <<
3209  ", pz=" << daughter->Pz() <<
3210  ", is not a daughter of " << it.current_index() <<
3211  "!" << std::endl;
3212  continue;
3213  }
3214 
3215  AliDebugStream(5) << "Found daughter " << daughterIndex <<
3216  " with pdg=" << daughterPdgCode <<
3217  ", px=" << daughter->Px() <<
3218  ", py=" << daughter->Py() <<
3219  ", pz=" << daughter->Pz() <<
3220  std::endl;
3221  if (daughterPdgCode == PdgCode) lastInPartonShower = kFALSE; // this parton is not the last parton in the shower
3222  if (TMath::Abs(daughterPdgCode) >= 111) hadronDaughter = kTRUE;
3223  }
3224  if (hadronDaughter) {
3225  AliDebugStream(5) << "This particle has at least a hadron among its daughters!" << std::endl;
3226  if (!lastInPartonShower) AliDebugStream(2) << "Odly, quark " << it.current_index() << " with PDG " << PdgCode << " (pt = " << part.first.Pt() << ", eta = " << part.first.Eta() << ") is not the last in the parton shower but at least a hadron found among its daughters?!" << std::endl;
3227  }
3228  else {
3229  AliDebugStream(5) << "This particle does not have hadrons among its daughters!" << std::endl;
3230  if (lastInPartonShower) AliDebugStream(2) << "Odly, quark " << it.current_index() << " with PDG " << PdgCode << " (pt = " << part.first.Pt() << ", eta = " << part.first.Eta() << ") is the last in the parton shower but no hadron found among its daughters?!" << std::endl;
3231  continue;
3232  }
3233 
3234  if (PdgCode == 4 || PdgCode == -4) {
3235  fHistManager.FillTH1("fHistCharmPt", part.first.Pt());
3236  fHistManager.FillTH1("fHistCharmEta", part.first.Eta());
3237  fHistManager.FillTH1("fHistCharmPhi", part.first.Phi_0_2pi());
3238  }
3239  else if (PdgCode == 5 || PdgCode == -5) {
3240  fHistManager.FillTH1("fHistBottomPt", part.first.Pt());
3241  fHistManager.FillTH1("fHistBottomEta", part.first.Eta());
3242  fHistManager.FillTH1("fHistBottomPhi", part.first.Phi_0_2pi());
3243  }
3244  nHQ++;
3245  }
3246  fHistManager.FillTH1("fHistNHeavyQuarks", nHQ);
3247  fHistManager.FillTH1("fHistHighestPartonPt",highestPartonPt);
3248  Int_t partonType = 0;
3249  Int_t absPdgHighParton = TMath::Abs(PdgHighParton);
3250  if (absPdgHighParton == 9 || absPdgHighParton == 21) {
3251  partonType = 7;
3252  }
3253  else {
3254  partonType = absPdgHighParton;
3255  }
3256  fHistManager.FillTH1("fHistHighestPartonType",partonType);
3257 }
3258 
3267 {
3268  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
3269 
3270  Double_t mass = part->Mass();
3271 
3272  // To make sure that the PDG mass value is not at the edge of a bin
3273  if (nbins % 2 == 0) {
3274  minMass = mass - range / 2 - range / nbins / 2;
3275  maxMass = mass + range / 2 - range / nbins / 2;
3276  }
3277  else {
3278  minMass = mass - range / 2;
3279  maxMass = mass + range / 2;
3280  }
3281 }
3282 
3289 {
3290  static TString label;
3291  label = "";
3292 
3293  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
3294  label = "NotSelTrigger";
3295  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
3296  return label.Data();
3297  }
3298  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
3299  label = "NoVertex";
3300  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
3301  return label.Data();
3302  }
3303  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
3304  label = "TooFewVtxContrib";
3305  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
3306  return label.Data();
3307  }
3308  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
3309  label = "ZVtxOutFid";
3310  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
3311  return label.Data();
3312  }
3313  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
3314  label = "Pileup";
3315  bitmap &= ~BIT(AliRDHFCuts::kPileup);
3316  return label.Data();
3317  }
3318  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
3319  label = "OutsideCentrality";
3320  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
3321  return label.Data();
3322  }
3323  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
3324  label = "PhysicsSelection";
3325  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
3326  return label.Data();
3327  }
3328  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
3329  label = "BadSPDVertex";
3330  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
3331  return label.Data();
3332  }
3333  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
3334  label = "ZVtxSPDOutFid";
3335  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
3336  return label.Data();
3337  }
3338  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
3339  label = "CentralityFlattening";
3340  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
3341  return label.Data();
3342  }
3343  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
3344  label = "BadTrackV0Correl";
3345  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
3346  return label.Data();
3347  }
3348 
3349  return label.Data();
3350 }
3351 
3358 {
3359  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
3360  PostData(eng.GetDataSlotNumber(), eng.GetTree());
3361  return eng.GetDataSlotNumber();
3362  }
3363  else {
3364  return -1;
3365  }
3366 }
3367 
3377 {
3378  // Get the pointer to the existing analysis manager via the static access method
3379  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
3380  if (!mgr) {
3381  ::Error("AddTaskDmesonJets", "No analysis manager to connect to.");
3382  return NULL;
3383  }
3384 
3385  // Check the analysis type using the event handlers connected to the analysis manager
3386  AliVEventHandler* handler = mgr->GetInputEventHandler();
3387  if (!handler) {
3388  ::Error("AddTaskEmcalJetSpectraQA", "This task requires an input event handler");
3389  return NULL;
3390  }
3391 
3392  EDataType_t dataType = kUnknownDataType;
3393 
3394  if (handler->InheritsFrom("AliESDInputHandler")) {
3395  dataType = kESD;
3396  }
3397  else if (handler->InheritsFrom("AliAODInputHandler")) {
3398  dataType = kAOD;
3399  }
3400 
3401  // Init the task and do settings
3402  if (ntracks == "usedefault") {
3403  if (dataType == kESD) {
3404  ntracks = "Tracks";
3405  }
3406  else if (dataType == kAOD) {
3407  ntracks = "tracks";
3408  }
3409  else {
3410  ntracks = "";
3411  }
3412  }
3413 
3414  if (nclusters == "usedefault") {
3415  if (dataType == kESD) {
3416  nclusters = "CaloClusters";
3417  }
3418  else if (dataType == kAOD) {
3419  nclusters = "caloClusters";
3420  }
3421  else {
3422  nclusters = "";
3423  }
3424  }
3425 
3426  if (nMCpart == "usedefault") {
3427  nMCpart = "mcparticles"; // Always needs AliAODMCParticle objects
3428  }
3429 
3430  TString name("AliAnalysisTaskDmesonJets");
3431  if (strcmp(suffix, "") != 0) {
3432  name += TString::Format("_%s", suffix.Data());
3433  }
3434 
3435  AliAnalysisTaskDmesonJets* jetTask = new AliAnalysisTaskDmesonJets(name, nMaxTrees);
3436 
3437  if (!ntracks.IsNull()) {
3438  AliHFTrackContainer* trackCont = new AliHFTrackContainer(ntracks);
3439  jetTask->AdoptParticleContainer(trackCont);
3440  }
3441 
3442  if (!nMCpart.IsNull()) {
3443  AliMCParticleContainer* partCont = new AliHFAODMCParticleContainer(nMCpart);
3444  partCont->SetEtaLimits(-1.5, 1.5);
3445  partCont->SetPtLimits(0, 1000);
3446  jetTask->AdoptParticleContainer(partCont);
3447  }
3448 
3449  jetTask->AddClusterContainer(nclusters.Data());
3450 
3451  // Final settings, pass to manager and set the containers
3452  mgr->AddTask(jetTask);
3453 
3454  // Create containers for input/output
3455  AliAnalysisDataContainer* cinput1 = mgr->GetCommonInputContainer();
3456  TString contname1(name);
3457  contname1 += "_histos";
3458  AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(contname1.Data(),
3459  TList::Class(), AliAnalysisManager::kOutputContainer,
3460  Form("%s", AliAnalysisManager::GetCommonFileName()));
3461 
3462  mgr->ConnectInput(jetTask, 0, cinput1);
3463  mgr->ConnectOutput(jetTask, 1, coutput1);
3464 
3465  for (Int_t i = 0; i < nMaxTrees; i++) {
3466  TString contname = TString::Format("%s_tree_%d", name.Data(), i);
3467  AliAnalysisDataContainer *coutput = mgr->CreateContainer(contname.Data(),
3468  TTree::Class(),AliAnalysisManager::kOutputContainer,
3469  Form("%s", AliAnalysisManager::GetCommonFileName()));
3470  mgr->ConnectOutput(jetTask, 2+i, coutput);
3471  }
3472  return jetTask;
3473 }
Double32_t fCorrZ
Z of the D meson after subtracting average background.
Lightweight class that encapsulates D meson jets.
Int_t pdg
void Print() const
Prints the content of this object in the standard output.
AliDmesonInfoSummary * fCurrentDmesonJetInfo
! Current D meson jet info
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.
UInt_t fRejectedOrigin
Bit mask with D meson origins that are rejected (used for MC analysis, i.e. signal-only, background-only and particle-level)
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)
TString fRhoName
Name of the object that holds the average background value.
Lightweight class that encapsulates D meson jets.
std::vector< AliTrackContainer * > fTrackContainers
! Track containers
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)
std::vector< AliClusterContainer * > fClusterContainers
! Cluster containers
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.
Double32_t fPartonPt
Transverse momentum of the parton.
friend bool operator<(const AnalysisEngine &lhs, const AnalysisEngine &rhs)
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()
Double32_t fInvMass
Invariant mass of the D0 meson candidate in GeV/c2.
Lightweight class that encapsulates D meson jets.
Double_t fInvMass2Prong
! 2-prong mass of the D* candidate (w/o the soft pion)
friend bool operator==(const AnalysisEngine &lhs, const AnalysisEngine &rhs)
Double_t fMinChargedPt
Minimum pt of the leading charged particle (or track)
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.
Double32_t fPtK
Transverse momentum of the kaon.
TRandom * fRandomGen
! Random number generator
Double32_t fCosPointing
Cosine of the pointing angle.
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.
Int_t fDataSlotNumber
! Data slot where the tree output is posted
Double32_t fPtPi
Transverse momentum of the pion.
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
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
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)
Double32_t fN
Number of jet constituents.
void SetR(Double_t r)
Definition: AliFJWrapper.h:127
Select tracks based on specific prescriptions of HF analysis.
Double32_t fPt
Transverse momentum of the jet in GeV/c.
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
Double32_t fR
Distance between D meson and jet axis.
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)
Double32_t fPt
Transverse momentum of the D meson in GeV/c.
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)
Double_t fMinNeutralPt
Minimum pt of the leading neutral particle (or cluster)
Bool_t fD0D0bar
! kTRUE if selected both as D0 and D0bar
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
Bool_t fD0Extended
Store extended information in the tree (only for D0 mesons)
AliRDHFCuts * LoadDMesonCutsFromFile(TString cutfname, TString cutsname)
std::list< AnalysisEngine >::iterator FindAnalysisEngine(const AnalysisEngine &eng)
std::vector< AliHFJetDefinition > fJetDefinitions
Jet definitions.
Double_t fMaxNeutralPt
Maximum pt of the leading neutral particle (or cluster)
Double_t fMaxMass
Max mass in histogram axis.
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 *="")
friend bool operator<(const AliHFJetDefinition &lhs, const AliHFJetDefinition &rhs)
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)
TClonesArray * fCandidateArray
! D meson candidate array
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)
Double32_t fCorrPt
Transverse momentum of the jet in GeV/c after subtracting average background.
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
friend bool operator==(const AliHFJetDefinition &lhs, const AliHFJetDefinition &rhs)
static FJJetAlgo ConvertToFJAlgo(EJetAlgo_t algo)
Double_t fMaxChargedPt
Maximum pt of the leading charged particle (or track)
decay
Definition: HFPtSpectrum.C:41
std::map< std::string, AliJetInfo > fJets
! list of jets
Look for the very first particle in the fragmentation tree.
AliAODMCParticle * fAncestor
! pointer to the ancestor particle in the shower tree of the D meson (only for particle level D meson...
Double32_t fDeltaInvMass
< Difference between the Kpipi and the Kpi invariant masses in GeV/c2
std::map< AliAODMCParticle *, Short_t > fPartons
! set of the partons in the shower that produced each D meson
Definition: External.C:220
Double_t minMass
Bool_t fRejectISR
! Reject initial state radiation
void RunParticleLevelAnalysis()
Run a particle level analysis.
std::map< std::string, AliParticleContainer * > fParticleCollArray
particle/track collection array
std::map< int, AliDmesonJetInfo > fDmesonJets
! Array containing the D meson jets
Double_t fTrackEfficiency
! Artificial tracking inefficiency (0...1) -> set automatically at ExecOnce by AliAnalysisTaskDmesonJ...
Bool_t fApplyKinematicCuts
Apply jet kinematic cuts.
Double32_t fSelectionType
Selection type: D0, D0bar, both.
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:290
const char * GetName() const
Generate a name for this jet definition.
AliEventNotFound(const std::string &class_name, const std::string &method_name)
UInt_t fAcceptedDecay
Bit mask with D meson decays that are accepted (only used for particle-level analysis) ...
virtual void Set(const AliDmesonJetInfo &source, std::string n)
unsigned short UShort_t
Definition: External.C:28
Bool_t fReconstructed
! Whether this D meson was reconstructed (only for particle level D mesons)
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)
Double_t InvMassDstarKpipi() const
const AliVEvent * fEvent
! pointer to the ESD/AOD event
TArrayI fPDGdaughters
List of the PDG code of the daughters.
const Int_t nbins
Double_t maxMass
TString fBranchName
AOD branch where the D meson candidate are found.
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.
AliHFAODMCParticleContainer * fMCContainer
! MC particle container
AliJetInfoSummary ** fCurrentJetInfo
! Current jet info
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:310
AliFJWrapper * fFastJetWrapper
! Fastjet wrapper
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
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)
Double_t fMinMass
Min mass in histogram axis.