AliPhysics  7f2a7c4 (7f2a7c4)
 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(),
688  fTrackContainers(),
689  fClusterContainers(),
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(),
732  fTrackContainers(),
733  fClusterContainers(),
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  fTrackContainers(source.fTrackContainers),
774  fClusterContainers(source.fClusterContainers),
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  break;
831  case kD0toKpiLikeSign:
832  fCandidatePDG = 421;
833  fCandidateName = "2ProngLikeSign";
834  fNDaughters = 2;
835  fPDGdaughters.Set(fNDaughters);
836  fPDGdaughters.Reset();
837  fPDGdaughters[0] = 211; // pi
838  fPDGdaughters[1] = 321; // K
839  fBranchName = "LikeSign2Prong";
840  break;
841  case kDstartoKpipi:
842  fCandidatePDG = 413;
843  fCandidateName = "DStar";
844  fNDaughters = 3;
845  fPDGdaughters.Set(fNDaughters);
846  fPDGdaughters.Reset();
847  fPDGdaughters[0] = 211; // pi soft
848  fPDGdaughters[1] = 211; // pi fromD0
849  fPDGdaughters[2] = 321; // K from D0
850  fBranchName = "Dstar";
851  fAcceptedDecay = kDecayDStartoKpipi;
852  break;
853  default:
854  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
855  }
856 
857  CalculateMassLimits(range, fCandidatePDG, fNMassBins, fMinMass, fMaxMass);
858 }
859 
864 {
865  if (fRDHFCuts) delete fRDHFCuts;
866  fRDHFCuts = cuts;
867 }
868 
873 {
874  if (!cuts) return;
875  if (fRDHFCuts) delete fRDHFCuts;
876  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
877 }
878 
883 {
884  static TString name;
885 
886  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
887 
888  return name.Data();
889 }
890 
895 {
896  static TString name;
897 
898  name = fCandidateName;
899  switch (fMCMode) {
900  case kBackgroundOnly:
901  name += "_kBackgroundOnly";
902  break;
903  case kSignalOnly:
904  name += "_kSignalOnly";
905  break;
906  case kMCTruth:
907  name += "_MCTruth";
908  break;
909  case kWrongPID:
910  name += "_WrongPID";
911  break;
912  default:
913  break;
914  }
915 
916  return name.Data();
917 }
918 
926 {
927  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
928 
929  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
930  fJetDefinitions.push_back(def);
931  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
932  "Total number of jet definitions is now %lu.",
933  def.GetName(), GetName(), fJetDefinitions.size());
934  // For detector level set maximum track pt to 100 GeV/c
935  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
936  }
937  else {
938  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
939  }
940 
941  return &(*it);
942 }
943 
955 {
956  AliHFJetDefinition def(type, r, algo, reco);
957 
958  return AddJetDefinition(def);
959 }
960 
966 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
967 {
968  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
969  while (it != fJetDefinitions.end() && (*it) != def) it++;
970  return it;
971 }
972 
979 {
980  if (lhs.fCandidateType > rhs.fCandidateType) return false;
981  else if (lhs.fCandidateType < rhs.fCandidateType) return true;
982  else {
983  if (lhs.fMCMode < rhs.fMCMode) return true;
984  else return false;
985  }
986 }
987 
994 {
995  if (lhs.fCandidateType != rhs.fCandidateType) return false;
996  if (lhs.fMCMode != rhs.fMCMode) return false;
997  return true;
998 }
999 
1008 {
1009  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) { // D0 candidate
1010  return ExtractD0Attributes(Dcand, DmesonJet, i);
1011  }
1012  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
1013  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
1014  }
1015  else {
1016  return kFALSE;
1017  }
1018 }
1019 
1028 {
1029  AliDebug(10,"Checking if D0 meson is selected");
1030  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
1031  if (isSelected == 0) return kFALSE;
1032 
1033  Int_t MCtruthPdgCode = 0;
1034 
1035  Double_t invMassD = 0;
1036 
1037  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
1038  // Checks also the origin, and if it matches the rejected origin mask, return false
1039  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly || fMCMode == kWrongPID) {
1040  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
1041  DmesonJet.fMCLabel = mcLab;
1042 
1043  // Retrieve the generated particle (if exists) and its PDG code
1044  if (mcLab >= 0) {
1045  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1046 
1047  if (aodMcPart) {
1048  // Check origin and return false if it matches the rejected origin mask
1049  if (fRejectedOrigin) {
1050  auto origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
1051  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1052  }
1053  MCtruthPdgCode = aodMcPart->PdgCode();
1054  }
1055  }
1056  }
1057 
1058  if (isSelected == 1) { // selected as a D0
1059  if (i != 0) return kFALSE; // only one mass hypothesis thanks to PID
1060 
1061  if (fMCMode == kNoMC ||
1062  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1063  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly) ||
1064  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kWrongPID)) {
1065  // 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)
1066  AliDebug(10,"Selected as D0");
1067  invMassD = Dcand->InvMassD0();
1068  }
1069  else { // conditions above not passed, so return FALSE
1070  return kFALSE;
1071  }
1072  }
1073  else if (isSelected == 2) { // selected as a D0bar
1074  if (i != 1) return kFALSE; // only one mass hypothesis thanks to PID
1075 
1076  if (fMCMode == kNoMC ||
1077  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1078  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly) ||
1079  (MCtruthPdgCode == fCandidatePDG && fMCMode == kWrongPID)) {
1080  // 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)
1081  AliDebug(10,"Selected as D0bar");
1082  invMassD = Dcand->InvMassD0bar();
1083  }
1084  else { // conditions above not passed, so return FALSE
1085  return kFALSE;
1086  }
1087  }
1088  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
1089  AliDebug(10,"Selected as either D0 or D0bar");
1090 
1091  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
1092  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1093  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kWrongPID))) {
1094  if (i != 0) return kFALSE;
1095  AliDebug(10, "MC truth is D0");
1096  invMassD = Dcand->InvMassD0();
1097  }
1098  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1099  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kWrongPID))) {
1100  if (i != 1) return kFALSE;
1101  AliDebug(10, "MC truth is D0bar");
1102  invMassD = Dcand->InvMassD0bar();
1103  }
1104  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
1105 
1106  // Only accept it if background-only OR background-and-signal was requested
1107  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
1108  // Select D0 or D0bar depending on the i-parameter
1109  if (i == 0) {
1110  AliDebug(10, "Returning invariant mass with D0 hypothesis");
1111  invMassD = Dcand->InvMassD0();
1112  }
1113  else if (i == 1) {
1114  AliDebug(10, "Returning invariant mass with D0bar hypothesis");
1115  invMassD = Dcand->InvMassD0bar();
1116  }
1117  else { // i > 1
1118  return kFALSE;
1119  }
1120  }
1121  else { // signal-only was requested but this is not a true D0
1122  return kFALSE;
1123  }
1124  }
1125  }
1126  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
1127  return kTRUE;
1128 }
1129 
1138 {
1139  AliDebug(10,"Checking if D* meson is selected");
1140  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoCascadeHF*>(DstarCand), AliRDHFCuts::kAll, fAodEvent);
1141  if (isSelected == 0) return kFALSE;
1142 
1143  if ((i == 1 && DstarCand->Charge()>0) || (i == 0 && DstarCand->Charge()<0) || i > 1) return kFALSE; // only one mass hypothesis for the D*
1144 
1145  Int_t MCtruthPdgCode = 0;
1146 
1147  Double_t invMassD = 0;
1148 
1149  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1150  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1151  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1152 
1153  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1154  AliDebug(10, Form("MC label is %d", mcLab));
1155  DmesonJet.fMCLabel = mcLab;
1156  if (mcLab >= 0) {
1157  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1158 
1159  if (aodMcPart) {
1160  if (fRejectedOrigin) {
1161  auto origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
1162  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1163  }
1164 
1165  MCtruthPdgCode = aodMcPart->PdgCode();
1166  AliDebug(10, Form("MC truth pdg code is %d",MCtruthPdgCode));
1167  }
1168  }
1169  }
1170 
1171  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1172  if (fMCMode == kNoMC ||
1173  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1174  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1175  // 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)
1176  invMassD = DstarCand->InvMassDstarKpipi();
1177  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1178  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1179  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1180  return kTRUE;
1181  }
1182  else { // conditions above not passed, so return FALSE
1183  return kFALSE;
1184  }
1185 }
1186 
1194 {
1195  if (!part) return kUnknownDecay;
1196  if (!mcArray) return kUnknownDecay;
1197 
1199 
1200  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1201 
1202  if (part->GetNDaughters() == 2) {
1203 
1204  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1205  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1206 
1207  if (!d1 || !d2) {
1208  return decay;
1209  }
1210 
1211  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1212  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1213 
1214  if (absPdgPart == 421) { // D0 -> K pi
1215  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1216  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1217  decay = kDecayD0toKpi;
1218  }
1219  }
1220 
1221  if (absPdgPart == 413) { // D* -> D0 pi
1222  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1223  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1224  if (D0decay == kDecayD0toKpi) {
1225  decay = kDecayDStartoKpipi;
1226  }
1227  }
1228  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1229  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1230  if (D0decay == kDecayD0toKpi) {
1231  decay = kDecayDStartoKpipi;
1232  }
1233  }
1234  }
1235  }
1236 
1237  return decay;
1238 }
1239 
1246 std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> AliAnalysisTaskDmesonJets::AnalysisEngine::CheckOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, Bool_t firstParton)
1247 {
1248  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1249 
1250  std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> result(kUnknownQuark, 0);
1251 
1252  if (!part) return result;
1253  if (!mcArray) return result;
1254 
1255  Int_t mother = part->GetMother();
1256  while (mother >= 0) {
1257  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1258  if (mcGranma) {
1259  Int_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
1260 
1261  if (abspdgGranma == 1) result = {kFromDown, mcGranma};
1262  if (abspdgGranma == 2) result = {kFromUp, mcGranma};
1263  if (abspdgGranma == 3) result = {kFromStrange, mcGranma};
1264  if (abspdgGranma == 4) result = {kFromCharm, mcGranma};
1265  if (abspdgGranma == 5) result = {kFromBottom, mcGranma};
1266  if (abspdgGranma == 6) result = {kFromTop, mcGranma};
1267  if (abspdgGranma == 9 || abspdgGranma == 21) result = {kFromGluon, mcGranma};
1268 
1269  // If looking for the very first parton in the hard scattering, it will continue the loop until it cannot find a mother particle
1270  if (result.first != kUnknownQuark && !firstParton) return result;
1271 
1272  mother = mcGranma->GetMother();
1273  }
1274  else {
1275  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1276  break;
1277  }
1278  }
1279 
1280  return result;
1281 }
1282 
1285 {
1286  for (auto& jetDef : fJetDefinitions) {
1287  jetDef.fJets.clear();
1288  }
1289 
1290  if (fMCMode == kMCTruth) {
1291  RunParticleLevelAnalysis();
1292  }
1293  else {
1294  RunDetectorLevelAnalysis();
1295  }
1296 }
1297 
1300 {
1301  // Fill the vertex info of the candidates
1302  // Needed for reduced delta AOD, where the vertex info has been deleted
1303  // to reduce the delta AOD file size
1305 
1306  const Int_t nD = fCandidateArray->GetEntriesFast();
1307 
1308  AliDmesonJetInfo DmesonJet;
1309 
1310  Int_t nAccCharm[3] = {0};
1311  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1312  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1313  if (!charmCand) continue;
1314  if(!(vHF.FillRecoCand(fAodEvent,charmCand))) continue;
1315 
1316  // region of interest + cuts
1317  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1318  Int_t nMassHypo = 0; // number of mass hypothesis accepted for this D meson
1319  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1320  DmesonJet.Reset();
1321  DmesonJet.fDmesonParticle = charmCand;
1322  DmesonJet.fSelectionType = im + 1;
1323  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1324  for (auto& def : fJetDefinitions) {
1325  if (!FindJet(charmCand, DmesonJet, def)) {
1326  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1327  def.GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1328  }
1329  }
1330  fDmesonJets[(icharm+1)*(1-(im*2))] = DmesonJet;
1331  nMassHypo++;
1332  nAccCharm[im]++;
1333  }
1334  }
1335  if (nMassHypo == 2) {
1336  nAccCharm[0]--;
1337  nAccCharm[1]--;
1338  nAccCharm[2] += 2;
1339  }
1340  if (nMassHypo == 2) { // both mass hypothesis accepted
1341  fDmesonJets[(icharm+1)].fSelectionType = 3;
1342  fDmesonJets[-(icharm+1)].fSelectionType = 3;
1343  }
1344  } // end of D cand loop
1345 
1346  TString hname;
1347 
1348  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1349  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1350  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1351  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1352 
1353  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1354  Int_t ntracks = 0;
1355  for (auto track_cont : fTrackContainers) ntracks += track_cont->GetNAcceptedTracks();
1356  fHistManager->FillTH2(hname, ntracks, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1357 
1358  hname = TString::Format("%s/fHistNDmesons", GetName());
1359  fHistManager->FillTH1(hname, nD);
1360 }
1361 
1372 {
1373  TString hname;
1374 
1376  fFastJetWrapper->SetR(jetDef.fRadius);
1379 
1380  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1381 
1382  if (jetDef.fJetType != AliJetContainer::kNeutralJet) {
1383  for (auto track_cont : fTrackContainers) {
1384  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
1385  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(Dcand);
1386  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1387  AddInputVectors(track_cont, 100, static_cast<TH2*>(fHistManager->FindObject(hname)), fTrackEfficiency);
1388 
1389  if (hftrack_cont) {
1390  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1391  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1392  const TObjArray& daughters = hftrack_cont->GetDaughterList();
1393  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1394  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1395  if (!hftrack_cont->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1396  }
1397  }
1398  }
1399  }
1400 
1401  if (jetDef.fJetType != AliJetContainer::kChargedJet) {
1402  for (auto clus_cont : fClusterContainers) {
1403  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1404  AddInputVectors(clus_cont, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1405  }
1406  }
1407 
1408  // run jet finder
1409  fFastJetWrapper->Run();
1410 
1411  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1412 
1413  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1414  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1415 
1416  Bool_t isDmesonJet = kFALSE;
1417 
1418  Double_t maxChPt = 0;
1419  Double_t maxNePt = 0;
1420  Double_t totalNeutralPt = 0;
1421 
1422  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1423  if (constituents[ic].user_index() == 0) {
1424  isDmesonJet = kTRUE;
1425  }
1426  else if (constituents[ic].user_index() >= 100) {
1427  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1428  }
1429  else if (constituents[ic].user_index() <= -100) {
1430  totalNeutralPt += constituents[ic].pt();
1431  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1432  }
1433  }
1434 
1435  if (isDmesonJet) {
1436  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1437  DmesonJet.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1438  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1439  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1440  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1441 
1442  return kTRUE;
1443  }
1444  }
1445 
1446  return kFALSE;
1447 }
1448 
1452 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist, Double_t eff)
1453 {
1454  auto itcont = cont->all_momentum();
1455  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1456  UInt_t rejectionReason = 0;
1457  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1458  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1459  continue;
1460  }
1461  if (fRandomGen && eff > 0 && eff < 1) {
1462  Double_t rnd = fRandomGen->Rndm();
1463  if (eff < rnd) {
1464  if (rejectHist) rejectHist->Fill(6, it->first.Pt());
1465  continue;
1466  }
1467  }
1468  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1469  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1470  }
1471 }
1472 
1475 {
1476  TString hname;
1477 
1478  fMCContainer->SetSpecialPDG(fCandidatePDG);
1479  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1480  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1481 
1482  if (!fMCContainer->IsSpecialPDGFound()) return;
1483 
1484  Int_t nAccCharm[3] = {0};
1485 
1486  for (auto &jetDef : fJetDefinitions) {
1487  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", GetName(), jetDef.GetName());
1488  TH1* histNDmesonsVsNconstituents = static_cast<TH1*>(fHistManager->FindObject(hname));
1489 
1490  switch (jetDef.fJetType) {
1492  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
1493  break;
1495  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
1496  break;
1498  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
1499  break;
1500  }
1501 
1503  fFastJetWrapper->SetR(jetDef.fRadius);
1506 
1507  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1508  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1509 
1510  fFastJetWrapper->Run();
1511 
1512  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1513 
1514  for (auto jet : jets_incl) {
1515  Int_t nDmesonsInJet = 0;
1516 
1517  for (auto constituent : jet.constituents()) {
1518  Int_t iPart = constituent.user_index() - 100;
1519  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
1520  if (!part) {
1521  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1522  continue;
1523  }
1524  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1525  nDmesonsInJet++;
1526  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1527  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1528  std::pair<int, AliDmesonJetInfo> element;
1529  element.first = iPart;
1530  dMesonJetIt = fDmesonJets.insert(element).first;
1531  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1532  (*dMesonJetIt).second.fDmesonParticle = part;
1533  (*dMesonJetIt).second.fSelectionType = part->PdgCode() > 0 ? 1 : 2;
1534 
1535  UShort_t p = 0;
1536  UInt_t rs = 0;
1537 
1538  auto firstParton = CheckOrigin(part, fMCContainer->GetArray(), kTRUE);
1539  p = 0;
1540  rs = firstParton.first;
1541  while (rs >>= 1) { p++; }
1542  (*dMesonJetIt).second.fFirstPartonType = p;
1543  (*dMesonJetIt).second.fFirstParton = firstParton.second;
1544 
1545  auto lastParton = CheckOrigin(part, fMCContainer->GetArray(), kFALSE);
1546  p = 0;
1547  rs = lastParton.first;
1548  while (rs >>= 1) { p++; }
1549  (*dMesonJetIt).second.fLastPartonType = p;
1550  (*dMesonJetIt).second.fLastParton = lastParton.second;
1551 
1552  if (part->PdgCode() > 0) {
1553  nAccCharm[0]++;
1554  }
1555  else {
1556  nAccCharm[1]++;
1557  }
1558  }
1559 
1560  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
1561  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = jet.constituents().size();
1562  } // if constituent is a D meson
1563  } // for each constituent
1564  if (nDmesonsInJet > 0) histNDmesonsVsNconstituents->Fill(jet.constituents().size(), nDmesonsInJet);
1565  } // for each jet
1566  } // for each jet definition
1567 
1568  if (fDmesonJets.size() != nAccCharm[0]+nAccCharm[1]) AliError(Form("I found %lu mesons (%d)?", fDmesonJets.size(), nAccCharm[0]+nAccCharm[1]));
1569  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1570  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1571  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1572  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1573 
1574  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1575  fHistManager->FillTH2(hname, fMCContainer->GetNAcceptedParticles(), nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1576 
1577  hname = TString::Format("%s/fHistNDmesons", GetName());
1578  fHistManager->FillTH1(hname, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]); // same as the number of accepted D mesons, since no selection is performed
1579 }
1580 
1585 {
1586  TString classname;
1587  if (fMCMode == kMCTruth) {
1588  classname = "AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary";
1589  fCurrentDmesonJetInfo = new AliDmesonMCInfoSummary();
1590  }
1591  else {
1592  switch (fCandidateType) {
1593  case kD0toKpi:
1594  case kD0toKpiLikeSign:
1595  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1596  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1597  break;
1598  case kDstartoKpipi:
1599  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1600  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1601  break;
1602  }
1603  }
1604  TString treeName = TString::Format("%s_%s", taskName, GetName());
1605  fTree = new TTree(treeName, treeName);
1606  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
1607  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1608  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1609  fCurrentJetInfo[i] = new AliJetInfoSummary();
1610  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1611  }
1612 
1613  return fTree;
1614 }
1615 
1620 {
1621  TString hname;
1622 
1623  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
1624 
1625  for (auto &jetDef : fJetDefinitions) {
1626 
1627  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
1628 
1629  Double_t radius = jetDef.fRadius;
1630 
1631  TString title[30] = {""};
1632  Int_t nbins[30] = {0 };
1633  Double_t min [30] = {0.};
1634  Double_t max [30] = {0.};
1635  Int_t dim = 0 ;
1636 
1637  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1638  nbins[dim] = nPtBins;
1639  min[dim] = 0;
1640  max[dim] = fMaxPt;
1641  dim++;
1642 
1643  if ((enabledAxis & kPositionD) != 0) {
1644  title[dim] = "#eta_{D}";
1645  nbins[dim] = 50;
1646  min[dim] = -1;
1647  max[dim] = 1;
1648  dim++;
1649 
1650  title[dim] = "#phi_{D} (rad)";
1651  nbins[dim] = 150;
1652  min[dim] = 0;
1653  max[dim] = TMath::TwoPi();
1654  dim++;
1655  }
1656 
1657  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1658  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1659  nbins[dim] = fNMassBins;
1660  min[dim] = fMinMass;
1661  max[dim] = fMaxMass;
1662  dim++;
1663  }
1664 
1665  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1666  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1667  nbins[dim] = fNMassBins;
1668  min[dim] = fMinMass;
1669  max[dim] = fMaxMass;
1670  dim++;
1671  }
1672 
1673  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1674  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1675  nbins[dim] = fNMassBins;
1676  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1677  dim++;
1678  }
1679 
1680  if (fCandidateType == kDstartoKpipi) {
1681  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1682  nbins[dim] = fNMassBins*6;
1683  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1684 
1685  // subtract mass of D0
1686  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1687  min[dim] -= D0mass;
1688  max[dim] -= D0mass;
1689 
1690  dim++;
1691  }
1692 
1693  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1694  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1695  nbins[dim] = 100;
1696  min[dim] = 0;
1697  max[dim] = 25;
1698  dim++;
1699  }
1700 
1701  title[dim] = "#it{z}_{D}";
1702  nbins[dim] = 110;
1703  min[dim] = 0;
1704  max[dim] = 1.10;
1705  dim++;
1706 
1707  if ((enabledAxis & kDeltaR) != 0) {
1708  title[dim] = "#Delta R_{D-jet}";
1709  nbins[dim] = 100;
1710  min[dim] = 0;
1711  max[dim] = radius * 1.5;
1712  dim++;
1713  }
1714 
1715  if ((enabledAxis & kDeltaEta) != 0) {
1716  title[dim] = "#eta_{D} - #eta_{jet}";
1717  nbins[dim] = 100;
1718  min[dim] = -radius * 1.2;
1719  max[dim] = radius * 1.2;
1720  dim++;
1721  }
1722 
1723  if ((enabledAxis & kDeltaPhi) != 0) {
1724  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1725  nbins[dim] = 100;
1726  min[dim] = -radius * 1.2;
1727  max[dim] = radius * 1.2;
1728  dim++;
1729  }
1730 
1731  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1732  nbins[dim] = nPtBins;
1733  min[dim] = 0;
1734  max[dim] = fMaxPt;
1735  dim++;
1736 
1737  if ((enabledAxis & kPositionJet) != 0) {
1738  title[dim] = "#eta_{jet}";
1739  nbins[dim] = 50;
1740  min[dim] = -1;
1741  max[dim] = 1;
1742  dim++;
1743 
1744  title[dim] = "#phi_{jet} (rad)";
1745  nbins[dim] = 150;
1746  min[dim] = 0;
1747  max[dim] = TMath::TwoPi();
1748  dim++;
1749  }
1750 
1751  if ((enabledAxis & kJetConstituents) != 0) {
1752  title[dim] = "No. of constituents";
1753  nbins[dim] = 50;
1754  min[dim] = -0.5;
1755  max[dim] = 49.5;
1756  dim++;
1757  }
1758 
1759  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1760  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1761  for (Int_t j = 0; j < dim; j++) {
1762  h->GetAxis(j)->SetTitle(title[j]);
1763  }
1764  }
1765 }
1766 
1771 {
1772  TString hname;
1773  fFirstPartons.clear();
1774  fLastPartons.clear();
1775  for (auto& dmeson_pair : fDmesonJets) {
1776  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
1777  Int_t accJets = 0;
1778  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1779  fCurrentJetInfo[ij]->Reset();
1780  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1781  if (!jet) continue;
1782  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1783  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1784  fHistManager->FillTH1(hname, jet->Pt());
1785  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1786  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1787  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1788  fHistManager->FillTH1(hname, jet->Eta());
1789  continue;
1790  }
1791  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
1792  accJets++;
1793  }
1794  if (accJets > 0) {
1795  fFirstPartons[dmeson_pair.second.fFirstParton] = dmeson_pair.second.fFirstPartonType;
1796  fLastPartons[dmeson_pair.second.fLastParton] = dmeson_pair.second.fLastPartonType;
1797 
1798  fTree->Fill();
1799  }
1800  else {
1801  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1802  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1803  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1804  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1805  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1806  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1807  if (fMCMode != kMCTruth) {
1808  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1809  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1810  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1811  }
1812  else if (fCandidateType == kDstartoKpipi) {
1813  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1814  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1815 
1816  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1817  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1818  }
1819  }
1820  }
1821  }
1822 
1823  hname = TString::Format("%s/fHistFirstPartonPt", GetName());
1824  TH1* histFirstPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
1825  hname = TString::Format("%s/fHistFirstPartonEta", GetName());
1826  TH1* histFirstPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
1827  hname = TString::Format("%s/fHistFirstPartonPhi", GetName());
1828  TH1* histFirstPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
1829  hname = TString::Format("%s/fHistFirstPartonType", GetName());
1830  TH1* histFirstPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
1831 
1832  for (auto parton : fFirstPartons) {
1833  if (!parton.first) continue;
1834  histFirstPartonPt->Fill(parton.first->Pt());
1835  histFirstPartonEta->Fill(parton.first->Eta());
1836  histFirstPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
1837  histFirstPartonType->Fill(parton.second);
1838  }
1839 
1840  hname = TString::Format("%s/fHistLastPartonPt", GetName());
1841  TH1* histLastPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
1842  hname = TString::Format("%s/fHistLastPartonEta", GetName());
1843  TH1* histLastPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
1844  hname = TString::Format("%s/fHistLastPartonPhi", GetName());
1845  TH1* histLastPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
1846  hname = TString::Format("%s/fHistLastPartonType", GetName());
1847  TH1* histLastPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
1848 
1849  for (auto parton : fLastPartons) {
1850  if (!parton.first) continue;
1851  histLastPartonPt->Fill(parton.first->Pt());
1852  histLastPartonEta->Fill(parton.first->Eta());
1853  histLastPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
1854  histLastPartonType->Fill(parton.second);
1855  }
1856 
1857  return kTRUE;
1858 }
1859 
1865 {
1866  TString hname;
1867  fFirstPartons.clear();
1868  fLastPartons.clear();
1869  for (auto& dmeson_pair : fDmesonJets) {
1870  Int_t accJets = 0;
1871  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1872  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1873  if (!jet) continue;
1874  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1875  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1876  fHistManager->FillTH1(hname, jet->Pt());
1877  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1878  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1879  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1880  fHistManager->FillTH1(hname, jet->Eta());
1881  continue;
1882  }
1883  accJets++;
1884  }
1885  if (accJets > 0) {
1886  fFirstPartons[dmeson_pair.second.fFirstParton] = dmeson_pair.second.fFirstPartonType;
1887  fLastPartons[dmeson_pair.second.fLastParton] = dmeson_pair.second.fLastPartonType;
1888  }
1889  else {
1890  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1891  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1892  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1893  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1894  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1895  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1896  if (fMCMode != kMCTruth) {
1897  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1898  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1899  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1900  }
1901  else if (fCandidateType == kDstartoKpipi) {
1902  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1903  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1904 
1905  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1906  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1907  }
1908  }
1909  }
1910  }
1911 
1912  hname = TString::Format("%s/fHistFirstPartonPt", GetName());
1913  TH1* histFirstPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
1914  hname = TString::Format("%s/fHistFirstPartonEta", GetName());
1915  TH1* histFirstPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
1916  hname = TString::Format("%s/fHistFirstPartonPhi", GetName());
1917  TH1* histFirstPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
1918  hname = TString::Format("%s/fHistFirstPartonType", GetName());
1919  TH1* histFirstPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
1920 
1921  for (auto parton : fFirstPartons) {
1922  if (!parton.first) continue;
1923  histFirstPartonPt->Fill(parton.first->Pt());
1924  histFirstPartonEta->Fill(parton.first->Eta());
1925  histFirstPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
1926  histFirstPartonType->Fill(parton.second);
1927  }
1928 
1929  hname = TString::Format("%s/fHistLastPartonPt", GetName());
1930  TH1* histLastPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
1931  hname = TString::Format("%s/fHistLastPartonEta", GetName());
1932  TH1* histLastPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
1933  hname = TString::Format("%s/fHistLastPartonPhi", GetName());
1934  TH1* histLastPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
1935  hname = TString::Format("%s/fHistLastPartonType", GetName());
1936  TH1* histLastPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
1937 
1938  for (auto parton : fLastPartons) {
1939  if (!parton.first) continue;
1940  histLastPartonPt->Fill(parton.first->Pt());
1941  histLastPartonEta->Fill(parton.first->Eta());
1942  histLastPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
1943  histLastPartonType->Fill(parton.second);
1944  }
1945 
1946  return kTRUE;
1947 }
1948 
1953 {
1954  TString hname;
1955 
1956  for (auto& dmeson_pair : fDmesonJets) {
1957  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
1958  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1959  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1960  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1961  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1962  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1963  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1964  }
1965  }
1966 
1967  for (auto &jetDef : fJetDefinitions) {
1968 
1969  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1970  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
1971 
1972  for (auto& dmeson_pair : fDmesonJets) {
1973  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
1974  if (!jet) continue;
1975  if (!jetDef.IsJetInAcceptance(*jet)) {
1976  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
1977  fHistManager->FillTH1(hname, jet->Pt());
1978  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
1979  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1980  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
1981  fHistManager->FillTH1(hname, jet->Eta());
1982  continue;
1983  }
1984  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
1985  }
1986  }
1987 
1988  return kTRUE;
1989 }
1990 
1997 {
1998  // Fill the THnSparse histogram.
1999 
2000  Double_t contents[30] = {0.};
2001 
2002  Double_t z = DmesonJet.GetZ(n);
2003  Double_t deltaPhi = 0;
2004  Double_t deltaEta = 0;
2005  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
2006 
2007  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
2008  if (it == DmesonJet.fJets.end()) return kFALSE;
2009 
2010  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
2011  TString title(h->GetAxis(i)->GetTitle());
2012  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
2013  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
2014  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
2015  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
2016  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
2017  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
2018  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
2019  else if (title=="#it{z}_{D}") contents[i] = z;
2020  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
2021  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
2022  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
2023  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
2024  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
2025  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
2026  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
2027  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
2028  }
2029 
2030  h->Fill(contents);
2031 
2032  return kTRUE;
2033 }
2034 
2035 // Definitions of class AliAnalysisTaskDmesonJets
2036 
2040 
2044  fAnalysisEngines(),
2045  fEnabledAxis(0),
2047  fHistManager(),
2048  fApplyKinematicCuts(kTRUE),
2049  fNOutputTrees(0),
2050  fTrackEfficiency(0),
2051  fMCContainer(0),
2052  fAodEvent(0),
2053  fFastJetWrapper(0)
2054 {
2055  SetMakeGeneralHistograms(kTRUE);
2056 }
2057 
2062  AliAnalysisTaskEmcalLight(name, kTRUE),
2063  fAnalysisEngines(),
2064  fEnabledAxis(k2ProngInvMass),
2065  fOutputType(kTreeOutput),
2066  fHistManager(name),
2067  fApplyKinematicCuts(kTRUE),
2068  fNOutputTrees(nOutputTrees),
2069  fTrackEfficiency(0),
2070  fMCContainer(0),
2071  fAodEvent(0),
2072  fFastJetWrapper(0)
2073 {
2074  SetMakeGeneralHistograms(kTRUE);
2075  for (Int_t i = 0; i < nOutputTrees; i++){
2076  DefineOutput(2+i, TTree::Class());
2077  }
2078 }
2079 
2082 {
2083  if (fFastJetWrapper) delete fFastJetWrapper;
2084 }
2085 
2093 {
2094  AliRDHFCuts* analysiscuts = 0;
2095  TFile* filecuts = TFile::Open(cutfname);
2096  if (!filecuts || filecuts->IsZombie()) {
2097  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Input file not found: will use std cuts.");
2098  filecuts = 0;
2099  }
2100 
2101  if (filecuts) analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
2102 
2103  if (!analysiscuts) {
2104  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Could not find analysis cuts '%s' in '%s'.", cutsname.Data(), cutfname.Data());
2105  if (filecuts) {
2106  filecuts->ls();
2107  }
2108  }
2109  else {
2110  ::Info("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Cuts '%s' loaded from file '%s'", cutsname.Data(), cutfname.Data());
2111  }
2112 
2113  return analysiscuts;
2114 }
2115 
2125 {
2127  return AddAnalysisEngine(type, MCmode, jetDef, cutfname);
2128 }
2129 
2139 {
2140  AliRDHFCuts* cuts = 0;
2141 
2142  if (!cutfname.IsNull()) {
2143  TString cutsname;
2144 
2145  switch (type) {
2146  case kD0toKpi:
2147  case kD0toKpiLikeSign:
2148  cutsname = "D0toKpiCuts";
2149  break;
2150  case kDstartoKpipi:
2151  cutsname = "DStartoKpipiCuts";
2152  break;
2153  default:
2154  return 0;
2155  }
2156 
2157  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
2158  if (cuts) cuts->PrintAll();
2159  }
2160 
2161  AnalysisEngine eng(type, MCmode, cuts);
2162 
2163  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
2164 
2165  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
2166  eng.AddJetDefinition(jetDef);
2167  it = fAnalysisEngines.insert(it, eng);
2168  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
2169  }
2170  else {
2171  AnalysisEngine* found_eng = &(*it);
2172  ::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());
2173  found_eng->AddJetDefinition(jetDef);
2174 
2175  if (cuts) {
2176  if (found_eng->fRDHFCuts != 0) ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
2177  found_eng->SetRDHFCuts(cuts);
2178  }
2179  }
2180 
2181  return &(*it);
2182 }
2183 
2184 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
2185 {
2186  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
2187  while (it != fAnalysisEngines.end() && (*it) < eng) it++;
2188  return it;
2189 }
2190 
2193 {
2194  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
2195 
2197 
2198  // Define histograms
2199  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
2200 
2201  TString hname;
2202  TString htitle;
2203  TH1* h = 0;
2204  Int_t treeSlot = 0;
2205 
2206  hname = "fHistCharmPt";
2207  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2208  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2209 
2210  hname = "fHistCharmEta";
2211  htitle = hname + ";#eta_{charm};counts";
2212  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2213 
2214  hname = "fHistCharmPhi";
2215  htitle = hname + ";#phi_{charm};counts";
2216  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2217 
2218  hname = "fHistCharmPt_Eta05";
2219  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2220  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2221 
2222  hname = "fHistBottomPt";
2223  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2224  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2225 
2226  hname = "fHistBottomEta";
2227  htitle = hname + ";#eta_{bottom};counts";
2228  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2229 
2230  hname = "fHistBottomPhi";
2231  htitle = hname + ";#phi_{bottom};counts";
2232  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2233 
2234  hname = "fHistBottomPt_Eta05";
2235  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2236  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2237 
2238  hname = "fHistHighestPartonPt";
2239  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2240  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2241 
2242  hname = "fHistHighestPartonType";
2243  htitle = hname + ";type;counts";
2244  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2245 
2246  hname = "fHistNHeavyQuarks";
2247  htitle = hname + ";number of heavy-quarks;counts";
2248  fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
2249 
2250  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
2251  for (auto &param : fAnalysisEngines) {
2252  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
2253 
2254  param.fHistManager = &fHistManager;
2255 
2256  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", param.GetName());
2257  htitle = hname + ";#it{N}_{tracks};#it{N}_{D};events";
2258  h = fHistManager.CreateTH2(hname, htitle, 251, -0.5, 250.5, 21, -0.5, 20.5);
2259 
2260  hname = TString::Format("%s/fHistNTotAcceptedDmesons", param.GetName());
2261  htitle = hname + ";;#it{N}_{D}";
2262  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
2263 
2264  hname = TString::Format("%s/fHistNDmesons", param.GetName());
2265  htitle = hname + ";#it{N}_{D};events";
2266  h = fHistManager.CreateTH1(hname, htitle, 501, -0.5, 500.5);
2267 
2268  hname = TString::Format("%s/fHistNEvents", param.GetName());
2269  htitle = hname + ";Event status;counts";
2270  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
2271  h->GetXaxis()->SetBinLabel(1, "Accepted");
2272  h->GetXaxis()->SetBinLabel(2, "Rejected");
2273 
2274  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
2275  htitle = hname + ";Rejection reason;counts";
2276  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
2277 
2278  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
2279  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
2280  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2281 
2282  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
2283  htitle = hname + ";#it{#eta}_{D};counts";
2284  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2285 
2286  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
2287  htitle = hname + ";#it{#phi}_{D};counts";
2288  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2289 
2290  if (param.fMCMode != kMCTruth) {
2291  if (param.fCandidateType == kD0toKpi || param.fCandidateType == kD0toKpiLikeSign) {
2292  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
2293  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2294  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
2295  }
2296  else if (param.fCandidateType == kDstartoKpipi) {
2297  Double_t min = 0;
2298  Double_t max = 0;
2299 
2300  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
2301  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2302  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
2303  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
2304 
2305  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2306  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
2307  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2308  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
2309  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
2310  }
2311  }
2312 
2313  if (param.fMCMode == kMCTruth) {
2314  hname = TString::Format("%s/fHistFirstPartonPt", param.GetName());
2315  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
2316  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2317 
2318  hname = TString::Format("%s/fHistFirstPartonEta", param.GetName());
2319  htitle = hname + ";#eta_{parton};counts";
2320  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2321 
2322  hname = TString::Format("%s/fHistFirstPartonPhi", param.GetName());
2323  htitle = hname + ";#phi_{parton};counts";
2324  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2325 
2326  hname = TString::Format("%s/fHistFirstPartonType", param.GetName());
2327  htitle = hname + ";type;counts";
2328  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2329 
2330  hname = TString::Format("%s/fHistLastPartonPt", param.GetName());
2331  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
2332  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2333 
2334  hname = TString::Format("%s/fHistLastPartonEta", param.GetName());
2335  htitle = hname + ";#eta_{parton};counts";
2336  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2337 
2338  hname = TString::Format("%s/fHistLastPartonPhi", param.GetName());
2339  htitle = hname + ";#phi_{parton};counts";
2340  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2341 
2342  hname = TString::Format("%s/fHistLastPartonType", param.GetName());
2343  htitle = hname + ";type;counts";
2344  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2345  }
2346 
2347  for (auto& jetDef : param.fJetDefinitions) {
2348  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef.GetName());
2349 
2350  if (param.fMCMode == kMCTruth) {
2351  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", param.GetName(), jetDef.GetName());
2352  htitle = hname + ";#it{N}_{constituents};#it{N}_{D};counts";
2353  h = fHistManager.CreateTH2(hname, htitle, 51, -0.5, 50.5, 10, 0.5, 10.5);
2354  }
2355 
2356  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef.GetName());
2357  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2358  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2359  SetRejectionReasonLabels(h->GetXaxis());
2360 
2361  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef.GetName());
2362  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2363  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2364  SetRejectionReasonLabels(h->GetXaxis());
2365 
2366  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef.GetName());
2367  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
2368  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2369  SetRejectionReasonLabels(h->GetXaxis());
2370 
2371  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef.GetName());
2372  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
2373  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
2374 
2375  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef.GetName());
2376  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
2377  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2378 
2379  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef.GetName());
2380  htitle = hname + ";#it{#eta}_{jet};counts";
2381  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2382 
2383  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef.GetName());
2384  htitle = hname + ";#it{#phi}_{jet};counts";
2385  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2386  }
2387  switch (fOutputType) {
2388  case kTreeOutput:
2389  param.BuildTree(GetName());
2390  if (treeSlot < fNOutputTrees) {
2391  param.AssignDataSlot(treeSlot+2);
2392  treeSlot++;
2394  }
2395  else {
2396  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
2397  }
2398  break;
2399  case kTHnOutput:
2400  param.BuildHnSparse(fEnabledAxis);
2401  break;
2402  case kNoOutput:
2403  break;
2404  }
2405  }
2406 
2408 
2409  PostData(1, fOutput);
2410 }
2411 
2415 {
2417 
2418  // Load the event
2419  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
2420 
2421  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
2422 
2423  fFastJetWrapper->SetAreaType(fastjet::active_area);
2425 
2426  if (!fAodEvent) {
2427  AliError(Form("This task need an AOD event (Task '%s'). Expect troubles...", GetName()));
2428  //return;
2429  }
2430 
2431  TRandom* rnd = 0;
2432  if (fTrackEfficiency > 0 && fTrackEfficiency < 1) rnd = new TRandom3(0);
2433 
2434  for (auto cont_it : fParticleCollArray) {
2435  AliHFAODMCParticleContainer* part_cont = dynamic_cast<AliHFAODMCParticleContainer*>(cont_it.second);
2436  if (part_cont) fMCContainer = part_cont;
2437  }
2438 
2439  for (auto &params : fAnalysisEngines) {
2440 
2441  params.fAodEvent = fAodEvent;
2442  params.fFastJetWrapper = fFastJetWrapper;
2443  params.fTrackEfficiency = fTrackEfficiency;
2444  params.fRandomGen = rnd;
2445 
2446  if (!params.fRDHFCuts) {
2447  if (params.fMCMode == kMCTruth) {
2448  ::Warning("AliAnalysisTaskDmesonJets::ExecOnce",
2449  "%s: RDHF cuts not provided for engine '%s'.",
2450  GetName(), params.GetName());
2451  }
2452  else {
2453  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2454  "%s: RDHF cuts not provided. Engine '%s' disabled.",
2455  GetName(), params.GetName());
2456  params.fInhibit = kTRUE;
2457  }
2458  }
2459 
2460  params.fMCContainer = fMCContainer;
2461 
2462  for (auto cont_it : fParticleCollArray) {
2463  AliTrackContainer* track_cont = dynamic_cast<AliTrackContainer*>(cont_it.second);
2464  if (track_cont) params.fTrackContainers.push_back(track_cont);
2465  }
2466 
2467  for (auto cont_it : fClusterCollArray) params.fClusterContainers.push_back(cont_it.second);
2468 
2469  if (fAodEvent) params.Init(fGeom, fAodEvent->GetRunNumber());
2470 
2471  if (params.fMCMode != kMCTruth && fAodEvent) {
2472  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
2473 
2474  if (params.fCandidateArray) {
2475  TString className;
2476  if (params.fCandidateType == kD0toKpi || params.fCandidateType == kD0toKpiLikeSign) {
2477  className = "AliAODRecoDecayHF2Prong";
2478  }
2479  else if (params.fCandidateType == kDstartoKpipi) {
2480  className = "AliAODRecoCascadeHF";
2481  }
2482  if (!params.fCandidateArray->GetClass()->InheritsFrom(className)) {
2483  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2484  "%s: Objects of type %s in %s are not inherited from %s! Task will be disabled!",
2485  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName(), className.Data());
2486  params.fCandidateArray = 0;
2487  params.fInhibit = kTRUE;
2488  }
2489  }
2490  else {
2491  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2492  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
2493  params.fBranchName.Data(), params.GetName());
2494  params.fInhibit = kTRUE;
2495  }
2496  }
2497 
2498  if (params.fMCMode != kNoMC) {
2499  if (!params.fMCContainer) {
2500  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2501  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
2502  params.GetName());
2503  params.fInhibit = kTRUE;
2504  }
2505  }
2506 
2507  if (params.fMCMode != kMCTruth) {
2508  if (params.fTrackContainers.size() == 0 && params.fClusterContainers.size() == 0) {
2509  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2510  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
2511  params.GetName());
2512  params.fInhibit = kTRUE;
2513  }
2514  }
2515  }
2516 }
2517 
2522 {
2523  TString hname;
2524 
2525  // Fix for temporary bug in ESDfilter
2526  // The AODs with null vertex pointer didn't pass the PhysSel
2527  // Now adding an entry in the histogram so as to check that this is actually cutting anything out
2528  if (fAodEvent && (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001)) {
2529  for (auto &eng : fAnalysisEngines) {
2530  if (eng.fInhibit) continue;
2531  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2532  fHistManager.FillTH1(hname, "ESDfilterBug");
2533  }
2534  return kFALSE;
2535  }
2536 
2537  if (fAodEvent) {
2538  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
2539  if (matchingAODdeltaAODlevel <= 0) {
2540  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
2541  for (auto &eng : fAnalysisEngines) {
2542  if (eng.fInhibit) continue;
2543  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2544  fHistManager.FillTH1(hname, "MismatchDeltaAOD");
2545  }
2546  return kFALSE;
2547  }
2548  }
2549 
2550  for (auto &eng : fAnalysisEngines) {
2551  eng.fDmesonJets.clear();
2552  if (eng.fInhibit) continue;
2553 
2554  //Event selection
2555  hname = TString::Format("%s/fHistNEvents", eng.GetName());
2556  if (fAodEvent) {
2557  Bool_t iseventselected = kTRUE;
2558  if (eng.fRDHFCuts) iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
2559  if (!iseventselected) {
2560  fHistManager.FillTH1(hname, "Rejected");
2561  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2562  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
2563  TString label;
2564  do {
2565  label = GetHFEventRejectionReasonLabel(bitmap);
2566  if (label.IsNull()) break;
2567  fHistManager.FillTH1(hname, label);
2568  } while (true);
2569 
2570  continue;
2571  }
2572  }
2573 
2574  fHistManager.FillTH1(hname, "Accepted");
2575 
2576  AliDebug(2, "Event selected");
2577 
2578  eng.RunAnalysis();
2579  }
2580  return kTRUE;
2581 }
2582 
2587 {
2588  for (auto &param : fAnalysisEngines) {
2589  if (param.fInhibit) continue;
2590 
2591  if (fOutputType == kTreeOutput) {
2592  param.FillTree(fApplyKinematicCuts);
2594  }
2595  else if (fOutputType == kTHnOutput) {
2596  param.FillHnSparse(fApplyKinematicCuts);
2597  }
2598  }
2600  return kTRUE;
2601 }
2602 
2605 {
2606  auto itcont = fMCContainer->all_momentum();
2607  Int_t nHQ = 0;
2608  Double_t highestPartonPt = 0;
2609  Int_t absPdgHighParton = 0;
2610  for (auto part : itcont) {
2611  Int_t absPdgCode = TMath::Abs(part.second->GetPdgCode());
2612 
2613  // Skip all particles that are not either quarks or gluons
2614  if (absPdgCode > 9 && absPdgCode != 21) continue;
2615 
2616  // Look for highest momentum parton
2617  if (highestPartonPt < part.first.Pt()) {
2618  highestPartonPt = part.first.Pt();
2619  absPdgHighParton = absPdgCode;
2620  }
2621  /*
2622  // Look for the mother PDG code
2623  Int_t motherIndex = part.second->GetMother();
2624  AliAODMCParticle *mother = 0;
2625  Int_t motherPdg = 0;
2626  Double_t motherPt = 0;
2627  if (motherIndex >= 0) {
2628  mother = fMCContainer->GetMCParticle(motherIndex);
2629  if (motherIndex) {
2630  motherPdg = TMath::Abs(mother->GetPdgCode());
2631  motherPt = mother->Pt();
2632  }
2633  }
2634  */
2635  if (absPdgCode != 4 && absPdgCode != 5) continue;
2636  Bool_t notLastInPartonShower = kFALSE;
2637  for (Int_t idaugh = 0; idaugh < 2; idaugh++){
2638  Int_t daughterIndex = part.second->GetDaughter(idaugh);
2639  if (daughterIndex < 0) {
2640  AliDebug(10, Form("Could not find daughter of heavy quark (pdg=%d, pt=%.3f)!", absPdgCode, part.first.Pt()));
2641  continue;
2642  }
2643  AliAODMCParticle *daughter = fMCContainer->GetMCParticle(daughterIndex);
2644  if (!daughter) {
2645  AliDebug(10, Form("Could not find daughter %d of heavy quark (pdg=%d, pt=%.3f)!", daughterIndex, absPdgCode, part.first.Pt()));
2646  continue;
2647  }
2648  Int_t daughterAbsPdgCode = TMath::Abs(daughter->GetPdgCode());
2649  if (daughterAbsPdgCode <= 9 || daughterAbsPdgCode == 21) notLastInPartonShower = kTRUE; // this parton is not the last parton in the shower
2650  AliDebug(10, Form("Found daughter with PDG=%d, pt=%.3f", daughterAbsPdgCode, daughter->Pt()));
2651  }
2652  if (notLastInPartonShower) continue;
2653 
2654  if (absPdgCode == 4) {
2655  fHistManager.FillTH1("fHistCharmPt", part.first.Pt());
2656  fHistManager.FillTH1("fHistCharmEta", part.first.Eta());
2657  fHistManager.FillTH1("fHistCharmPhi", part.first.Phi_0_2pi());
2658  if (TMath::Abs(part.first.Eta()) < 0.5) fHistManager.FillTH1("fHistCharmPt_Eta05", part.first.Pt());
2659  }
2660  else if (absPdgCode == 5) {
2661  fHistManager.FillTH1("fHistBottomPt", part.first.Pt());
2662  fHistManager.FillTH1("fHistBottomEta", part.first.Eta());
2663  fHistManager.FillTH1("fHistBottomPhi", part.first.Phi_0_2pi());
2664  if (TMath::Abs(part.first.Eta()) < 0.5) fHistManager.FillTH1("fHistBottomPt_Eta05", part.first.Pt());
2665  }
2666  nHQ++;
2667  }
2668  fHistManager.FillTH1("fHistNHeavyQuarks", nHQ);
2669  fHistManager.FillTH1("fHistHighestPartonPt",highestPartonPt);
2670  Int_t partonType = 0;
2671  if (absPdgHighParton == 9 || absPdgHighParton == 21) {
2672  partonType = 7;
2673  }
2674  else {
2675  partonType = absPdgHighParton;
2676  }
2677  fHistManager.FillTH1("fHistHighestPartonType",partonType);
2678 }
2679 
2688 {
2689  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
2690 
2691  Double_t mass = part->Mass();
2692 
2693  // To make sure that the PDG mass value is not at the edge of a bin
2694  if (nbins % 2 == 0) {
2695  minMass = mass - range / 2 - range / nbins / 2;
2696  maxMass = mass + range / 2 - range / nbins / 2;
2697  }
2698  else {
2699  minMass = mass - range / 2;
2700  maxMass = mass + range / 2;
2701  }
2702 }
2703 
2710 {
2711  static TString label;
2712  label = "";
2713 
2714  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
2715  label = "NotSelTrigger";
2716  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
2717  return label.Data();
2718  }
2719  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
2720  label = "NoVertex";
2721  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
2722  return label.Data();
2723  }
2724  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
2725  label = "TooFewVtxContrib";
2726  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
2727  return label.Data();
2728  }
2729  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
2730  label = "ZVtxOutFid";
2731  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
2732  return label.Data();
2733  }
2734  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
2735  label = "Pileup";
2736  bitmap &= ~BIT(AliRDHFCuts::kPileup);
2737  return label.Data();
2738  }
2739  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
2740  label = "OutsideCentrality";
2741  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
2742  return label.Data();
2743  }
2744  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
2745  label = "PhysicsSelection";
2746  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
2747  return label.Data();
2748  }
2749  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
2750  label = "BadSPDVertex";
2751  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
2752  return label.Data();
2753  }
2754  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
2755  label = "ZVtxSPDOutFid";
2756  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
2757  return label.Data();
2758  }
2759  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
2760  label = "CentralityFlattening";
2761  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
2762  return label.Data();
2763  }
2764  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
2765  label = "BadTrackV0Correl";
2766  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
2767  return label.Data();
2768  }
2769 
2770  return label.Data();
2771 }
2772 
2779 {
2780  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
2781  PostData(eng.GetDataSlotNumber(), eng.GetTree());
2782  return eng.GetDataSlotNumber();
2783  }
2784  else {
2785  return -1;
2786  }
2787 }
2788 
2798 {
2799  // Get the pointer to the existing analysis manager via the static access method
2800  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
2801  if (!mgr) {
2802  ::Error("AddTaskDmesonJets", "No analysis manager to connect to.");
2803  return NULL;
2804  }
2805 
2806  // Check the analysis type using the event handlers connected to the analysis manager
2807  AliVEventHandler* handler = mgr->GetInputEventHandler();
2808  if (!handler) {
2809  ::Error("AddTaskEmcalJetSpectraQA", "This task requires an input event handler");
2810  return NULL;
2811  }
2812 
2813  EDataType_t dataType = kUnknownDataType;
2814 
2815  if (handler->InheritsFrom("AliESDInputHandler")) {
2816  dataType = kESD;
2817  }
2818  else if (handler->InheritsFrom("AliAODInputHandler")) {
2819  dataType = kAOD;
2820  }
2821 
2822  // Init the task and do settings
2823  if (ntracks == "usedefault") {
2824  if (dataType == kESD) {
2825  ntracks = "Tracks";
2826  }
2827  else if (dataType == kAOD) {
2828  ntracks = "tracks";
2829  }
2830  else {
2831  ntracks = "";
2832  }
2833  }
2834 
2835  if (nclusters == "usedefault") {
2836  if (dataType == kESD) {
2837  nclusters = "CaloClusters";
2838  }
2839  else if (dataType == kAOD) {
2840  nclusters = "caloClusters";
2841  }
2842  else {
2843  nclusters = "";
2844  }
2845  }
2846 
2847  if (nMCpart == "usedefault") {
2848  nMCpart = "mcparticles"; // Always needs AliAODMCParticle objects
2849  }
2850 
2851  TString name("AliAnalysisTaskDmesonJets");
2852  if (strcmp(suffix, "") != 0) {
2853  name += TString::Format("_%s", suffix.Data());
2854  }
2855 
2856  AliAnalysisTaskDmesonJets* jetTask = new AliAnalysisTaskDmesonJets(name, nMaxTrees);
2857 
2858  if (!ntracks.IsNull()) {
2859  AliHFTrackContainer* trackCont = new AliHFTrackContainer(ntracks);
2860  jetTask->AdoptParticleContainer(trackCont);
2861  }
2862 
2863  if (!nMCpart.IsNull()) {
2864  AliMCParticleContainer* partCont = new AliHFAODMCParticleContainer(nMCpart);
2865  partCont->SetEtaLimits(-1.5, 1.5);
2866  partCont->SetPtLimits(0, 1000);
2867  jetTask->AdoptParticleContainer(partCont);
2868  }
2869 
2870  jetTask->AddClusterContainer(nclusters.Data());
2871 
2872  // Final settings, pass to manager and set the containers
2873  mgr->AddTask(jetTask);
2874 
2875  // Create containers for input/output
2876  AliAnalysisDataContainer* cinput1 = mgr->GetCommonInputContainer();
2877  TString contname1(name);
2878  contname1 += "_histos";
2879  AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(contname1.Data(),
2880  TList::Class(), AliAnalysisManager::kOutputContainer,
2881  Form("%s", AliAnalysisManager::GetCommonFileName()));
2882 
2883  mgr->ConnectInput(jetTask, 0, cinput1);
2884  mgr->ConnectOutput(jetTask, 1, coutput1);
2885 
2886  for (Int_t i = 0; i < nMaxTrees; i++) {
2887  TString contname = TString::Format("%s_tree_%d", name.Data(), i);
2888  AliAnalysisDataContainer *coutput = mgr->CreateContainer(contname.Data(),
2889  TTree::Class(),AliAnalysisManager::kOutputContainer,
2890  Form("%s", AliAnalysisManager::GetCommonFileName()));
2891  mgr->ConnectOutput(jetTask, 2+i, coutput);
2892  }
2893  return jetTask;
2894 }
2895 
Int_t pdg
AliClusterContainer * AddClusterContainer(std::string branchName, std::string contName="")
void Print() const
Prints the content of this object in the standard output.
std::list< AnalysisEngine > fAnalysisEngines
Array of analysis parameters.
AliEMCALGeometry * fGeom
!emcal geometry
void SetDMesonCandidate(AliAODRecoDecay *c)
Analysis task for D meson jets.
virtual void Set(const AliDmesonJetInfo &source)
Invariant mass of the D0 meson candidate in GeV/c2.
virtual Int_t Run()
double Double_t
Definition: External.C:58
Bool_t IsJetInAcceptance(const AliJetInfo &jet) const
Double_t fSoftPionPt
! Transverse momentum of the soft pion of the D* candidate
EDataType_t
Switch for the data type.
virtual AliAODMCParticle * GetMCParticle(Int_t i=-1) const
const char * title
Definition: MakeQAPdf.C:26
AliRDHFCuts * fRDHFCuts
D meson candidates cuts.
virtual void Set(const AliDmesonJetInfo &source)
Lightweight class that encapsulates D meson jets.
Double_t fTrackEfficiency
Artificial tracking inefficiency (0...1)
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
AliHFAODMCParticleContainer * fMCContainer
! MC particle container
ERecoScheme_t fRecoScheme
Jet recombination scheme (pt scheme, E scheme, ...)
const TObjArray & GetDaughterList() const
Double_t fMaxNeutralPt
Transverse momentum of the leading neutral particle (or cluster)
void RunDetectorLevelAnalysis()
Run a detector level analysis.
Double_t mass
Container with name, TClonesArray and cuts for particles.
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
virtual void UserCreateOutputObjects()
Creates the output containers.
void FillPartonLevelHistograms()
Fill histograms with parton-level information.
static Int_t CheckMatchingAODdeltaAODevents()
Lightweight class that encapsulates D meson jets.
Double_t fInvMass2Prong
! 2-prong mass of the D* candidate (w/o the soft pion)
static void CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t &minMass, Double_t &maxMass)
AliAODEvent * fAodEvent
! AOD event
Bool_t FindJet(AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, AliHFJetDefinition &jetDef)
UInt_t fEnabledAxis
Use bit defined in EAxis_t to enable axis in the THnSparse.
AliAODMCParticle * fFirstParton
! pointer to the first parton in the shower tree of the D meson (only for particle level D mesons) ...
EMCMode_t fMCMode
MC mode: No MC (data and MC detector level), background-only (MC), signal-only (MC), MC truth (particle level)
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
Short_t fLastPartonType
! type of the last parton in the shower tree (only for particle level D mesons)
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:103
AnalysisEngine()
This is the default constructor, used for ROOT I/O purposes.
void Init(const AliEMCALGeometry *const geom, Int_t runNumber)
Initialize the analysis engine.
static FJRecoScheme ConvertToFJRecoScheme(ERecoScheme_t reco)
AliVParticle * fDmesonParticle
! pointer to the particle object
static TString GenerateJetName(EJetType_t jetType, EJetAlgo_t jetAlgo, ERecoScheme_t recoScheme, Double_t radius, AliParticleContainer *partCont, AliClusterContainer *clusCont, TString tag)
void SetR(Double_t r)
Definition: AliFJWrapper.h:108
Select tracks based on specific prescriptions of HF analysis.
ECandidateType_t fCandidateType
Candidate type.
Double_t GetDistance(std::string n, Double_t &deta, Double_t &dphi) const
Int_t fMCLabel
! MC label, i.e. index of the generator level D meson (only for detector level D meson candidates) ...
const Int_t nPtBins
virtual void Set(const AliDmesonJetInfo &source)
void AddInputVectors(AliEmcalContainer *cont, Int_t offset, TH2 *rejectHist=0, Double_t eff=0.)
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:34
EJetType_t fJetType
Jet type (charged, full, neutral)
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
Create a new TH2 within the container.
virtual void Set(const AliDmesonJetInfo &source)
AliTLorentzVector fD
! 4-momentum of the D meson candidate
Double_t GetDistance(const AliJetInfo &jet, Double_t &deta, Double_t &dphi) const
TObject * FindObject(const char *name) const
Find an object inside the container.
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:658
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:102
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
AliAODTrack * GetBachelor() const
AliRDHFCuts * LoadDMesonCutsFromFile(TString cutfname, TString cutsname)
bool operator<(const AliAnalysisTaskDmesonJets::AliHFJetDefinition &lhs, const AliAnalysisTaskDmesonJets::AliHFJetDefinition &rhs)
std::list< AnalysisEngine >::iterator FindAnalysisEngine(const AnalysisEngine &eng)
std::vector< AliHFJetDefinition > fJetDefinitions
Jet definitions.
EJetAlgo_t fJetAlgo
Jet algorithm (kt, anti-kt,...)
static const char * GetHFEventRejectionReasonLabel(UInt_t &bitmap)
Lightweight class that encapsulates D*.
virtual void Clear(const Option_t *="")
Bool_t ExtractDstarAttributes(const AliAODRecoCascadeHF *DstarCand, AliDmesonJetInfo &DmesonJet, UInt_t i)
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
const AliJetInfo * GetJet(std::string n) const
Short_t fFirstPartonType
! type of the first parton in the shower tree (only for particle level D mesons)
const AliMCParticleIterableMomentumContainer all_momentum() const
virtual ~AliAnalysisTaskDmesonJets()
This is the standard destructor.
Byte_t fSelectionType
! for D0: 0=not selected, 1=D0, 2=D0bar, 3=both
static std::pair< AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle * > CheckOrigin(const AliAODMCParticle *part, TClonesArray *mcArray, Bool_t firstParton=kFALSE)
Bool_t ExtractD0Attributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
Select MC particles based on specific prescriptions of HF analysis.
static EMesonDecayChannel_t CheckDecayChannel(const AliAODMCParticle *part, TClonesArray *mcArray)
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:106
AliDmesonJetInfo & operator=(const AliDmesonJetInfo &source)
Double_t fMaxChargedPt
Transverse momentum of the leading charged particle (or track)
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
static FJJetAlgo ConvertToFJAlgo(EJetAlgo_t algo)
decay
Definition: HFPtSpectrum.C:41
std::map< std::string, AliJetInfo > fJets
! list of jets
virtual void PrintAll() const
Definition: External.C:220
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Double_t minMass
void RunParticleLevelAnalysis()
Run a particle level analysis.
std::map< std::string, AliParticleContainer * > fParticleCollArray
particle/track collection array
Bool_t fApplyKinematicCuts
Apply jet kinematic cuts.
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:104
AliTLorentzVector fMomentum
4-momentum of the jet
std::map< std::string, AliClusterContainer * > fClusterCollArray
cluster collection array
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
Create a new THnSparse within the container.
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)