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