AliPhysics  ff1d528 (ff1d528)
 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 #include "AliAnalysisVertexingHF.h"
42 
43 // Aliroot EMCal jet framework
44 #include "AliEmcalJetTask.h"
45 #include "AliEmcalJet.h"
46 #include "AliJetContainer.h"
47 #include "AliParticleContainer.h"
48 #include "AliEmcalParticle.h"
49 #include "AliFJWrapper.h"
50 #include "AliRhoParameter.h"
51 
53 
54 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfo
55 
59 
67 {
68  dphi = TVector2::Phi_mpi_pi(fMomentum.Phi() - jet.Phi());;
69  deta = fMomentum.Eta() - jet.Eta();
70  return TMath::Sqrt(dphi*dphi + deta*deta);
71 }
72 
78 {
79  Double_t deta = 0;
80  Double_t dphi = 0;
81  return GetDistance(jet, deta, dphi);
82 }
83 
84 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonJetInfo
85 
89 
92  fDmesonParticle(0),
93  fD(),
94  fSoftPionPt(0),
95  fInvMass2Prong(0),
96  fJets(),
97  fMCLabel(-1),
98  fReconstructed(kFALSE),
99  fFirstParton(0),
100  fFirstPartonType(0),
101  fLastParton(0),
102  fLastPartonType(0),
103  fSelectionType(0)
104 {
105 }
106 
111  fDmesonParticle(source.fDmesonParticle),
112  fD(source.fD),
113  fSoftPionPt(source.fSoftPionPt),
114  fInvMass2Prong(source.fInvMass2Prong),
115  fJets(source.fJets),
116  fMCLabel(source.fMCLabel),
117  fReconstructed(source.fReconstructed),
118  fFirstParton(source.fFirstParton),
119  fFirstPartonType(source.fFirstPartonType),
120  fLastParton(source.fLastParton),
121  fLastPartonType(source.fLastPartonType),
122  fSelectionType(source.fSelectionType)
123 {
124 }
125 
130 {
131  new (this) AliDmesonJetInfo(source);
132  return *this;
133 }
134 
137 {
138  fD.SetPtEtaPhiE(0,0,0,0);
139  fSoftPionPt = 0;
140  fInvMass2Prong = 0;
141  fDmesonParticle = 0;
142  fMCLabel = -1;
143  fReconstructed = kFALSE;
144  fFirstParton = 0;
145  fFirstPartonType = 0;
146  fLastParton = 0;
147  fLastPartonType = 0;
148  for (auto &jet : fJets) {
149  jet.second.fMomentum.SetPtEtaPhiE(0,0,0,0);
150  jet.second.fNConstituents = 0;
151  jet.second.fNEF = 0;
152  jet.second.fMaxChargedPt = 0;
153  jet.second.fMaxNeutralPt = 0;
154  }
155 }
156 
159 {
160  Printf("Printing D Meson Jet object.");
161  Printf("D Meson: pT = %.3f, eta = %.3f, phi = %.3f, inv. mass = %.3f", fD.Pt(), fD.Eta(), fD.Phi_0_2pi(), fD.M());
162  Printf("Soft pion pT: %.3f. 2-Prong Invariant mass = %.3f", fSoftPionPt, fInvMass2Prong);
163  for (auto &jet : fJets) {
164  Printf("Jet %s: pT = %.3f, eta = %.3f, phi = %.3f", jet.first.c_str(), jet.second.Pt(), jet.second.Eta(), jet.second.Phi_0_2pi());
165  Printf("Jet N Consituents = %d", jet.second.fNConstituents);
166  }
167 }
168 
173 {
174  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
175  if (it == fJets.end()) return 0;
176 
177  Double_t z = 0;
178 
179  if ((*it).second.Pt() > 0) {
180  TVector3 dvect = fD.Vect();
181  TVector3 jvect = (*it).second.fMomentum.Vect();
182 
183  Double_t jetMom = jvect * jvect;
184 
185  if (jetMom < 1e-6) {
186  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetZ", "Zero jet momentum!");
187  z = 0.999;
188  }
189  else {
190  z = (dvect * jvect) / jetMom;
191  }
192 
193  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
194  }
195 
196  return z;
197 }
198 
206 {
207  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
208  if (it == fJets.end()) return 0;
209 
210  return GetDistance((*it).second, deta, dphi);
211 }
212 
218 {
219  Double_t deta = 0;
220  Double_t dphi = 0;
221  return GetDistance(n, deta, dphi);
222 }
223 
231 {
232  dphi = TVector2::Phi_mpi_pi(fD.Phi() - jet.Phi());;
233  deta = fD.Eta() - jet.Eta();
234  return TMath::Sqrt(dphi*dphi + deta*deta);
235 }
236 
242 {
243  Double_t deta = 0;
244  Double_t dphi = 0;
245  return GetDistance(jet, deta, dphi);
246 }
247 
253 {
254  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
255  if (it == fJets.end()) {
256  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
257  n.c_str());
258  return 0;
259  }
260  return &((*it).second);
261 }
262 
268 {
269  std::map<std::string, AliJetInfo>::iterator it = fJets.find(n);
270  if (it == fJets.end()) {
271  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
272  n.c_str());
273  return 0;
274  }
275  return &((*it).second);
276 }
277 
278 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoSummary
279 
283 
289  fPt(0),
290  fEta(0),
291  fPhi(0),
292  fR(0),
293  fZ(0),
294  fN(0)
295 {
296  Set(source, n);
297 }
298 
301 {
302  fPt = 0;
303  fEta = 0;
304  fPhi = 0;
305  fR = 0;
306  fZ = 0;
307  fN = 0;
308 }
309 
315 {
316  std::map<std::string, AliJetInfo>::const_iterator it = source.fJets.find(n);
317  if (it == source.fJets.end()) return;
318 
319  Set((*it).second);
320 
321  fR = source.GetDistance(n);
322  fZ = source.GetZ(n);
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::AliJetInfoPbPbSummary
339 
343 
350  fCorrPt(0),
351  fArea(0)
352 {
353  Set(source, n);
354 }
355 
358 {
360  fCorrPt = 0;
361  fArea = 0;
362 }
363 
368 {
369  AliJetInfoSummary::Set(source);
370  fArea = source.fArea;
371  fCorrPt = source.fCorrPt;
372 }
373 
374 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonInfoSummary
375 
379 
384  fPt(0),
385  fEta(0),
386  fPhi(0)
387 {
388  Set(source);
389 }
390 
395 {
396  fPt = source.fD.Pt();
397  fEta = source.fD.Eta();
398  fPhi = source.fD.Phi_0_2pi();
399 }
400 
403 {
404  fPt = 0;
405  fEta = 0;
406  fPhi = 0;
407 }
408 
409 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary
410 
414 
419  AliDmesonInfoSummary(source),
420  fFirstPartonType(0),
421  fFirstPartonPt(0),
422  fLastPartonType(0),
423  fLastPartonPt(0)
424 {
425  Set(source);
426 }
427 
432 {
434 
435  fFirstPartonType = source.fFirstPartonType;
436  fLastPartonType = source.fLastPartonType;
437 
438  if (source.fFirstParton) {
439  fFirstPartonPt = source.fFirstParton->Pt();
440  }
441  else {
442  fFirstPartonPt = 0.;
443  }
444 
445  if (source.fLastParton) {
446  fLastPartonPt = source.fLastParton->Pt();
447  }
448  else {
449  fLastPartonPt = 0.;
450  }
451 }
452 
455 {
457  fFirstPartonType = 0,
458  fFirstPartonPt = 0.;
459  fLastPartonType = 0,
460  fLastPartonPt = 0.;
461 }
462 
463 // Definitions of class AliAnalysisTaskDmesonJets::AliD0InfoSummary
464 
468 
473  AliDmesonInfoSummary(source),
474  fInvMass(source.fD.M()),
475  fSelectionType(0)
476 {
477 }
478 
483 {
484  fInvMass = source.fD.M();
485  fSelectionType = source.fSelectionType;
487 }
488 
491 {
493  fSelectionType = 0;
494  fInvMass = 0;
495 }
496 
497 // Definitions of class AliAnalysisTaskDmesonJets::AliDStarInfoSummary
498 
502 
507  AliDmesonInfoSummary(source),
508  f2ProngInvMass(source.fInvMass2Prong),
509  fDeltaInvMass(source.fD.M() - source.fInvMass2Prong)
510 {
511 }
512 
517 {
518  f2ProngInvMass = source.fInvMass2Prong;
519  fDeltaInvMass = source.fD.M() - source.fInvMass2Prong;
521 }
522 
525 {
527 
528  f2ProngInvMass = 0;
529  fDeltaInvMass = 0;
530 }
531 
532 // Definitions of class AliAnalysisTaskDmesonJets::AliJetDefinition
533 
537 
540  TObject(),
541  fJetType(AliJetContainer::kChargedJet),
542  fRadius(0),
543  fJetAlgo(AliJetContainer::antikt_algorithm),
544  fRecoScheme(AliJetContainer::pt_scheme),
545  fMinJetPt(0.),
546  fMaxJetPt(500.),
547  fMinJetPhi(0.),
548  fMaxJetPhi(0.),
549  fMinJetEta(-1.),
550  fMaxJetEta(1.),
551  fMinChargedPt(0.),
552  fMaxChargedPt(0.),
553  fMinNeutralPt(0.),
554  fMaxNeutralPt(0.),
555  fRhoName(),
556  fRho(0)
557 {
558 }
559 
567  TObject(),
568  fJetType(type),
569  fRadius(r),
570  fJetAlgo(algo),
571  fRecoScheme(reco),
572  fMinJetPt(0.),
573  fMaxJetPt(500.),
574  fMinJetPhi(0.),
575  fMaxJetPhi(0.),
576  fMinJetEta(-1.),
577  fMaxJetEta(1.),
578  fMinChargedPt(0.),
579  fMaxChargedPt(0.),
580  fMinNeutralPt(0.),
581  fMaxNeutralPt(0.),
582  fRhoName(),
583  fRho(0)
584 {
585 }
586 
594  TObject(),
595  fJetType(type),
596  fRadius(r),
597  fJetAlgo(algo),
598  fRecoScheme(reco),
599  fMinJetPt(0.),
600  fMaxJetPt(500.),
601  fMinJetPhi(0.),
602  fMaxJetPhi(0.),
603  fMinJetEta(-1.),
604  fMaxJetEta(1.),
605  fMinChargedPt(0.),
606  fMaxChargedPt(0.),
607  fMinNeutralPt(0.),
608  fMaxNeutralPt(0.),
609  fRhoName(rhoName),
610  fRho(0)
611 {
612 }
613 
618  TObject(),
619  fJetType(source.fJetType),
620  fRadius(source.fRadius),
621  fJetAlgo(source.fJetAlgo),
622  fRecoScheme(source.fRecoScheme),
623  fMinJetPt(source.fMinJetPt),
624  fMaxJetPt(source.fMaxJetPt),
625  fMinJetPhi(source.fMinJetPhi),
626  fMaxJetPhi(source.fMaxJetPhi),
627  fMinJetEta(source.fMinJetEta),
628  fMaxJetEta(source.fMaxJetEta),
629  fMinChargedPt(source.fMinChargedPt),
630  fMaxChargedPt(source.fMaxChargedPt),
631  fMinNeutralPt(source.fMinNeutralPt),
632  fMaxNeutralPt(source.fMaxNeutralPt),
633  fRhoName(source.fRhoName),
634  fRho(0)
635 {
636 }
637 
642 {
643  new (this) AliHFJetDefinition(source);
644  return *this;
645 }
646 
649 {
650  static TString name;
651 
652  name = AliJetContainer::GenerateJetName(fJetType, fJetAlgo, fRecoScheme, fRadius, 0, 0, "Jet");
653 
654  return name.Data();
655 }
656 
662 {
663  if (fMinJetEta < fMaxJetEta && (jet.Eta() < fMinJetEta || jet.Eta() > fMaxJetEta)) return kFALSE;
664  if (fMinJetPhi < fMaxJetPhi && (jet.Phi() < fMinJetPhi || jet.Phi() > fMaxJetPhi)) return kFALSE;
665  if (jet.Pt() > fMaxJetPt || jet.Pt() < fMinJetPt) return kFALSE;
666  if (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt) return kFALSE;
667  if (jet.fMaxNeutralPt < fMinNeutralPt || jet.fMaxNeutralPt > fMaxNeutralPt) return kFALSE;
668 
669  return kTRUE;
670 }
671 
677 {
678  const AliJetInfo* jet = dMesonJet.GetJet(n);
679  if (!jet) return kFALSE;
680  return IsJetInAcceptance((*jet));
681 }
682 
689 {
690  if (lhs.fJetType > rhs.fJetType) return false;
691  else if (lhs.fJetType < rhs.fJetType) return true;
692  else {
693  if (lhs.fRadius > rhs.fRadius) return false;
694  else if (lhs.fRadius < rhs.fRadius) return true;
695  else {
696  if (lhs.fJetAlgo > rhs.fJetAlgo) return false;
697  else if (lhs.fJetAlgo < rhs.fJetAlgo) return true;
698  else {
699  if (lhs.fRecoScheme < rhs.fRecoScheme) return true;
700  else return false;
701  }
702  }
703  }
704 }
705 
712 {
713  if (lhs.fJetType != rhs.fJetType) return false;
714  if (lhs.fRadius != rhs.fRadius) return false;
715  if (lhs.fJetAlgo != rhs.fJetAlgo) return false;
716  if (lhs.fRecoScheme != rhs.fRecoScheme) return false;
717  return true;
718 }
719 
720 // Definitions of class AliAnalysisTaskDmesonJets::AnalysisEngine
721 
725 
728  TObject(),
729  fFirstPartons(),
730  fLastPartons(),
731  fCandidateType(kD0toKpi),
732  fCandidateName(),
733  fCandidatePDG(0),
734  fNDaughters(0),
735  fPDGdaughters(),
736  fBranchName(),
737  fMCMode(kNoMC),
738  fNMassBins(0),
739  fMinMass(0),
740  fMaxMass(0),
741  fRDHFCuts(0),
742  fRejectedOrigin(0),
743  fAcceptedDecay(0),
744  fInhibit(kFALSE),
745  fJetDefinitions(),
746  fPtBinWidth(0.5),
747  fMaxPt(100),
748  fRandomGen(0),
749  fTrackEfficiency(0),
750  fDataSlotNumber(-1),
751  fTree(0),
752  fCurrentDmesonJetInfo(0),
753  fCurrentJetInfo(0),
754  fCandidateArray(0),
755  fMCContainer(),
756  fTrackContainers(),
757  fClusterContainers(),
758  fAodEvent(0),
759  fFastJetWrapper(0),
760  fHistManager(0)
761 {
762 }
763 
772  TObject(),
773  fFirstPartons(),
774  fLastPartons(),
775  fCandidateType(type),
776  fCandidateName(),
777  fCandidatePDG(0),
778  fNDaughters(0),
779  fPDGdaughters(),
780  fBranchName(),
781  fMCMode(MCmode),
782  fNMassBins(nMassBins),
783  fMinMass(0),
784  fMaxMass(0),
785  fRDHFCuts(cuts),
786  fRejectedOrigin(0),
787  fAcceptedDecay(kAnyDecay),
788  fInhibit(kFALSE),
789  fJetDefinitions(),
790  fPtBinWidth(0.5),
791  fMaxPt(100),
792  fRandomGen(0),
793  fTrackEfficiency(0),
794  fDataSlotNumber(-1),
795  fTree(0),
796  fCurrentDmesonJetInfo(0),
797  fCurrentJetInfo(0),
798  fCandidateArray(0),
799  fMCContainer(),
800  fTrackContainers(),
801  fClusterContainers(),
802  fAodEvent(0),
803  fFastJetWrapper(0),
804  fHistManager(0)
805 {
806  SetCandidateProperties(range);
807 }
808 
813  TObject(source),
814  fFirstPartons(source.fFirstPartons),
815  fLastPartons(source.fLastPartons),
816  fCandidateType(source.fCandidateType),
817  fCandidateName(source.fCandidateName),
818  fCandidatePDG(source.fCandidatePDG),
819  fNDaughters(source.fNDaughters),
820  fPDGdaughters(source.fPDGdaughters),
821  fBranchName(source.fBranchName),
822  fMCMode(source.fMCMode),
823  fNMassBins(source.fNMassBins),
824  fMinMass(source.fMinMass),
825  fMaxMass(source.fMaxMass),
826  fRDHFCuts(),
827  fRejectedOrigin(source.fRejectedOrigin),
828  fAcceptedDecay(source.fAcceptedDecay),
829  fInhibit(source.fInhibit),
830  fJetDefinitions(source.fJetDefinitions),
831  fPtBinWidth(source.fPtBinWidth),
832  fMaxPt(source.fMaxPt),
833  fRandomGen(source.fRandomGen),
835  fDataSlotNumber(-1),
836  fTree(0),
837  fCurrentDmesonJetInfo(0),
838  fCurrentJetInfo(0),
839  fCandidateArray(source.fCandidateArray),
840  fMCContainer(source.fMCContainer),
841  fTrackContainers(source.fTrackContainers),
842  fClusterContainers(source.fClusterContainers),
843  fAodEvent(source.fAodEvent),
845  fHistManager(source.fHistManager)
846 {
847  SetRDHFCuts(source.fRDHFCuts);
848 }
849 
850 // Destructor
852 {
853  delete fRDHFCuts;
854 }
855 
860 {
861  new (this) AnalysisEngine(source);
862  return *this;
863 }
864 
869 {
870  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
871  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
872  }
873 
874  return kFALSE;
875 }
876 
878 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const /*geom*/, Int_t /*runNumber*/)
879 {
880 }
881 
886 {
887  switch (fCandidateType) {
888  case kD0toKpi:
889  fCandidatePDG = 421;
890  fCandidateName = "D0";
891  fNDaughters = 2;
892  fPDGdaughters.Set(fNDaughters);
893  fPDGdaughters.Reset();
894  fPDGdaughters[0] = 211; // pi
895  fPDGdaughters[1] = 321; // K
896  fBranchName = "D0toKpi";
897  fAcceptedDecay = kDecayD0toKpi;
898  break;
899  case kD0toKpiLikeSign:
900  fCandidatePDG = 421;
901  fCandidateName = "2ProngLikeSign";
902  fNDaughters = 2;
903  fPDGdaughters.Set(fNDaughters);
904  fPDGdaughters.Reset();
905  fPDGdaughters[0] = 211; // pi
906  fPDGdaughters[1] = 321; // K
907  fBranchName = "LikeSign2Prong";
908  break;
909  case kDstartoKpipi:
910  fCandidatePDG = 413;
911  fCandidateName = "DStar";
912  fNDaughters = 3;
913  fPDGdaughters.Set(fNDaughters);
914  fPDGdaughters.Reset();
915  fPDGdaughters[0] = 211; // pi soft
916  fPDGdaughters[1] = 211; // pi fromD0
917  fPDGdaughters[2] = 321; // K from D0
918  fBranchName = "Dstar";
919  fAcceptedDecay = kDecayDStartoKpipi;
920  break;
921  default:
922  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
923  }
924 
925  CalculateMassLimits(range, fCandidatePDG, fNMassBins, fMinMass, fMaxMass);
926 }
927 
932 {
933  if (fRDHFCuts) delete fRDHFCuts;
934  fRDHFCuts = cuts;
935 }
936 
941 {
942  if (!cuts) return;
943  if (fRDHFCuts) delete fRDHFCuts;
944  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
945 }
946 
951 {
952  static TString name;
953 
954  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
955 
956  return name.Data();
957 }
958 
963 {
964  static TString name;
965 
966  name = fCandidateName;
967  switch (fMCMode) {
968  case kBackgroundOnly:
969  name += "_kBackgroundOnly";
970  break;
971  case kSignalOnly:
972  name += "_kSignalOnly";
973  break;
974  case kMCTruth:
975  name += "_MCTruth";
976  break;
977  case kWrongPID:
978  name += "_WrongPID";
979  break;
980  default:
981  break;
982  }
983 
984  return name.Data();
985 }
986 
994 {
995  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
996 
997  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
998  fJetDefinitions.push_back(def);
999  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
1000  "Total number of jet definitions is now %lu.",
1001  def.GetName(), GetName(), fJetDefinitions.size());
1002  // For detector level set maximum track pt to 100 GeV/c
1003  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
1004  }
1005  else {
1006  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
1007  }
1008 
1009  return &(*it);
1010 }
1011 
1023 {
1024  AliHFJetDefinition def(type, r, algo, reco);
1025 
1026  return AddJetDefinition(def);
1027 }
1028 
1034 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
1035 {
1036  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
1037  while (it != fJetDefinitions.end() && (*it) != def) it++;
1038  return it;
1039 }
1040 
1047 {
1048  if (lhs.fCandidateType > rhs.fCandidateType) return false;
1049  else if (lhs.fCandidateType < rhs.fCandidateType) return true;
1050  else {
1051  if (lhs.fMCMode < rhs.fMCMode) return true;
1052  else return false;
1053  }
1054 }
1055 
1062 {
1063  if (lhs.fCandidateType != rhs.fCandidateType) return false;
1064  if (lhs.fMCMode != rhs.fMCMode) return false;
1065  return true;
1066 }
1067 
1076 {
1077  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) { // D0 candidate
1078  return ExtractD0Attributes(Dcand, DmesonJet, i);
1079  }
1080  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
1081  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
1082  }
1083  else {
1084  return kFALSE;
1085  }
1086 }
1087 
1096 {
1097  AliDebug(10,"Checking if D0 meson is selected");
1098  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
1099  if (isSelected == 0) return kFALSE;
1100 
1101  Int_t MCtruthPdgCode = 0;
1102 
1103  Double_t invMassD = 0;
1104 
1105  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
1106  // Checks also the origin, and if it matches the rejected origin mask, return false
1107  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly || fMCMode == kWrongPID) {
1108  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
1109  DmesonJet.fMCLabel = mcLab;
1110 
1111  // Retrieve the generated particle (if exists) and its PDG code
1112  if (mcLab >= 0) {
1113  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1114 
1115  if (aodMcPart) {
1116  // Check origin and return false if it matches the rejected origin mask
1117  if (fRejectedOrigin) {
1118  auto origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
1119  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1120  }
1121  MCtruthPdgCode = aodMcPart->PdgCode();
1122  }
1123  }
1124  }
1125 
1126  if (isSelected == 1) { // selected as a D0
1127  if (i != 0) return kFALSE; // only one mass hypothesis thanks to PID
1128 
1129  if (fMCMode == kNoMC ||
1130  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1131  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly) ||
1132  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kWrongPID)) {
1133  // 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)
1134  AliDebug(10,"Selected as D0");
1135  invMassD = Dcand->InvMassD0();
1136  }
1137  else { // conditions above not passed, so return FALSE
1138  return kFALSE;
1139  }
1140  }
1141  else if (isSelected == 2) { // selected as a D0bar
1142  if (i != 1) return kFALSE; // only one mass hypothesis thanks to PID
1143 
1144  if (fMCMode == kNoMC ||
1145  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1146  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly) ||
1147  (MCtruthPdgCode == fCandidatePDG && fMCMode == kWrongPID)) {
1148  // 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)
1149  AliDebug(10,"Selected as D0bar");
1150  invMassD = Dcand->InvMassD0bar();
1151  }
1152  else { // conditions above not passed, so return FALSE
1153  return kFALSE;
1154  }
1155  }
1156  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
1157  AliDebug(10,"Selected as either D0 or D0bar");
1158 
1159  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
1160  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1161  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kWrongPID))) {
1162  if (i != 0) return kFALSE;
1163  AliDebug(10, "MC truth is D0");
1164  invMassD = Dcand->InvMassD0();
1165  }
1166  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1167  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kWrongPID))) {
1168  if (i != 1) return kFALSE;
1169  AliDebug(10, "MC truth is D0bar");
1170  invMassD = Dcand->InvMassD0bar();
1171  }
1172  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
1173 
1174  // Only accept it if background-only OR background-and-signal was requested
1175  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
1176  // Select D0 or D0bar depending on the i-parameter
1177  if (i == 0) {
1178  AliDebug(10, "Returning invariant mass with D0 hypothesis");
1179  invMassD = Dcand->InvMassD0();
1180  }
1181  else if (i == 1) {
1182  AliDebug(10, "Returning invariant mass with D0bar hypothesis");
1183  invMassD = Dcand->InvMassD0bar();
1184  }
1185  else { // i > 1
1186  return kFALSE;
1187  }
1188  }
1189  else { // signal-only was requested but this is not a true D0
1190  return kFALSE;
1191  }
1192  }
1193  }
1194  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
1195  return kTRUE;
1196 }
1197 
1206 {
1207  AliDebug(10,"Checking if D* meson is selected");
1208  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoCascadeHF*>(DstarCand), AliRDHFCuts::kAll, fAodEvent);
1209  if (isSelected == 0) return kFALSE;
1210 
1211  if ((i == 1 && DstarCand->Charge()>0) || (i == 0 && DstarCand->Charge()<0) || i > 1) return kFALSE; // only one mass hypothesis for the D*
1212 
1213  Int_t MCtruthPdgCode = 0;
1214 
1215  Double_t invMassD = 0;
1216 
1217  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1218  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1219  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1220 
1221  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1222  AliDebug(10, Form("MC label is %d", mcLab));
1223  DmesonJet.fMCLabel = mcLab;
1224  if (mcLab >= 0) {
1225  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1226 
1227  if (aodMcPart) {
1228  if (fRejectedOrigin) {
1229  auto origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
1230  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1231  }
1232 
1233  MCtruthPdgCode = aodMcPart->PdgCode();
1234  AliDebug(10, Form("MC truth pdg code is %d",MCtruthPdgCode));
1235  }
1236  }
1237  }
1238 
1239  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1240  if (fMCMode == kNoMC ||
1241  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1242  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1243  // 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)
1244  invMassD = DstarCand->InvMassDstarKpipi();
1245  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1246  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1247  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1248  return kTRUE;
1249  }
1250  else { // conditions above not passed, so return FALSE
1251  return kFALSE;
1252  }
1253 }
1254 
1262 {
1263  if (!part) return kUnknownDecay;
1264  if (!mcArray) return kUnknownDecay;
1265 
1267 
1268  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1269 
1270  if (part->GetNDaughters() == 2) {
1271 
1272  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1273  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1274 
1275  if (!d1 || !d2) {
1276  return decay;
1277  }
1278 
1279  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1280  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1281 
1282  if (absPdgPart == 421) { // D0 -> K pi
1283  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1284  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1285  decay = kDecayD0toKpi;
1286  }
1287  }
1288 
1289  if (absPdgPart == 413) { // D* -> D0 pi
1290  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1291  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1292  if (D0decay == kDecayD0toKpi) {
1293  decay = kDecayDStartoKpipi;
1294  }
1295  }
1296  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1297  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1298  if (D0decay == kDecayD0toKpi) {
1299  decay = kDecayDStartoKpipi;
1300  }
1301  }
1302  }
1303  }
1304 
1305  return decay;
1306 }
1307 
1314 std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> AliAnalysisTaskDmesonJets::AnalysisEngine::CheckOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, Bool_t firstParton)
1315 {
1316  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1317 
1318  std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> result(kUnknownQuark, 0);
1319 
1320  if (!part) return result;
1321  if (!mcArray) return result;
1322 
1323  Int_t mother = part->GetMother();
1324  while (mother >= 0) {
1325  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1326  if (mcGranma) {
1327  Int_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
1328 
1329  if (abspdgGranma == 1) result = {kFromDown, mcGranma};
1330  if (abspdgGranma == 2) result = {kFromUp, mcGranma};
1331  if (abspdgGranma == 3) result = {kFromStrange, mcGranma};
1332  if (abspdgGranma == 4) result = {kFromCharm, mcGranma};
1333  if (abspdgGranma == 5) result = {kFromBottom, mcGranma};
1334  if (abspdgGranma == 6) result = {kFromTop, mcGranma};
1335  if (abspdgGranma == 9 || abspdgGranma == 21) result = {kFromGluon, mcGranma};
1336 
1337  // If looking for the very first parton in the hard scattering, it will continue the loop until it cannot find a mother particle
1338  if (result.first != kUnknownQuark && !firstParton) return result;
1339 
1340  mother = mcGranma->GetMother();
1341  }
1342  else {
1343  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1344  break;
1345  }
1346  }
1347 
1348  return result;
1349 }
1350 
1353 {
1354  for (auto& jetDef : fJetDefinitions) {
1355  jetDef.fJets.clear();
1356  }
1357 
1358  if (fMCMode == kMCTruth) {
1359  RunParticleLevelAnalysis();
1360  }
1361  else {
1362  RunDetectorLevelAnalysis();
1363  }
1364 }
1365 
1368 {
1369  // Fill the vertex info of the candidates
1370  // Needed for reduced delta AOD, where the vertex info has been deleted
1371  // to reduce the delta AOD file size
1373 
1374  const Int_t nD = fCandidateArray->GetEntriesFast();
1375 
1376  AliDmesonJetInfo DmesonJet;
1377 
1378  Int_t nAccCharm[3] = {0};
1379  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1380  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1381  if (!charmCand) continue;
1382  if(!(vHF.FillRecoCand(fAodEvent,charmCand))) continue;
1383 
1384  // region of interest + cuts
1385  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1386  Int_t nMassHypo = 0; // number of mass hypothesis accepted for this D meson
1387  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1388  DmesonJet.Reset();
1389  DmesonJet.fDmesonParticle = charmCand;
1390  DmesonJet.fSelectionType = im + 1;
1391  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1392  for (auto& def : fJetDefinitions) {
1393  if (!FindJet(charmCand, DmesonJet, def)) {
1394  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1395  def.GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1396  }
1397  }
1398  fDmesonJets[(icharm+1)*(1-(im*2))] = DmesonJet;
1399  nMassHypo++;
1400  nAccCharm[im]++;
1401  }
1402  }
1403  if (nMassHypo == 2) {
1404  nAccCharm[0]--;
1405  nAccCharm[1]--;
1406  nAccCharm[2] += 2;
1407  }
1408  if (nMassHypo == 2) { // both mass hypothesis accepted
1409  fDmesonJets[(icharm+1)].fSelectionType = 3;
1410  fDmesonJets[-(icharm+1)].fSelectionType = 3;
1411  }
1412  } // end of D cand loop
1413 
1414  TString hname;
1415 
1416  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1417  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1418  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1419  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1420 
1421  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1422  Int_t ntracks = 0;
1423  for (auto track_cont : fTrackContainers) ntracks += track_cont->GetNAcceptedTracks();
1424  fHistManager->FillTH2(hname, ntracks, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1425 
1426  hname = TString::Format("%s/fHistNDmesons", GetName());
1427  fHistManager->FillTH1(hname, nD);
1428 }
1429 
1440 {
1441  TString hname;
1442 
1443  Double_t rho = 0;
1444  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
1445 
1447  fFastJetWrapper->SetR(jetDef.fRadius);
1450 
1451  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1452 
1453  if (jetDef.fJetType != AliJetContainer::kNeutralJet) {
1454  for (auto track_cont : fTrackContainers) {
1455  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
1456  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(Dcand);
1457  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1458  AddInputVectors(track_cont, 100, static_cast<TH2*>(fHistManager->FindObject(hname)), fTrackEfficiency);
1459 
1460  if (hftrack_cont) {
1461  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1462  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1463  const TObjArray& daughters = hftrack_cont->GetDaughterList();
1464  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1465  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1466  if (!hftrack_cont->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1467  }
1468  }
1469  }
1470  }
1471 
1472  if (jetDef.fJetType != AliJetContainer::kChargedJet) {
1473  for (auto clus_cont : fClusterContainers) {
1474  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1475  AddInputVectors(clus_cont, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1476  }
1477  }
1478 
1479  // run jet finder
1480  fFastJetWrapper->Run();
1481 
1482  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1483 
1484  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1485  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1486 
1487  Bool_t isDmesonJet = kFALSE;
1488 
1489  Double_t maxChPt = 0;
1490  Double_t maxNePt = 0;
1491  Double_t totalNeutralPt = 0;
1492 
1493  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1494  if (constituents[ic].user_index() == 0) {
1495  isDmesonJet = kTRUE;
1496  }
1497  else if (constituents[ic].user_index() >= 100) {
1498  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1499  }
1500  else if (constituents[ic].user_index() <= -100) {
1501  totalNeutralPt += constituents[ic].pt();
1502  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1503  }
1504  }
1505 
1506  if (isDmesonJet) {
1507  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1508  DmesonJet.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1509  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1510  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1511  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1512  DmesonJet.fJets[jetDef.GetName()].fArea = jets_incl[ijet].area();
1513  DmesonJet.fJets[jetDef.GetName()].fCorrPt = DmesonJet.fJets[jetDef.GetName()].fMomentum.Pt() - jets_incl[ijet].area() * rho;
1514 
1515 
1516  return kTRUE;
1517  }
1518  }
1519 
1520  return kFALSE;
1521 }
1522 
1526 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist, Double_t eff)
1527 {
1528  auto itcont = cont->all_momentum();
1529  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1530  UInt_t rejectionReason = 0;
1531  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1532  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1533  continue;
1534  }
1535  if (fRandomGen && eff > 0 && eff < 1) {
1536  Double_t rnd = fRandomGen->Rndm();
1537  if (eff < rnd) {
1538  if (rejectHist) rejectHist->Fill(6, it->first.Pt());
1539  continue;
1540  }
1541  }
1542  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1543  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1544  }
1545 }
1546 
1549 {
1550  TString hname;
1551 
1552  fMCContainer->SetSpecialPDG(fCandidatePDG);
1553  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1554  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1555 
1556  if (!fMCContainer->IsSpecialPDGFound()) return;
1557 
1558  Int_t nAccCharm[3] = {0};
1559 
1560  for (auto &jetDef : fJetDefinitions) {
1561  Double_t rho = 0;
1562  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
1563  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", GetName(), jetDef.GetName());
1564  TH1* histNDmesonsVsNconstituents = static_cast<TH1*>(fHistManager->FindObject(hname));
1565 
1566  switch (jetDef.fJetType) {
1568  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
1569  break;
1571  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
1572  break;
1574  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
1575  break;
1576  }
1577 
1579  fFastJetWrapper->SetR(jetDef.fRadius);
1582 
1583  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1584  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1585 
1586  fFastJetWrapper->Run();
1587 
1588  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1589 
1590  for (auto jet : jets_incl) {
1591  Int_t nDmesonsInJet = 0;
1592 
1593  for (auto constituent : jet.constituents()) {
1594  Int_t iPart = constituent.user_index() - 100;
1595  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
1596  if (!part) {
1597  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1598  continue;
1599  }
1600  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1601  nDmesonsInJet++;
1602  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1603  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1604  std::pair<int, AliDmesonJetInfo> element;
1605  element.first = iPart;
1606  dMesonJetIt = fDmesonJets.insert(element).first;
1607  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1608  (*dMesonJetIt).second.fDmesonParticle = part;
1609  (*dMesonJetIt).second.fSelectionType = part->PdgCode() > 0 ? 1 : 2;
1610 
1611  UShort_t p = 0;
1612  UInt_t rs = 0;
1613 
1614  auto firstParton = CheckOrigin(part, fMCContainer->GetArray(), kTRUE);
1615  p = 0;
1616  rs = firstParton.first;
1617  while (rs >>= 1) { p++; }
1618  (*dMesonJetIt).second.fFirstPartonType = p;
1619  (*dMesonJetIt).second.fFirstParton = firstParton.second;
1620 
1621  auto lastParton = CheckOrigin(part, fMCContainer->GetArray(), kFALSE);
1622  p = 0;
1623  rs = lastParton.first;
1624  while (rs >>= 1) { p++; }
1625  (*dMesonJetIt).second.fLastPartonType = p;
1626  (*dMesonJetIt).second.fLastParton = lastParton.second;
1627 
1628  if (part->PdgCode() > 0) {
1629  nAccCharm[0]++;
1630  }
1631  else {
1632  nAccCharm[1]++;
1633  }
1634  }
1635 
1636  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
1637  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = jet.constituents().size();
1638  (*dMesonJetIt).second.fJets[jetDef.GetName()].fArea = jet.area();
1639  (*dMesonJetIt).second.fJets[jetDef.GetName()].fCorrPt = (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.Pt() - jet.area() * rho;
1640  } // if constituent is a D meson
1641  } // for each constituent
1642  if (nDmesonsInJet > 0) histNDmesonsVsNconstituents->Fill(jet.constituents().size(), nDmesonsInJet);
1643  } // for each jet
1644  } // for each jet definition
1645 
1646  if (fDmesonJets.size() != nAccCharm[0]+nAccCharm[1]) AliError(Form("I found %lu mesons (%d)?", fDmesonJets.size(), nAccCharm[0]+nAccCharm[1]));
1647  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1648  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1649  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1650  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1651 
1652  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1653  fHistManager->FillTH2(hname, fMCContainer->GetNAcceptedParticles(), nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1654 
1655  hname = TString::Format("%s/fHistNDmesons", GetName());
1656  fHistManager->FillTH1(hname, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]); // same as the number of accepted D mesons, since no selection is performed
1657 }
1658 
1663 {
1664  TString classname;
1665  if (fMCMode == kMCTruth) {
1666  classname = "AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary";
1667  fCurrentDmesonJetInfo = new AliDmesonMCInfoSummary();
1668  }
1669  else {
1670  switch (fCandidateType) {
1671  case kD0toKpi:
1672  case kD0toKpiLikeSign:
1673  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1674  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1675  break;
1676  case kDstartoKpipi:
1677  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1678  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1679  break;
1680  }
1681  }
1682  TString treeName = TString::Format("%s_%s", taskName, GetName());
1683  fTree = new TTree(treeName, treeName);
1684  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
1685  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1686  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1687  if (fJetDefinitions[i].fRhoName.IsNull()) {
1688  fCurrentJetInfo[i] = new AliJetInfoSummary();
1689  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1690  }
1691  else {
1692  fCurrentJetInfo[i] = new AliJetInfoPbPbSummary();
1693  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoPbPbSummary", &fCurrentJetInfo[i]);
1694  }
1695  }
1696 
1697  return fTree;
1698 }
1699 
1704 {
1705  TString hname;
1706 
1707  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
1708 
1709  for (auto &jetDef : fJetDefinitions) {
1710 
1711  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
1712 
1713  Double_t radius = jetDef.fRadius;
1714 
1715  TString title[30] = {""};
1716  Int_t nbins[30] = {0 };
1717  Double_t min [30] = {0.};
1718  Double_t max [30] = {0.};
1719  Int_t dim = 0 ;
1720 
1721  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1722  nbins[dim] = nPtBins;
1723  min[dim] = 0;
1724  max[dim] = fMaxPt;
1725  dim++;
1726 
1727  if ((enabledAxis & kPositionD) != 0) {
1728  title[dim] = "#eta_{D}";
1729  nbins[dim] = 50;
1730  min[dim] = -1;
1731  max[dim] = 1;
1732  dim++;
1733 
1734  title[dim] = "#phi_{D} (rad)";
1735  nbins[dim] = 150;
1736  min[dim] = 0;
1737  max[dim] = TMath::TwoPi();
1738  dim++;
1739  }
1740 
1741  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1742  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1743  nbins[dim] = fNMassBins;
1744  min[dim] = fMinMass;
1745  max[dim] = fMaxMass;
1746  dim++;
1747  }
1748 
1749  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1750  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1751  nbins[dim] = fNMassBins;
1752  min[dim] = fMinMass;
1753  max[dim] = fMaxMass;
1754  dim++;
1755  }
1756 
1757  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1758  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1759  nbins[dim] = fNMassBins;
1760  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1761  dim++;
1762  }
1763 
1764  if (fCandidateType == kDstartoKpipi) {
1765  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1766  nbins[dim] = fNMassBins*6;
1767  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1768 
1769  // subtract mass of D0
1770  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1771  min[dim] -= D0mass;
1772  max[dim] -= D0mass;
1773 
1774  dim++;
1775  }
1776 
1777  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1778  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1779  nbins[dim] = 100;
1780  min[dim] = 0;
1781  max[dim] = 25;
1782  dim++;
1783  }
1784 
1785  title[dim] = "#it{z}_{D}";
1786  nbins[dim] = 110;
1787  min[dim] = 0;
1788  max[dim] = 1.10;
1789  dim++;
1790 
1791  if ((enabledAxis & kDeltaR) != 0) {
1792  title[dim] = "#Delta R_{D-jet}";
1793  nbins[dim] = 100;
1794  min[dim] = 0;
1795  max[dim] = radius * 1.5;
1796  dim++;
1797  }
1798 
1799  if ((enabledAxis & kDeltaEta) != 0) {
1800  title[dim] = "#eta_{D} - #eta_{jet}";
1801  nbins[dim] = 100;
1802  min[dim] = -radius * 1.2;
1803  max[dim] = radius * 1.2;
1804  dim++;
1805  }
1806 
1807  if ((enabledAxis & kDeltaPhi) != 0) {
1808  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1809  nbins[dim] = 100;
1810  min[dim] = -radius * 1.2;
1811  max[dim] = radius * 1.2;
1812  dim++;
1813  }
1814 
1815  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1816  nbins[dim] = nPtBins;
1817  min[dim] = 0;
1818  max[dim] = fMaxPt;
1819  dim++;
1820 
1821  if ((enabledAxis & kPositionJet) != 0) {
1822  title[dim] = "#eta_{jet}";
1823  nbins[dim] = 50;
1824  min[dim] = -1;
1825  max[dim] = 1;
1826  dim++;
1827 
1828  title[dim] = "#phi_{jet} (rad)";
1829  nbins[dim] = 150;
1830  min[dim] = 0;
1831  max[dim] = TMath::TwoPi();
1832  dim++;
1833  }
1834 
1835  if ((enabledAxis & kJetConstituents) != 0) {
1836  title[dim] = "No. of constituents";
1837  nbins[dim] = 50;
1838  min[dim] = -0.5;
1839  max[dim] = 49.5;
1840  dim++;
1841  }
1842 
1843  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1844  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1845  for (Int_t j = 0; j < dim; j++) {
1846  h->GetAxis(j)->SetTitle(title[j]);
1847  }
1848  }
1849 }
1850 
1855 {
1856  TString hname;
1857  fFirstPartons.clear();
1858  fLastPartons.clear();
1859  for (auto& dmeson_pair : fDmesonJets) {
1860  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
1861  Int_t accJets = 0;
1862  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1863  fCurrentJetInfo[ij]->Reset();
1864  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1865  if (!jet) continue;
1866  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1867  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1868  fHistManager->FillTH1(hname, jet->Pt());
1869  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1870  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1871  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1872  fHistManager->FillTH1(hname, jet->Eta());
1873  continue;
1874  }
1875  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
1876  accJets++;
1877  }
1878  if (accJets > 0) {
1879  fFirstPartons[dmeson_pair.second.fFirstParton] = dmeson_pair.second.fFirstPartonType;
1880  fLastPartons[dmeson_pair.second.fLastParton] = dmeson_pair.second.fLastPartonType;
1881 
1882  fTree->Fill();
1883  }
1884  else {
1885  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1886  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1887  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1888  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1889  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1890  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1891  if (fMCMode != kMCTruth) {
1892  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1893  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1894  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1895  }
1896  else if (fCandidateType == kDstartoKpipi) {
1897  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1898  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1899 
1900  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1901  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1902  }
1903  }
1904  }
1905  }
1906 
1907  hname = TString::Format("%s/fHistFirstPartonPt", GetName());
1908  TH1* histFirstPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
1909  hname = TString::Format("%s/fHistFirstPartonEta", GetName());
1910  TH1* histFirstPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
1911  hname = TString::Format("%s/fHistFirstPartonPhi", GetName());
1912  TH1* histFirstPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
1913  hname = TString::Format("%s/fHistFirstPartonType", GetName());
1914  TH1* histFirstPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
1915 
1916  for (auto parton : fFirstPartons) {
1917  if (!parton.first) continue;
1918  histFirstPartonPt->Fill(parton.first->Pt());
1919  histFirstPartonEta->Fill(parton.first->Eta());
1920  histFirstPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
1921  histFirstPartonType->Fill(parton.second);
1922  }
1923 
1924  hname = TString::Format("%s/fHistLastPartonPt", GetName());
1925  TH1* histLastPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
1926  hname = TString::Format("%s/fHistLastPartonEta", GetName());
1927  TH1* histLastPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
1928  hname = TString::Format("%s/fHistLastPartonPhi", GetName());
1929  TH1* histLastPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
1930  hname = TString::Format("%s/fHistLastPartonType", GetName());
1931  TH1* histLastPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
1932 
1933  for (auto parton : fLastPartons) {
1934  if (!parton.first) continue;
1935  histLastPartonPt->Fill(parton.first->Pt());
1936  histLastPartonEta->Fill(parton.first->Eta());
1937  histLastPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
1938  histLastPartonType->Fill(parton.second);
1939  }
1940 
1941  return kTRUE;
1942 }
1943 
1949 {
1950  TString hname;
1951  fFirstPartons.clear();
1952  fLastPartons.clear();
1953  for (auto& dmeson_pair : fDmesonJets) {
1954  Int_t accJets = 0;
1955  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1956  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1957  if (!jet) continue;
1958  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1959  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1960  fHistManager->FillTH1(hname, jet->Pt());
1961  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1962  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1963  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1964  fHistManager->FillTH1(hname, jet->Eta());
1965  continue;
1966  }
1967  accJets++;
1968  }
1969  if (accJets > 0) {
1970  fFirstPartons[dmeson_pair.second.fFirstParton] = dmeson_pair.second.fFirstPartonType;
1971  fLastPartons[dmeson_pair.second.fLastParton] = dmeson_pair.second.fLastPartonType;
1972  }
1973  else {
1974  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1975  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1976  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1977  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1978  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1979  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1980  if (fMCMode != kMCTruth) {
1981  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1982  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1983  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1984  }
1985  else if (fCandidateType == kDstartoKpipi) {
1986  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1987  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1988 
1989  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1990  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1991  }
1992  }
1993  }
1994  }
1995 
1996  hname = TString::Format("%s/fHistFirstPartonPt", GetName());
1997  TH1* histFirstPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
1998  hname = TString::Format("%s/fHistFirstPartonEta", GetName());
1999  TH1* histFirstPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
2000  hname = TString::Format("%s/fHistFirstPartonPhi", GetName());
2001  TH1* histFirstPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
2002  hname = TString::Format("%s/fHistFirstPartonType", GetName());
2003  TH1* histFirstPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
2004 
2005  for (auto parton : fFirstPartons) {
2006  if (!parton.first) continue;
2007  histFirstPartonPt->Fill(parton.first->Pt());
2008  histFirstPartonEta->Fill(parton.first->Eta());
2009  histFirstPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
2010  histFirstPartonType->Fill(parton.second);
2011  }
2012 
2013  hname = TString::Format("%s/fHistLastPartonPt", GetName());
2014  TH1* histLastPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
2015  hname = TString::Format("%s/fHistLastPartonEta", GetName());
2016  TH1* histLastPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
2017  hname = TString::Format("%s/fHistLastPartonPhi", GetName());
2018  TH1* histLastPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
2019  hname = TString::Format("%s/fHistLastPartonType", GetName());
2020  TH1* histLastPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
2021 
2022  for (auto parton : fLastPartons) {
2023  if (!parton.first) continue;
2024  histLastPartonPt->Fill(parton.first->Pt());
2025  histLastPartonEta->Fill(parton.first->Eta());
2026  histLastPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
2027  histLastPartonType->Fill(parton.second);
2028  }
2029 
2030  return kTRUE;
2031 }
2032 
2037 {
2038  TString hname;
2039 
2040  for (auto& dmeson_pair : fDmesonJets) {
2041  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
2042  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2043  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2044  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2045  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2046  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2047  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2048  }
2049  }
2050 
2051  for (auto &jetDef : fJetDefinitions) {
2052 
2053  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
2054  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
2055 
2056  for (auto& dmeson_pair : fDmesonJets) {
2057  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
2058  if (!jet) continue;
2059  if (!jetDef.IsJetInAcceptance(*jet)) {
2060  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
2061  fHistManager->FillTH1(hname, jet->Pt());
2062  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
2063  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2064  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
2065  fHistManager->FillTH1(hname, jet->Eta());
2066  continue;
2067  }
2068  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
2069  }
2070  }
2071 
2072  return kTRUE;
2073 }
2074 
2081 {
2082  // Fill the THnSparse histogram.
2083 
2084  Double_t contents[30] = {0.};
2085 
2086  Double_t z = DmesonJet.GetZ(n);
2087  Double_t deltaPhi = 0;
2088  Double_t deltaEta = 0;
2089  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
2090 
2091  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
2092  if (it == DmesonJet.fJets.end()) return kFALSE;
2093 
2094  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
2095  TString title(h->GetAxis(i)->GetTitle());
2096  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
2097  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
2098  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
2099  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
2100  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
2101  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
2102  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
2103  else if (title=="#it{z}_{D}") contents[i] = z;
2104  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
2105  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
2106  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
2107  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
2108  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
2109  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
2110  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
2111  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
2112  }
2113 
2114  h->Fill(contents);
2115 
2116  return kTRUE;
2117 }
2118 
2119 // Definitions of class AliAnalysisTaskDmesonJets
2120 
2124 
2128  fAnalysisEngines(),
2129  fEnabledAxis(0),
2131  fHistManager(),
2132  fApplyKinematicCuts(kTRUE),
2133  fNOutputTrees(0),
2134  fTrackEfficiency(0),
2135  fMCContainer(0),
2136  fAodEvent(0),
2137  fFastJetWrapper(0)
2138 {
2139  SetMakeGeneralHistograms(kTRUE);
2140 }
2141 
2146  AliAnalysisTaskEmcalLight(name, kTRUE),
2147  fAnalysisEngines(),
2148  fEnabledAxis(k2ProngInvMass),
2149  fOutputType(kTreeOutput),
2150  fHistManager(name),
2151  fApplyKinematicCuts(kTRUE),
2152  fNOutputTrees(nOutputTrees),
2153  fTrackEfficiency(0),
2154  fMCContainer(0),
2155  fAodEvent(0),
2156  fFastJetWrapper(0)
2157 {
2158  SetMakeGeneralHistograms(kTRUE);
2159  for (Int_t i = 0; i < nOutputTrees; i++){
2160  DefineOutput(2+i, TTree::Class());
2161  }
2162 }
2163 
2166 {
2167  if (fFastJetWrapper) delete fFastJetWrapper;
2168 }
2169 
2177 {
2178  AliRDHFCuts* analysiscuts = 0;
2179  TFile* filecuts = TFile::Open(cutfname);
2180  if (!filecuts || filecuts->IsZombie()) {
2181  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Input file not found: will use std cuts.");
2182  filecuts = 0;
2183  }
2184 
2185  if (filecuts) analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
2186 
2187  if (!analysiscuts) {
2188  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Could not find analysis cuts '%s' in '%s'.", cutsname.Data(), cutfname.Data());
2189  if (filecuts) {
2190  filecuts->ls();
2191  }
2192  }
2193  else {
2194  ::Info("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Cuts '%s' loaded from file '%s'", cutsname.Data(), cutfname.Data());
2195  }
2196 
2197  return analysiscuts;
2198 }
2199 
2209 {
2211  return AddAnalysisEngine(type, cutfname, MCmode, jetDef, rhoName);
2212 }
2213 
2223 {
2224  AliRDHFCuts* cuts = 0;
2225 
2226  if (!cutfname.IsNull()) {
2227  TString cutsname;
2228 
2229  switch (type) {
2230  case kD0toKpi:
2231  case kD0toKpiLikeSign:
2232  cutsname = "D0toKpiCuts";
2233  break;
2234  case kDstartoKpipi:
2235  cutsname = "DStartoKpipiCuts";
2236  break;
2237  default:
2238  return 0;
2239  }
2240 
2241  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
2242  if (cuts) cuts->PrintAll();
2243  }
2244 
2245  AnalysisEngine eng(type, MCmode, cuts);
2246 
2247  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
2248 
2249  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
2250  eng.AddJetDefinition(jetDef);
2251  it = fAnalysisEngines.insert(it, eng);
2252  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
2253  }
2254  else {
2255  AnalysisEngine* found_eng = &(*it);
2256  ::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());
2257  found_eng->AddJetDefinition(jetDef);
2258 
2259  if (cuts) {
2260  if (found_eng->fRDHFCuts != 0) ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
2261  found_eng->SetRDHFCuts(cuts);
2262  }
2263  }
2264 
2265  return &(*it);
2266 }
2267 
2268 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
2269 {
2270  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
2271  while (it != fAnalysisEngines.end() && (*it) < eng) it++;
2272  return it;
2273 }
2274 
2277 {
2278  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
2279 
2281 
2282  // Define histograms
2283  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
2284 
2285  TString hname;
2286  TString htitle;
2287  TH1* h = 0;
2288  Int_t treeSlot = 0;
2289 
2290  hname = "fHistCharmPt";
2291  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2292  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2293 
2294  hname = "fHistCharmEta";
2295  htitle = hname + ";#eta_{charm};counts";
2296  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2297 
2298  hname = "fHistCharmPhi";
2299  htitle = hname + ";#phi_{charm};counts";
2300  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2301 
2302  hname = "fHistCharmPt_Eta05";
2303  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2304  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2305 
2306  hname = "fHistBottomPt";
2307  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2308  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2309 
2310  hname = "fHistBottomEta";
2311  htitle = hname + ";#eta_{bottom};counts";
2312  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2313 
2314  hname = "fHistBottomPhi";
2315  htitle = hname + ";#phi_{bottom};counts";
2316  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2317 
2318  hname = "fHistBottomPt_Eta05";
2319  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2320  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2321 
2322  hname = "fHistHighestPartonPt";
2323  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2324  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2325 
2326  hname = "fHistHighestPartonType";
2327  htitle = hname + ";type;counts";
2328  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2329 
2330  hname = "fHistNHeavyQuarks";
2331  htitle = hname + ";number of heavy-quarks;counts";
2332  fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
2333 
2334  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
2335  for (auto &param : fAnalysisEngines) {
2336  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
2337 
2338  param.fHistManager = &fHistManager;
2339 
2340  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", param.GetName());
2341  htitle = hname + ";#it{N}_{tracks};#it{N}_{D};events";
2342  h = fHistManager.CreateTH2(hname, htitle, 251, -0.5, 250.5, 21, -0.5, 20.5);
2343 
2344  hname = TString::Format("%s/fHistNTotAcceptedDmesons", param.GetName());
2345  htitle = hname + ";;#it{N}_{D}";
2346  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
2347 
2348  hname = TString::Format("%s/fHistNDmesons", param.GetName());
2349  htitle = hname + ";#it{N}_{D};events";
2350  h = fHistManager.CreateTH1(hname, htitle, 501, -0.5, 500.5);
2351 
2352  hname = TString::Format("%s/fHistNEvents", param.GetName());
2353  htitle = hname + ";Event status;counts";
2354  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
2355  h->GetXaxis()->SetBinLabel(1, "Accepted");
2356  h->GetXaxis()->SetBinLabel(2, "Rejected");
2357 
2358  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
2359  htitle = hname + ";Rejection reason;counts";
2360  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
2361 
2362  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
2363  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
2364  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2365 
2366  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
2367  htitle = hname + ";#it{#eta}_{D};counts";
2368  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2369 
2370  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
2371  htitle = hname + ";#it{#phi}_{D};counts";
2372  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2373 
2374  if (param.fMCMode != kMCTruth) {
2375  if (param.fCandidateType == kD0toKpi || param.fCandidateType == kD0toKpiLikeSign) {
2376  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
2377  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2378  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
2379  }
2380  else if (param.fCandidateType == kDstartoKpipi) {
2381  Double_t min = 0;
2382  Double_t max = 0;
2383 
2384  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
2385  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2386  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
2387  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
2388 
2389  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2390  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
2391  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2392  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
2393  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
2394  }
2395  }
2396 
2397  if (param.fMCMode == kMCTruth) {
2398  hname = TString::Format("%s/fHistFirstPartonPt", param.GetName());
2399  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
2400  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2401 
2402  hname = TString::Format("%s/fHistFirstPartonEta", param.GetName());
2403  htitle = hname + ";#eta_{parton};counts";
2404  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2405 
2406  hname = TString::Format("%s/fHistFirstPartonPhi", param.GetName());
2407  htitle = hname + ";#phi_{parton};counts";
2408  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2409 
2410  hname = TString::Format("%s/fHistFirstPartonType", param.GetName());
2411  htitle = hname + ";type;counts";
2412  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2413 
2414  hname = TString::Format("%s/fHistLastPartonPt", param.GetName());
2415  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
2416  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2417 
2418  hname = TString::Format("%s/fHistLastPartonEta", param.GetName());
2419  htitle = hname + ";#eta_{parton};counts";
2420  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2421 
2422  hname = TString::Format("%s/fHistLastPartonPhi", param.GetName());
2423  htitle = hname + ";#phi_{parton};counts";
2424  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2425 
2426  hname = TString::Format("%s/fHistLastPartonType", param.GetName());
2427  htitle = hname + ";type;counts";
2428  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2429  }
2430 
2431  for (auto& jetDef : param.fJetDefinitions) {
2432  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef.GetName());
2433 
2434  if (param.fMCMode == kMCTruth) {
2435  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", param.GetName(), jetDef.GetName());
2436  htitle = hname + ";#it{N}_{constituents};#it{N}_{D};counts";
2437  h = fHistManager.CreateTH2(hname, htitle, 51, -0.5, 50.5, 10, 0.5, 10.5);
2438  }
2439 
2440  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef.GetName());
2441  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2442  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2443  SetRejectionReasonLabels(h->GetXaxis());
2444 
2445  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef.GetName());
2446  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2447  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2448  SetRejectionReasonLabels(h->GetXaxis());
2449 
2450  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef.GetName());
2451  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
2452  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2453  SetRejectionReasonLabels(h->GetXaxis());
2454 
2455  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef.GetName());
2456  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
2457  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
2458 
2459  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef.GetName());
2460  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
2461  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2462 
2463  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef.GetName());
2464  htitle = hname + ";#it{#eta}_{jet};counts";
2465  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2466 
2467  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef.GetName());
2468  htitle = hname + ";#it{#phi}_{jet};counts";
2469  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2470  }
2471  switch (fOutputType) {
2472  case kTreeOutput:
2473  param.BuildTree(GetName());
2474  if (treeSlot < fNOutputTrees) {
2475  param.AssignDataSlot(treeSlot+2);
2476  treeSlot++;
2478  }
2479  else {
2480  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
2481  }
2482  break;
2483  case kTHnOutput:
2484  param.BuildHnSparse(fEnabledAxis);
2485  break;
2486  case kNoOutput:
2487  break;
2488  }
2489  }
2490 
2492 
2493  PostData(1, fOutput);
2494 }
2495 
2499 {
2501 
2502  // Load the event
2503  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
2504 
2505  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
2506 
2507  // TODO: make this settable
2508  fFastJetWrapper->SetAreaType(fastjet::active_area_explicit_ghosts);
2509  fFastJetWrapper->SetGhostArea(0.005);
2510 
2511  if (!fAodEvent) {
2512  AliError(Form("This task need an AOD event (Task '%s'). Expect troubles...", GetName()));
2513  //return;
2514  }
2515 
2516  TRandom* rnd = 0;
2517  if (fTrackEfficiency > 0 && fTrackEfficiency < 1) rnd = new TRandom3(0);
2518 
2519  for (auto cont_it : fParticleCollArray) {
2520  AliHFAODMCParticleContainer* part_cont = dynamic_cast<AliHFAODMCParticleContainer*>(cont_it.second);
2521  if (part_cont) fMCContainer = part_cont;
2522  }
2523 
2524  for (auto &params : fAnalysisEngines) {
2525 
2526  params.fAodEvent = fAodEvent;
2527  params.fFastJetWrapper = fFastJetWrapper;
2528  params.fTrackEfficiency = fTrackEfficiency;
2529  params.fRandomGen = rnd;
2530 
2531  for (auto &jetdef: params.fJetDefinitions) {
2532  if (!jetdef.fRhoName.IsNull()) {
2533  jetdef.fRho = dynamic_cast<AliRhoParameter*>(fInputEvent->FindListObject(jetdef.fRhoName));
2534  if (!jetdef.fRho) {
2535  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2536  "%s: Could not find rho object '%s' for engine '%s'",
2537  jetdef.fRhoName.Data(), GetName(), params.GetName());
2538  }
2539  }
2540  }
2541 
2542  if (!params.fRDHFCuts) {
2543  if (params.fMCMode == kMCTruth) {
2544  ::Warning("AliAnalysisTaskDmesonJets::ExecOnce",
2545  "%s: RDHF cuts not provided for engine '%s'.",
2546  GetName(), params.GetName());
2547  }
2548  else {
2549  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2550  "%s: RDHF cuts not provided. Engine '%s' disabled.",
2551  GetName(), params.GetName());
2552  params.fInhibit = kTRUE;
2553  }
2554  }
2555 
2556  params.fMCContainer = fMCContainer;
2557 
2558  for (auto cont_it : fParticleCollArray) {
2559  AliTrackContainer* track_cont = dynamic_cast<AliTrackContainer*>(cont_it.second);
2560  if (track_cont) params.fTrackContainers.push_back(track_cont);
2561  }
2562 
2563  for (auto cont_it : fClusterCollArray) params.fClusterContainers.push_back(cont_it.second);
2564 
2565  if (fAodEvent) params.Init(fGeom, fAodEvent->GetRunNumber());
2566 
2567  if (params.fMCMode != kMCTruth && fAodEvent) {
2568  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
2569 
2570  if (params.fCandidateArray) {
2571  TString className;
2572  if (params.fCandidateType == kD0toKpi || params.fCandidateType == kD0toKpiLikeSign) {
2573  className = "AliAODRecoDecayHF2Prong";
2574  }
2575  else if (params.fCandidateType == kDstartoKpipi) {
2576  className = "AliAODRecoCascadeHF";
2577  }
2578  if (!params.fCandidateArray->GetClass()->InheritsFrom(className)) {
2579  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2580  "%s: Objects of type %s in %s are not inherited from %s! Task will be disabled!",
2581  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName(), className.Data());
2582  params.fCandidateArray = 0;
2583  params.fInhibit = kTRUE;
2584  }
2585  }
2586  else {
2587  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2588  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
2589  params.fBranchName.Data(), params.GetName());
2590  params.fInhibit = kTRUE;
2591  }
2592  }
2593 
2594  if (params.fMCMode != kNoMC) {
2595  if (!params.fMCContainer) {
2596  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2597  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
2598  params.GetName());
2599  params.fInhibit = kTRUE;
2600  }
2601  }
2602 
2603  if (params.fMCMode != kMCTruth) {
2604  if (params.fTrackContainers.size() == 0 && params.fClusterContainers.size() == 0) {
2605  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2606  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
2607  params.GetName());
2608  params.fInhibit = kTRUE;
2609  }
2610  }
2611  }
2612 }
2613 
2618 {
2619  TString hname;
2620 
2621  // Fix for temporary bug in ESDfilter
2622  // The AODs with null vertex pointer didn't pass the PhysSel
2623  // Now adding an entry in the histogram so as to check that this is actually cutting anything out
2624  if (fAodEvent && (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001)) {
2625  for (auto &eng : fAnalysisEngines) {
2626  if (eng.fInhibit) continue;
2627  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2628  fHistManager.FillTH1(hname, "ESDfilterBug");
2629  }
2630  return kFALSE;
2631  }
2632 
2633  if (fAodEvent) {
2634  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
2635  if (matchingAODdeltaAODlevel <= 0) {
2636  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
2637  for (auto &eng : fAnalysisEngines) {
2638  if (eng.fInhibit) continue;
2639  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2640  fHistManager.FillTH1(hname, "MismatchDeltaAOD");
2641  }
2642  return kFALSE;
2643  }
2644  }
2645 
2646  for (auto &eng : fAnalysisEngines) {
2647  eng.fDmesonJets.clear();
2648  if (eng.fInhibit) continue;
2649 
2650  //Event selection
2651  hname = TString::Format("%s/fHistNEvents", eng.GetName());
2652  if (fAodEvent) {
2653  Bool_t iseventselected = kTRUE;
2654  if (eng.fRDHFCuts) iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
2655  if (!iseventselected) {
2656  fHistManager.FillTH1(hname, "Rejected");
2657  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2658  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
2659  TString label;
2660  do {
2661  label = GetHFEventRejectionReasonLabel(bitmap);
2662  if (label.IsNull()) break;
2663  fHistManager.FillTH1(hname, label);
2664  } while (true);
2665 
2666  continue;
2667  }
2668  }
2669 
2670  fHistManager.FillTH1(hname, "Accepted");
2671 
2672  AliDebug(2, "Event selected");
2673 
2674  eng.RunAnalysis();
2675  }
2676  return kTRUE;
2677 }
2678 
2683 {
2684  for (auto &param : fAnalysisEngines) {
2685  if (param.fInhibit) continue;
2686 
2687  if (fOutputType == kTreeOutput) {
2688  param.FillTree(fApplyKinematicCuts);
2690  }
2691  else if (fOutputType == kTHnOutput) {
2692  param.FillHnSparse(fApplyKinematicCuts);
2693  }
2694  }
2696  return kTRUE;
2697 }
2698 
2701 {
2702  auto itcont = fMCContainer->all_momentum();
2703  Int_t nHQ = 0;
2704  Double_t highestPartonPt = 0;
2705  Int_t absPdgHighParton = 0;
2706  for (auto part : itcont) {
2707  Int_t absPdgCode = TMath::Abs(part.second->GetPdgCode());
2708 
2709  // Skip all particles that are not either quarks or gluons
2710  if (absPdgCode > 9 && absPdgCode != 21) continue;
2711 
2712  // Look for highest momentum parton
2713  if (highestPartonPt < part.first.Pt()) {
2714  highestPartonPt = part.first.Pt();
2715  absPdgHighParton = absPdgCode;
2716  }
2717  /*
2718  // Look for the mother PDG code
2719  Int_t motherIndex = part.second->GetMother();
2720  AliAODMCParticle *mother = 0;
2721  Int_t motherPdg = 0;
2722  Double_t motherPt = 0;
2723  if (motherIndex >= 0) {
2724  mother = fMCContainer->GetMCParticle(motherIndex);
2725  if (motherIndex) {
2726  motherPdg = TMath::Abs(mother->GetPdgCode());
2727  motherPt = mother->Pt();
2728  }
2729  }
2730  */
2731  if (absPdgCode != 4 && absPdgCode != 5) continue;
2732  Bool_t notLastInPartonShower = kFALSE;
2733  for (Int_t idaugh = 0; idaugh < 2; idaugh++){
2734  Int_t daughterIndex = part.second->GetDaughter(idaugh);
2735  if (daughterIndex < 0) {
2736  AliDebug(10, Form("Could not find daughter of heavy quark (pdg=%d, pt=%.3f)!", absPdgCode, part.first.Pt()));
2737  continue;
2738  }
2739  AliAODMCParticle *daughter = fMCContainer->GetMCParticle(daughterIndex);
2740  if (!daughter) {
2741  AliDebug(10, Form("Could not find daughter %d of heavy quark (pdg=%d, pt=%.3f)!", daughterIndex, absPdgCode, part.first.Pt()));
2742  continue;
2743  }
2744  Int_t daughterAbsPdgCode = TMath::Abs(daughter->GetPdgCode());
2745  if (daughterAbsPdgCode <= 9 || daughterAbsPdgCode == 21) notLastInPartonShower = kTRUE; // this parton is not the last parton in the shower
2746  AliDebug(10, Form("Found daughter with PDG=%d, pt=%.3f", daughterAbsPdgCode, daughter->Pt()));
2747  }
2748  if (notLastInPartonShower) continue;
2749 
2750  if (absPdgCode == 4) {
2751  fHistManager.FillTH1("fHistCharmPt", part.first.Pt());
2752  fHistManager.FillTH1("fHistCharmEta", part.first.Eta());
2753  fHistManager.FillTH1("fHistCharmPhi", part.first.Phi_0_2pi());
2754  if (TMath::Abs(part.first.Eta()) < 0.5) fHistManager.FillTH1("fHistCharmPt_Eta05", part.first.Pt());
2755  }
2756  else if (absPdgCode == 5) {
2757  fHistManager.FillTH1("fHistBottomPt", part.first.Pt());
2758  fHistManager.FillTH1("fHistBottomEta", part.first.Eta());
2759  fHistManager.FillTH1("fHistBottomPhi", part.first.Phi_0_2pi());
2760  if (TMath::Abs(part.first.Eta()) < 0.5) fHistManager.FillTH1("fHistBottomPt_Eta05", part.first.Pt());
2761  }
2762  nHQ++;
2763  }
2764  fHistManager.FillTH1("fHistNHeavyQuarks", nHQ);
2765  fHistManager.FillTH1("fHistHighestPartonPt",highestPartonPt);
2766  Int_t partonType = 0;
2767  if (absPdgHighParton == 9 || absPdgHighParton == 21) {
2768  partonType = 7;
2769  }
2770  else {
2771  partonType = absPdgHighParton;
2772  }
2773  fHistManager.FillTH1("fHistHighestPartonType",partonType);
2774 }
2775 
2784 {
2785  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
2786 
2787  Double_t mass = part->Mass();
2788 
2789  // To make sure that the PDG mass value is not at the edge of a bin
2790  if (nbins % 2 == 0) {
2791  minMass = mass - range / 2 - range / nbins / 2;
2792  maxMass = mass + range / 2 - range / nbins / 2;
2793  }
2794  else {
2795  minMass = mass - range / 2;
2796  maxMass = mass + range / 2;
2797  }
2798 }
2799 
2806 {
2807  static TString label;
2808  label = "";
2809 
2810  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
2811  label = "NotSelTrigger";
2812  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
2813  return label.Data();
2814  }
2815  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
2816  label = "NoVertex";
2817  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
2818  return label.Data();
2819  }
2820  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
2821  label = "TooFewVtxContrib";
2822  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
2823  return label.Data();
2824  }
2825  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
2826  label = "ZVtxOutFid";
2827  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
2828  return label.Data();
2829  }
2830  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
2831  label = "Pileup";
2832  bitmap &= ~BIT(AliRDHFCuts::kPileup);
2833  return label.Data();
2834  }
2835  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
2836  label = "OutsideCentrality";
2837  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
2838  return label.Data();
2839  }
2840  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
2841  label = "PhysicsSelection";
2842  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
2843  return label.Data();
2844  }
2845  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
2846  label = "BadSPDVertex";
2847  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
2848  return label.Data();
2849  }
2850  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
2851  label = "ZVtxSPDOutFid";
2852  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
2853  return label.Data();
2854  }
2855  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
2856  label = "CentralityFlattening";
2857  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
2858  return label.Data();
2859  }
2860  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
2861  label = "BadTrackV0Correl";
2862  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
2863  return label.Data();
2864  }
2865 
2866  return label.Data();
2867 }
2868 
2875 {
2876  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
2877  PostData(eng.GetDataSlotNumber(), eng.GetTree());
2878  return eng.GetDataSlotNumber();
2879  }
2880  else {
2881  return -1;
2882  }
2883 }
2884 
2894 {
2895  // Get the pointer to the existing analysis manager via the static access method
2896  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
2897  if (!mgr) {
2898  ::Error("AddTaskDmesonJets", "No analysis manager to connect to.");
2899  return NULL;
2900  }
2901 
2902  // Check the analysis type using the event handlers connected to the analysis manager
2903  AliVEventHandler* handler = mgr->GetInputEventHandler();
2904  if (!handler) {
2905  ::Error("AddTaskEmcalJetSpectraQA", "This task requires an input event handler");
2906  return NULL;
2907  }
2908 
2909  EDataType_t dataType = kUnknownDataType;
2910 
2911  if (handler->InheritsFrom("AliESDInputHandler")) {
2912  dataType = kESD;
2913  }
2914  else if (handler->InheritsFrom("AliAODInputHandler")) {
2915  dataType = kAOD;
2916  }
2917 
2918  // Init the task and do settings
2919  if (ntracks == "usedefault") {
2920  if (dataType == kESD) {
2921  ntracks = "Tracks";
2922  }
2923  else if (dataType == kAOD) {
2924  ntracks = "tracks";
2925  }
2926  else {
2927  ntracks = "";
2928  }
2929  }
2930 
2931  if (nclusters == "usedefault") {
2932  if (dataType == kESD) {
2933  nclusters = "CaloClusters";
2934  }
2935  else if (dataType == kAOD) {
2936  nclusters = "caloClusters";
2937  }
2938  else {
2939  nclusters = "";
2940  }
2941  }
2942 
2943  if (nMCpart == "usedefault") {
2944  nMCpart = "mcparticles"; // Always needs AliAODMCParticle objects
2945  }
2946 
2947  TString name("AliAnalysisTaskDmesonJets");
2948  if (strcmp(suffix, "") != 0) {
2949  name += TString::Format("_%s", suffix.Data());
2950  }
2951 
2952  AliAnalysisTaskDmesonJets* jetTask = new AliAnalysisTaskDmesonJets(name, nMaxTrees);
2953 
2954  if (!ntracks.IsNull()) {
2955  AliHFTrackContainer* trackCont = new AliHFTrackContainer(ntracks);
2956  jetTask->AdoptParticleContainer(trackCont);
2957  }
2958 
2959  if (!nMCpart.IsNull()) {
2960  AliMCParticleContainer* partCont = new AliHFAODMCParticleContainer(nMCpart);
2961  partCont->SetEtaLimits(-1.5, 1.5);
2962  partCont->SetPtLimits(0, 1000);
2963  jetTask->AdoptParticleContainer(partCont);
2964  }
2965 
2966  jetTask->AddClusterContainer(nclusters.Data());
2967 
2968  // Final settings, pass to manager and set the containers
2969  mgr->AddTask(jetTask);
2970 
2971  // Create containers for input/output
2972  AliAnalysisDataContainer* cinput1 = mgr->GetCommonInputContainer();
2973  TString contname1(name);
2974  contname1 += "_histos";
2975  AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(contname1.Data(),
2976  TList::Class(), AliAnalysisManager::kOutputContainer,
2977  Form("%s", AliAnalysisManager::GetCommonFileName()));
2978 
2979  mgr->ConnectInput(jetTask, 0, cinput1);
2980  mgr->ConnectOutput(jetTask, 1, coutput1);
2981 
2982  for (Int_t i = 0; i < nMaxTrees; i++) {
2983  TString contname = TString::Format("%s_tree_%d", name.Data(), i);
2984  AliAnalysisDataContainer *coutput = mgr->CreateContainer(contname.Data(),
2985  TTree::Class(),AliAnalysisManager::kOutputContainer,
2986  Form("%s", AliAnalysisManager::GetCommonFileName()));
2987  mgr->ConnectOutput(jetTask, 2+i, coutput);
2988  }
2989  return jetTask;
2990 }
2991 
Int_t pdg
AliClusterContainer * AddClusterContainer(std::string branchName, std::string contName="")
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
void SetDMesonCandidate(AliAODRecoDecay *c)
Analysis task for D meson jets.
virtual void Set(const AliDmesonJetInfo &source)
Invariant mass of the D0 meson candidate in GeV/c2.
virtual Int_t Run()
double Double_t
Definition: External.C:58
Bool_t IsJetInAcceptance(const AliJetInfo &jet) const
Double_t fSoftPionPt
! Transverse momentum of the soft pion of the D* candidate
EDataType_t
Switch for the data type.
virtual AliAODMCParticle * GetMCParticle(Int_t i=-1) const
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)
virtual void Set(const AliDmesonJetInfo &source, std::string n)
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
AliHFAODMCParticleContainer * fMCContainer
! MC particle container
ERecoScheme_t fRecoScheme
Jet recombination scheme (pt scheme, E scheme, ...)
const TObjArray & GetDaughterList() const
Double_t fMaxNeutralPt
Transverse momentum of the leading neutral particle (or cluster)
void RunDetectorLevelAnalysis()
Run a detector level analysis.
Double_t mass
Container with name, TClonesArray and cuts for particles.
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
virtual void UserCreateOutputObjects()
Creates the output containers.
void FillPartonLevelHistograms()
Fill histograms with parton-level information.
static Int_t CheckMatchingAODdeltaAODevents()
Lightweight class that encapsulates D meson jets.
Double_t fInvMass2Prong
! 2-prong mass of the D* candidate (w/o the soft pion)
static void CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t &minMass, Double_t &maxMass)
AliAODEvent * fAodEvent
! AOD event
Bool_t FindJet(AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, AliHFJetDefinition &jetDef)
UInt_t fEnabledAxis
Use bit defined in EAxis_t to enable axis in the THnSparse.
AliAODMCParticle * fFirstParton
! pointer to the first parton in the shower tree of the D meson (only for particle level D mesons) ...
AliRhoParameter * fRho
Object that holds the average background value.
EMCMode_t fMCMode
MC mode: No MC (data and MC detector level), background-only (MC), signal-only (MC), MC truth (particle level)
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
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:103
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:108
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:34
EJetType_t fJetType
Jet type (charged, full, neutral)
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
Create a new TH2 within the container.
virtual void Set(const AliDmesonJetInfo &source)
AliTLorentzVector fD
! 4-momentum of the D meson candidate
Double_t GetDistance(const AliJetInfo &jet, Double_t &deta, Double_t &dphi) const
TObject * FindObject(const char *name) const
Find an object inside the container.
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
AliAnalysisTaskDmesonJets()
This is the default constructor, used for ROOT I/O purposes.
virtual void Reset()
Reset the current object.
AnalysisEngine & operator=(const AnalysisEngine &source)
Lightweight class that encapsulates D meson jets for PbPb analysis.
Struct that encapsulates analysis parameters.
std::vector< AliHFJetDefinition >::iterator FindJetDefinition(const AliHFJetDefinition &eng)
Base task in the EMCAL framework (lighter version of AliAnalysisTaskEmcal)
Double_t Phi_0_2pi() const
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:102
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
AliAODTrack * GetBachelor() const
AliRDHFCuts * LoadDMesonCutsFromFile(TString cutfname, TString cutsname)
bool operator<(const AliAnalysisTaskDmesonJets::AliHFJetDefinition &lhs, const AliAnalysisTaskDmesonJets::AliHFJetDefinition &rhs)
std::list< AnalysisEngine >::iterator FindAnalysisEngine(const AnalysisEngine &eng)
std::vector< AliHFJetDefinition > fJetDefinitions
Jet definitions.
EJetAlgo_t fJetAlgo
Jet algorithm (kt, anti-kt,...)
static const char * GetHFEventRejectionReasonLabel(UInt_t &bitmap)
Lightweight class that encapsulates D*.
virtual void Clear(const Option_t *="")
Bool_t ExtractDstarAttributes(const AliAODRecoCascadeHF *DstarCand, AliDmesonJetInfo &DmesonJet, UInt_t i)
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
const AliJetInfo * GetJet(std::string n) const
Short_t 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)
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:106
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
virtual void PrintAll() const
Definition: External.C:220
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Double_t minMass
void RunParticleLevelAnalysis()
Run a particle level analysis.
std::map< std::string, AliParticleContainer * > fParticleCollArray
particle/track collection array
Bool_t fApplyKinematicCuts
Apply jet kinematic cuts.
AnalysisEngine * AddAnalysisEngine(ECandidateType_t type, TString cutfname, EMCMode_t bkgMode, EJetType_t jettype, Double_t jetradius, TString rhoName="")
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
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:104
AliTLorentzVector fMomentum
4-momentum of the jet
std::map< std::string, AliClusterContainer * > fClusterCollArray
cluster collection array
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
Create a new THnSparse within the container.
Double_t fCorrPt
Transverse momentum of the jet after subtracting the average background.
AliHFJetDefinition * AddJetDefinition(EJetType_t type, Double_t r, EJetAlgo_t algo, ERecoScheme_t reco)
Container for MC-true particles within the EMCAL framework.
Lightweight class that encapsulates D0.
Bool_t IsAnyJetInAcceptance(const AliDmesonJetInfo &dMesonJet) const
void AdoptParticleContainer(AliParticleContainer *cont)
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)