AliPhysics  4ea6a45 (4ea6a45)
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 #include <array>
19 
20 // Root
21 #include <TClonesArray.h>
22 #include <TDatabasePDG.h>
23 #include <TParticlePDG.h>
24 #include <TVector3.h>
25 #include <THnSparse.h>
26 #include <TParticle.h>
27 #include <TMath.h>
28 #include <THashList.h>
29 #include <TFile.h>
30 #include <TRandom3.h>
31 
32 // Aliroot general
33 #include "AliLog.h"
34 #include "AliEMCALGeometry.h"
35 #include "AliAnalysisManager.h"
36 #include "AliVEventHandler.h"
37 
38 // Aliroot HF
39 #include "AliAODRecoCascadeHF.h"
41 #include "AliRDHFCutsD0toKpi.h"
44 #include "AliHFTrackContainer.h"
45 #include "AliAnalysisVertexingHF.h"
46 
47 // Aliroot EMCal jet framework
48 #include "AliEmcalJetTask.h"
49 #include "AliEmcalJet.h"
50 #include "AliJetContainer.h"
51 #include "AliParticleContainer.h"
52 #include "AliEmcalParticle.h"
53 #include "AliFJWrapper.h"
54 #include "AliRhoParameter.h"
55 
57 
58 AliAnalysisTaskDmesonJets::AliEventNotFound::AliEventNotFound(const std::string& class_name, const std::string& method_name) :
59  std::exception(),
60  fClassName(class_name),
61  fAccessMethodName(method_name)
62 {
63  std::stringstream what_str;
64  what_str << "ALICE event not found in class '" << fClassName << "' using method '" << method_name << "'.";
65  fWhat = what_str.str();
66 }
67 
68 #if !(defined(__CINT__) || defined(__MAKECINT__))
70 {
71  return fWhat.c_str();
72 }
73 #endif
74 
75 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfo
76 
80 
88 {
89  dphi = TVector2::Phi_mpi_pi(fMomentum.Phi() - jet.Phi());;
90  deta = fMomentum.Eta() - jet.Eta();
91  return TMath::Sqrt(dphi*dphi + deta*deta);
92 }
93 
99 {
100  Double_t deta = 0;
101  Double_t dphi = 0;
102  return GetDistance(jet, deta, dphi);
103 }
104 
105 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonJetInfo
106 
110 
113  fDmesonParticle(0),
114  fD(),
115  fSoftPionPt(0),
116  fInvMass2Prong(0),
117  fJets(),
118  fMCLabel(-1),
119  fReconstructed(kFALSE),
120  fParton(0),
121  fPartonType(0),
122  fAncestor(0),
123  fD0D0bar(kFALSE),
124  fSelectionType(0),
125  fEvent(nullptr)
126 {
127 }
128 
134  fD(source.fD),
135  fSoftPionPt(source.fSoftPionPt),
137  fJets(source.fJets),
138  fMCLabel(source.fMCLabel),
140  fParton(source.fParton),
141  fPartonType(source.fPartonType),
142  fAncestor(source.fAncestor),
143  fD0D0bar(source.fD0D0bar),
145  fEvent(source.fEvent)
146 {
147 }
148 
153 {
154  new (this) AliDmesonJetInfo(source);
155  return *this;
156 }
157 
160 {
161  fD.SetPtEtaPhiE(0,0,0,0);
162  fSoftPionPt = 0;
163  fInvMass2Prong = 0;
164  fDmesonParticle = 0;
165  fMCLabel = -1;
166  fReconstructed = kFALSE;
167  fParton = 0;
168  fPartonType = 0;
169  fAncestor = 0;
170  fD0D0bar = kFALSE;
171  for (auto &jet : fJets) {
172  jet.second.fMomentum.SetPtEtaPhiE(0,0,0,0);
173  jet.second.fNConstituents = 0;
174  jet.second.fNEF = 0;
175  jet.second.fMaxChargedPt = 0;
176  jet.second.fMaxNeutralPt = 0;
177  }
178 }
179 
182 {
183  Printf("Printing D Meson Jet object.");
184  Printf("D Meson: pT = %.3f, eta = %.3f, phi = %.3f, inv. mass = %.3f", fD.Pt(), fD.Eta(), fD.Phi_0_2pi(), fD.M());
185  Printf("Soft pion pT: %.3f. 2-Prong Invariant mass = %.3f", fSoftPionPt, fInvMass2Prong);
186  for (auto &jet : fJets) {
187  Printf("Jet %s: pT = %.3f, eta = %.3f, phi = %.3f", jet.first.c_str(), jet.second.Pt(), jet.second.Eta(), jet.second.Phi_0_2pi());
188  Printf("Jet N Consituents = %d", jet.second.fNConstituents);
189  }
190 }
191 
196 {
197  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
198  if (it == fJets.end()) return 0;
199 
200  Double_t z = 0;
201 
202  if ((*it).second.Pt() > 0) {
203  TVector3 dvect = fD.Vect();
204  TVector3 jvect = (*it).second.fMomentum.Vect();
205 
206  Double_t jetMom = jvect * jvect;
207 
208  if (jetMom < 1e-6) {
209  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetZ", "Zero jet momentum!");
210  z = 0.999;
211  }
212  else {
213  z = (dvect * jvect) / jetMom;
214  }
215 
216  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
217  }
218 
219  return z;
220 }
221 
226 {
227  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
228  if (it == fJets.end()) return 0;
229 
230  Double_t z = 0;
231 
232  if ((*it).second.Pt() > 0) {
233  TVector3 dvect = fD.Vect();
234  TVector3 jvect = (*it).second.fMomentum.Vect();
235  // If the corr pt is < 0, assign 0.
236  Double_t corrpt = (*it).second.fCorrPt > 0 ? (*it).second.fCorrPt : 0.;
237  jvect.SetPerp(corrpt);
238 
239  Double_t jetMom = jvect * jvect;
240 
241  if (jetMom < 1e-6) {
242  z = 1.0;
243  }
244  else {
245  z = (dvect * jvect) / jetMom;
246  }
247  }
248 
249  return z;
250 }
251 
259 {
260  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
261  if (it == fJets.end()) return 0;
262 
263  return GetDistance((*it).second, deta, dphi);
264 }
265 
271 {
272  Double_t deta = 0;
273  Double_t dphi = 0;
274  return GetDistance(n, deta, dphi);
275 }
276 
284 {
285  dphi = TVector2::Phi_mpi_pi(fD.Phi() - jet.Phi());;
286  deta = fD.Eta() - jet.Eta();
287  return TMath::Sqrt(dphi*dphi + deta*deta);
288 }
289 
295 {
296  Double_t deta = 0;
297  Double_t dphi = 0;
298  return GetDistance(jet, deta, dphi);
299 }
300 
306 {
307  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
308  if (it == fJets.end()) {
309  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
310  n.c_str());
311  return 0;
312  }
313  return &((*it).second);
314 }
315 
321 {
322  std::map<std::string, AliJetInfo>::iterator it = fJets.find(n);
323  if (it == fJets.end()) {
324  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
325  n.c_str());
326  return 0;
327  }
328  return &((*it).second);
329 }
330 
331 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoSummary
332 
336 
342  fPt(0),
343  fEta(0),
344  fPhi(0),
345  fR(0),
346  fZ(0),
347  fN(0)
348 {
349  Set(source, n);
350 }
351 
354 {
355  fPt = 0;
356  fEta = 0;
357  fPhi = 0;
358  fR = 0;
359  fZ = 0;
360  fN = 0;
361 }
362 
368 {
369  std::map<std::string, AliJetInfo>::const_iterator it = source.fJets.find(n);
370  if (it == source.fJets.end()) return;
371 
372  Set((*it).second);
373 
374  fR = source.GetDistance(n);
375  fZ = source.GetZ(n);
376 }
377 
382 {
383  fPt = source.Pt();
384  fEta = source.Eta();
385  fPhi = source.Phi_0_2pi();
386  fN = source.GetNConstituents();
387  fR = 0;
388  fZ = 0;
389 }
390 
391 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoPbPbSummary
392 
396 
403  fCorrPt(0),
404  fCorrZ(0),
405  fArea(0)
406 {
407  Set(source, n);
408 }
409 
412 {
414  fCorrPt = 0;
415  fCorrZ = 0;
416  fArea = 0;
417 }
418 
423 {
424  AliJetInfoSummary::Set(source);
425  fArea = source.fArea;
426  fCorrPt = source.fCorrPt;
427 }
428 
434 {
435  AliJetInfoSummary::Set(source, n);
436  fCorrZ = source.GetCorrZ(n);
437 }
438 
439 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonInfoSummary
440 
444 
449  fPt(0),
450  fEta(0),
451  fPhi(0)
452 {
453  Set(source);
454 }
455 
460 {
461  fPt = source.fD.Pt();
462  fEta = source.fD.Eta();
463  fPhi = source.fD.Phi_0_2pi();
464 }
465 
468 {
469  fPt = 0;
470  fEta = 0;
471  fPhi = 0;
472 }
473 
474 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary
475 
479 
484  AliDmesonInfoSummary(source),
485  fPartonType(0),
486  fPartonPt(0),
487  fAncestorPDG(0)
488 {
489  Set(source);
490 }
491 
496 {
498 
499  fPartonType = source.fPartonType;
500 
501  if (source.fParton) {
502  fPartonPt = source.fParton->Pt();
503  }
504  else {
505  fPartonPt = 0.;
506  }
507 
508  if (source.fAncestor) {
509  fAncestorPDG = (UShort_t)((UInt_t)(TMath::Abs(source.fAncestor->GetPdgCode())));
510  }
511 }
512 
515 {
517  fPartonType = 0,
518  fPartonPt = 0.;
519  fAncestorPDG = 0;
520 }
521 
522 // Definitions of class AliAnalysisTaskDmesonJets::AliD0InfoSummary
523 
527 
532  AliDmesonInfoSummary(source),
533  fInvMass(source.fD.M()),
534  fSelectionType(0)
535 {
536 }
537 
542 {
543  fInvMass = source.fD.M();
546 }
547 
550 {
552  fSelectionType = 0;
553  fInvMass = 0;
554 }
555 
556 // Definitions of class AliAnalysisTaskDmesonJets::AliD0ExtendedInfoSummary
557 
561 
566  AliD0InfoSummary(source),
567  fDCA(0),
568  fCosThetaStar(0),
569  fd0K(0),
570  fd0Pi(0),
571  fd0d0(0),
572  fCosPointing(0),
573  fMaxNormd0(0)
574 {
575  Set(source);
576 }
577 
582 {
583  AliD0InfoSummary::Set(source);
584 
585  AliAODRecoDecayHF2Prong* recoDecay = dynamic_cast<AliAODRecoDecayHF2Prong*>(source.fDmesonParticle);
586  if (recoDecay) {
587  fDCA = recoDecay->GetDCA();
588  if (source.fSelectionType == 1) { // D0
589  fCosThetaStar = recoDecay->CosThetaStarD0();
590  fPtK = recoDecay->PtProng(0);
591  fPtPi = recoDecay->PtProng(1);
592  fd0K = recoDecay->Getd0Prong(0);
593  fd0Pi = recoDecay->Getd0Prong(1);
594  }
595  else { //D0bar
596  fCosThetaStar = recoDecay->CosThetaStarD0bar();
597  fPtK = recoDecay->PtProng(1);
598  fPtPi = recoDecay->PtProng(0);
599  fd0K = recoDecay->Getd0Prong(1);
600  fd0Pi = recoDecay->Getd0Prong(0);
601  }
602 
603  fMaxNormd0 = 0.;
604  // Based on Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent* aod)
605  // Line 480 and following
606  if (source.fEvent) {
607  for (Int_t ipr=0; ipr < 2; ipr++) {
608  Double_t diffIP = 0., errdiffIP = 0.;
609  recoDecay->Getd0MeasMinusExpProng(ipr, source.fEvent->GetMagneticField(), diffIP, errdiffIP);
610  Double_t normdd0 = 0.;
611  if (errdiffIP > 0.) {
612  normdd0 = diffIP / errdiffIP;
613  }
614  else {
615  if (diffIP == 0) {
616  normdd0 = 0;
617  }
618  else {
619  normdd0 = diffIP > 0 ? 9999. : -9999.;
620  }
621  }
622  if (TMath::Abs(normdd0) > TMath::Abs(fMaxNormd0)) {
623  fMaxNormd0 = normdd0;
624  }
625  }
626  }
627  else {
628  throw AliAnalysisTaskDmesonJets::AliEventNotFound("AliAnalysisTaskDmesonJets::AliDmesonJetInfo", "fEvent");
629  }
630 
631  fd0d0 = recoDecay->Prodd0d0();
632  fCosPointing = recoDecay->CosPointingAngle();
633  }
634 }
635 
638 {
640  fDCA = 0;
641  fCosThetaStar = 0;
642  fd0K = 0;
643  fd0Pi = 0;
644  fd0d0 = 0;
645  fCosPointing = 0;
646  fMaxNormd0 = 0;
647 }
648 
649 // Definitions of class AliAnalysisTaskDmesonJets::AliDStarInfoSummary
650 
654 
659  AliDmesonInfoSummary(source),
660  f2ProngInvMass(source.fInvMass2Prong),
661  fDeltaInvMass(source.fD.M() - source.fInvMass2Prong)
662 {
663 }
664 
669 {
671  fDeltaInvMass = source.fD.M() - source.fInvMass2Prong;
673 }
674 
677 {
679 
680  f2ProngInvMass = 0;
681  fDeltaInvMass = 0;
682 }
683 
684 // Definitions of class AliAnalysisTaskDmesonJets::AliJetDefinition
685 
689 
692  TObject(),
693  fJetType(AliJetContainer::kChargedJet),
694  fRadius(0),
695  fJetAlgo(AliJetContainer::antikt_algorithm),
696  fRecoScheme(AliJetContainer::pt_scheme),
697  fMinJetPt(0.),
698  fMaxJetPt(500.),
699  fMinJetPhi(0.),
700  fMaxJetPhi(0.),
701  fMinJetEta(-1.),
702  fMaxJetEta(1.),
703  fMinChargedPt(0.),
704  fMaxChargedPt(0.),
705  fMinNeutralPt(0.),
706  fMaxNeutralPt(0.),
707  fRhoName(),
708  fRho(0)
709 {
710 }
711 
719  TObject(),
720  fJetType(type),
721  fRadius(r),
722  fJetAlgo(algo),
723  fRecoScheme(reco),
724  fMinJetPt(0.),
725  fMaxJetPt(500.),
726  fMinJetPhi(0.),
727  fMaxJetPhi(0.),
728  fMinJetEta(-1.),
729  fMaxJetEta(1.),
730  fMinChargedPt(0.),
731  fMaxChargedPt(0.),
732  fMinNeutralPt(0.),
733  fMaxNeutralPt(0.),
734  fRhoName(),
735  fRho(0)
736 {
737 }
738 
746  TObject(),
747  fJetType(type),
748  fRadius(r),
749  fJetAlgo(algo),
750  fRecoScheme(reco),
751  fMinJetPt(0.),
752  fMaxJetPt(0.),
753  fMinJetPhi(0.),
754  fMaxJetPhi(0.),
755  fMinJetEta(0.),
756  fMaxJetEta(0.),
757  fMinChargedPt(0.),
758  fMaxChargedPt(0.),
759  fMinNeutralPt(0.),
760  fMaxNeutralPt(0.),
761  fRhoName(rhoName),
762  fRho(0)
763 {
764 }
765 
770  TObject(),
771  fJetType(source.fJetType),
772  fRadius(source.fRadius),
773  fJetAlgo(source.fJetAlgo),
774  fRecoScheme(source.fRecoScheme),
775  fMinJetPt(source.fMinJetPt),
776  fMaxJetPt(source.fMaxJetPt),
777  fMinJetPhi(source.fMinJetPhi),
778  fMaxJetPhi(source.fMaxJetPhi),
779  fMinJetEta(source.fMinJetEta),
780  fMaxJetEta(source.fMaxJetEta),
785  fRhoName(source.fRhoName),
786  fRho(0)
787 {
788 }
789 
794 {
795  new (this) AliHFJetDefinition(source);
796  return *this;
797 }
798 
801 {
802  static TString name;
803 
805 
806  return name.Data();
807 }
808 
814 {
815  if (fMinJetEta < fMaxJetEta && (jet.Eta() < fMinJetEta || jet.Eta() > fMaxJetEta)) return kFALSE;
816  if (fMinJetPhi < fMaxJetPhi && (jet.Phi() < fMinJetPhi || jet.Phi() > fMaxJetPhi)) return kFALSE;
817  if (fMinJetPt < fMaxJetPt && (jet.Pt() > fMaxJetPt || jet.Pt() < fMinJetPt)) return kFALSE;
818  if (fMinChargedPt < fMaxChargedPt && (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt)) return kFALSE;
819  if (fMinNeutralPt < fMaxNeutralPt && (jet.fMaxNeutralPt < fMinNeutralPt || jet.fMaxNeutralPt > fMaxNeutralPt)) return kFALSE;
820 
821  return kTRUE;
822 }
823 
829 {
830  const AliJetInfo* jet = dMesonJet.GetJet(n);
831  if (!jet) return kFALSE;
832  return IsJetInAcceptance((*jet));
833 }
834 
841 {
842  if (lhs.fJetType > rhs.fJetType) return false;
843  else if (lhs.fJetType < rhs.fJetType) return true;
844  else {
845  if (lhs.fRadius > rhs.fRadius) return false;
846  else if (lhs.fRadius < rhs.fRadius) return true;
847  else {
848  if (lhs.fJetAlgo > rhs.fJetAlgo) return false;
849  else if (lhs.fJetAlgo < rhs.fJetAlgo) return true;
850  else {
851  if (lhs.fRecoScheme < rhs.fRecoScheme) return true;
852  else return false;
853  }
854  }
855  }
856 }
857 
864 {
865  if (lhs.fJetType != rhs.fJetType) return false;
866  if (lhs.fRadius != rhs.fRadius) return false;
867  if (lhs.fJetAlgo != rhs.fJetAlgo) return false;
868  if (lhs.fRecoScheme != rhs.fRecoScheme) return false;
869  return true;
870 }
871 
872 // Definitions of class AliAnalysisTaskDmesonJets::AnalysisEngine
873 
877 
880  TObject(),
881  fPartons(),
882  fCandidateType(kD0toKpi),
883  fCandidateName(),
884  fCandidatePDG(0),
885  fNDaughters(0),
886  fPDGdaughters(),
887  fBranchName(),
888  fMCMode(kNoMC),
889  fNMassBins(0),
890  fMinMass(0),
891  fMaxMass(0),
892  fRDHFCuts(0),
893  fRejectedOrigin(0),
894  fAcceptedDecay(0),
895  fInhibit(kFALSE),
896  fJetDefinitions(),
897  fPtBinWidth(0.5),
898  fMaxPt(100),
899  fD0Extended(kFALSE),
900  fRandomGen(0),
901  fTrackEfficiency(0),
902  fRejectISR(kFALSE),
903  fDataSlotNumber(-1),
904  fTree(0),
905  fCurrentDmesonJetInfo(0),
906  fCurrentJetInfo(0),
907  fCandidateArray(0),
908  fMCContainer(),
909  fTrackContainers(),
910  fClusterContainers(),
911  fAodEvent(0),
912  fFastJetWrapper(0),
913  fHistManager(0),
914  fCent(-1)
915 {
916 }
917 
926  TObject(),
927  fPartons(),
928  fCandidateType(type),
929  fCandidateName(),
930  fCandidatePDG(0),
931  fNDaughters(0),
932  fPDGdaughters(),
933  fBranchName(),
934  fMCMode(MCmode),
935  fNMassBins(nMassBins),
936  fMinMass(0),
937  fMaxMass(0),
938  fRDHFCuts(cuts),
939  fRejectedOrigin(0),
941  fInhibit(kFALSE),
942  fJetDefinitions(),
943  fPtBinWidth(0.5),
944  fMaxPt(100),
945  fD0Extended(kFALSE),
946  fRandomGen(0),
947  fTrackEfficiency(0),
948  fDataSlotNumber(-1),
949  fTree(0),
951  fCurrentJetInfo(0),
952  fCandidateArray(0),
953  fMCContainer(),
956  fAodEvent(0),
957  fFastJetWrapper(0),
958  fHistManager(0),
959  fCent(-1)
960 {
961  SetCandidateProperties(range);
962 }
963 
968  TObject(source),
969  fPartons(source.fPartons),
973  fNDaughters(source.fNDaughters),
975  fBranchName(source.fBranchName),
976  fMCMode(source.fMCMode),
977  fNMassBins(source.fNMassBins),
978  fMinMass(source.fMinMass),
979  fMaxMass(source.fMaxMass),
980  fRDHFCuts(),
983  fInhibit(source.fInhibit),
985  fPtBinWidth(source.fPtBinWidth),
986  fMaxPt(source.fMaxPt),
987  fD0Extended(source.fD0Extended),
988  fRandomGen(source.fRandomGen),
990  fDataSlotNumber(-1),
991  fTree(0),
993  fCurrentJetInfo(0),
995  fMCContainer(source.fMCContainer),
998  fAodEvent(source.fAodEvent),
1000  fHistManager(source.fHistManager),
1001  fCent(-1)
1002 {
1003  SetRDHFCuts(source.fRDHFCuts);
1004 }
1005 
1006 // Destructor
1008 {
1009  delete fRDHFCuts;
1010 }
1011 
1016 {
1017  new (this) AnalysisEngine(source);
1018  return *this;
1019 }
1020 
1025 {
1026  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
1027  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
1028  }
1029 
1030  return kFALSE;
1031 }
1032 
1034 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const /*geom*/, Int_t /*runNumber*/)
1035 {
1036 }
1037 
1042 {
1043  switch (fCandidateType) {
1044  case kD0toKpi:
1045  fCandidatePDG = 421;
1046  fCandidateName = "D0";
1047  fNDaughters = 2;
1049  fPDGdaughters.Reset();
1050  fPDGdaughters[0] = 211; // pi
1051  fPDGdaughters[1] = 321; // K
1052  fBranchName = "D0toKpi";
1054  break;
1055  case kD0toKpiLikeSign:
1056  fCandidatePDG = 421;
1057  fCandidateName = "2ProngLikeSign";
1058  fNDaughters = 2;
1060  fPDGdaughters.Reset();
1061  fPDGdaughters[0] = 211; // pi
1062  fPDGdaughters[1] = 321; // K
1063  fBranchName = "LikeSign2Prong";
1064  break;
1065  case kDstartoKpipi:
1066  fCandidatePDG = 413;
1067  fCandidateName = "DStar";
1068  fNDaughters = 3;
1070  fPDGdaughters.Reset();
1071  fPDGdaughters[0] = 211; // pi soft
1072  fPDGdaughters[1] = 211; // pi fromD0
1073  fPDGdaughters[2] = 321; // K from D0
1074  fBranchName = "Dstar";
1076  break;
1077  default:
1078  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
1079  }
1080 
1082 }
1083 
1088 {
1089  if (fRDHFCuts) delete fRDHFCuts;
1090  fRDHFCuts = cuts;
1091 }
1092 
1097 {
1098  if (!cuts) return;
1099  if (fRDHFCuts) delete fRDHFCuts;
1100  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
1101 }
1102 
1107 {
1108  static TString name;
1109 
1110  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
1111 
1112  return name.Data();
1113 }
1114 
1119 {
1121  switch (fMCMode) {
1122  case kBackgroundOnly:
1123  fName += "_BackgroundOnly";
1124  break;
1125  case kSignalOnly:
1126  fName += "_SignalOnly";
1127  break;
1128  case kMCTruth:
1129  fName += "_MCTruth";
1130  break;
1131  case kD0Reflection:
1132  fName += "_D0Reflection";
1133  break;
1134  case kOnlyWrongPIDAccepted:
1135  fName += "_OnlyWrongPIDAccepted";
1136  break;
1137  default:
1138  break;
1139  }
1140 
1141  if (fRDHFCuts) fName += TString::Format("_%s", fRDHFCuts->GetName());
1142 
1143  return fName.Data();
1144 }
1145 
1153 {
1154  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
1155 
1156  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
1157  fJetDefinitions.push_back(def);
1158  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
1159  "Total number of jet definitions is now %lu.",
1160  def.GetName(), GetName(), fJetDefinitions.size());
1161  // For detector level set maximum track pt to 100 GeV/c
1162  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
1163  }
1164  else {
1165  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
1166  }
1167 
1168  return &(*it);
1169 }
1170 
1182 {
1183  AliHFJetDefinition def(type, r, algo, reco);
1184 
1185  return AddJetDefinition(def);
1186 }
1187 
1193 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
1194 {
1195  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
1196  while (it != fJetDefinitions.end() && (*it) != def) it++;
1197  return it;
1198 }
1199 
1204 {
1205  for (auto &jetdef : fJetDefinitions) jetdef.SetJetPhiRange(min, max);
1206 }
1207 
1212 {
1213  for (auto &jetdef : fJetDefinitions) jetdef.SetJetEtaRange(min, max);
1214 }
1215 
1220 {
1221  for (auto &jetdef : fJetDefinitions) jetdef.SetJetPtRange(min, max);
1222 }
1223 
1228 {
1229  for (auto &jetdef : fJetDefinitions) jetdef.SetChargedPtRange(min, max);
1230 }
1231 
1236 {
1237  for (auto &jetdef : fJetDefinitions) jetdef.SetNeutralPtRange(min, max);
1238 }
1239 
1246 {
1247  if (lhs.fCandidateType < rhs.fCandidateType) {
1248  return true;
1249  }
1250  else if (lhs.fCandidateType > rhs.fCandidateType) {
1251  return false;
1252  }
1253  else if (lhs.fMCMode < rhs.fMCMode) {
1254  return true;
1255  }
1256  else if (lhs.fMCMode > rhs.fMCMode) {
1257  return false;
1258  }
1259  else if (lhs.fRDHFCuts && !rhs.fRDHFCuts) {
1260  return true;
1261  }
1262  else if (lhs.fRDHFCuts && rhs.fRDHFCuts && strcmp(lhs.fRDHFCuts->GetName(), rhs.fRDHFCuts->GetName()) < 0) {
1263  return true;
1264  }
1265  else {
1266  return false;
1267  }
1268 }
1269 
1276 {
1277  if (lhs.fCandidateType != rhs.fCandidateType) return false;
1278  if (lhs.fMCMode != rhs.fMCMode) return false;
1279  if (lhs.fRDHFCuts == nullptr && rhs.fRDHFCuts != nullptr) return false;
1280  if (lhs.fRDHFCuts != nullptr && rhs.fRDHFCuts == nullptr) return false;
1281  if (lhs.fRDHFCuts && rhs.fRDHFCuts && strcmp(lhs.fRDHFCuts->GetName(), rhs.fRDHFCuts->GetName()) != 0) return false;
1282  return true;
1283 }
1284 
1293 {
1294  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) { // D0 candidate
1295  return ExtractD0Attributes(Dcand, DmesonJet, i);
1296  }
1297  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
1298  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
1299  }
1300  else {
1301  return kFALSE;
1302  }
1303 }
1304 
1313 {
1314  AliDebug(10,"Checking if D0 meson is selected");
1315  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
1316  if (isSelected == 0) return kFALSE;
1317 
1318  Int_t MCtruthPdgCode = 0;
1319 
1320  Double_t invMassD = 0;
1321 
1322  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
1323  // Checks also the origin, and if it matches the rejected origin mask, return false
1324  if (fMCMode != kNoMC) {
1325  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
1326  DmesonJet.fMCLabel = mcLab;
1327 
1328  // Retrieve the generated particle (if exists) and its PDG code
1329  if (mcLab >= 0) {
1330  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1331 
1332  if (aodMcPart) {
1333  // Check origin and return false if it matches the rejected origin mask
1334  if (fRejectedOrigin) {
1335  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
1336  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1337  }
1338  MCtruthPdgCode = aodMcPart->PdgCode();
1339  }
1340  }
1341  }
1342 
1343  if (isSelected == 1) { // selected as a D0
1344  if (i != 0) return kFALSE; // only one mass hypothesis thanks to PID
1345 
1346  if (fMCMode == kNoMC ||
1347  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1348  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly) ||
1349  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kD0Reflection || fMCMode == kOnlyWrongPIDAccepted))) {
1350  // 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)
1351  AliDebug(10,"Selected as D0");
1352  invMassD = Dcand->InvMassD0();
1353  }
1354  else { // conditions above not passed, so return FALSE
1355  return kFALSE;
1356  }
1357  }
1358  else if (isSelected == 2) { // selected as a D0bar
1359  if (i != 1) return kFALSE; // only one mass hypothesis thanks to PID
1360 
1361  if (fMCMode == kNoMC ||
1362  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1363  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly) ||
1364  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kD0Reflection || fMCMode == kOnlyWrongPIDAccepted))) {
1365  // 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)
1366  AliDebug(10,"Selected as D0bar");
1367  invMassD = Dcand->InvMassD0bar();
1368  }
1369  else { // conditions above not passed, so return FALSE
1370  return kFALSE;
1371  }
1372  }
1373  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
1374  AliDebug(10,"Selected as either D0 or D0bar");
1375 
1376  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
1377  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1378  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kD0Reflection))) {
1379  if (i != 0) return kFALSE;
1380  AliDebug(10, "MC truth is D0");
1381  invMassD = Dcand->InvMassD0();
1382  }
1383  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1384  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kD0Reflection))) {
1385  if (i != 1) return kFALSE;
1386  AliDebug(10, "MC truth is D0bar");
1387  invMassD = Dcand->InvMassD0bar();
1388  }
1389  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
1390 
1391  // Only accept it if background-only OR background-and-signal was requested
1392  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
1393  // Select D0 or D0bar depending on the i-parameter
1394  if (i == 0) {
1395  AliDebug(10, "Returning invariant mass with D0 hypothesis");
1396  invMassD = Dcand->InvMassD0();
1397  }
1398  else if (i == 1) {
1399  AliDebug(10, "Returning invariant mass with D0bar hypothesis");
1400  invMassD = Dcand->InvMassD0bar();
1401  }
1402  else { // i > 1
1403  return kFALSE;
1404  }
1405  }
1406  else { // signal-only was requested but this is not a true D0
1407  return kFALSE;
1408  }
1409  }
1410  }
1411  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
1412  return kTRUE;
1413 }
1414 
1423 {
1424  AliDebug(10,"Checking if D* meson is selected");
1425  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoCascadeHF*>(DstarCand), AliRDHFCuts::kAll, fAodEvent);
1426  if (isSelected == 0) return kFALSE;
1427 
1428  if ((i == 1 && DstarCand->Charge()>0) || (i == 0 && DstarCand->Charge()<0) || i > 1) return kFALSE; // only one mass hypothesis for the D*
1429 
1430  Int_t MCtruthPdgCode = 0;
1431 
1432  Double_t invMassD = 0;
1433 
1434  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1435  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1436  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1437 
1438  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1439  AliDebug(10, Form("MC label is %d", mcLab));
1440  DmesonJet.fMCLabel = mcLab;
1441  if (mcLab >= 0) {
1442  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1443 
1444  if (aodMcPart) {
1445  if (fRejectedOrigin) {
1446  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
1447  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1448  }
1449 
1450  MCtruthPdgCode = aodMcPart->PdgCode();
1451  AliDebug(10, Form("MC truth pdg code is %d",MCtruthPdgCode));
1452  }
1453  }
1454  }
1455 
1456  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1457  if (fMCMode == kNoMC ||
1458  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1459  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1460  // 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)
1461  invMassD = DstarCand->InvMassDstarKpipi();
1462  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1463  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1464  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1465  return kTRUE;
1466  }
1467  else { // conditions above not passed, so return FALSE
1468  return kFALSE;
1469  }
1470 }
1471 
1479 {
1480  if (!part) return kUnknownDecay;
1481  if (!mcArray) return kUnknownDecay;
1482 
1484 
1485  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1486 
1487  if (part->GetNDaughters() == 2) {
1488 
1489  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1490  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1491 
1492  if (!d1 || !d2) {
1493  return decay;
1494  }
1495 
1496  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1497  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1498 
1499  if (absPdgPart == 421) { // D0 -> K pi
1500  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1501  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1502  decay = kDecayD0toKpi;
1503  }
1504  }
1505 
1506  if (absPdgPart == 413) { // D* -> D0 pi
1507  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1508  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1509  if (D0decay == kDecayD0toKpi) {
1510  decay = kDecayDStartoKpipi;
1511  }
1512  }
1513  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1514  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1515  if (D0decay == kDecayD0toKpi) {
1516  decay = kDecayDStartoKpipi;
1517  }
1518  }
1519  }
1520  }
1521 
1522  return decay;
1523 }
1524 
1531 std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> AliAnalysisTaskDmesonJets::AnalysisEngine::IsPromptCharm(const AliAODMCParticle* part, TClonesArray* mcArray)
1532 {
1533  std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> result(kUnknownQuark, 0);
1534 
1535  if (!part) return result;
1536  if (!mcArray) return result;
1537 
1538  static std::set<UInt_t> partons = { 4, 5 };
1539 
1540  AliAODMCParticle* parton = FindParticleOrigin(part, mcArray, kFindLast, partons);
1541  if (parton) {
1542  result.second = parton;
1543  UInt_t absPdgParton = TMath::Abs(parton->GetPdgCode());
1544  if (absPdgParton == 4) result.first = kFromCharm;
1545  else if (absPdgParton == 5) result.first = kFromBottom;
1546  }
1547 
1548  return result;
1549 }
1550 
1558 
1559 AliAODMCParticle* AliAnalysisTaskDmesonJets::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode)
1560 {
1561  static std::set<UInt_t> pdgSet;
1562 
1563  return FindParticleOrigin(part, mcArray, mode, pdgSet);
1564 }
1565 
1574 
1575 AliAODMCParticle* AliAnalysisTaskDmesonJets::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode, const std::set<UInt_t>& pdgSet)
1576 {
1577  AliAODMCParticle* result = nullptr;
1578 
1579  Int_t mother = part->GetMother();
1580  while (mother >= 0) {
1581  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1582  if (mcGranma) {
1583  UInt_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
1584 
1585  // If the current particle is one of the particle types that is being searched assign it to the result pointer
1586  if (pdgSet.empty() || pdgSet.count(abspdgGranma) > 0) {
1587  result = mcGranma;
1588  // If the last particle in the fragmentation tree (first when going reverse) was requested then stop the loop
1589  if (mode == kFindLast) break;
1590  }
1591  if (mother == mcGranma->GetMother()) { // avoid infinite loop!
1592  AliWarningClassStream() << "Particle " << mother << " (PDG=" << mcGranma->PdgCode() << ") is the mother of itself!?" << std::endl;
1593  break;
1594  }
1595  mother = mcGranma->GetMother();
1596  }
1597  else {
1598  AliErrorClassStream() << "Could not retrieve mother particle " << mother << "!" << std::endl;
1599  break;
1600  }
1601  }
1602 
1603  return result;
1604 }
1605 
1608 {
1609  for (auto& jetDef : fJetDefinitions) {
1610  jetDef.fJets.clear();
1611  }
1612 
1613  if (fMCMode == kMCTruth) {
1615  }
1616  else {
1618  }
1619 }
1620 
1623 {
1624  // Fill the vertex info of the candidates
1625  // Needed for reduced delta AOD, where the vertex info has been deleted
1626  // to reduce the delta AOD file size
1628 
1629  const Int_t nD = fCandidateArray->GetEntriesFast();
1630 
1631  AliDmesonJetInfo DmesonJet;
1632  DmesonJet.fEvent = this->fAodEvent;
1633 
1634  std::map<AliHFJetDefinition*,Double_t> maxJetPt;
1635  for (auto& def : fJetDefinitions) maxJetPt[&def] = 0;
1636  Double_t maxDPt = 0;
1637 
1638  std::array<int, 3> nAccCharm = {0};
1639  std::array<std::array<int, 3>, 5> nAccCharmPt = {{{0}}};
1640 
1641  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1642  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1643  if (!charmCand) continue;
1644  if(!(vHF.FillRecoCand(fAodEvent,charmCand))) continue;
1645 
1646  // region of interest + cuts
1647  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1648  Int_t nMassHypo = 0; // number of mass hypothesis accepted for this D meson
1649  if (charmCand->Pt() > maxDPt) maxDPt = charmCand->Pt();
1650  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1651  DmesonJet.Reset();
1652  DmesonJet.fDmesonParticle = charmCand;
1653  DmesonJet.fSelectionType = im + 1;
1654  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1655  for (auto& def : fJetDefinitions) {
1656  if (FindJet(charmCand, DmesonJet, def)) {
1657  Double_t jetPt = DmesonJet.fJets[def.GetName()].fMomentum.Pt();
1658  if (jetPt > maxJetPt[&def]) maxJetPt[&def] = jetPt;
1659  }
1660  else {
1661  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1662  def.GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1663  }
1664  }
1665  fDmesonJets[(icharm+1)*(1-(im*2))] = DmesonJet;
1666  nMassHypo++;
1667  nAccCharm[im]++;
1668 
1669  for (int i = 0; i < nAccCharmPt.size(); i++) {
1670  if (charmCand->Pt() < i) break;
1671  nAccCharmPt[i][im]++;
1672  }
1673  }
1674  }
1675  if (nMassHypo == 2) { // both mass hypothesis accepted
1676  nAccCharm[0]--;
1677  nAccCharm[1]--;
1678  nAccCharm[2]++;
1679 
1680  for (int i = 0; i < nAccCharmPt.size(); i++) {
1681  if (charmCand->Pt() < i) break;
1682  nAccCharmPt[i][0]--;
1683  nAccCharmPt[i][1]--;
1684  nAccCharmPt[i][2]++;
1685  }
1686 
1687  fDmesonJets[(icharm+1)].fD0D0bar = kTRUE;
1688  fDmesonJets[-(icharm+1)].fD0D0bar = kTRUE;
1689  }
1690  } // end of D cand loop
1691 
1692  TString hname;
1693 
1694  Int_t ntracks = 0;
1695 
1696  for (auto track_cont : fTrackContainers) {
1697  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
1698  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(nullptr);
1699  ntracks += track_cont->GetNAcceptEntries();
1700  }
1701 
1702  for (auto& def : fJetDefinitions) {
1703  if (!def.fRho) continue;
1704  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
1705  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
1706 
1707  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
1708  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
1709 
1710  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
1711  fHistManager->FillTH2(hname, fCent, def.fRho->GetVal());
1712 
1713  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
1714  fHistManager->FillTH2(hname, fCent, maxJetPt[&def]);
1715 
1716  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
1717  fHistManager->FillTH2(hname, fCent, maxDPt);
1718 
1719  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
1720  fHistManager->FillTH2(hname, ntracks, def.fRho->GetVal());
1721 
1722  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
1723  fHistManager->FillTH2(hname, ntracks, maxJetPt[&def]);
1724 
1725  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
1726  fHistManager->FillTH2(hname, ntracks, maxDPt);
1727  }
1728 
1729  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1730  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1731  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1732  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1733 
1734  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1735  fHistManager->FillTH2(hname, ntracks, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1736 
1737  for (int i = 0; i < nAccCharmPt.size(); i++) {
1738  hname = TString::Format("%s/fHistNTotAcceptedDmesonsPt%d", GetName(), i);
1739  fHistManager->FillTH1(hname, "D", nAccCharmPt[i][0]);
1740  fHistManager->FillTH1(hname, "Anti-D", nAccCharmPt[i][1]);
1741  fHistManager->FillTH1(hname, "Both", nAccCharmPt[i][2]);
1742 
1743  hname = TString::Format("%s/fHistNAcceptedDmesonsPt%d", GetName(), i);
1744  fHistManager->FillTH1(hname, nAccCharmPt[i][0]+nAccCharmPt[i][1]+nAccCharmPt[i][2]);
1745  }
1746 
1747  hname = TString::Format("%s/fHistNDmesons", GetName());
1748  fHistManager->FillTH1(hname, nD);
1749 }
1750 
1761 {
1762  TString hname;
1763 
1764  Double_t rho = 0;
1765  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
1766 
1768  fFastJetWrapper->SetR(jetDef.fRadius);
1771 
1772  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1773 
1774  if (jetDef.fJetType != AliJetContainer::kNeutralJet) {
1775  for (auto track_cont : fTrackContainers) {
1776  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
1777  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(Dcand);
1778  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1779  AddInputVectors(track_cont, 100, static_cast<TH2*>(fHistManager->FindObject(hname)), fTrackEfficiency);
1780 
1781  if (hftrack_cont) {
1782  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1783  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1784  const TObjArray& daughters = hftrack_cont->GetDaughterList();
1785  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1786  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1787  if (!hftrack_cont->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1788  }
1789  }
1790  }
1791  }
1792 
1793  if (jetDef.fJetType != AliJetContainer::kChargedJet) {
1794  for (auto clus_cont : fClusterContainers) {
1795  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1796  AddInputVectors(clus_cont, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1797  }
1798  }
1799 
1800  // run jet finder
1801  fFastJetWrapper->Run();
1802 
1803  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1804 
1805  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1806  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1807 
1808  Bool_t isDmesonJet = kFALSE;
1809 
1810  Double_t maxChPt = 0;
1811  Double_t maxNePt = 0;
1812  Double_t totalNeutralPt = 0;
1813  Int_t nConst = 1;
1814 
1815  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1816  if (constituents[ic].user_index() == 0) {
1817  isDmesonJet = kTRUE;
1818  }
1819  else if (constituents[ic].user_index() >= 100) {
1820  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1821  nConst++;
1822  }
1823  else if (constituents[ic].user_index() <= -100) {
1824  totalNeutralPt += constituents[ic].pt();
1825  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1826  nConst++;
1827  }
1828  }
1829 
1830  if (isDmesonJet) {
1831  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1832  DmesonJet.fJets[jetDef.GetName()].fNConstituents = nConst;
1833  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1834  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1835  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1836  DmesonJet.fJets[jetDef.GetName()].fArea = jets_incl[ijet].area();
1837  DmesonJet.fJets[jetDef.GetName()].fCorrPt = DmesonJet.fJets[jetDef.GetName()].fMomentum.Pt() - jets_incl[ijet].area() * rho;
1838 
1839  return kTRUE;
1840  }
1841  }
1842 
1843  return kFALSE;
1844 }
1845 
1849 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist, Double_t eff)
1850 {
1851  auto itcont = cont->all_momentum();
1852  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1853  UInt_t rejectionReason = 0;
1854  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1855  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1856  continue;
1857  }
1858  if (fRandomGen && eff > 0 && eff < 1) {
1859  Double_t rnd = fRandomGen->Rndm();
1860  if (eff < rnd) {
1861  if (rejectHist) rejectHist->Fill(6, it->first.Pt());
1862  continue;
1863  }
1864  }
1865  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1866  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1867  }
1868 }
1869 
1872 {
1873  TString hname;
1874 
1879 
1880  if (!fMCContainer->IsSpecialPDGFound()) return;
1881 
1882  std::array<int,2> nAccCharm = {0};
1883  std::array<std::array<int, 2>, 5> nAccCharmPt = {{{0}}};
1884 
1885  std::map<AliHFJetDefinition*, Double_t> maxJetPt;
1886  Double_t maxDPt = 0;
1887 
1888  for (auto &jetDef : fJetDefinitions) {
1889  maxJetPt[&jetDef] = 0;
1890  Double_t rho = 0;
1891  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
1892  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", GetName(), jetDef.GetName());
1893  TH1* histNDmesonsVsNconstituents = static_cast<TH1*>(fHistManager->FindObject(hname));
1894 
1895  switch (jetDef.fJetType) {
1897  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
1898  break;
1900  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
1901  break;
1903  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
1904  break;
1905  }
1906 
1908  fFastJetWrapper->SetR(jetDef.fRadius);
1911 
1912  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1913  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1914 
1915  fFastJetWrapper->Run();
1916 
1917  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1918 
1919  for (auto jet : jets_incl) {
1920  Int_t nDmesonsInJet = 0;
1921 
1922  for (auto constituent : jet.constituents()) {
1923  Int_t iPart = constituent.user_index() - 100;
1924  if (constituent.perp() < 1e-6) continue; // reject ghost particles
1925  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
1926  if (!part) {
1927  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1928  continue;
1929  }
1930  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1931  nDmesonsInJet++;
1932  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1933  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1934  if (part->Pt() > maxDPt) maxDPt = part->Pt();
1935  std::pair<int, AliDmesonJetInfo> element;
1936  element.first = iPart;
1937  dMesonJetIt = fDmesonJets.insert(element).first;
1938  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1939  (*dMesonJetIt).second.fDmesonParticle = part;
1940  (*dMesonJetIt).second.fSelectionType = part->PdgCode() > 0 ? 1 : 2;
1941 
1942  UShort_t p = 0;
1943  UInt_t rs = 0;
1944 
1945  auto origin = IsPromptCharm(part, fMCContainer->GetArray());
1946  p = 0;
1947  rs = origin.first;
1948  while (rs >>= 1) { p++; }
1949  (*dMesonJetIt).second.fPartonType = p;
1950  (*dMesonJetIt).second.fParton = origin.second;
1951 
1952  (*dMesonJetIt).second.fAncestor = FindParticleOrigin(part, fMCContainer->GetArray(), kFindFirst);
1953 
1954  Int_t im = -1;
1955  if (part->PdgCode() > 0) { // D0
1956  im = 0;
1957  }
1958  else { // D0bar
1959  im = 1;
1960  }
1961 
1962  nAccCharm[im]++;
1963  for (int i = 0; i < nAccCharmPt.size(); i++) {
1964  if (part->Pt() < i) break;
1965  nAccCharmPt[i][im]++;
1966  }
1967  }
1968 
1969  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
1970  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = jet.constituents().size();
1971  (*dMesonJetIt).second.fJets[jetDef.GetName()].fArea = jet.area();
1972  (*dMesonJetIt).second.fJets[jetDef.GetName()].fCorrPt = (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.Pt() - jet.area() * rho;
1973  if (jet.perp() > maxJetPt[&jetDef]) maxJetPt[&jetDef] = jet.perp();
1974  } // if constituent is a D meson
1975  } // for each constituent
1976  if (nDmesonsInJet > 0) histNDmesonsVsNconstituents->Fill(jet.constituents().size(), nDmesonsInJet);
1977  } // for each jet
1978  } // for each jet definition
1979 
1981 
1982  for (auto& def : fJetDefinitions) {
1983  if (!def.fRho) continue;
1984  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
1985  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
1986 
1987  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
1988  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
1989 
1990  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
1991  fHistManager->FillTH2(hname, fCent, def.fRho->GetVal());
1992 
1993  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
1994  fHistManager->FillTH2(hname, fCent, maxJetPt[&def]);
1995 
1996  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
1997  fHistManager->FillTH2(hname, fCent, maxDPt);
1998 
1999  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
2000  fHistManager->FillTH2(hname, npart, def.fRho->GetVal());
2001 
2002  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
2003  fHistManager->FillTH2(hname, npart, maxJetPt[&def]);
2004 
2005  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
2006  fHistManager->FillTH2(hname, npart, maxDPt);
2007  }
2008 
2009  if (fDmesonJets.size() != nAccCharm[0]+nAccCharm[1]) AliError(Form("I found %lu mesons (%d)?", fDmesonJets.size(), nAccCharm[0]+nAccCharm[1]));
2010  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
2011  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
2012  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
2013 
2014  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
2015  fHistManager->FillTH2(hname, npart, nAccCharm[0]+nAccCharm[1]);
2016 
2017  for (int i = 0; i < nAccCharmPt.size(); i++) {
2018  hname = TString::Format("%s/fHistNTotAcceptedDmesonsPt%d", GetName(), i);
2019  fHistManager->FillTH1(hname, "D", nAccCharmPt[i][0]);
2020  fHistManager->FillTH1(hname, "Anti-D", nAccCharmPt[i][1]);
2021 
2022  hname = TString::Format("%s/fHistNAcceptedDmesonsPt%d", GetName(), i);
2023  fHistManager->FillTH1(hname, nAccCharmPt[i][0]+nAccCharmPt[i][1]);
2024  }
2025 
2026  hname = TString::Format("%s/fHistNDmesons", GetName());
2027  fHistManager->FillTH1(hname, nAccCharm[0]+nAccCharm[1]); // same as the number of accepted D mesons, since no selection is performed
2028 }
2029 
2034 {
2035  TString classname;
2036  if (fMCMode == kMCTruth) {
2037  classname = "AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary";
2039  }
2040  else {
2041  switch (fCandidateType) {
2042  case kD0toKpi:
2043  case kD0toKpiLikeSign:
2044  if (fD0Extended) {
2045  classname = "AliAnalysisTaskDmesonJets::AliD0ExtendedInfoSummary";
2047  }
2048  else {
2049  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
2051  }
2052  break;
2053  case kDstartoKpipi:
2054  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
2056  break;
2057  }
2058  }
2059  TString treeName = TString::Format("%s_%s", taskName, GetName());
2060  fTree = new TTree(treeName, treeName);
2061  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
2063  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
2064  if (fJetDefinitions[i].fRhoName.IsNull()) {
2066  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
2067  }
2068  else {
2070  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoPbPbSummary", &fCurrentJetInfo[i]);
2071  }
2072  }
2073 
2074  return fTree;
2075 }
2076 
2081 {
2082  TString hname;
2083 
2084  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
2085 
2086  for (auto &jetDef : fJetDefinitions) {
2087 
2088  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
2089 
2090  Double_t radius = jetDef.fRadius;
2091 
2092  TString title[30] = {""};
2093  Int_t nbins[30] = {0 };
2094  Double_t min [30] = {0.};
2095  Double_t max [30] = {0.};
2096  Int_t dim = 0 ;
2097 
2098  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
2099  nbins[dim] = nPtBins;
2100  min[dim] = 0;
2101  max[dim] = fMaxPt;
2102  dim++;
2103 
2104  if ((enabledAxis & kPositionD) != 0) {
2105  title[dim] = "#eta_{D}";
2106  nbins[dim] = 50;
2107  min[dim] = -1;
2108  max[dim] = 1;
2109  dim++;
2110 
2111  title[dim] = "#phi_{D} (rad)";
2112  nbins[dim] = 150;
2113  min[dim] = 0;
2114  max[dim] = TMath::TwoPi();
2115  dim++;
2116  }
2117 
2118  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
2119  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
2120  nbins[dim] = fNMassBins;
2121  min[dim] = fMinMass;
2122  max[dim] = fMaxMass;
2123  dim++;
2124  }
2125 
2127  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
2128  nbins[dim] = fNMassBins;
2129  min[dim] = fMinMass;
2130  max[dim] = fMaxMass;
2131  dim++;
2132  }
2133 
2134  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
2135  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
2136  nbins[dim] = fNMassBins;
2137  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
2138  dim++;
2139  }
2140 
2141  if (fCandidateType == kDstartoKpipi) {
2142  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
2143  nbins[dim] = fNMassBins*6;
2144  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
2145 
2146  // subtract mass of D0
2147  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2148  min[dim] -= D0mass;
2149  max[dim] -= D0mass;
2150 
2151  dim++;
2152  }
2153 
2154  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
2155  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
2156  nbins[dim] = 100;
2157  min[dim] = 0;
2158  max[dim] = 25;
2159  dim++;
2160  }
2161 
2162  title[dim] = "#it{z}_{D}";
2163  nbins[dim] = 110;
2164  min[dim] = 0;
2165  max[dim] = 1.10;
2166  dim++;
2167 
2168  if ((enabledAxis & kDeltaR) != 0) {
2169  title[dim] = "#Delta R_{D-jet}";
2170  nbins[dim] = 100;
2171  min[dim] = 0;
2172  max[dim] = radius * 1.5;
2173  dim++;
2174  }
2175 
2176  if ((enabledAxis & kDeltaEta) != 0) {
2177  title[dim] = "#eta_{D} - #eta_{jet}";
2178  nbins[dim] = 100;
2179  min[dim] = -radius * 1.2;
2180  max[dim] = radius * 1.2;
2181  dim++;
2182  }
2183 
2184  if ((enabledAxis & kDeltaPhi) != 0) {
2185  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
2186  nbins[dim] = 100;
2187  min[dim] = -radius * 1.2;
2188  max[dim] = radius * 1.2;
2189  dim++;
2190  }
2191 
2192  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
2193  nbins[dim] = nPtBins;
2194  min[dim] = 0;
2195  max[dim] = fMaxPt;
2196  dim++;
2197 
2198  if ((enabledAxis & kPositionJet) != 0) {
2199  title[dim] = "#eta_{jet}";
2200  nbins[dim] = 50;
2201  min[dim] = -1;
2202  max[dim] = 1;
2203  dim++;
2204 
2205  title[dim] = "#phi_{jet} (rad)";
2206  nbins[dim] = 150;
2207  min[dim] = 0;
2208  max[dim] = TMath::TwoPi();
2209  dim++;
2210  }
2211 
2212  if ((enabledAxis & kJetConstituents) != 0) {
2213  title[dim] = "No. of constituents";
2214  nbins[dim] = 50;
2215  min[dim] = -0.5;
2216  max[dim] = 49.5;
2217  dim++;
2218  }
2219 
2220  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
2221  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
2222  for (Int_t j = 0; j < dim; j++) {
2223  h->GetAxis(j)->SetTitle(title[j]);
2224  }
2225  }
2226 }
2227 
2232 {
2233  TString hname;
2234  fPartons.clear();
2235 
2236  TH1* histAncestor = nullptr;
2237  TH1* histPrompt = nullptr;
2238 
2239  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2240  hname = TString::Format("%s/fHistPrompt", GetName());
2241  histPrompt = static_cast<TH1*>(fHistManager->FindObject(hname));
2242 
2243  hname = TString::Format("%s/fHistAncestor", GetName());
2244  histAncestor = static_cast<TH1*>(fHistManager->FindObject(hname));
2245  }
2246 
2247  for (auto& dmeson_pair : fDmesonJets) {
2248  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
2249  Int_t accJets = 0;
2250  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
2251  fCurrentJetInfo[ij]->Reset();
2252  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
2253  if (!jet) continue;
2254  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
2255  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
2256  fHistManager->FillTH1(hname, jet->Pt());
2257  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
2258  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2259  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
2260  fHistManager->FillTH1(hname, jet->Eta());
2261  continue;
2262  }
2263  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
2264  accJets++;
2265  }
2266  if (accJets > 0) {
2267  if (histPrompt) {
2268  if (dmeson_pair.second.fParton) {
2269  fPartons[dmeson_pair.second.fParton] = dmeson_pair.second.fPartonType;
2270  UInt_t absPdgParton = TMath::Abs(dmeson_pair.second.fParton->GetPdgCode());
2271  if (absPdgParton == 4) {
2272  histPrompt->Fill("Prompt", 1);
2273  }
2274  else if (absPdgParton == 5) {
2275  histPrompt->Fill("Non-Prompt", 1);
2276  }
2277  else {
2278  histPrompt->Fill("Unknown", 1);
2279  }
2280  }
2281  else {
2282  histPrompt->Fill("Unknown", 1);
2283  }
2284  }
2285 
2286  if (histAncestor) {
2287  if (dmeson_pair.second.fAncestor) {
2288  UInt_t absPdgAncestor = TMath::Abs(dmeson_pair.second.fAncestor->GetPdgCode());
2289  if (absPdgAncestor == 4) {
2290  histAncestor->Fill("Charm", 1);
2291  }
2292  else if (absPdgAncestor == 5) {
2293  histAncestor->Fill("Bottom", 1);
2294  }
2295  else if (absPdgAncestor == 2212) {
2296  histAncestor->Fill("Proton", 1);
2297  }
2298  else {
2299  histAncestor->Fill("Unknown", 1);
2300  }
2301  }
2302  else {
2303  histAncestor->Fill("Unknown", 1);
2304  }
2305  }
2306 
2307  fTree->Fill();
2308  }
2309  else {
2310  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2311  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2312  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2313  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2314  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2315  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2316  if (fMCMode != kMCTruth) {
2318  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
2319  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
2320  }
2321  else if (fCandidateType == kDstartoKpipi) {
2322  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
2323  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
2324 
2325  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
2326  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
2327  }
2328  }
2329  }
2330  }
2331 
2332  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2333  hname = TString::Format("%s/fHistPartonPt", GetName());
2334  TH1* histPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
2335  hname = TString::Format("%s/fHistPartonEta", GetName());
2336  TH1* histPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
2337  hname = TString::Format("%s/fHistPartonPhi", GetName());
2338  TH1* histPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
2339  hname = TString::Format("%s/fHistPartonType", GetName());
2340  TH1* histPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
2341 
2342  for (auto parton : fPartons) {
2343  if (!parton.first) continue;
2344  histPartonPt->Fill(parton.first->Pt());
2345  histPartonEta->Fill(parton.first->Eta());
2346  histPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
2347  histPartonType->Fill(parton.second);
2348  }
2349  }
2350 
2351  return kTRUE;
2352 }
2353 
2359 {
2360  TString hname;
2361 
2362  TH1* histAncestor = nullptr;
2363  TH1* histPrompt = nullptr;
2364 
2365  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2366  hname = TString::Format("%s/fHistPrompt", GetName());
2367  histPrompt = static_cast<TH1*>(fHistManager->FindObject(hname));
2368 
2369  hname = TString::Format("%s/fHistAncestor", GetName());
2370  histAncestor = static_cast<TH1*>(fHistManager->FindObject(hname));
2371  }
2372 
2373  fPartons.clear();
2374  for (auto& dmeson_pair : fDmesonJets) {
2375  Int_t accJets = 0;
2376  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
2377  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
2378  if (!jet) continue;
2379  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
2380  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
2381  fHistManager->FillTH1(hname, jet->Pt());
2382  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
2383  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2384  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
2385  fHistManager->FillTH1(hname, jet->Eta());
2386  continue;
2387  }
2388  accJets++;
2389  }
2390  if (accJets > 0) {
2391  if (histPrompt) {
2392  if (dmeson_pair.second.fParton) {
2393  fPartons[dmeson_pair.second.fParton] = dmeson_pair.second.fPartonType;
2394  UInt_t absPdgParton = TMath::Abs(dmeson_pair.second.fParton->GetPdgCode());
2395  if (absPdgParton == 4) {
2396  histPrompt->Fill("Prompt", 1);
2397  }
2398  else if (absPdgParton == 5) {
2399  histPrompt->Fill("Non-Prompt", 1);
2400  }
2401  else {
2402  histPrompt->Fill("Unknown", 1);
2403  }
2404  }
2405  else {
2406  histPrompt->Fill("Unknown", 1);
2407  }
2408  }
2409 
2410  if (histAncestor) {
2411  if (dmeson_pair.second.fAncestor) {
2412  UInt_t absPdgAncestor = TMath::Abs(dmeson_pair.second.fAncestor->GetPdgCode());
2413  if (absPdgAncestor == 4) {
2414  histAncestor->Fill("Charm", 1);
2415  }
2416  else if (absPdgAncestor == 5) {
2417  histAncestor->Fill("Bottom", 1);
2418  }
2419  else if (absPdgAncestor == 2212) {
2420  histAncestor->Fill("Proton", 1);
2421  }
2422  else {
2423  histAncestor->Fill("Unknown", 1);
2424  }
2425  }
2426  else {
2427  histAncestor->Fill("Unknown", 1);
2428  }
2429  }
2430  }
2431  else {
2432  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2433  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2434  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2435  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2436  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2437  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2438  if (fMCMode != kMCTruth) {
2440  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
2441  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
2442  }
2443  else if (fCandidateType == kDstartoKpipi) {
2444  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
2445  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
2446 
2447  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
2448  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
2449  }
2450  }
2451  }
2452  }
2453 
2454  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
2455  hname = TString::Format("%s/fHistPartonPt", GetName());
2456  TH1* histPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
2457  hname = TString::Format("%s/fHistPartonEta", GetName());
2458  TH1* histPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
2459  hname = TString::Format("%s/fHistPartonPhi", GetName());
2460  TH1* histPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
2461  hname = TString::Format("%s/fHistPartonType", GetName());
2462  TH1* histPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
2463 
2464  for (auto parton : fPartons) {
2465  if (!parton.first) continue;
2466  histPartonPt->Fill(parton.first->Pt());
2467  histPartonEta->Fill(parton.first->Eta());
2468  histPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
2469  histPartonType->Fill(parton.second);
2470  }
2471  }
2472 
2473  return kTRUE;
2474 }
2475 
2480 {
2481  TString hname;
2482 
2483  for (auto& dmeson_pair : fDmesonJets) {
2484  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
2485  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2486  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2487  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2488  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2489  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2490  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2491  }
2492  }
2493 
2494  for (auto &jetDef : fJetDefinitions) {
2495 
2496  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
2497  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
2498 
2499  for (auto& dmeson_pair : fDmesonJets) {
2500  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
2501  if (!jet) continue;
2502  if (!jetDef.IsJetInAcceptance(*jet)) {
2503  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
2504  fHistManager->FillTH1(hname, jet->Pt());
2505  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
2506  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2507  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
2508  fHistManager->FillTH1(hname, jet->Eta());
2509  continue;
2510  }
2511  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
2512  }
2513  }
2514 
2515  return kTRUE;
2516 }
2517 
2524 {
2525  // Fill the THnSparse histogram.
2526 
2527  Double_t contents[30] = {0.};
2528 
2529  Double_t z = DmesonJet.GetZ(n);
2530  Double_t deltaPhi = 0;
2531  Double_t deltaEta = 0;
2532  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
2533 
2534  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
2535  if (it == DmesonJet.fJets.end()) return kFALSE;
2536 
2537  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
2538  TString title(h->GetAxis(i)->GetTitle());
2539  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
2540  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
2541  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
2542  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
2543  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
2544  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
2545  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
2546  else if (title=="#it{z}_{D}") contents[i] = z;
2547  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
2548  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
2549  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
2550  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
2551  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
2552  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
2553  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
2554  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
2555  }
2556 
2557  h->Fill(contents);
2558 
2559  return kTRUE;
2560 }
2561 
2562 // Definitions of class AliAnalysisTaskDmesonJets
2563 
2565 ClassImp(AliAnalysisTaskDmesonJets);
2567 
2571  fAnalysisEngines(),
2572  fEnabledAxis(0),
2574  fHistManager(),
2575  fApplyKinematicCuts(kTRUE),
2576  fNOutputTrees(0),
2577  fTrackEfficiency(0),
2578  fRejectISR(kFALSE),
2579  fJetAreaType(fastjet::active_area),
2580  fJetGhostArea(0.005),
2581  fMCContainer(0),
2582  fAodEvent(0),
2583  fFastJetWrapper(0)
2584 {
2585  SetMakeGeneralHistograms(kTRUE);
2586 }
2587 
2592  AliAnalysisTaskEmcalLight(name, kTRUE),
2593  fAnalysisEngines(),
2596  fHistManager(name),
2597  fApplyKinematicCuts(kTRUE),
2598  fNOutputTrees(nOutputTrees),
2599  fTrackEfficiency(0),
2600  fRejectISR(kFALSE),
2601  fJetAreaType(fastjet::active_area),
2602  fJetGhostArea(0.005),
2603  fMCContainer(0),
2604  fAodEvent(0),
2605  fFastJetWrapper(0)
2606 {
2607  SetMakeGeneralHistograms(kTRUE);
2608  for (Int_t i = 0; i < nOutputTrees; i++){
2609  DefineOutput(2+i, TTree::Class());
2610  }
2611 }
2612 
2615 {
2616  if (fFastJetWrapper) delete fFastJetWrapper;
2617 }
2618 
2626 {
2627  AliRDHFCuts* analysiscuts = 0;
2628  TFile* filecuts = TFile::Open(cutfname);
2629  if (!filecuts || filecuts->IsZombie()) {
2630  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Input file not found: will use std cuts.");
2631  filecuts = 0;
2632  }
2633 
2634  if (filecuts) analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
2635 
2636  if (!analysiscuts) {
2637  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Could not find analysis cuts '%s' in '%s'.", cutsname.Data(), cutfname.Data());
2638  if (filecuts) {
2639  filecuts->ls();
2640  }
2641  }
2642  else {
2643  ::Info("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Cuts '%s' loaded from file '%s'", cutsname.Data(), cutfname.Data());
2644  }
2645 
2646  return analysiscuts;
2647 }
2648 
2661 {
2663  return AddAnalysisEngine(type, cutfname, cuttype, MCmode, jetDef, rhoName);
2664 }
2665 
2677 {
2678  AliRDHFCuts* cuts = 0;
2679 
2680  if (!cutfname.IsNull()) {
2681  TString cutsname;
2682 
2683  switch (type) {
2684  case kD0toKpi:
2685  case kD0toKpiLikeSign:
2686  cutsname = "D0toKpiCuts";
2687  break;
2688  case kDstartoKpipi:
2689  cutsname = "DStartoKpipiCuts";
2690  break;
2691  default:
2692  return 0;
2693  }
2694 
2695  if (!cuttype.IsNull()) {
2696  cutsname += TString::Format("_%s", cuttype.Data());
2697  }
2698 
2699  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
2700  if (cuts) cuts->PrintAll();
2701  }
2702 
2703  AnalysisEngine eng(type, MCmode, cuts);
2704 
2705  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
2706 
2707  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
2708  eng.AddJetDefinition(jetDef);
2709  it = fAnalysisEngines.insert(it, eng);
2710  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
2711  }
2712  else {
2713  AnalysisEngine* found_eng = &(*it);
2714  ::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());
2715  found_eng->AddJetDefinition(jetDef);
2716 
2717  if (cuts) {
2718  if (found_eng->fRDHFCuts != 0) ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
2719  found_eng->SetRDHFCuts(cuts);
2720  }
2721  }
2722 
2723  return &(*it);
2724 }
2725 
2726 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
2727 {
2728  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
2729  while (it != fAnalysisEngines.end() && (*it) != eng) it++;
2730  return it;
2731 }
2732 
2735 {
2736  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
2737 
2739 
2740  // Define histograms
2741  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
2742 
2743  TString hname;
2744  TString htitle;
2745  TH1* h = 0;
2746  Int_t treeSlot = 0;
2747 
2748  Int_t maxTracks = 6000;
2749  Double_t maxRho = 500;
2750  if (fForceBeamType == kpp) {
2751  maxRho = 50;
2752  maxTracks = 200;
2753  }
2754  else if (fForceBeamType == kpA) {
2755  maxRho = 200;
2756  maxTracks = 500;
2757  }
2758 
2759  hname = "fHistCharmPt";
2760  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2761  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2762 
2763  hname = "fHistCharmEta";
2764  htitle = hname + ";#eta_{charm};counts";
2765  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2766 
2767  hname = "fHistCharmPhi";
2768  htitle = hname + ";#phi_{charm};counts";
2769  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2770 
2771  hname = "fHistBottomPt";
2772  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2773  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2774 
2775  hname = "fHistBottomEta";
2776  htitle = hname + ";#eta_{bottom};counts";
2777  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2778 
2779  hname = "fHistBottomPhi";
2780  htitle = hname + ";#phi_{bottom};counts";
2781  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2782 
2783  hname = "fHistHighestPartonPt";
2784  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2785  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2786 
2787  hname = "fHistHighestPartonType";
2788  htitle = hname + ";type;counts";
2789  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2790 
2791  hname = "fHistNHeavyQuarks";
2792  htitle = hname + ";number of heavy-quarks;counts";
2793  fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
2794 
2795  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
2796  for (auto &param : fAnalysisEngines) {
2797  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
2798 
2799  param.fHistManager = &fHistManager;
2800 
2801  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", param.GetName());
2802  htitle = hname + ";#it{N}_{tracks};#it{N}_{D};events";
2803  h = fHistManager.CreateTH2(hname, htitle, 251, -0.5, 250.5, 21, -0.5, 20.5);
2804 
2805  hname = TString::Format("%s/fHistNTotAcceptedDmesons", param.GetName());
2806  htitle = hname + ";;#it{N}_{D}";
2807  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
2808 
2809  for (int i = 0 ; i < 5; i++) {
2810  hname = TString::Format("%s/fHistNAcceptedDmesonsPt%d", param.GetName(), i);
2811  htitle = hname + ";#it{N}_{D};events";
2812  h = fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
2813 
2814  hname = TString::Format("%s/fHistNTotAcceptedDmesonsPt%d", param.GetName(), i);
2815  htitle = hname + ";;#it{N}_{D}";
2816  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
2817  }
2818 
2819  hname = TString::Format("%s/fHistNDmesons", param.GetName());
2820  htitle = hname + ";#it{N}_{D};events";
2821  h = fHistManager.CreateTH1(hname, htitle, 501, -0.5, 500.5);
2822 
2823  hname = TString::Format("%s/fHistNEvents", param.GetName());
2824  htitle = hname + ";Event status;counts";
2825  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
2826  h->GetXaxis()->SetBinLabel(1, "Accepted");
2827  h->GetXaxis()->SetBinLabel(2, "Rejected");
2828 
2829  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
2830  htitle = hname + ";Rejection reason;counts";
2831  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
2832 
2833  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
2834  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
2835  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2836 
2837  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
2838  htitle = hname + ";#it{#eta}_{D};counts";
2839  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2840 
2841  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
2842  htitle = hname + ";#it{#phi}_{D};counts";
2843  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2844 
2845  if (param.fMCMode != kMCTruth) {
2846  if (param.fCandidateType == kD0toKpi || param.fCandidateType == kD0toKpiLikeSign) {
2847  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
2848  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2849  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
2850  }
2851  else if (param.fCandidateType == kDstartoKpipi) {
2852  Double_t min = 0;
2853  Double_t max = 0;
2854 
2855  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
2856  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2857  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
2858  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
2859 
2860  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2861  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
2862  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2863  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
2864  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
2865  }
2866  }
2867 
2868  if (param.fMCMode == kMCTruth) {
2869  hname = TString::Format("%s/fHistPartonPt", param.GetName());
2870  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
2871  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2872 
2873  hname = TString::Format("%s/fHistPartonEta", param.GetName());
2874  htitle = hname + ";#eta_{parton};counts";
2875  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2876 
2877  hname = TString::Format("%s/fHistPartonPhi", param.GetName());
2878  htitle = hname + ";#phi_{parton};counts";
2879  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2880 
2881  hname = TString::Format("%s/fHistPartonType", param.GetName());
2882  htitle = hname + ";type;counts";
2883  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2884 
2885  hname = TString::Format("%s/fHistPrompt", param.GetName());
2886  htitle = hname + ";Type;counts";
2887  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
2888  h->GetXaxis()->SetBinLabel(1, "Unknown");
2889  h->GetXaxis()->SetBinLabel(2, "Prompt");
2890  h->GetXaxis()->SetBinLabel(3, "Non-Prompt");
2891 
2892  hname = TString::Format("%s/fHistAncestor", param.GetName());
2893  htitle = hname + ";Ancestor;counts";
2894  h = fHistManager.CreateTH1(hname, htitle, 4, 0, 4);
2895  h->GetXaxis()->SetBinLabel(1, "Unknown");
2896  h->GetXaxis()->SetBinLabel(2, "Charm");
2897  h->GetXaxis()->SetBinLabel(3, "Bottom");
2898  h->GetXaxis()->SetBinLabel(4, "Proton");
2899  }
2900 
2901  for (auto& jetDef : param.fJetDefinitions) {
2902  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef.GetName());
2903 
2904  if (param.fMCMode == kMCTruth) {
2905  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", param.GetName(), jetDef.GetName());
2906  htitle = hname + ";#it{N}_{constituents};#it{N}_{D};counts";
2907  h = fHistManager.CreateTH2(hname, htitle, 51, -0.5, 50.5, 10, 0.5, 10.5);
2908  }
2909 
2910  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef.GetName());
2911  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2912  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2913  SetRejectionReasonLabels(h->GetXaxis());
2914 
2915  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef.GetName());
2916  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2917  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2918  SetRejectionReasonLabels(h->GetXaxis());
2919 
2920  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef.GetName());
2921  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
2922  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2923  SetRejectionReasonLabels(h->GetXaxis());
2924 
2925  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef.GetName());
2926  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
2927  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
2928 
2929  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef.GetName());
2930  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
2931  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2932 
2933  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef.GetName());
2934  htitle = hname + ";#it{#eta}_{jet};counts";
2935  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2936 
2937  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef.GetName());
2938  htitle = hname + ";#it{#phi}_{jet};counts";
2939  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2940 
2941  if (!jetDef.fRhoName.IsNull()) {
2942  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", param.GetName(), jetDef.GetName());
2943  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});#rho (GeV/#it{c} #times rad^{-1});counts";
2944  fHistManager.CreateTH2(hname, htitle, 300, 0, 150, 1000, 0, maxRho);
2945 
2946  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", param.GetName(), jetDef.GetName());
2947  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});#rho (GeV/#it{c} #times rad^{-1});counts";
2948  fHistManager.CreateTH2(hname, htitle, 300, 0, 150, 1000, 0, maxRho);
2949 
2950  hname = TString::Format("%s/%s/fHistRhoVsCent", param.GetName(), jetDef.GetName());
2951  htitle = hname + ";Centrality (%);#rho (GeV/#it{c} #times rad^{-1});counts";
2952  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 1000, 0, maxRho);
2953 
2954  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", param.GetName(), jetDef.GetName());
2955  htitle = hname + ";Centrality (%);#it{p}_{T,jet} (GeV/#it{c});counts";
2956  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
2957 
2958  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", param.GetName(), jetDef.GetName());
2959  htitle = hname + ";Centrality (%);#it{p}_{T,D} (GeV/#it{c});counts";
2960  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
2961 
2962  hname = TString::Format("%s/%s/fHistRhoVsNTracks", param.GetName(), jetDef.GetName());
2963  htitle = hname + ";no. of tracks;#rho (GeV/#it{c} #times rad^{-1});counts";
2964  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 1000, 0, maxRho);
2965 
2966  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", param.GetName(), jetDef.GetName());
2967  htitle = hname + ";no. of tracks;#it{p}_{T,jet} (GeV/#it{c});counts";
2968  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
2969 
2970  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", param.GetName(), jetDef.GetName());
2971  htitle = hname + ";no. of tracks;#it{p}_{T,D} (GeV/#it{c});counts";
2972  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
2973  }
2974  }
2975  switch (fOutputType) {
2976  case kTreeOutput:
2977  param.BuildTree(GetName());
2978  if (treeSlot < fNOutputTrees) {
2979  param.AssignDataSlot(treeSlot+2);
2980  treeSlot++;
2982  }
2983  else {
2984  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
2985  }
2986  break;
2987  case kTHnOutput:
2988  param.BuildHnSparse(fEnabledAxis);
2989  break;
2990  case kNoOutput:
2991  break;
2992  }
2993  }
2994 
2996 
2997  PostData(1, fOutput);
2998 }
2999 
3003 {
3005 
3006  // Load the event
3007  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
3008 
3009  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
3010 
3011  fFastJetWrapper->SetAreaType((fastjet::AreaType)fJetAreaType);
3013 
3014  if (!fAodEvent) {
3015  AliError(Form("This task need an AOD event (Task '%s'). Expect troubles...", GetName()));
3016  //return;
3017  }
3018 
3019  TRandom* rnd = 0;
3020  if (fTrackEfficiency > 0 && fTrackEfficiency < 1) rnd = new TRandom3(0);
3021 
3022  for (auto cont_it : fParticleCollArray) {
3023  AliHFAODMCParticleContainer* part_cont = dynamic_cast<AliHFAODMCParticleContainer*>(cont_it.second);
3024  if (part_cont) fMCContainer = part_cont;
3025  }
3026 
3027  for (auto &params : fAnalysisEngines) {
3028 
3029  params.fAodEvent = fAodEvent;
3030  params.fFastJetWrapper = fFastJetWrapper;
3031  params.fTrackEfficiency = fTrackEfficiency;
3032  params.fRejectISR = fRejectISR;
3033  params.fRandomGen = rnd;
3034 
3035  for (auto &jetdef: params.fJetDefinitions) {
3036  if (!jetdef.fRhoName.IsNull()) {
3037  jetdef.fRho = dynamic_cast<AliRhoParameter*>(fInputEvent->FindListObject(jetdef.fRhoName));
3038  if (!jetdef.fRho) {
3039  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3040  "%s: Could not find rho object '%s' for engine '%s'",
3041  GetName(), jetdef.fRhoName.Data(), params.GetName());
3042  }
3043  }
3044  }
3045 
3046  if (!params.fRDHFCuts) {
3047  if (params.fMCMode == kMCTruth) {
3048  ::Warning("AliAnalysisTaskDmesonJets::ExecOnce",
3049  "%s: RDHF cuts not provided for engine '%s'.",
3050  GetName(), params.GetName());
3051  }
3052  else {
3053  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3054  "%s: RDHF cuts not provided. Engine '%s' disabled.",
3055  GetName(), params.GetName());
3056  params.fInhibit = kTRUE;
3057  }
3058  }
3059 
3060  params.fMCContainer = fMCContainer;
3061 
3062  for (auto cont_it : fParticleCollArray) {
3063  AliTrackContainer* track_cont = dynamic_cast<AliTrackContainer*>(cont_it.second);
3064  if (track_cont) params.fTrackContainers.push_back(track_cont);
3065  }
3066 
3067  for (auto cont_it : fClusterCollArray) params.fClusterContainers.push_back(cont_it.second);
3068 
3069  if (fAodEvent) params.Init(fGeom, fAodEvent->GetRunNumber());
3070 
3071  if (params.fMCMode != kMCTruth && fAodEvent) {
3072  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
3073 
3074  if (params.fCandidateArray) {
3075  TString className;
3076  if (params.fCandidateType == kD0toKpi || params.fCandidateType == kD0toKpiLikeSign) {
3077  className = "AliAODRecoDecayHF2Prong";
3078  }
3079  else if (params.fCandidateType == kDstartoKpipi) {
3080  className = "AliAODRecoCascadeHF";
3081  }
3082  if (!params.fCandidateArray->GetClass()->InheritsFrom(className)) {
3083  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3084  "%s: Objects of type %s in %s are not inherited from %s! Task will be disabled!",
3085  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName(), className.Data()); // @suppress("Ambiguous problem")
3086  params.fCandidateArray = 0;
3087  params.fInhibit = kTRUE;
3088  }
3089  }
3090  else {
3091  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3092  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
3093  params.fBranchName.Data(), params.GetName());
3094  params.fInhibit = kTRUE;
3095  }
3096  }
3097 
3098  if (params.fMCMode != kNoMC) {
3099  if (!params.fMCContainer) {
3100  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3101  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
3102  params.GetName());
3103  params.fInhibit = kTRUE;
3104  }
3105  }
3106 
3107  if (params.fMCMode != kMCTruth) {
3108  if (params.fTrackContainers.size() == 0 && params.fClusterContainers.size() == 0) {
3109  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3110  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
3111  params.GetName());
3112  params.fInhibit = kTRUE;
3113  }
3114  }
3115  }
3116 }
3117 
3122 {
3123  TString hname;
3124 
3125  // Fix for temporary bug in ESDfilter
3126  // The AODs with null vertex pointer didn't pass the PhysSel
3127  // Now adding an entry in the histogram so as to check that this is actually cutting anything out
3128  if (fAodEvent && (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001)) {
3129  for (auto &eng : fAnalysisEngines) {
3130  if (eng.fInhibit) continue;
3131  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
3132  fHistManager.FillTH1(hname, "ESDfilterBug");
3133  }
3134  return kFALSE;
3135  }
3136 
3137  if (fAodEvent) {
3138  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
3139  if (matchingAODdeltaAODlevel <= 0) {
3140  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
3141  for (auto &eng : fAnalysisEngines) {
3142  if (eng.fInhibit) continue;
3143  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
3144  fHistManager.FillTH1(hname, "MismatchDeltaAOD");
3145  }
3146  return kFALSE;
3147  }
3148  }
3149 
3150  for (auto &eng : fAnalysisEngines) {
3151  eng.fDmesonJets.clear();
3152  if (eng.fInhibit) continue;
3153 
3154  eng.fCent = fCent;
3155 
3156  //Event selection
3157  hname = TString::Format("%s/fHistNEvents", eng.GetName());
3158  if (fAodEvent) {
3159  Bool_t iseventselected = kTRUE;
3160  if (eng.fRDHFCuts) iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
3161  if (!iseventselected) {
3162  fHistManager.FillTH1(hname, "Rejected");
3163  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
3164  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
3165  TString label;
3166  do {
3167  label = GetHFEventRejectionReasonLabel(bitmap);
3168  if (label.IsNull()) break;
3169  fHistManager.FillTH1(hname, label);
3170  } while (true);
3171 
3172  continue;
3173  }
3174  }
3175 
3176  fHistManager.FillTH1(hname, "Accepted");
3177 
3178  AliDebug(2, "Event selected");
3179 
3180  eng.RunAnalysis();
3181  }
3182  return kTRUE;
3183 }
3184 
3189 {
3190  for (auto &param : fAnalysisEngines) {
3191  if (param.fInhibit) continue;
3192 
3193  if (fOutputType == kTreeOutput) {
3194  param.FillTree(fApplyKinematicCuts);
3196  }
3197  else if (fOutputType == kTHnOutput) {
3198  param.FillHnSparse(fApplyKinematicCuts);
3199  }
3200  }
3202  return kTRUE;
3203 }
3204 
3207 {
3208  auto itcont = fMCContainer->all_momentum();
3209  Int_t nHQ = 0;
3210  Double_t highestPartonPt = 0;
3211  Int_t PdgHighParton = 0;
3212  for (auto it = itcont.begin(); it != itcont.end(); it++) {
3213  auto part = *it;
3214  if (part.first.Pt() == 0) continue;
3215 
3216  Int_t PdgCode = part.second->GetPdgCode();
3217 
3218  // Skip all particles that are not either quarks or gluons
3219  if ((PdgCode < -9 || PdgCode > 9) && PdgCode != 21 && PdgCode != -21) continue;
3220 
3221  AliDebugStream(5) << "Parton " << it.current_index() <<
3222  " with pdg=" << PdgCode <<
3223  ", px=" << part.first.Px() <<
3224  ", py=" << part.first.Py() <<
3225  ", pz=" << part.first.Pz() <<
3226  ", n daughters = " << part.second->GetNDaughters() <<
3227  std::endl;
3228 
3229  // Skip partons that do not have any children
3230  // Unclear how this can happen, it would seem that this parton were not fragmented by the generator
3231  if (part.second->GetNDaughters() == 0) continue;
3232 
3233  // Look for highest momentum parton
3234  if (highestPartonPt < part.first.Pt()) {
3235  highestPartonPt = part.first.Pt();
3236  PdgHighParton = PdgCode;
3237  }
3238 
3239  // Skip partons that are not HF
3240  if (PdgCode != 4 && PdgCode != 5 && PdgCode != -4 && PdgCode != -5) continue;
3241 
3242  Bool_t lastInPartonShower = kTRUE;
3243  Bool_t hadronDaughter = kFALSE;
3244  for (Int_t daughterIndex = part.second->GetFirstDaughter(); daughterIndex <= part.second->GetLastDaughter(); daughterIndex++){
3245  if (daughterIndex < 0) {
3246  AliDebugStream(5) << "Could not find daughter index!" << std::endl;
3247  continue;
3248  }
3249  AliAODMCParticle *daughter = fMCContainer->GetMCParticle(daughterIndex);
3250  if (!daughter) {
3251  AliDebugStream(5) << "Could not find particle with index " << daughterIndex << "!" << std::endl;
3252  continue;
3253  }
3254  Int_t daughterPdgCode = daughter->GetPdgCode();
3255  if (daughter->GetMother() != it.current_index()) {
3256  AliDebugStream(5) << "Particle " << daughterIndex << " with pdg=" << daughterPdgCode <<
3257  ", px=" << daughter->Px() <<
3258  ", py=" << daughter->Py() <<
3259  ", pz=" << daughter->Pz() <<
3260  ", is not a daughter of " << it.current_index() <<
3261  "!" << std::endl;
3262  continue;
3263  }
3264 
3265  AliDebugStream(5) << "Found daughter " << daughterIndex <<
3266  " with pdg=" << daughterPdgCode <<
3267  ", px=" << daughter->Px() <<
3268  ", py=" << daughter->Py() <<
3269  ", pz=" << daughter->Pz() <<
3270  std::endl;
3271  if (daughterPdgCode == PdgCode) lastInPartonShower = kFALSE; // this parton is not the last parton in the shower
3272  if (TMath::Abs(daughterPdgCode) >= 111) hadronDaughter = kTRUE;
3273  }
3274  if (hadronDaughter) {
3275  AliDebugStream(5) << "This particle has at least a hadron among its daughters!" << std::endl;
3276  if (!lastInPartonShower) AliDebugStream(2) << "Odly, quark " << it.current_index() << " with PDG " << PdgCode << " (pt = " << part.first.Pt() << ", eta = " << part.first.Eta() << ") is not the last in the parton shower but at least a hadron found among its daughters?!" << std::endl;
3277  }
3278  else {
3279  AliDebugStream(5) << "This particle does not have hadrons among its daughters!" << std::endl;
3280  if (lastInPartonShower) AliDebugStream(2) << "Odly, quark " << it.current_index() << " with PDG " << PdgCode << " (pt = " << part.first.Pt() << ", eta = " << part.first.Eta() << ") is the last in the parton shower but no hadron found among its daughters?!" << std::endl;
3281  continue;
3282  }
3283 
3284  if (PdgCode == 4 || PdgCode == -4) {
3285  fHistManager.FillTH1("fHistCharmPt", part.first.Pt());
3286  fHistManager.FillTH1("fHistCharmEta", part.first.Eta());
3287  fHistManager.FillTH1("fHistCharmPhi", part.first.Phi_0_2pi());
3288  }
3289  else if (PdgCode == 5 || PdgCode == -5) {
3290  fHistManager.FillTH1("fHistBottomPt", part.first.Pt());
3291  fHistManager.FillTH1("fHistBottomEta", part.first.Eta());
3292  fHistManager.FillTH1("fHistBottomPhi", part.first.Phi_0_2pi());
3293  }
3294  nHQ++;
3295  }
3296  fHistManager.FillTH1("fHistNHeavyQuarks", nHQ);
3297  fHistManager.FillTH1("fHistHighestPartonPt",highestPartonPt);
3298  Int_t partonType = 0;
3299  Int_t absPdgHighParton = TMath::Abs(PdgHighParton);
3300  if (absPdgHighParton == 9 || absPdgHighParton == 21) {
3301  partonType = 7;
3302  }
3303  else {
3304  partonType = absPdgHighParton;
3305  }
3306  fHistManager.FillTH1("fHistHighestPartonType",partonType);
3307 }
3308 
3317 {
3318  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
3319 
3320  Double_t mass = part->Mass();
3321 
3322  // To make sure that the PDG mass value is not at the edge of a bin
3323  if (nbins % 2 == 0) {
3324  minMass = mass - range / 2 - range / nbins / 2;
3325  maxMass = mass + range / 2 - range / nbins / 2;
3326  }
3327  else {
3328  minMass = mass - range / 2;
3329  maxMass = mass + range / 2;
3330  }
3331 }
3332 
3339 {
3340  static TString label;
3341  label = "";
3342 
3343  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
3344  label = "NotSelTrigger";
3345  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
3346  return label.Data();
3347  }
3348  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
3349  label = "NoVertex";
3350  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
3351  return label.Data();
3352  }
3353  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
3354  label = "TooFewVtxContrib";
3355  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
3356  return label.Data();
3357  }
3358  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
3359  label = "ZVtxOutFid";
3360  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
3361  return label.Data();
3362  }
3363  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
3364  label = "Pileup";
3365  bitmap &= ~BIT(AliRDHFCuts::kPileup);
3366  return label.Data();
3367  }
3368  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
3369  label = "OutsideCentrality";
3370  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
3371  return label.Data();
3372  }
3373  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
3374  label = "PhysicsSelection";
3375  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
3376  return label.Data();
3377  }
3378  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
3379  label = "BadSPDVertex";
3380  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
3381  return label.Data();
3382  }
3383  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
3384  label = "ZVtxSPDOutFid";
3385  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
3386  return label.Data();
3387  }
3388  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
3389  label = "CentralityFlattening";
3390  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
3391  return label.Data();
3392  }
3393  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
3394  label = "BadTrackV0Correl";
3395  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
3396  return label.Data();
3397  }
3398 
3399  return label.Data();
3400 }
3401 
3408 {
3409  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
3410  PostData(eng.GetDataSlotNumber(), eng.GetTree());
3411  return eng.GetDataSlotNumber();
3412  }
3413  else {
3414  return -1;
3415  }
3416 }
3417 
3427 {
3428  // Get the pointer to the existing analysis manager via the static access method
3429  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
3430  if (!mgr) {
3431  ::Error("AddTaskDmesonJets", "No analysis manager to connect to.");
3432  return NULL;
3433  }
3434 
3435  // Check the analysis type using the event handlers connected to the analysis manager
3436  AliVEventHandler* handler = mgr->GetInputEventHandler();
3437  if (!handler) {
3438  ::Error("AddTaskEmcalJetSpectraQA", "This task requires an input event handler");
3439  return NULL;
3440  }
3441 
3442  EDataType_t dataType = kUnknownDataType;
3443 
3444  if (handler->InheritsFrom("AliESDInputHandler")) {
3445  dataType = kESD;
3446  }
3447  else if (handler->InheritsFrom("AliAODInputHandler")) {
3448  dataType = kAOD;
3449  }
3450 
3451  // Init the task and do settings
3452  if (ntracks == "usedefault") {
3453  if (dataType == kESD) {
3454  ntracks = "Tracks";
3455  }
3456  else if (dataType == kAOD) {
3457  ntracks = "tracks";
3458  }
3459  else {
3460  ntracks = "";
3461  }
3462  }
3463 
3464  if (nclusters == "usedefault") {
3465  if (dataType == kESD) {
3466  nclusters = "CaloClusters";
3467  }
3468  else if (dataType == kAOD) {
3469  nclusters = "caloClusters";
3470  }
3471  else {
3472  nclusters = "";
3473  }
3474  }
3475 
3476  if (nMCpart == "usedefault") {
3477  nMCpart = "mcparticles"; // Always needs AliAODMCParticle objects
3478  }
3479 
3480  TString name("AliAnalysisTaskDmesonJets");
3481  if (strcmp(suffix, "") != 0) {
3482  name += TString::Format("_%s", suffix.Data());
3483  }
3484 
3485  AliAnalysisTaskDmesonJets* jetTask = new AliAnalysisTaskDmesonJets(name, nMaxTrees);
3486 
3487  if (!ntracks.IsNull()) {
3488  AliHFTrackContainer* trackCont = new AliHFTrackContainer(ntracks);
3489  jetTask->AdoptParticleContainer(trackCont);
3490  }
3491 
3492  if (!nMCpart.IsNull()) {
3493  AliMCParticleContainer* partCont = new AliHFAODMCParticleContainer(nMCpart);
3494  partCont->SetEtaLimits(-1.5, 1.5);
3495  partCont->SetPtLimits(0, 1000);
3496  jetTask->AdoptParticleContainer(partCont);
3497  }
3498 
3499  jetTask->AddClusterContainer(nclusters.Data());
3500 
3501  // Final settings, pass to manager and set the containers
3502  mgr->AddTask(jetTask);
3503 
3504  // Create containers for input/output
3505  AliAnalysisDataContainer* cinput1 = mgr->GetCommonInputContainer();
3506  TString contname1(name);
3507  contname1 += "_histos";
3508  AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(contname1.Data(),
3509  TList::Class(), AliAnalysisManager::kOutputContainer,
3510  Form("%s", AliAnalysisManager::GetCommonFileName()));
3511 
3512  mgr->ConnectInput(jetTask, 0, cinput1);
3513  mgr->ConnectOutput(jetTask, 1, coutput1);
3514 
3515  for (Int_t i = 0; i < nMaxTrees; i++) {
3516  TString contname = TString::Format("%s_tree_%d", name.Data(), i);
3517  AliAnalysisDataContainer *coutput = mgr->CreateContainer(contname.Data(),
3518  TTree::Class(),AliAnalysisManager::kOutputContainer,
3519  Form("%s", AliAnalysisManager::GetCommonFileName()));
3520  mgr->ConnectOutput(jetTask, 2+i, coutput);
3521  }
3522  return jetTask;
3523 }
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.
Bool_t fInhibit
!inhibit execution of the task
virtual Int_t Run()
double Double_t
Definition: External.C:58
Bool_t IsJetInAcceptance(const AliJetInfo &jet) const
Double_t fSoftPionPt
! Transverse momentum of the soft pion of the D* candidate
EDataType_t
Switch for the data type.
virtual AliAODMCParticle * GetMCParticle(Int_t i=-1) const
void Getd0MeasMinusExpProng(Int_t ip, Double_t magf, Double_t &d0diff, Double_t &errd0diff) const
const char * title
Definition: MakeQAPdf.C:27
AliRDHFCuts * fRDHFCuts
D meson candidates cuts.
virtual void Set(const AliDmesonJetInfo &source)
TString fRhoName
Name of the object that holds the average background value.
Lightweight class that encapsulates D meson jets.
std::vector< AliTrackContainer * > fTrackContainers
! Track containers
AnalysisEngine * AddAnalysisEngine(ECandidateType_t type, TString cutfname, TString cuttype, EMCMode_t bkgMode, EJetType_t jettype, Double_t jetradius, TString rhoName="")
Double_t fTrackEfficiency
Artificial tracking inefficiency (0...1)
std::vector< AliClusterContainer * > fClusterContainers
! Cluster containers
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
AliHFAODMCParticleContainer * fMCContainer
! MC particle container
ERecoScheme_t fRecoScheme
Jet recombination scheme (pt scheme, E scheme, ...)
const TObjArray & GetDaughterList() const
Double_t fMaxNeutralPt
Transverse momentum of the leading neutral particle (or cluster)
void RunDetectorLevelAnalysis()
Run a detector level analysis.
Double_t mass
Container with name, TClonesArray and cuts for particles.
Double32_t fPartonPt
Transverse momentum of the parton.
friend bool operator<(const AnalysisEngine &lhs, const AnalysisEngine &rhs)
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
virtual void UserCreateOutputObjects()
Creates the output containers.
void FillPartonLevelHistograms()
Fill histograms with parton-level information.
static Int_t CheckMatchingAODdeltaAODevents()
Double32_t fInvMass
Invariant mass of the D0 meson candidate in GeV/c2.
Lightweight class that encapsulates D meson jets.
Double_t fInvMass2Prong
! 2-prong mass of the D* candidate (w/o the soft pion)
friend bool operator==(const AnalysisEngine &lhs, const AnalysisEngine &rhs)
Double_t fMinChargedPt
Minimum pt of the leading charged particle (or track)
static void CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t &minMass, Double_t &maxMass)
AliAODEvent * fAodEvent
! AOD event
Bool_t FindJet(AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, AliHFJetDefinition &jetDef)
UInt_t fEnabledAxis
Use bit defined in EAxis_t to enable axis in the THnSparse.
Double32_t fPtK
Transverse momentum of the kaon.
TRandom * fRandomGen
! Random number generator
Double32_t fCosPointing
Cosine of the pointing angle.
AliRhoParameter * fRho
Object that holds the average background value.
EMCMode_t fMCMode
MC mode: No MC (data and MC detector level), background-only (MC), signal-only (MC), MC truth (particle level)
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
Double_t fJetGhostArea
Area of the ghost particles.
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:122
AnalysisEngine()
This is the default constructor, used for ROOT I/O purposes.
Int_t fDataSlotNumber
! Data slot where the tree output is posted
Double32_t fPtPi
Transverse momentum of the pion.
void Init(const AliEMCALGeometry *const geom, Int_t runNumber)
Initialize the analysis engine.
static FJRecoScheme ConvertToFJRecoScheme(ERecoScheme_t reco)
AliVParticle * fDmesonParticle
! pointer to the particle object
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
static TString GenerateJetName(EJetType_t jetType, EJetAlgo_t jetAlgo, ERecoScheme_t recoScheme, Double_t radius, AliParticleContainer *partCont, AliClusterContainer *clusCont, TString tag)
static std::pair< AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle * > IsPromptCharm(const AliAODMCParticle *part, TClonesArray *mcArray)
virtual void Set(const AliDmesonJetInfo &source)
Double32_t fN
Number of jet constituents.
void SetR(Double_t r)
Definition: AliFJWrapper.h:127
Select tracks based on specific prescriptions of HF analysis.
Double32_t fPt
Transverse momentum of the jet in GeV/c.
ECandidateType_t fCandidateType
Candidate type.
Double_t GetDistance(std::string n, Double_t &deta, Double_t &dphi) const
Int_t fMCLabel
! MC label, i.e. index of the generator level D meson (only for detector level D meson candidates) ...
EBeamType_t fForceBeamType
forced beam type
const Int_t nPtBins
Double32_t fR
Distance between D meson and jet axis.
virtual void Set(const AliDmesonJetInfo &source)
void AddInputVectors(AliEmcalContainer *cont, Int_t offset, TH2 *rejectHist=0, Double_t eff=0.)
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:34
EJetType_t fJetType
Jet type (charged, full, neutral)
Double32_t fPt
Transverse momentum of the D meson in GeV/c.
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
Create a new TH2 within the container.
virtual void Set(const AliDmesonJetInfo &source)
AliTLorentzVector fD
! 4-momentum of the D meson candidate
Double_t GetDistance(const AliJetInfo &jet, Double_t &deta, Double_t &dphi) const
TObject * FindObject(const char *name) const
Find an object inside the container.
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
Double_t CosThetaStarD0bar() const
angle of K
AliAnalysisTaskDmesonJets()
This is the default constructor, used for ROOT I/O purposes.
virtual void Reset()
Reset the current object.
AnalysisEngine & operator=(const AnalysisEngine &source)
Lightweight class that encapsulates D meson jets for PbPb analysis.
Struct that encapsulates analysis parameters.
std::vector< AliHFJetDefinition >::iterator FindJetDefinition(const AliHFJetDefinition &eng)
Base task in the EMCAL framework (lighter version of AliAnalysisTaskEmcal)
Double_t Phi_0_2pi() const
static AliAODMCParticle * FindParticleOrigin(const AliAODMCParticle *part, TClonesArray *mcArray, EFindParticleOriginMode_t mode, const std::set< UInt_t > &pdgSet)
Double_t fMinNeutralPt
Minimum pt of the leading neutral particle (or cluster)
Bool_t fD0D0bar
! kTRUE if selected both as D0 and D0bar
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:121
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
AliAODTrack * GetBachelor() const
Bool_t fD0Extended
Store extended information in the tree (only for D0 mesons)
AliRDHFCuts * LoadDMesonCutsFromFile(TString cutfname, TString cutsname)
std::list< AnalysisEngine >::iterator FindAnalysisEngine(const AnalysisEngine &eng)
std::vector< AliHFJetDefinition > fJetDefinitions
Jet definitions.
Double_t fMaxNeutralPt
Maximum pt of the leading neutral particle (or cluster)
Double_t fMaxMass
Max mass in histogram axis.
EJetAlgo_t fJetAlgo
Jet algorithm (kt, anti-kt,...)
static const char * GetHFEventRejectionReasonLabel(UInt_t &bitmap)
Lightweight class that encapsulates D*.
virtual void Clear(const Option_t *="")
friend bool operator<(const AliHFJetDefinition &lhs, const AliHFJetDefinition &rhs)
Int_t mode
Definition: anaM.C:41
Bool_t ExtractDstarAttributes(const AliAODRecoCascadeHF *DstarCand, AliDmesonJetInfo &DmesonJet, UInt_t i)
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
const AliJetInfo * GetJet(std::string n) const
Short_t fPartonType
! type of the parton in the shower tree (only for particle level D mesons)
TClonesArray * fCandidateArray
! D meson candidate array
const AliMCParticleIterableMomentumContainer all_momentum() const
virtual ~AliAnalysisTaskDmesonJets()
This is the standard destructor.
Byte_t fSelectionType
! for D0: 0=not selected, 1=D0, 2=D0bar
Bool_t ExtractD0Attributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
Select MC particles based on specific prescriptions of HF analysis.
static EMesonDecayChannel_t CheckDecayChannel(const AliAODMCParticle *part, TClonesArray *mcArray)
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:125
AliDmesonJetInfo & operator=(const AliDmesonJetInfo &source)
Double_t fMaxChargedPt
Transverse momentum of the leading charged particle (or track)
Double32_t fCorrPt
Transverse momentum of the jet in GeV/c after subtracting average background.
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
friend bool operator==(const AliHFJetDefinition &lhs, const AliHFJetDefinition &rhs)
static FJJetAlgo ConvertToFJAlgo(EJetAlgo_t algo)
Double_t fMaxChargedPt
Maximum pt of the leading charged particle (or track)
decay
Definition: HFPtSpectrum.C:41
std::map< std::string, AliJetInfo > fJets
! list of jets
Look for the very first particle in the fragmentation tree.
AliAODMCParticle * fAncestor
! pointer to the ancestor particle in the shower tree of the D meson (only for particle level D meson...
Double32_t fDeltaInvMass
< Difference between the Kpipi and the Kpi invariant masses in GeV/c2
std::map< AliAODMCParticle *, Short_t > fPartons
! set of the partons in the shower that produced each D meson
Definition: External.C:220
Double_t minMass
Bool_t fRejectISR
! Reject initial state radiation
void RunParticleLevelAnalysis()
Run a particle level analysis.
std::map< std::string, AliParticleContainer * > fParticleCollArray
particle/track collection array
std::map< int, AliDmesonJetInfo > fDmesonJets
! Array containing the D meson jets
Double_t fTrackEfficiency
! Artificial tracking inefficiency (0...1) -> set automatically at ExecOnce by AliAnalysisTaskDmesonJ...
Bool_t fApplyKinematicCuts
Apply jet kinematic cuts.
Double32_t fSelectionType
Selection type: D0, D0bar, both.
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:290
const char * GetName() const
Generate a name for this jet definition.
AliEventNotFound(const std::string &class_name, const std::string &method_name)
UInt_t fAcceptedDecay
Bit mask with D meson decays that are accepted (only used for particle-level analysis) ...
virtual void Set(const AliDmesonJetInfo &source, std::string n)
unsigned short UShort_t
Definition: External.C:28
Bool_t fReconstructed
! Whether this D meson was reconstructed (only for particle level D mesons)
Bool_t ExtractRecoDecayAttributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
Double_t InvMassD0() const
Bool_t fRejectISR
Reject initial state radiation.
EOutputType_t fOutputType
Output type: none, TTree or THnSparse.
virtual void AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
Double_t InvMassDstarKpipi() const
const AliVEvent * fEvent
! pointer to the ESD/AOD event
TArrayI fPDGdaughters
List of the PDG code of the daughters.
const Int_t nbins
Double_t maxMass
TString fBranchName
AOD branch where the D meson candidate are found.
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
Int_t fNOutputTrees
Maximum number of output trees.
AliHFJetDefinition()
This is the default constructor, used for ROOT I/O purposes.
AliHFJetDefinition & operator=(const AliHFJetDefinition &source)
Int_t GetNAcceptedParticles() const
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:123
std::map< std::string, AliClusterContainer * > fClusterCollArray
cluster collection array
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
Create a new THnSparse within the container.
Double_t fCorrPt
Transverse momentum of the jet after subtracting the average background.
AliHFAODMCParticleContainer * fMCContainer
! MC particle container
AliJetInfoSummary ** fCurrentJetInfo
! Current jet info
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:310
AliFJWrapper * fFastJetWrapper
! Fastjet wrapper
AliHFJetDefinition * AddJetDefinition(EJetType_t type, Double_t r, EJetAlgo_t algo, ERecoScheme_t reco)
Container for MC-true particles within the EMCAL framework.
Lightweight class that encapsulates D0.
Bool_t IsAnyJetInAcceptance(const AliDmesonJetInfo &dMesonJet) const
Double_t fCent
!event centrality
void SetCharge(EChargeCut_t c)
Container for jet within the EMCAL jet framework.
Definition: External.C:196
Int_t PostDataFromAnalysisEngine(const AnalysisEngine &eng)
AliFJWrapper * fFastJetWrapper
! Fastjet wrapper
THistManager fHistManager
Histogram manager.
void Reset()
Reset all fields to their default values.
static AliAnalysisTaskDmesonJets * AddTaskDmesonJets(TString ntracks="usedefault", TString nclusters="usedefault", TString nMCpart="", Int_t nMaxTrees=2, TString suffix="")
virtual void Set(const AliDmesonJetInfo &source, std::string n)
Double_t fMinMass
Min mass in histogram axis.