AliPhysics  97344c9 (97344c9)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskDmesonJets.cxx
Go to the documentation of this file.
1 /*************************************************************************
2 * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15 
16 // Root
17 #include <TClonesArray.h>
18 #include <TDatabasePDG.h>
19 #include <TParticlePDG.h>
20 #include <TVector3.h>
21 #include <THnSparse.h>
22 #include <TParticle.h>
23 #include <TMath.h>
24 #include <THashList.h>
25 #include <TFile.h>
26 #include <TRandom3.h>
27 
28 // Aliroot general
29 #include "AliLog.h"
30 #include "AliEMCALGeometry.h"
31 #include "AliAnalysisManager.h"
32 #include "AliVEventHandler.h"
33 
34 // Aliroot HF
35 #include "AliAODRecoCascadeHF.h"
37 #include "AliRDHFCutsD0toKpi.h"
40 #include "AliHFTrackContainer.h"
41 
42 // Aliroot EMCal jet framework
43 #include "AliEmcalJetTask.h"
44 #include "AliEmcalJet.h"
45 #include "AliJetContainer.h"
46 #include "AliParticleContainer.h"
47 #include "AliEmcalParticle.h"
48 #include "AliFJWrapper.h"
49 
51 
52 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfo
53 
57 
65 {
66  dphi = TVector2::Phi_mpi_pi(fMomentum.Phi() - jet.Phi());;
67  deta = fMomentum.Eta() - jet.Eta();
68  return TMath::Sqrt(dphi*dphi + deta*deta);
69 }
70 
76 {
77  Double_t deta = 0;
78  Double_t dphi = 0;
79  return GetDistance(jet, deta, dphi);
80 }
81 
82 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonJetInfo
83 
87 
90  fDmesonParticle(0),
91  fD(),
92  fSoftPionPt(0),
93  fInvMass2Prong(0),
94  fJets(),
95  fMCLabel(-1),
96  fReconstructed(kFALSE),
97  fFirstParton(0),
98  fFirstPartonType(0),
99  fLastParton(0),
100  fLastPartonType(0),
101  fSelectionType(0)
102 {
103 }
104 
109  fDmesonParticle(source.fDmesonParticle),
110  fD(source.fD),
111  fSoftPionPt(source.fSoftPionPt),
112  fInvMass2Prong(source.fInvMass2Prong),
113  fJets(source.fJets),
114  fMCLabel(source.fMCLabel),
115  fReconstructed(source.fReconstructed),
116  fFirstParton(source.fFirstParton),
117  fFirstPartonType(source.fFirstPartonType),
118  fLastParton(source.fLastParton),
119  fLastPartonType(source.fLastPartonType),
120  fSelectionType(source.fSelectionType)
121 {
122 }
123 
128 {
129  new (this) AliDmesonJetInfo(source);
130  return *this;
131 }
132 
135 {
136  fD.SetPtEtaPhiE(0,0,0,0);
137  fSoftPionPt = 0;
138  fInvMass2Prong = 0;
139  fDmesonParticle = 0;
140  fMCLabel = -1;
141  fReconstructed = kFALSE;
142  fFirstParton = 0;
143  fFirstPartonType = 0;
144  fLastParton = 0;
145  fLastPartonType = 0;
146  for (auto &jet : fJets) {
147  jet.second.fMomentum.SetPtEtaPhiE(0,0,0,0);
148  jet.second.fNConstituents = 0;
149  jet.second.fNEF = 0;
150  jet.second.fMaxChargedPt = 0;
151  jet.second.fMaxNeutralPt = 0;
152  }
153 }
154 
157 {
158  Printf("Printing D Meson Jet object.");
159  Printf("D Meson: pT = %.3f, eta = %.3f, phi = %.3f, inv. mass = %.3f", fD.Pt(), fD.Eta(), fD.Phi_0_2pi(), fD.M());
160  Printf("Soft pion pT: %.3f. 2-Prong Invariant mass = %.3f", fSoftPionPt, fInvMass2Prong);
161  for (auto &jet : fJets) {
162  Printf("Jet %s: pT = %.3f, eta = %.3f, phi = %.3f", jet.first.c_str(), jet.second.Pt(), jet.second.Eta(), jet.second.Phi_0_2pi());
163  Printf("Jet N Consituents = %d", jet.second.fNConstituents);
164  }
165 }
166 
171 {
172  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
173  if (it == fJets.end()) return 0;
174 
175  Double_t z = 0;
176 
177  if ((*it).second.Pt() > 0) {
178  TVector3 dvect = fD.Vect();
179  TVector3 jvect = (*it).second.fMomentum.Vect();
180 
181  Double_t jetMom = jvect * jvect;
182 
183  if (jetMom < 1e-6) {
184  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetZ", "Zero jet momentum!");
185  z = 0.999;
186  }
187  else {
188  z = (dvect * jvect) / jetMom;
189  }
190 
191  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
192  }
193 
194  return z;
195 }
196 
204 {
205  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
206  if (it == fJets.end()) return 0;
207 
208  return GetDistance((*it).second, deta, dphi);
209 }
210 
216 {
217  Double_t deta = 0;
218  Double_t dphi = 0;
219  return GetDistance(n, deta, dphi);
220 }
221 
229 {
230  dphi = TVector2::Phi_mpi_pi(fD.Phi() - jet.Phi());;
231  deta = fD.Eta() - jet.Eta();
232  return TMath::Sqrt(dphi*dphi + deta*deta);
233 }
234 
240 {
241  Double_t deta = 0;
242  Double_t dphi = 0;
243  return GetDistance(jet, deta, dphi);
244 }
245 
251 {
252  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
253  if (it == fJets.end()) {
254  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
255  n.c_str());
256  return 0;
257  }
258  return &((*it).second);
259 }
260 
266 {
267  std::map<std::string, AliJetInfo>::iterator it = fJets.find(n);
268  if (it == fJets.end()) {
269  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
270  n.c_str());
271  return 0;
272  }
273  return &((*it).second);
274 }
275 
276 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoSummary
277 
281 
287  fPt(0),
288  fEta(0),
289  fPhi(0),
290  fR(0),
291  fZ(0),
292  fN(0)
293 {
294  Set(source, n);
295 }
296 
299 {
300  fPt = 0;
301  fEta = 0;
302  fPhi = 0;
303  fR = 0;
304  fZ = 0;
305  fN = 0;
306 }
307 
313 {
314  std::map<std::string, AliJetInfo>::const_iterator it = source.fJets.find(n);
315  if (it == source.fJets.end()) return;
316 
317  fPt = (*it).second.Pt();
318  fEta = (*it).second.Eta();
319  fPhi = (*it).second.Phi_0_2pi();
320  fR = source.GetDistance(n);
321  fZ = source.GetZ(n);
322  fN = (*it).second.GetNConstituents();
323 }
324 
329 {
330  fPt = source.Pt();
331  fEta = source.Eta();
332  fPhi = source.Phi_0_2pi();
333  fN = source.GetNConstituents();
334  fR = 0;
335  fZ = 0;
336 }
337 
338 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonInfoSummary
339 
343 
348  fPt(0),
349  fEta(0),
350  fPhi(0)
351 {
352  Set(source);
353 }
354 
359 {
360  fPt = source.fD.Pt();
361  fEta = source.fD.Eta();
362  fPhi = source.fD.Phi_0_2pi();
363 }
364 
367 {
368  fPt = 0;
369  fEta = 0;
370  fPhi = 0;
371 }
372 
373 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary
374 
378 
383  AliDmesonInfoSummary(source),
384  fFirstPartonType(0),
385  fFirstPartonPt(0),
386  fLastPartonType(0),
387  fLastPartonPt(0)
388 {
389  Set(source);
390 }
391 
396 {
398 
399  fFirstPartonType = source.fFirstPartonType;
400  fLastPartonType = source.fLastPartonType;
401 
402  if (source.fFirstParton) {
403  fFirstPartonPt = source.fFirstParton->Pt();
404  }
405  else {
406  fFirstPartonPt = 0.;
407  }
408 
409  if (source.fLastParton) {
410  fLastPartonPt = source.fLastParton->Pt();
411  }
412  else {
413  fLastPartonPt = 0.;
414  }
415 }
416 
419 {
421  fFirstPartonType = 0,
422  fFirstPartonPt = 0.;
423  fLastPartonType = 0,
424  fLastPartonPt = 0.;
425 }
426 
427 // Definitions of class AliAnalysisTaskDmesonJets::AliD0InfoSummary
428 
432 
437  AliDmesonInfoSummary(source),
438  fInvMass(source.fD.M()),
439  fSelectionType(0)
440 {
441 }
442 
447 {
448  fInvMass = source.fD.M();
449  fSelectionType = source.fSelectionType;
451 }
452 
455 {
457  fSelectionType = 0;
458  fInvMass = 0;
459 }
460 
461 // Definitions of class AliAnalysisTaskDmesonJets::AliDStarInfoSummary
462 
466 
471  AliDmesonInfoSummary(source),
472  f2ProngInvMass(source.fInvMass2Prong),
473  fDeltaInvMass(source.fD.M() - source.fInvMass2Prong)
474 {
475 }
476 
481 {
482  f2ProngInvMass = source.fInvMass2Prong;
483  fDeltaInvMass = source.fD.M() - source.fInvMass2Prong;
485 }
486 
489 {
491 
492  f2ProngInvMass = 0;
493  fDeltaInvMass = 0;
494 }
495 
496 // Definitions of class AliAnalysisTaskDmesonJets::AliJetDefinition
497 
501 
504  TObject(),
505  fJetType(AliJetContainer::kChargedJet),
506  fRadius(0),
507  fJetAlgo(AliJetContainer::antikt_algorithm),
508  fRecoScheme(AliJetContainer::pt_scheme),
509  fMinJetPt(0.),
510  fMaxJetPt(500.),
511  fMinJetPhi(0.),
512  fMaxJetPhi(0.),
513  fMinJetEta(-1.),
514  fMaxJetEta(1.),
515  fMinChargedPt(0.),
516  fMaxChargedPt(0.),
517  fMinNeutralPt(0.),
518  fMaxNeutralPt(0.)
519 {
520 }
521 
529  TObject(),
530  fJetType(type),
531  fRadius(r),
532  fJetAlgo(algo),
533  fRecoScheme(reco),
534  fMinJetPt(0.),
535  fMaxJetPt(500.),
536  fMinJetPhi(0.),
537  fMaxJetPhi(0.),
538  fMinJetEta(-1.),
539  fMaxJetEta(1.),
540  fMinChargedPt(0.),
541  fMaxChargedPt(0.),
542  fMinNeutralPt(0.),
543  fMaxNeutralPt(0.)
544 {
545 }
546 
551  TObject(),
552  fJetType(source.fJetType),
553  fRadius(source.fRadius),
554  fJetAlgo(source.fJetAlgo),
555  fRecoScheme(source.fRecoScheme),
556  fMinJetPt(source.fMinJetPt),
557  fMaxJetPt(source.fMaxJetPt),
558  fMinJetPhi(source.fMinJetPhi),
559  fMaxJetPhi(source.fMaxJetPhi),
560  fMinJetEta(source.fMinJetEta),
561  fMaxJetEta(source.fMaxJetEta),
562  fMinChargedPt(source.fMinChargedPt),
563  fMaxChargedPt(source.fMaxChargedPt),
564  fMinNeutralPt(source.fMinNeutralPt),
565  fMaxNeutralPt(source.fMaxNeutralPt)
566 {
567 }
568 
573 {
574  new (this) AliHFJetDefinition(source);
575  return *this;
576 }
577 
580 {
581  static TString name;
582 
583  name = AliJetContainer::GenerateJetName(fJetType, fJetAlgo, fRecoScheme, fRadius, 0, 0, "Jet");
584 
585  return name.Data();
586 }
587 
593 {
594  if (fMinJetEta < fMaxJetEta && (jet.Eta() < fMinJetEta || jet.Eta() > fMaxJetEta)) return kFALSE;
595  if (fMinJetPhi < fMaxJetPhi && (jet.Phi() < fMinJetPhi || jet.Phi() > fMaxJetPhi)) return kFALSE;
596  if (jet.Pt() > fMaxJetPt || jet.Pt() < fMinJetPt) return kFALSE;
597  if (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt) return kFALSE;
598  if (jet.fMaxNeutralPt < fMinNeutralPt || jet.fMaxNeutralPt > fMaxNeutralPt) return kFALSE;
599 
600  return kTRUE;
601 }
602 
608 {
609  const AliJetInfo* jet = dMesonJet.GetJet(n);
610  if (!jet) return kFALSE;
611  return IsJetInAcceptance((*jet));
612 }
613 
620 {
621  if (lhs.fJetType > rhs.fJetType) return false;
622  else if (lhs.fJetType < rhs.fJetType) return true;
623  else {
624  if (lhs.fRadius > rhs.fRadius) return false;
625  else if (lhs.fRadius < rhs.fRadius) return true;
626  else {
627  if (lhs.fJetAlgo > rhs.fJetAlgo) return false;
628  else if (lhs.fJetAlgo < rhs.fJetAlgo) return true;
629  else {
630  if (lhs.fRecoScheme < rhs.fRecoScheme) return true;
631  else return false;
632  }
633  }
634  }
635 }
636 
643 {
644  if (lhs.fJetType != rhs.fJetType) return false;
645  if (lhs.fRadius != rhs.fRadius) return false;
646  if (lhs.fJetAlgo != rhs.fJetAlgo) return false;
647  if (lhs.fRecoScheme != rhs.fRecoScheme) return false;
648  return true;
649 }
650 
651 // Definitions of class AliAnalysisTaskDmesonJets::AnalysisEngine
652 
656 
659  TObject(),
660  fFirstPartons(),
661  fLastPartons(),
662  fCandidateType(kD0toKpi),
663  fCandidateName(),
664  fCandidatePDG(0),
665  fNDaughters(0),
666  fPDGdaughters(),
667  fBranchName(),
668  fMCMode(kNoMC),
669  fNMassBins(0),
670  fMinMass(0),
671  fMaxMass(0),
672  fRDHFCuts(0),
673  fRejectedOrigin(0),
674  fAcceptedDecay(0),
675  fInhibit(kFALSE),
676  fJetDefinitions(),
677  fPtBinWidth(0.5),
678  fMaxPt(100),
679  fRandomGen(0),
680  fTrackEfficiency(0),
681  fDataSlotNumber(-1),
682  fTree(0),
683  fCurrentDmesonJetInfo(0),
684  fCurrentJetInfo(0),
685  fCandidateArray(0),
686  fMCContainer(0),
687  fTrackContainer(0),
688  fClusterContainer(0),
689  fAodEvent(0),
690  fFastJetWrapper(0),
691  fHistManager(0)
692 {
693 }
694 
703  TObject(),
704  fFirstPartons(),
705  fLastPartons(),
706  fCandidateType(type),
707  fCandidateName(),
708  fCandidatePDG(0),
709  fNDaughters(0),
710  fPDGdaughters(),
711  fBranchName(),
712  fMCMode(MCmode),
713  fNMassBins(nMassBins),
714  fMinMass(0),
715  fMaxMass(0),
716  fRDHFCuts(cuts),
717  fRejectedOrigin(0),
718  fAcceptedDecay(kAnyDecay),
719  fInhibit(kFALSE),
720  fJetDefinitions(),
721  fPtBinWidth(0.5),
722  fMaxPt(100),
723  fRandomGen(0),
724  fTrackEfficiency(0),
725  fDataSlotNumber(-1),
726  fTree(0),
727  fCurrentDmesonJetInfo(0),
728  fCurrentJetInfo(0),
729  fCandidateArray(0),
730  fMCContainer(0),
731  fTrackContainer(0),
732  fClusterContainer(0),
733  fAodEvent(0),
734  fFastJetWrapper(0),
735  fHistManager(0)
736 {
737  SetCandidateProperties(range);
738 }
739 
744  TObject(source),
745  fFirstPartons(source.fFirstPartons),
746  fLastPartons(source.fLastPartons),
747  fCandidateType(source.fCandidateType),
748  fCandidateName(source.fCandidateName),
749  fCandidatePDG(source.fCandidatePDG),
750  fNDaughters(source.fNDaughters),
751  fPDGdaughters(source.fPDGdaughters),
752  fBranchName(source.fBranchName),
753  fMCMode(source.fMCMode),
754  fNMassBins(source.fNMassBins),
755  fMinMass(source.fMinMass),
756  fMaxMass(source.fMaxMass),
757  fRDHFCuts(),
758  fRejectedOrigin(source.fRejectedOrigin),
759  fAcceptedDecay(source.fAcceptedDecay),
760  fInhibit(source.fInhibit),
761  fJetDefinitions(source.fJetDefinitions),
762  fPtBinWidth(source.fPtBinWidth),
763  fMaxPt(source.fMaxPt),
764  fRandomGen(source.fRandomGen),
766  fDataSlotNumber(-1),
767  fTree(0),
768  fCurrentDmesonJetInfo(0),
769  fCurrentJetInfo(0),
770  fCandidateArray(source.fCandidateArray),
771  fMCContainer(source.fMCContainer),
772  fTrackContainer(source.fTrackContainer),
773  fClusterContainer(source.fClusterContainer),
774  fAodEvent(source.fAodEvent),
776  fHistManager(source.fHistManager)
777 {
778  SetRDHFCuts(source.fRDHFCuts);
779 }
780 
781 // Destructor
783 {
784  delete fRDHFCuts;
785 }
786 
791 {
792  new (this) AnalysisEngine(source);
793  return *this;
794 }
795 
800 {
801  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
802  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
803  }
804 
805  return kFALSE;
806 }
807 
809 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const /*geom*/, Int_t /*runNumber*/)
810 {
811 }
812 
817 {
818  switch (fCandidateType) {
819  case kD0toKpi:
820  fCandidatePDG = 421;
821  fCandidateName = "D0";
822  fNDaughters = 2;
823  fPDGdaughters.Set(fNDaughters);
824  fPDGdaughters.Reset();
825  fPDGdaughters[0] = 211; // pi
826  fPDGdaughters[1] = 321; // K
827  fBranchName = "D0toKpi";
828  fAcceptedDecay = kDecayD0toKpi;
829  if (!fRDHFCuts) {
830  fRDHFCuts = new AliRDHFCutsD0toKpi();
831  fRDHFCuts->SetStandardCutsPP2010();
832  fRDHFCuts->GetPidHF()->SetOldPid(kFALSE);
833  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
834  fRDHFCuts->SetTriggerClass("","");
835  }
836  break;
837  case kD0toKpiLikeSign:
838  fCandidatePDG = 421;
839  fCandidateName = "2ProngLikeSign";
840  fNDaughters = 2;
841  fPDGdaughters.Set(fNDaughters);
842  fPDGdaughters.Reset();
843  fPDGdaughters[0] = 211; // pi
844  fPDGdaughters[1] = 321; // K
845  fBranchName = "LikeSign2Prong";
846  if (!fRDHFCuts) {
847  fRDHFCuts = new AliRDHFCutsD0toKpi();
848  fRDHFCuts->SetStandardCutsPP2010();
849  fRDHFCuts->GetPidHF()->SetOldPid(kFALSE);
850  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
851  fRDHFCuts->SetTriggerClass("","");
852  }
853  break;
854  case kDstartoKpipi:
855  fCandidatePDG = 413;
856  fCandidateName = "DStar";
857  fNDaughters = 3;
858  fPDGdaughters.Set(fNDaughters);
859  fPDGdaughters.Reset();
860  fPDGdaughters[0] = 211; // pi soft
861  fPDGdaughters[1] = 211; // pi fromD0
862  fPDGdaughters[2] = 321; // K from D0
863  fBranchName = "Dstar";
864  fAcceptedDecay = kDecayDStartoKpipi;
865  if (!fRDHFCuts) {
866  fRDHFCuts = new AliRDHFCutsDStartoKpipi();
867  fRDHFCuts->SetStandardCutsPP2010();
868  fRDHFCuts->GetPidHF()->SetOldPid(kFALSE);
869  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
870  fRDHFCuts->SetTriggerClass("","");
871  }
872  break;
873  default:
874  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
875  }
876 
877  CalculateMassLimits(range, fCandidatePDG, fNMassBins, fMinMass, fMaxMass);
878 }
879 
884 {
885  if (fRDHFCuts) delete fRDHFCuts;
886  fRDHFCuts = cuts;
887 }
888 
893 {
894  if (!cuts) return;
895  if (fRDHFCuts) delete fRDHFCuts;
896  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
897 }
898 
903 {
904  static TString name;
905 
906  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
907 
908  return name.Data();
909 }
910 
915 {
916  static TString name;
917 
918  name = fCandidateName;
919  switch (fMCMode) {
920  case kBackgroundOnly:
921  name += "_kBackgroundOnly";
922  break;
923  case kSignalOnly:
924  name += "_kSignalOnly";
925  break;
926  case kMCTruth:
927  name += "_MCTruth";
928  break;
929  case kWrongPID:
930  name += "_WrongPID";
931  break;
932  default:
933  break;
934  }
935 
936  return name.Data();
937 }
938 
946 {
947  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
948 
949  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
950  fJetDefinitions.push_back(def);
951  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
952  "Total number of jet definitions is now %lu.",
953  def.GetName(), GetName(), fJetDefinitions.size());
954  // For detector level set maximum track pt to 100 GeV/c
955  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
956  }
957  else {
958  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
959  }
960 
961  return &(*it);
962 }
963 
975 {
976  AliHFJetDefinition def(type, r, algo, reco);
977 
978  return AddJetDefinition(def);
979 }
980 
986 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
987 {
988  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
989  while (it != fJetDefinitions.end() && (*it) != def) it++;
990  return it;
991 }
992 
999 {
1000  if (lhs.fCandidateType > rhs.fCandidateType) return false;
1001  else if (lhs.fCandidateType < rhs.fCandidateType) return true;
1002  else {
1003  if (lhs.fMCMode < rhs.fMCMode) return true;
1004  else return false;
1005  }
1006 }
1007 
1014 {
1015  if (lhs.fCandidateType != rhs.fCandidateType) return false;
1016  if (lhs.fMCMode != rhs.fMCMode) return false;
1017  return true;
1018 }
1019 
1028 {
1029  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) { // D0 candidate
1030  return ExtractD0Attributes(Dcand, DmesonJet, i);
1031  }
1032  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
1033  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
1034  }
1035  else {
1036  return kFALSE;
1037  }
1038 }
1039 
1048 {
1049  AliDebug(10,"Checking if D0 meson is selected");
1050  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
1051  if (isSelected == 0) return kFALSE;
1052 
1053  Int_t MCtruthPdgCode = 0;
1054 
1055  Double_t invMassD = 0;
1056 
1057  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
1058  // Checks also the origin, and if it matches the rejected origin mask, return false
1059  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly || fMCMode == kWrongPID) {
1060  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
1061  DmesonJet.fMCLabel = mcLab;
1062 
1063  // Retrieve the generated particle (if exists) and its PDG code
1064  if (mcLab >= 0) {
1065  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1066 
1067  if (aodMcPart) {
1068  // Check origin and return false if it matches the rejected origin mask
1069  if (fRejectedOrigin) {
1070  auto origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
1071  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1072  }
1073  MCtruthPdgCode = aodMcPart->PdgCode();
1074  }
1075  }
1076  }
1077 
1078  if (isSelected == 1) { // selected as a D0
1079  if (i != 0) return kFALSE; // only one mass hypothesis thanks to PID
1080 
1081  if (fMCMode == kNoMC ||
1082  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1083  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly) ||
1084  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kWrongPID)) {
1085  // 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)
1086  AliDebug(10,"Selected as D0");
1087  invMassD = Dcand->InvMassD0();
1088  }
1089  else { // conditions above not passed, so return FALSE
1090  return kFALSE;
1091  }
1092  }
1093  else if (isSelected == 2) { // selected as a D0bar
1094  if (i != 1) return kFALSE; // only one mass hypothesis thanks to PID
1095 
1096  if (fMCMode == kNoMC ||
1097  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1098  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly) ||
1099  (MCtruthPdgCode == fCandidatePDG && fMCMode == kWrongPID)) {
1100  // 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)
1101  AliDebug(10,"Selected as D0bar");
1102  invMassD = Dcand->InvMassD0bar();
1103  }
1104  else { // conditions above not passed, so return FALSE
1105  return kFALSE;
1106  }
1107  }
1108  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
1109  AliDebug(10,"Selected as either D0 or D0bar");
1110 
1111  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
1112  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1113  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kWrongPID))) {
1114  if (i != 0) return kFALSE;
1115  AliDebug(10, "MC truth is D0");
1116  invMassD = Dcand->InvMassD0();
1117  }
1118  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1119  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kWrongPID))) {
1120  if (i != 1) return kFALSE;
1121  AliDebug(10, "MC truth is D0bar");
1122  invMassD = Dcand->InvMassD0bar();
1123  }
1124  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
1125 
1126  // Only accept it if background-only OR background-and-signal was requested
1127  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
1128  // Select D0 or D0bar depending on the i-parameter
1129  if (i == 0) {
1130  AliDebug(10, "Returning invariant mass with D0 hypothesis");
1131  invMassD = Dcand->InvMassD0();
1132  }
1133  else if (i == 1) {
1134  AliDebug(10, "Returning invariant mass with D0bar hypothesis");
1135  invMassD = Dcand->InvMassD0bar();
1136  }
1137  else { // i > 1
1138  return kFALSE;
1139  }
1140  }
1141  else { // signal-only was requested but this is not a true D0
1142  return kFALSE;
1143  }
1144  }
1145  }
1146  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
1147  return kTRUE;
1148 }
1149 
1158 {
1159  AliDebug(10,"Checking if D* meson is selected");
1160  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoCascadeHF*>(DstarCand), AliRDHFCuts::kAll, fAodEvent);
1161  if (isSelected == 0) return kFALSE;
1162 
1163  if ((i == 1 && DstarCand->Charge()>0) || (i == 0 && DstarCand->Charge()<0) || i > 1) return kFALSE; // only one mass hypothesis for the D*
1164 
1165  Int_t MCtruthPdgCode = 0;
1166 
1167  Double_t invMassD = 0;
1168 
1169  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1170  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1171  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1172 
1173  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1174  AliDebug(10, Form("MC label is %d", mcLab));
1175  DmesonJet.fMCLabel = mcLab;
1176  if (mcLab >= 0) {
1177  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1178 
1179  if (aodMcPart) {
1180  if (fRejectedOrigin) {
1181  auto origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
1182  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1183  }
1184 
1185  MCtruthPdgCode = aodMcPart->PdgCode();
1186  AliDebug(10, Form("MC truth pdg code is %d",MCtruthPdgCode));
1187  }
1188  }
1189  }
1190 
1191  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1192  if (fMCMode == kNoMC ||
1193  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1194  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1195  // 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)
1196  invMassD = DstarCand->InvMassDstarKpipi();
1197  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1198  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1199  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1200  return kTRUE;
1201  }
1202  else { // conditions above not passed, so return FALSE
1203  return kFALSE;
1204  }
1205 }
1206 
1214 {
1215  if (!part) return kUnknownDecay;
1216  if (!mcArray) return kUnknownDecay;
1217 
1219 
1220  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1221 
1222  if (part->GetNDaughters() == 2) {
1223 
1224  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1225  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1226 
1227  if (!d1 || !d2) {
1228  return decay;
1229  }
1230 
1231  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1232  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1233 
1234  if (absPdgPart == 421) { // D0 -> K pi
1235  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1236  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1237  decay = kDecayD0toKpi;
1238  }
1239  }
1240 
1241  if (absPdgPart == 413) { // D* -> D0 pi
1242  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1243  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1244  if (D0decay == kDecayD0toKpi) {
1245  decay = kDecayDStartoKpipi;
1246  }
1247  }
1248  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1249  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1250  if (D0decay == kDecayD0toKpi) {
1251  decay = kDecayDStartoKpipi;
1252  }
1253  }
1254  }
1255  }
1256 
1257  return decay;
1258 }
1259 
1266 std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> AliAnalysisTaskDmesonJets::AnalysisEngine::CheckOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, Bool_t firstParton)
1267 {
1268  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1269 
1270  std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> result(kUnknownQuark, 0);
1271 
1272  if (!part) return result;
1273  if (!mcArray) return result;
1274 
1275  Int_t mother = part->GetMother();
1276  while (mother >= 0) {
1277  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1278  if (mcGranma) {
1279  Int_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
1280 
1281  if (abspdgGranma == 1) result = {kFromDown, mcGranma};
1282  if (abspdgGranma == 2) result = {kFromUp, mcGranma};
1283  if (abspdgGranma == 3) result = {kFromStrange, mcGranma};
1284  if (abspdgGranma == 4) result = {kFromCharm, mcGranma};
1285  if (abspdgGranma == 5) result = {kFromBottom, mcGranma};
1286  if (abspdgGranma == 6) result = {kFromTop, mcGranma};
1287  if (abspdgGranma == 9 || abspdgGranma == 21) result = {kFromGluon, mcGranma};
1288 
1289  // If looking for the very first parton in the hard scattering, it will continue the loop until it cannot find a mother particle
1290  if (result.first != kUnknownQuark && !firstParton) return result;
1291 
1292  mother = mcGranma->GetMother();
1293  }
1294  else {
1295  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1296  break;
1297  }
1298  }
1299 
1300  return result;
1301 }
1302 
1305 {
1306  for (auto& jetDef : fJetDefinitions) {
1307  jetDef.fJets.clear();
1308  }
1309 
1310  if (fMCMode == kMCTruth) {
1311  RunParticleLevelAnalysis();
1312  }
1313  else {
1314  RunDetectorLevelAnalysis();
1315  }
1316 }
1317 
1323 {
1325  fFastJetWrapper->SetR(jetDef.fRadius);
1328 
1329  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1330  fTrackContainer->SetDMesonCandidate(0);
1331  AddInputVectors(fTrackContainer, 100);
1332  }
1333 
1334  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1335  AddInputVectors(fClusterContainer, -100);
1336  }
1337 
1338  // run jet finder
1339  fFastJetWrapper->Run();
1340 
1341  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1342 
1343  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1344  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1345 
1346  Double_t maxChPt = 0;
1347  Double_t maxNePt = 0;
1348  Double_t totalNeutralPt = 0;
1349 
1350  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1351  if (constituents[ic].user_index() >= 100) {
1352  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1353  }
1354  else if (constituents[ic].user_index() <= -100) {
1355  totalNeutralPt += constituents[ic].pt();
1356  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1357  }
1358  }
1359 
1360  jetDef.fJets.push_back(
1361  AliJetInfo(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E(),
1362  constituents.size(), totalNeutralPt / jets_incl[ijet].pt(), maxChPt, maxNePt)
1363  );
1364  }
1365 }
1366 
1375 {
1376  if (jetDef.fJets.size() == 0) FindJets(jetDef);
1377 
1378  Double_t d_closest = 999;
1379  AliJetInfo* jet_closest = 0;
1380  const AliJetInfo* truth_jet = 0;
1381  try {
1382  truth_jet = &(dmeson.fJets.at(jetDef.GetName()));
1383  }
1384  catch(...) {
1385  return jet_distance_pair(0, 999);
1386  }
1387 
1388  for (auto& jet : jetDef.fJets) {
1389  Double_t d = truth_jet->GetDistance(jet);
1390 
1391  if (d > dMax) continue;
1392  if (d < d_closest) {
1393  d_closest = d;
1394  jet_closest = &jet;
1395  }
1396  }
1397 
1398  if (jet_closest && applyKinCuts) {
1399  if (!jetDef.IsJetInAcceptance(*jet_closest)) jet_closest = 0;
1400  }
1401 
1402  if (jet_closest) {
1403  AliDebug(2, Form("Found closest jet (pt=%.3f, eta=%.3f, phi=%.3f) to d meson (pt=%.3f, eta=%.3f, phi=%.3f) with d = %.3f",
1404  jet_closest->Pt(), jet_closest->Eta(), jet_closest->Phi_0_2pi(),
1405  dmeson.fD.Pt(), dmeson.fD.Eta(), dmeson.fD.Phi_0_2pi(),
1406  d_closest));
1407  }
1408 
1409  return jet_distance_pair(jet_closest, d_closest);
1410 }
1411 
1414 {
1415  const Int_t nD = fCandidateArray->GetEntriesFast();
1416 
1417  AliDmesonJetInfo DmesonJet;
1418 
1419  Int_t nAccCharm[3] = {0};
1420  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1421  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1422  if (!charmCand) continue;
1423 
1424  // region of interest + cuts
1425  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1426  Int_t nMassHypo = 0; // number of mass hypothesis accepted for this D meson
1427  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1428  DmesonJet.Reset();
1429  DmesonJet.fDmesonParticle = charmCand;
1430  DmesonJet.fSelectionType = im + 1;
1431  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1432  for (auto& def : fJetDefinitions) {
1433  if (!FindJet(charmCand, DmesonJet, def)) {
1434  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1435  def.GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1436  }
1437  }
1438  fDmesonJets[(icharm+1)*(1-(im*2))] = DmesonJet;
1439  nMassHypo++;
1440  nAccCharm[im]++;
1441  }
1442  }
1443  if (nMassHypo == 2) {
1444  nAccCharm[0]--;
1445  nAccCharm[1]--;
1446  nAccCharm[2] += 2;
1447  }
1448  if (nMassHypo == 2) { // both mass hypothesis accepted
1449  fDmesonJets[(icharm+1)].fSelectionType = 3;
1450  fDmesonJets[-(icharm+1)].fSelectionType = 3;
1451  }
1452  } // end of D cand loop
1453 
1454  TString hname;
1455 
1456  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1457  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1458  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1459  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1460 
1461  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1462  fHistManager->FillTH2(hname, fTrackContainer->GetNAcceptedTracks(), nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1463 
1464  hname = TString::Format("%s/fHistNDmesons", GetName());
1465  fHistManager->FillTH1(hname, nD);
1466 }
1467 
1478 {
1479  TString hname;
1480 
1482  fFastJetWrapper->SetR(jetDef.fRadius);
1485 
1486  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1487 
1488  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1489  fTrackContainer->SetDMesonCandidate(Dcand);
1490  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1491  AddInputVectors(fTrackContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)), fTrackEfficiency);
1492 
1493  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1494  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1495  const TObjArray& daughters = fTrackContainer->GetDaughterList();
1496  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1497  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1498  if (!fTrackContainer->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1499  }
1500  }
1501 
1502  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1503  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1504  AddInputVectors(fClusterContainer, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1505  }
1506 
1507  // run jet finder
1508  fFastJetWrapper->Run();
1509 
1510  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1511 
1512  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1513  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1514 
1515  Bool_t isDmesonJet = kFALSE;
1516 
1517  Double_t maxChPt = 0;
1518  Double_t maxNePt = 0;
1519  Double_t totalNeutralPt = 0;
1520 
1521  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1522  if (constituents[ic].user_index() == 0) {
1523  isDmesonJet = kTRUE;
1524  }
1525  else if (constituents[ic].user_index() >= 100) {
1526  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1527  }
1528  else if (constituents[ic].user_index() <= -100) {
1529  totalNeutralPt += constituents[ic].pt();
1530  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1531  }
1532  }
1533 
1534  if (isDmesonJet) {
1535  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1536  DmesonJet.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1537  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1538  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1539  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1540 
1541  return kTRUE;
1542  }
1543  }
1544 
1545  return kFALSE;
1546 }
1547 
1551 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist, Double_t eff)
1552 {
1553  auto itcont = cont->all_momentum();
1554  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1555  UInt_t rejectionReason = 0;
1556  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1557  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1558  continue;
1559  }
1560  if (fRandomGen && eff > 0 && eff < 1) {
1561  Double_t rnd = fRandomGen->Rndm();
1562  if (eff < rnd) {
1563  if (rejectHist) rejectHist->Fill(6, it->first.Pt());
1564  continue;
1565  }
1566  }
1567  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1568  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1569  }
1570 }
1571 
1574 {
1575  TString hname;
1576 
1577  fMCContainer->SetSpecialPDG(fCandidatePDG);
1578  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1579  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1580 
1581  if (!fMCContainer->IsSpecialPDGFound()) return;
1582 
1583  Int_t nAccCharm[3] = {0};
1584 
1585  for (auto &jetDef : fJetDefinitions) {
1586  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", GetName(), jetDef.GetName());
1587  TH1* histNDmesonsVsNconstituents = static_cast<TH1*>(fHistManager->FindObject(hname));
1588 
1589  switch (jetDef.fJetType) {
1591  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
1592  break;
1594  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
1595  break;
1597  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
1598  break;
1599  }
1600 
1602  fFastJetWrapper->SetR(jetDef.fRadius);
1605 
1606  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1607  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1608 
1609  fFastJetWrapper->Run();
1610 
1611  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1612 
1613  for (auto jet : jets_incl) {
1614  Int_t nDmesonsInJet = 0;
1615 
1616  for (auto constituent : jet.constituents()) {
1617  Int_t iPart = constituent.user_index() - 100;
1618  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
1619  if (!part) {
1620  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1621  continue;
1622  }
1623  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1624  nDmesonsInJet++;
1625  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1626  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1627  std::pair<int, AliDmesonJetInfo> element;
1628  element.first = iPart;
1629  dMesonJetIt = fDmesonJets.insert(element).first;
1630  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1631  (*dMesonJetIt).second.fDmesonParticle = part;
1632  (*dMesonJetIt).second.fSelectionType = part->PdgCode() > 0 ? 1 : 2;
1633 
1634  UShort_t p = 0;
1635  UInt_t rs = 0;
1636 
1637  auto firstParton = CheckOrigin(part, fMCContainer->GetArray(), kTRUE);
1638  p = 0;
1639  rs = firstParton.first;
1640  while (rs >>= 1) { p++; }
1641  (*dMesonJetIt).second.fFirstPartonType = p;
1642  (*dMesonJetIt).second.fFirstParton = firstParton.second;
1643 
1644  auto lastParton = CheckOrigin(part, fMCContainer->GetArray(), kFALSE);
1645  p = 0;
1646  rs = lastParton.first;
1647  while (rs >>= 1) { p++; }
1648  (*dMesonJetIt).second.fLastPartonType = p;
1649  (*dMesonJetIt).second.fLastParton = lastParton.second;
1650 
1651  if (part->PdgCode() > 0) {
1652  nAccCharm[0]++;
1653  }
1654  else {
1655  nAccCharm[1]++;
1656  }
1657  }
1658 
1659  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
1660  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = jet.constituents().size();
1661  } // if constituent is a D meson
1662  } // for each constituent
1663  if (nDmesonsInJet > 0) histNDmesonsVsNconstituents->Fill(jet.constituents().size(), nDmesonsInJet);
1664  } // for each jet
1665  } // for each jet definition
1666 
1667  if (fDmesonJets.size() != nAccCharm[0]+nAccCharm[1]) AliError(Form("I found %lu mesons (%d)?", fDmesonJets.size(), nAccCharm[0]+nAccCharm[1]));
1668  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1669  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1670  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1671  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1672 
1673  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1674  fHistManager->FillTH2(hname, fMCContainer->GetNAcceptedParticles(), nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1675 
1676  hname = TString::Format("%s/fHistNDmesons", GetName());
1677  fHistManager->FillTH1(hname, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]); // same as the number of accepted D mesons, since no selection is performed
1678 }
1679 
1684 {
1685  TString classname;
1686  if (fMCMode == kMCTruth) {
1687  classname = "AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary";
1688  fCurrentDmesonJetInfo = new AliDmesonMCInfoSummary();
1689  }
1690  else {
1691  switch (fCandidateType) {
1692  case kD0toKpi:
1693  case kD0toKpiLikeSign:
1694  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1695  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1696  break;
1697  case kDstartoKpipi:
1698  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1699  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1700  break;
1701  }
1702  }
1703  TString treeName = TString::Format("%s_%s", taskName, GetName());
1704  fTree = new TTree(treeName, treeName);
1705  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
1706  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1707  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1708  fCurrentJetInfo[i] = new AliJetInfoSummary();
1709  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1710  }
1711 
1712  return fTree;
1713 }
1714 
1719 {
1720  TString hname;
1721 
1722  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
1723 
1724  for (auto &jetDef : fJetDefinitions) {
1725 
1726  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
1727 
1728  Double_t radius = jetDef.fRadius;
1729 
1730  TString title[30] = {""};
1731  Int_t nbins[30] = {0 };
1732  Double_t min [30] = {0.};
1733  Double_t max [30] = {0.};
1734  Int_t dim = 0 ;
1735 
1736  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1737  nbins[dim] = nPtBins;
1738  min[dim] = 0;
1739  max[dim] = fMaxPt;
1740  dim++;
1741 
1742  if ((enabledAxis & kPositionD) != 0) {
1743  title[dim] = "#eta_{D}";
1744  nbins[dim] = 50;
1745  min[dim] = -1;
1746  max[dim] = 1;
1747  dim++;
1748 
1749  title[dim] = "#phi_{D} (rad)";
1750  nbins[dim] = 150;
1751  min[dim] = 0;
1752  max[dim] = TMath::TwoPi();
1753  dim++;
1754  }
1755 
1756  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1757  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1758  nbins[dim] = fNMassBins;
1759  min[dim] = fMinMass;
1760  max[dim] = fMaxMass;
1761  dim++;
1762  }
1763 
1764  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1765  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1766  nbins[dim] = fNMassBins;
1767  min[dim] = fMinMass;
1768  max[dim] = fMaxMass;
1769  dim++;
1770  }
1771 
1772  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1773  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1774  nbins[dim] = fNMassBins;
1775  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1776  dim++;
1777  }
1778 
1779  if (fCandidateType == kDstartoKpipi) {
1780  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1781  nbins[dim] = fNMassBins*6;
1782  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1783 
1784  // subtract mass of D0
1785  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1786  min[dim] -= D0mass;
1787  max[dim] -= D0mass;
1788 
1789  dim++;
1790  }
1791 
1792  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1793  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1794  nbins[dim] = 100;
1795  min[dim] = 0;
1796  max[dim] = 25;
1797  dim++;
1798  }
1799 
1800  title[dim] = "#it{z}_{D}";
1801  nbins[dim] = 110;
1802  min[dim] = 0;
1803  max[dim] = 1.10;
1804  dim++;
1805 
1806  if ((enabledAxis & kDeltaR) != 0) {
1807  title[dim] = "#Delta R_{D-jet}";
1808  nbins[dim] = 100;
1809  min[dim] = 0;
1810  max[dim] = radius * 1.5;
1811  dim++;
1812  }
1813 
1814  if ((enabledAxis & kDeltaEta) != 0) {
1815  title[dim] = "#eta_{D} - #eta_{jet}";
1816  nbins[dim] = 100;
1817  min[dim] = -radius * 1.2;
1818  max[dim] = radius * 1.2;
1819  dim++;
1820  }
1821 
1822  if ((enabledAxis & kDeltaPhi) != 0) {
1823  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1824  nbins[dim] = 100;
1825  min[dim] = -radius * 1.2;
1826  max[dim] = radius * 1.2;
1827  dim++;
1828  }
1829 
1830  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1831  nbins[dim] = nPtBins;
1832  min[dim] = 0;
1833  max[dim] = fMaxPt;
1834  dim++;
1835 
1836  if ((enabledAxis & kPositionJet) != 0) {
1837  title[dim] = "#eta_{jet}";
1838  nbins[dim] = 50;
1839  min[dim] = -1;
1840  max[dim] = 1;
1841  dim++;
1842 
1843  title[dim] = "#phi_{jet} (rad)";
1844  nbins[dim] = 150;
1845  min[dim] = 0;
1846  max[dim] = TMath::TwoPi();
1847  dim++;
1848  }
1849 
1850  if ((enabledAxis & kJetConstituents) != 0) {
1851  title[dim] = "No. of constituents";
1852  nbins[dim] = 50;
1853  min[dim] = -0.5;
1854  max[dim] = 49.5;
1855  dim++;
1856  }
1857 
1858  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1859  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1860  for (Int_t j = 0; j < dim; j++) {
1861  h->GetAxis(j)->SetTitle(title[j]);
1862  }
1863  }
1864 }
1865 
1870 {
1871  TString hname;
1872  fFirstPartons.clear();
1873  fLastPartons.clear();
1874  for (auto& dmeson_pair : fDmesonJets) {
1875  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
1876  Int_t accJets = 0;
1877  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1878  fCurrentJetInfo[ij]->Reset();
1879  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1880  if (!jet) continue;
1881  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1882  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1883  fHistManager->FillTH1(hname, jet->Pt());
1884  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1885  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1886  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1887  fHistManager->FillTH1(hname, jet->Eta());
1888  continue;
1889  }
1890  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
1891  accJets++;
1892  }
1893  if (accJets > 0) {
1894  fFirstPartons[dmeson_pair.second.fFirstParton] = dmeson_pair.second.fFirstPartonType;
1895  fLastPartons[dmeson_pair.second.fLastParton] = dmeson_pair.second.fLastPartonType;
1896 
1897  fTree->Fill();
1898  }
1899  else {
1900  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1901  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1902  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1903  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1904  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1905  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1906  if (fMCMode != kMCTruth) {
1907  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1908  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1909  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1910  }
1911  else if (fCandidateType == kDstartoKpipi) {
1912  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1913  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1914 
1915  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1916  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1917  }
1918  }
1919  }
1920  }
1921 
1922  hname = TString::Format("%s/fHistFirstPartonPt", GetName());
1923  TH1* histFirstPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
1924  hname = TString::Format("%s/fHistFirstPartonEta", GetName());
1925  TH1* histFirstPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
1926  hname = TString::Format("%s/fHistFirstPartonPhi", GetName());
1927  TH1* histFirstPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
1928  hname = TString::Format("%s/fHistFirstPartonType", GetName());
1929  TH1* histFirstPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
1930 
1931  for (auto parton : fFirstPartons) {
1932  if (!parton.first) continue;
1933  histFirstPartonPt->Fill(parton.first->Pt());
1934  histFirstPartonEta->Fill(parton.first->Eta());
1935  histFirstPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
1936  histFirstPartonType->Fill(parton.second);
1937  }
1938 
1939  hname = TString::Format("%s/fHistLastPartonPt", GetName());
1940  TH1* histLastPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
1941  hname = TString::Format("%s/fHistLastPartonEta", GetName());
1942  TH1* histLastPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
1943  hname = TString::Format("%s/fHistLastPartonPhi", GetName());
1944  TH1* histLastPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
1945  hname = TString::Format("%s/fHistLastPartonType", GetName());
1946  TH1* histLastPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
1947 
1948  for (auto parton : fLastPartons) {
1949  if (!parton.first) continue;
1950  histLastPartonPt->Fill(parton.first->Pt());
1951  histLastPartonEta->Fill(parton.first->Eta());
1952  histLastPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
1953  histLastPartonType->Fill(parton.second);
1954  }
1955 
1956  return kTRUE;
1957 }
1958 
1964 {
1965  TString hname;
1966  fFirstPartons.clear();
1967  fLastPartons.clear();
1968  for (auto& dmeson_pair : fDmesonJets) {
1969  Int_t accJets = 0;
1970  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1971  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1972  if (!jet) continue;
1973  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1974  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1975  fHistManager->FillTH1(hname, jet->Pt());
1976  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1977  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1978  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1979  fHistManager->FillTH1(hname, jet->Eta());
1980  continue;
1981  }
1982  accJets++;
1983  }
1984  if (accJets > 0) {
1985  fFirstPartons[dmeson_pair.second.fFirstParton] = dmeson_pair.second.fFirstPartonType;
1986  fLastPartons[dmeson_pair.second.fLastParton] = dmeson_pair.second.fLastPartonType;
1987  }
1988  else {
1989  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1990  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1991  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1992  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1993  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1994  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1995  if (fMCMode != kMCTruth) {
1996  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1997  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1998  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1999  }
2000  else if (fCandidateType == kDstartoKpipi) {
2001  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
2002  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
2003 
2004  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
2005  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
2006  }
2007  }
2008  }
2009  }
2010 
2011  hname = TString::Format("%s/fHistFirstPartonPt", GetName());
2012  TH1* histFirstPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
2013  hname = TString::Format("%s/fHistFirstPartonEta", GetName());
2014  TH1* histFirstPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
2015  hname = TString::Format("%s/fHistFirstPartonPhi", GetName());
2016  TH1* histFirstPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
2017  hname = TString::Format("%s/fHistFirstPartonType", GetName());
2018  TH1* histFirstPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
2019 
2020  for (auto parton : fFirstPartons) {
2021  if (!parton.first) continue;
2022  histFirstPartonPt->Fill(parton.first->Pt());
2023  histFirstPartonEta->Fill(parton.first->Eta());
2024  histFirstPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
2025  histFirstPartonType->Fill(parton.second);
2026  }
2027 
2028  hname = TString::Format("%s/fHistLastPartonPt", GetName());
2029  TH1* histLastPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
2030  hname = TString::Format("%s/fHistLastPartonEta", GetName());
2031  TH1* histLastPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
2032  hname = TString::Format("%s/fHistLastPartonPhi", GetName());
2033  TH1* histLastPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
2034  hname = TString::Format("%s/fHistLastPartonType", GetName());
2035  TH1* histLastPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
2036 
2037  for (auto parton : fLastPartons) {
2038  if (!parton.first) continue;
2039  histLastPartonPt->Fill(parton.first->Pt());
2040  histLastPartonEta->Fill(parton.first->Eta());
2041  histLastPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
2042  histLastPartonType->Fill(parton.second);
2043  }
2044 
2045  return kTRUE;
2046 }
2047 
2052 {
2053  TString hname;
2054 
2055  for (auto& dmeson_pair : fDmesonJets) {
2056  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
2057  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2058  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2059  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2060  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2061  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2062  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2063  }
2064  }
2065 
2066  for (auto &jetDef : fJetDefinitions) {
2067 
2068  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
2069  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
2070 
2071  for (auto& dmeson_pair : fDmesonJets) {
2072  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
2073  if (!jet) continue;
2074  if (!jetDef.IsJetInAcceptance(*jet)) {
2075  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
2076  fHistManager->FillTH1(hname, jet->Pt());
2077  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
2078  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2079  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
2080  fHistManager->FillTH1(hname, jet->Eta());
2081  continue;
2082  }
2083  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
2084  }
2085  }
2086 
2087  return kTRUE;
2088 }
2089 
2096 {
2097  // Fill the THnSparse histogram.
2098 
2099  Double_t contents[30] = {0.};
2100 
2101  Double_t z = DmesonJet.GetZ(n);
2102  Double_t deltaPhi = 0;
2103  Double_t deltaEta = 0;
2104  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
2105 
2106  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
2107  if (it == DmesonJet.fJets.end()) return kFALSE;
2108 
2109  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
2110  TString title(h->GetAxis(i)->GetTitle());
2111  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
2112  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
2113  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
2114  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
2115  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
2116  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
2117  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
2118  else if (title=="#it{z}_{D}") contents[i] = z;
2119  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
2120  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
2121  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
2122  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
2123  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
2124  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
2125  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
2126  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
2127  }
2128 
2129  h->Fill(contents);
2130 
2131  return kTRUE;
2132 }
2133 
2134 // Definitions of class AliAnalysisTaskDmesonJets
2135 
2139 
2143  fAnalysisEngines(),
2144  fEnabledAxis(0),
2146  fHistManager(),
2147  fApplyKinematicCuts(kTRUE),
2148  fNOutputTrees(0),
2149  fTrackEfficiency(0),
2150  fAodEvent(0),
2151  fFastJetWrapper(0),
2152  fMCContainer(0)
2153 {
2154  SetMakeGeneralHistograms(kTRUE);
2155 }
2156 
2161  AliAnalysisTaskEmcalLight(name, kTRUE),
2162  fAnalysisEngines(),
2163  fEnabledAxis(k2ProngInvMass),
2164  fOutputType(kTreeOutput),
2165  fHistManager(name),
2166  fApplyKinematicCuts(kTRUE),
2167  fNOutputTrees(nOutputTrees),
2168  fTrackEfficiency(0),
2169  fAodEvent(0),
2170  fFastJetWrapper(0),
2171  fMCContainer(0)
2172 {
2173  SetMakeGeneralHistograms(kTRUE);
2174  for (Int_t i = 0; i < nOutputTrees; i++){
2175  DefineOutput(2+i, TTree::Class());
2176  }
2177 }
2178 
2181 {
2182  if (fFastJetWrapper) delete fFastJetWrapper;
2183 }
2184 
2192 {
2193  AliRDHFCuts* analysiscuts = 0;
2194  TFile* filecuts = TFile::Open(cutfname);
2195  if (!filecuts || filecuts->IsZombie()) {
2196  ::Warning("AddTaskDmesonJets", "Input file not found: will use std cuts.");
2197  filecuts = 0;
2198  }
2199 
2200  if (filecuts) {
2201  analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
2202  if (!analysiscuts) {
2203  ::Warning("AddTaskDmesonJetCorr", "Could not find analysis cuts '%s' in '%s'. Using std cuts.", cutsname.Data(), cutfname.Data());
2204  }
2205  }
2206 
2207  return analysiscuts;
2208 }
2209 
2219 {
2221  return AddAnalysisEngine(type, MCmode, jetDef, cutfname);
2222 }
2223 
2233 {
2234  AliRDHFCuts* cuts = 0;
2235 
2236  if (!cutfname.IsNull()) {
2237  TString cutsname;
2238 
2239  switch (type) {
2240  case kD0toKpi:
2241  case kD0toKpiLikeSign:
2242  cutsname = "D0toKpiCuts";
2243  break;
2244  case kDstartoKpipi:
2245  cutsname = "DStartoKpipiCuts";
2246  break;
2247  default:
2248  return 0;
2249  }
2250 
2251  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
2252  }
2253 
2254  AnalysisEngine eng(type, MCmode, cuts);
2255 
2256  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
2257 
2258  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
2259  eng.AddJetDefinition(jetDef);
2260  it = fAnalysisEngines.insert(it, eng);
2261  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
2262  }
2263  else {
2264  AnalysisEngine* found_eng = &(*it);
2265  ::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());
2266  found_eng->AddJetDefinition(jetDef);
2267 
2268  if (cuts) {
2269  if (found_eng->fRDHFCuts != 0) ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
2270  found_eng->SetRDHFCuts(cuts);
2271  }
2272  }
2273 
2274  return &(*it);
2275 }
2276 
2277 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
2278 {
2279  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
2280  while (it != fAnalysisEngines.end() && (*it) < eng) it++;
2281  return it;
2282 }
2283 
2286 {
2287  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
2288 
2290 
2291  // Define histograms
2292  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
2293 
2294  TString hname;
2295  TString htitle;
2296  TH1* h = 0;
2297  Int_t treeSlot = 0;
2298 
2299  hname = "fHistCharmPt";
2300  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2301  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2302 
2303  hname = "fHistCharmEta";
2304  htitle = hname + ";#eta_{charm};counts";
2305  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2306 
2307  hname = "fHistCharmPhi";
2308  htitle = hname + ";#phi_{charm};counts";
2309  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2310 
2311  hname = "fHistCharmPt_Eta05";
2312  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2313  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2314 
2315  hname = "fHistBottomPt";
2316  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2317  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2318 
2319  hname = "fHistBottomEta";
2320  htitle = hname + ";#eta_{bottom};counts";
2321  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2322 
2323  hname = "fHistBottomPhi";
2324  htitle = hname + ";#phi_{bottom};counts";
2325  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2326 
2327  hname = "fHistBottomPt_Eta05";
2328  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2329  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2330 
2331  hname = "fHistHighestPartonPt";
2332  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2333  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2334 
2335  hname = "fHistHighestPartonType";
2336  htitle = hname + ";type;counts";
2337  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2338 
2339  hname = "fHistNHeavyQuarks";
2340  htitle = hname + ";number of heavy-quarks;counts";
2341  fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
2342 
2343  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
2344  for (auto &param : fAnalysisEngines) {
2345  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
2346 
2347  param.fHistManager = &fHistManager;
2348 
2349  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", param.GetName());
2350  htitle = hname + ";#it{N}_{tracks};#it{N}_{D};events";
2351  h = fHistManager.CreateTH2(hname, htitle, 251, -0.5, 250.5, 21, -0.5, 20.5);
2352 
2353  hname = TString::Format("%s/fHistNTotAcceptedDmesons", param.GetName());
2354  htitle = hname + ";;#it{N}_{D}";
2355  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
2356 
2357  hname = TString::Format("%s/fHistNDmesons", param.GetName());
2358  htitle = hname + ";#it{N}_{D};events";
2359  h = fHistManager.CreateTH1(hname, htitle, 501, -0.5, 500.5);
2360 
2361  hname = TString::Format("%s/fHistNEvents", param.GetName());
2362  htitle = hname + ";Event status;counts";
2363  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
2364  h->GetXaxis()->SetBinLabel(1, "Accepted");
2365  h->GetXaxis()->SetBinLabel(2, "Rejected");
2366 
2367  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
2368  htitle = hname + ";Rejection reason;counts";
2369  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
2370 
2371  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
2372  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
2373  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2374 
2375  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
2376  htitle = hname + ";#it{#eta}_{D};counts";
2377  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2378 
2379  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
2380  htitle = hname + ";#it{#phi}_{D};counts";
2381  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2382 
2383  if (param.fMCMode != kMCTruth) {
2384  if (param.fCandidateType == kD0toKpi || param.fCandidateType == kD0toKpiLikeSign) {
2385  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
2386  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2387  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
2388  }
2389  else if (param.fCandidateType == kDstartoKpipi) {
2390  Double_t min = 0;
2391  Double_t max = 0;
2392 
2393  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
2394  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2395  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
2396  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
2397 
2398  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2399  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
2400  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2401  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
2402  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
2403  }
2404  }
2405 
2406  if (param.fMCMode == kMCTruth) {
2407  hname = TString::Format("%s/fHistFirstPartonPt", param.GetName());
2408  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
2409  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2410 
2411  hname = TString::Format("%s/fHistFirstPartonEta", param.GetName());
2412  htitle = hname + ";#eta_{parton};counts";
2413  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2414 
2415  hname = TString::Format("%s/fHistFirstPartonPhi", param.GetName());
2416  htitle = hname + ";#phi_{parton};counts";
2417  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2418 
2419  hname = TString::Format("%s/fHistFirstPartonType", param.GetName());
2420  htitle = hname + ";type;counts";
2421  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2422 
2423  hname = TString::Format("%s/fHistLastPartonPt", param.GetName());
2424  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
2425  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2426 
2427  hname = TString::Format("%s/fHistLastPartonEta", param.GetName());
2428  htitle = hname + ";#eta_{parton};counts";
2429  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2430 
2431  hname = TString::Format("%s/fHistLastPartonPhi", param.GetName());
2432  htitle = hname + ";#phi_{parton};counts";
2433  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2434 
2435  hname = TString::Format("%s/fHistLastPartonType", param.GetName());
2436  htitle = hname + ";type;counts";
2437  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2438  }
2439 
2440  for (auto& jetDef : param.fJetDefinitions) {
2441  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef.GetName());
2442 
2443  if (param.fMCMode == kMCTruth) {
2444  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", param.GetName(), jetDef.GetName());
2445  htitle = hname + ";#it{N}_{constituents};#it{N}_{D};counts";
2446  h = fHistManager.CreateTH2(hname, htitle, 51, -0.5, 50.5, 10, 0.5, 10.5);
2447  }
2448 
2449  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef.GetName());
2450  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2451  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2452  SetRejectionReasonLabels(h->GetXaxis());
2453 
2454  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef.GetName());
2455  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2456  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2457  SetRejectionReasonLabels(h->GetXaxis());
2458 
2459  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef.GetName());
2460  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
2461  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2462  SetRejectionReasonLabels(h->GetXaxis());
2463 
2464  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef.GetName());
2465  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
2466  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
2467 
2468  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef.GetName());
2469  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
2470  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2471 
2472  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef.GetName());
2473  htitle = hname + ";#it{#eta}_{jet};counts";
2474  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2475 
2476  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef.GetName());
2477  htitle = hname + ";#it{#phi}_{jet};counts";
2478  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2479  }
2480  switch (fOutputType) {
2481  case kTreeOutput:
2482  param.BuildTree(GetName());
2483  if (treeSlot < fNOutputTrees) {
2484  param.AssignDataSlot(treeSlot+2);
2485  treeSlot++;
2487  }
2488  else {
2489  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
2490  }
2491  break;
2492  case kTHnOutput:
2493  param.BuildHnSparse(fEnabledAxis);
2494  break;
2495  case kNoOutput:
2496  break;
2497  }
2498  }
2499 
2501 
2502  PostData(1, fOutput);
2503 }
2504 
2508 {
2510 
2511  // Load the event
2512  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
2513 
2514  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
2515 
2516  fFastJetWrapper->SetAreaType(fastjet::active_area);
2518 
2519  if (!fAodEvent) {
2520  AliError(Form("This task need an AOD event (Task '%s'). Expect troubles...", GetName()));
2521  //return;
2522  }
2523 
2526 
2527  TRandom* rnd = 0;
2528  if (fTrackEfficiency > 0 && fTrackEfficiency < 1) rnd = new TRandom3(0);
2529 
2530  for (auto &params : fAnalysisEngines) {
2531 
2532  params.fAodEvent = fAodEvent;
2533  params.fFastJetWrapper = fFastJetWrapper;
2534  params.fTrackEfficiency = fTrackEfficiency;
2535  params.fRandomGen = rnd;
2536 
2537  if (fAodEvent) params.Init(fGeom, fAodEvent->GetRunNumber());
2538 
2539  if (params.fMCMode != kMCTruth && fAodEvent) {
2540  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
2541 
2542  if (params.fCandidateArray) {
2543  TString className;
2544  if (params.fCandidateType == kD0toKpi || params.fCandidateType == kD0toKpiLikeSign) {
2545  className = "AliAODRecoDecayHF2Prong";
2546  }
2547  else if (params.fCandidateType == kDstartoKpipi) {
2548  className = "AliAODRecoCascadeHF";
2549  }
2550  if (!params.fCandidateArray->GetClass()->InheritsFrom(className)) {
2551  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2552  "%s: Objects of type %s in %s are not inherited from %s! Task will be disabled!",
2553  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName(), className.Data());
2554  params.fCandidateArray = 0;
2555  params.fInhibit = kTRUE;
2556  }
2557  }
2558  else {
2559  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2560  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
2561  params.fBranchName.Data(), params.GetName());
2562  params.fInhibit = kTRUE;
2563  }
2564  }
2565 
2566  if (params.fMCMode != kNoMC) {
2567  if (!fMCContainer) {
2568  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2569  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
2570  params.GetName());
2571  params.fInhibit = kTRUE;
2572  }
2573  params.fMCContainer = fMCContainer;
2574  }
2575 
2576  if (params.fMCMode != kMCTruth) {
2577  params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(0));
2578  if (!params.fTrackContainer) params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(1));
2579 
2580  params.fClusterContainer = GetClusterContainer(0);
2581 
2582  if (!params.fTrackContainer && !params.fClusterContainer) {
2583  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2584  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
2585  params.GetName());
2586  params.fInhibit = kTRUE;
2587  }
2588  }
2589  }
2590 }
2591 
2596 {
2597  TString hname;
2598 
2599  // Fix for temporary bug in ESDfilter
2600  // The AODs with null vertex pointer didn't pass the PhysSel
2601  // Now adding an entry in the histogram so as to check that this is actually cutting anything out
2602  if (fAodEvent && (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001)) {
2603  for (auto &eng : fAnalysisEngines) {
2604  if (eng.fInhibit) continue;
2605  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2606  fHistManager.FillTH1(hname, "ESDfilterBug");
2607  }
2608  return kFALSE;
2609  }
2610 
2611  for (auto &eng : fAnalysisEngines) {
2612  eng.fDmesonJets.clear();
2613  if (eng.fInhibit) continue;
2614 
2615  //Event selection
2616  hname = TString::Format("%s/fHistNEvents", eng.GetName());
2617  if (fAodEvent) {
2618  Bool_t iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
2619  if (!iseventselected) {
2620  fHistManager.FillTH1(hname, "Rejected");
2621  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2622  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
2623  TString label;
2624  do {
2625  label = GetHFEventRejectionReasonLabel(bitmap);
2626  if (label.IsNull()) break;
2627  fHistManager.FillTH1(hname, label);
2628  } while (true);
2629 
2630  continue;
2631  }
2632  }
2633 
2634  fHistManager.FillTH1(hname, "Accepted");
2635 
2636  AliDebug(2, "Event selected");
2637 
2638  eng.RunAnalysis();
2639  }
2640  return kTRUE;
2641 }
2642 
2647 {
2648  for (auto &param : fAnalysisEngines) {
2649  if (param.fInhibit) continue;
2650 
2651  if (fOutputType == kTreeOutput) {
2652  param.FillTree(fApplyKinematicCuts);
2654  }
2655  else if (fOutputType == kTHnOutput) {
2656  param.FillHnSparse(fApplyKinematicCuts);
2657  }
2658  }
2660  return kTRUE;
2661 }
2662 
2665 {
2666  auto itcont = fMCContainer->all_momentum();
2667  Int_t nHQ = 0;
2668  Double_t highestPartonPt = 0;
2669  Int_t absPdgHighParton = 0;
2670  for (auto part : itcont) {
2671  Int_t absPdgCode = TMath::Abs(part.second->GetPdgCode());
2672 
2673  // Skip all particles that are not either quarks or gluons
2674  if (absPdgCode > 9 && absPdgCode != 21) continue;
2675 
2676  // Look for highest momentum parton
2677  if (highestPartonPt < part.first.Pt()) {
2678  highestPartonPt = part.first.Pt();
2679  absPdgHighParton = absPdgCode;
2680  }
2681  /*
2682  // Look for the mother PDG code
2683  Int_t motherIndex = part.second->GetMother();
2684  AliAODMCParticle *mother = 0;
2685  Int_t motherPdg = 0;
2686  Double_t motherPt = 0;
2687  if (motherIndex >= 0) {
2688  mother = fMCContainer->GetMCParticle(motherIndex);
2689  if (motherIndex) {
2690  motherPdg = TMath::Abs(mother->GetPdgCode());
2691  motherPt = mother->Pt();
2692  }
2693  }
2694  */
2695  if (absPdgCode != 4 && absPdgCode != 5) continue;
2696  Bool_t notLastInPartonShower = kFALSE;
2697  for (Int_t idaugh = 0; idaugh < 2; idaugh++){
2698  Int_t daughterIndex = part.second->GetDaughter(idaugh);
2699  if (daughterIndex < 0) {
2700  AliDebug(10, Form("Could not find daughter of heavy quark (pdg=%d, pt=%.3f)!", absPdgCode, part.first.Pt()));
2701  continue;
2702  }
2703  AliAODMCParticle *daughter = fMCContainer->GetMCParticle(daughterIndex);
2704  if (!daughter) {
2705  AliDebug(10, Form("Could not find daughter %d of heavy quark (pdg=%d, pt=%.3f)!", daughterIndex, absPdgCode, part.first.Pt()));
2706  continue;
2707  }
2708  Int_t daughterAbsPdgCode = TMath::Abs(daughter->GetPdgCode());
2709  if (daughterAbsPdgCode <= 9 || daughterAbsPdgCode == 21) notLastInPartonShower = kTRUE; // this parton is not the last parton in the shower
2710  AliDebug(10, Form("Found daughter with PDG=%d, pt=%.3f", daughterAbsPdgCode, daughter->Pt()));
2711  }
2712  if (notLastInPartonShower) continue;
2713 
2714  if (absPdgCode == 4) {
2715  fHistManager.FillTH1("fHistCharmPt", part.first.Pt());
2716  fHistManager.FillTH1("fHistCharmEta", part.first.Eta());
2717  fHistManager.FillTH1("fHistCharmPhi", part.first.Phi_0_2pi());
2718  if (TMath::Abs(part.first.Eta()) < 0.5) fHistManager.FillTH1("fHistCharmPt_Eta05", part.first.Pt());
2719  }
2720  else if (absPdgCode == 5) {
2721  fHistManager.FillTH1("fHistBottomPt", part.first.Pt());
2722  fHistManager.FillTH1("fHistBottomEta", part.first.Eta());
2723  fHistManager.FillTH1("fHistBottomPhi", part.first.Phi_0_2pi());
2724  if (TMath::Abs(part.first.Eta()) < 0.5) fHistManager.FillTH1("fHistBottomPt_Eta05", part.first.Pt());
2725  }
2726  nHQ++;
2727  }
2728  fHistManager.FillTH1("fHistNHeavyQuarks", nHQ);
2729  fHistManager.FillTH1("fHistHighestPartonPt",highestPartonPt);
2730  Int_t partonType = 0;
2731  if (absPdgHighParton == 9 || absPdgHighParton == 21) {
2732  partonType = 7;
2733  }
2734  else {
2735  partonType = absPdgHighParton;
2736  }
2737  fHistManager.FillTH1("fHistHighestPartonType",partonType);
2738 }
2739 
2748 {
2749  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
2750 
2751  Double_t mass = part->Mass();
2752 
2753  // To make sure that the PDG mass value is not at the edge of a bin
2754  if (nbins % 2 == 0) {
2755  minMass = mass - range / 2 - range / nbins / 2;
2756  maxMass = mass + range / 2 - range / nbins / 2;
2757  }
2758  else {
2759  minMass = mass - range / 2;
2760  maxMass = mass + range / 2;
2761  }
2762 }
2763 
2770 {
2771  static TString label;
2772  label = "";
2773 
2774  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
2775  label = "NotSelTrigger";
2776  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
2777  return label.Data();
2778  }
2779  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
2780  label = "NoVertex";
2781  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
2782  return label.Data();
2783  }
2784  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
2785  label = "TooFewVtxContrib";
2786  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
2787  return label.Data();
2788  }
2789  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
2790  label = "ZVtxOutFid";
2791  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
2792  return label.Data();
2793  }
2794  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
2795  label = "Pileup";
2796  bitmap &= ~BIT(AliRDHFCuts::kPileup);
2797  return label.Data();
2798  }
2799  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
2800  label = "OutsideCentrality";
2801  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
2802  return label.Data();
2803  }
2804  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
2805  label = "PhysicsSelection";
2806  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
2807  return label.Data();
2808  }
2809  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
2810  label = "BadSPDVertex";
2811  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
2812  return label.Data();
2813  }
2814  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
2815  label = "ZVtxSPDOutFid";
2816  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
2817  return label.Data();
2818  }
2819  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
2820  label = "CentralityFlattening";
2821  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
2822  return label.Data();
2823  }
2824  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
2825  label = "BadTrackV0Correl";
2826  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
2827  return label.Data();
2828  }
2829 
2830  return label.Data();
2831 }
2832 
2839 {
2840  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
2841  PostData(eng.GetDataSlotNumber(), eng.GetTree());
2842  return eng.GetDataSlotNumber();
2843  }
2844  else {
2845  return -1;
2846  }
2847 }
2848 
2858 {
2859  // Get the pointer to the existing analysis manager via the static access method
2860  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
2861  if (!mgr) {
2862  ::Error("AddTaskDmesonJets", "No analysis manager to connect to.");
2863  return NULL;
2864  }
2865 
2866  // Check the analysis type using the event handlers connected to the analysis manager
2867  AliVEventHandler* handler = mgr->GetInputEventHandler();
2868  if (!handler) {
2869  ::Error("AddTaskEmcalJetSpectraQA", "This task requires an input event handler");
2870  return NULL;
2871  }
2872 
2873  EDataType_t dataType = kUnknownDataType;
2874 
2875  if (handler->InheritsFrom("AliESDInputHandler")) {
2876  dataType = kESD;
2877  }
2878  else if (handler->InheritsFrom("AliAODInputHandler")) {
2879  dataType = kAOD;
2880  }
2881 
2882  // Init the task and do settings
2883  if (ntracks == "usedefault") {
2884  if (dataType == kESD) {
2885  ntracks = "Tracks";
2886  }
2887  else if (dataType == kAOD) {
2888  ntracks = "tracks";
2889  }
2890  else {
2891  ntracks = "";
2892  }
2893  }
2894 
2895  if (nclusters == "usedefault") {
2896  if (dataType == kESD) {
2897  nclusters = "CaloClusters";
2898  }
2899  else if (dataType == kAOD) {
2900  nclusters = "caloClusters";
2901  }
2902  else {
2903  nclusters = "";
2904  }
2905  }
2906 
2907  if (nMCpart == "usedefault") {
2908  nMCpart = "mcparticles"; // Always needs AliAODMCParticle objects
2909  }
2910 
2911  TString name("AliAnalysisTaskDmesonJets");
2912  if (strcmp(suffix, "") != 0) {
2913  name += TString::Format("_%s", suffix.Data());
2914  }
2915 
2916  AliAnalysisTaskDmesonJets* jetTask = new AliAnalysisTaskDmesonJets(name, nMaxTrees);
2917 
2918  if (!ntracks.IsNull()) {
2919  AliHFTrackContainer* trackCont = new AliHFTrackContainer(ntracks);
2920  jetTask->AdoptParticleContainer(trackCont);
2921  }
2922 
2923  if (!nMCpart.IsNull()) {
2924  AliMCParticleContainer* partCont = new AliHFAODMCParticleContainer(nMCpart);
2925  partCont->SetEtaLimits(-1.5, 1.5);
2926  partCont->SetPtLimits(0, 1000);
2927  jetTask->AdoptParticleContainer(partCont);
2928  }
2929 
2930  jetTask->AddClusterContainer(nclusters);
2931 
2932  // Final settings, pass to manager and set the containers
2933  mgr->AddTask(jetTask);
2934 
2935  // Create containers for input/output
2936  AliAnalysisDataContainer* cinput1 = mgr->GetCommonInputContainer();
2937  TString contname1(name);
2938  contname1 += "_histos";
2939  AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(contname1.Data(),
2940  TList::Class(), AliAnalysisManager::kOutputContainer,
2941  Form("%s", AliAnalysisManager::GetCommonFileName()));
2942 
2943  mgr->ConnectInput(jetTask, 0, cinput1);
2944  mgr->ConnectOutput(jetTask, 1, coutput1);
2945 
2946  for (Int_t i = 0; i < nMaxTrees; i++) {
2947  TString contname = TString::Format("%s_tree_%d", name.Data(), i);
2948  AliAnalysisDataContainer *coutput = mgr->CreateContainer(contname.Data(),
2949  TTree::Class(),AliAnalysisManager::kOutputContainer,
2950  Form("%s", AliAnalysisManager::GetCommonFileName()));
2951  mgr->ConnectOutput(jetTask, 2+i, coutput);
2952  }
2953  return jetTask;
2954 }
2955 
Int_t pdg
void Print() const
Prints the content of this object in the standard output.
std::list< AnalysisEngine > fAnalysisEngines
Array of analysis parameters.
AliEMCALGeometry * fGeom
!emcal geometry
Analysis task for D meson jets.
virtual void Set(const AliDmesonJetInfo &source)
Invariant mass of the D0 meson candidate in GeV/c2.
virtual Int_t Run()
double Double_t
Definition: External.C:58
Bool_t IsJetInAcceptance(const AliJetInfo &jet) const
Double_t fSoftPionPt
! Transverse momentum of the soft pion of the D* candidate
EDataType_t
Switch for the data type.
virtual AliAODMCParticle * GetMCParticle(Int_t i=-1) const
const char * title
Definition: MakeQAPdf.C:26
AliRDHFCuts * fRDHFCuts
D meson candidates cuts.
virtual void Set(const AliDmesonJetInfo &source)
Lightweight class that encapsulates D meson jets.
Double_t fTrackEfficiency
Artificial tracking inefficiency (0...1)
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
AliHFAODMCParticleContainer * fMCContainer
! MC particle container
ERecoScheme_t fRecoScheme
Jet recombination scheme (pt scheme, E scheme, ...)
Double_t fMaxNeutralPt
Transverse momentum of the leading neutral particle (or cluster)
void RunDetectorLevelAnalysis()
Run a detector level analysis.
Double_t mass
jet_distance_pair FindJetMacthedToGeneratedDMeson(const AliDmesonJetInfo &dmeson, AliHFJetDefinition &jetDef, Double_t dMax, Bool_t applyKinCuts)
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
virtual void UserCreateOutputObjects()
Creates the output containers.
void FillPartonLevelHistograms()
Fill histograms with parton-level information.
Lightweight class that encapsulates D meson jets.
Double_t fInvMass2Prong
! 2-prong mass of the D* candidate (w/o the soft pion)
static void CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t &minMass, Double_t &maxMass)
AliAODEvent * fAodEvent
! AOD event
Bool_t FindJet(AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, AliHFJetDefinition &jetDef)
UInt_t fEnabledAxis
Use bit defined in EAxis_t to enable axis in the THnSparse.
AliAODMCParticle * fFirstParton
! pointer to the first parton in the shower tree of the D meson (only for particle level D mesons) ...
EMCMode_t fMCMode
MC mode: No MC (data and MC detector level), background-only (MC), signal-only (MC), MC truth (particle level)
Short_t fLastPartonType
! type of the last parton in the shower tree (only for particle level D mesons)
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:96
AnalysisEngine()
This is the default constructor, used for ROOT I/O purposes.
void Init(const AliEMCALGeometry *const geom, Int_t runNumber)
Initialize the analysis engine.
static FJRecoScheme ConvertToFJRecoScheme(ERecoScheme_t reco)
AliVParticle * fDmesonParticle
! pointer to the particle object
static TString GenerateJetName(EJetType_t jetType, EJetAlgo_t jetAlgo, ERecoScheme_t recoScheme, Double_t radius, AliParticleContainer *partCont, AliClusterContainer *clusCont, TString tag)
void SetR(Double_t r)
Definition: AliFJWrapper.h:101
Select tracks based on specific prescriptions of HF analysis.
ECandidateType_t fCandidateType
Candidate type.
Double_t GetDistance(std::string n, Double_t &deta, Double_t &dphi) const
Int_t fMCLabel
! MC label, i.e. index of the generator level D meson (only for detector level D meson candidates) ...
const Int_t nPtBins
virtual void Set(const AliDmesonJetInfo &source)
void AddInputVectors(AliEmcalContainer *cont, Int_t offset, TH2 *rejectHist=0, Double_t eff=0.)
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:33
AliClusterContainer * AddClusterContainer(const char *n)
EJetType_t fJetType
Jet type (charged, full, neutral)
std::pair< AliJetInfo *, Double_t > jet_distance_pair
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
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
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Definition: THistManager.h:504
AliAnalysisTaskDmesonJets()
This is the default constructor, used for ROOT I/O purposes.
virtual void Reset()
Reset the current object.
AnalysisEngine & operator=(const AnalysisEngine &source)
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
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:95
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
AliClusterContainer * GetClusterContainer(Int_t i=0) const
AliAODTrack * GetBachelor() const
AliRDHFCuts * LoadDMesonCutsFromFile(TString cutfname, TString cutsname)
bool operator<(const AliAnalysisTaskDmesonJets::AliHFJetDefinition &lhs, const AliAnalysisTaskDmesonJets::AliHFJetDefinition &rhs)
std::list< AnalysisEngine >::iterator FindAnalysisEngine(const AnalysisEngine &eng)
std::vector< AliHFJetDefinition > fJetDefinitions
Jet definitions.
EJetAlgo_t fJetAlgo
Jet algorithm (kt, anti-kt,...)
static const char * GetHFEventRejectionReasonLabel(UInt_t &bitmap)
Lightweight class that encapsulates D*.
virtual void Clear(const Option_t *="")
Bool_t ExtractDstarAttributes(const AliAODRecoCascadeHF *DstarCand, AliDmesonJetInfo &DmesonJet, UInt_t i)
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
const AliJetInfo * GetJet(std::string n) const
Short_t fFirstPartonType
! type of the first parton in the shower tree (only for particle level D mesons)
const AliMCParticleIterableMomentumContainer all_momentum() const
virtual ~AliAnalysisTaskDmesonJets()
This is the standard destructor.
Byte_t fSelectionType
! for D0: 0=not selected, 1=D0, 2=D0bar, 3=both
static std::pair< AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle * > CheckOrigin(const AliAODMCParticle *part, TClonesArray *mcArray, Bool_t firstParton=kFALSE)
Bool_t ExtractD0Attributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
std::vector< AliJetInfo > fJets
! Inclusive jets reconstructed in the current event (includes D meson candidate daughters, if any)
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:99
AliDmesonJetInfo & operator=(const AliDmesonJetInfo &source)
Double_t fMaxChargedPt
Transverse momentum of the leading charged particle (or track)
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
static FJJetAlgo ConvertToFJAlgo(EJetAlgo_t algo)
decay
Definition: HFPtSpectrum.C:41
std::map< std::string, AliJetInfo > fJets
! list of jets
Definition: External.C:220
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Double_t minMass
void RunParticleLevelAnalysis()
Run a particle level analysis.
Bool_t fApplyKinematicCuts
Apply jet kinematic cuts.
const char * GetName() const
Generate a name for this jet definition.
unsigned short UShort_t
Definition: External.C:28
Bool_t ExtractRecoDecayAttributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
Double_t InvMassD0() const
EOutputType_t fOutputType
Output type: none, TTree or THnSparse.
virtual void AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
bool operator==(const AliAnalysisTaskDmesonJets::AliHFJetDefinition &lhs, const AliAnalysisTaskDmesonJets::AliHFJetDefinition &rhs)
Double_t InvMassDstarKpipi() const
AnalysisEngine * AddAnalysisEngine(ECandidateType_t type, EMCMode_t bkgMode, EJetType_t jettype, Double_t jetradius, TString cutfname="")
const Int_t nbins
Double_t maxMass
bool Bool_t
Definition: External.C:53
Int_t fNOutputTrees
Maximum number of output trees.
AliHFJetDefinition()
This is the default constructor, used for ROOT I/O purposes.
AliHFJetDefinition & operator=(const AliHFJetDefinition &source)
Int_t GetNAcceptedParticles() const
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:97
AliTLorentzVector fMomentum
4-momentum of the jet
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
AliHFJetDefinition * AddJetDefinition(EJetType_t type, Double_t r, EJetAlgo_t algo, ERecoScheme_t reco)
Container for MC-true particles within the EMCAL framework.
Lightweight class that encapsulates D0.
Bool_t IsAnyJetInAcceptance(const AliDmesonJetInfo &dMesonJet) const
void AdoptParticleContainer(AliParticleContainer *cont)
void SetCharge(EChargeCut_t c)
Container for jet within the EMCAL jet framework.
Definition: External.C:196
AliAODMCParticle * fLastParton
! pointer to the last parton in the shower tree of the D meson (only for particle level D mesons) ...
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)