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