AliPhysics  d2444a6 (d2444a6)
 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 
27 // Aliroot general
28 #include "AliLog.h"
29 #include "AliEMCALGeometry.h"
30 
31 // Aliroot HF
32 #include "AliAODRecoCascadeHF.h"
34 #include "AliRDHFCutsD0toKpi.h"
37 #include "AliHFTrackContainer.h"
38 
39 // Aliroot EMCal jet framework
40 #include "AliEmcalJetTask.h"
41 #include "AliEmcalJet.h"
42 #include "AliJetContainer.h"
43 #include "AliParticleContainer.h"
44 #include "AliEmcalParticle.h"
45 #include "AliFJWrapper.h"
46 
48 
49 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfo
50 
54 
61 Double_t AliAnalysisTaskDmesonJets::AliJetInfo::GetDistance(const AliJetInfo& jet, Double_t& deta, Double_t& dphi) const
62 {
63  dphi = TVector2::Phi_mpi_pi(fMomentum.Phi() - jet.Phi());;
64  deta = fMomentum.Eta() - jet.Eta();
65  return TMath::Sqrt(dphi*dphi + deta*deta);
66 }
67 
73 {
74  Double_t deta = 0;
75  Double_t dphi = 0;
76  return GetDistance(jet, deta, dphi);
77 }
78 
79 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonJetInfo
80 
84 
87 {
88  fD.SetPtEtaPhiE(0,0,0,0);
89  fSoftPionPt = 0;
90  fInvMass2Prong = 0;
91  fDmesonParticle = 0;
92  fMCLabel = -1;
93  fReconstructed = kFALSE;
94  for (auto &jet : fJets) {
95  jet.second.fMomentum.SetPtEtaPhiE(0,0,0,0);
96  jet.second.fNConstituents = 0;
97  jet.second.fNEF = 0;
98  jet.second.fMaxChargedPt = 0;
99  jet.second.fMaxNeutralPt = 0;
100  }
101 }
102 
105 {
106  Printf("Printing D Meson Jet object.");
107  Printf("D Meson: pT = %.3f, eta = %.3f, phi = %.3f, inv. mass = %.3f", fD.Pt(), fD.Eta(), fD.Phi_0_2pi(), fD.M());
108  Printf("Soft pion pT: %.3f. 2-Prong Invariant mass = %.3f", fSoftPionPt, fInvMass2Prong);
109  for (auto &jet : fJets) {
110  Printf("Jet %s: pT = %.3f, eta = %.3f, phi = %.3f", jet.first.c_str(), jet.second.Pt(), jet.second.Eta(), jet.second.Phi_0_2pi());
111  Printf("Jet N Consituents = %d", jet.second.fNConstituents);
112  }
113 }
114 
119 {
120  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
121  if (it == fJets.end()) return 0;
122 
123  Double_t z = 0;
124 
125  if ((*it).second.Pt() > 0) {
126  TVector3 dvect = fD.Vect();
127  TVector3 jvect = (*it).second.fMomentum.Vect();
128 
129  Double_t jetMom = jvect * jvect;
130 
131  if (jetMom < 1e-6) {
132  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetZ", "Zero jet momentum!");
133  z = 0.999;
134  }
135  else {
136  z = (dvect * jvect) / jetMom;
137  }
138 
139  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
140  }
141 
142  return z;
143 }
144 
151 Double_t AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetDistance(std::string n, Double_t& deta, Double_t& dphi) const
152 {
153  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
154  if (it == fJets.end()) return 0;
155 
156  return GetDistance((*it).second, deta, dphi);
157 }
158 
164 {
165  Double_t deta = 0;
166  Double_t dphi = 0;
167  return GetDistance(n, deta, dphi);
168 }
169 
176 Double_t AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetDistance(const AliJetInfo& jet, Double_t& deta, Double_t& dphi) const
177 {
178  dphi = TVector2::Phi_mpi_pi(fD.Phi() - jet.Phi());;
179  deta = fD.Eta() - jet.Eta();
180  return TMath::Sqrt(dphi*dphi + deta*deta);
181 }
182 
188 {
189  Double_t deta = 0;
190  Double_t dphi = 0;
191  return GetDistance(jet, deta, dphi);
192 }
193 
199 {
200  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
201  if (it == fJets.end()) {
202  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
203  n.c_str());
204  return 0;
205  }
206  return &((*it).second);
207 }
208 
214 {
215  std::map<std::string, AliJetInfo>::iterator it = fJets.find(n);
216  if (it == fJets.end()) {
217  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
218  n.c_str());
219  return 0;
220  }
221  return &((*it).second);
222 }
223 
224 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoSummary
225 
229 
235  fPt(0),
236  fEta(0),
237  fPhi(0),
238  fR(0),
239  fZ(0)
240 {
241  Set(source, n);
242 }
243 
246 {
247  fPt = 0;
248  fEta = 0;
249  fPhi = 0;
250  fR = 0;
251  fZ = 0;
252 }
253 
259 {
260  std::map<std::string, AliJetInfo>::const_iterator it = source.fJets.find(n);
261  if (it == source.fJets.end()) return;
262 
263  fPt = (*it).second.Pt();
264  fEta = (*it).second.Eta();
265  fPhi = (*it).second.Phi_0_2pi();
266  fR = source.GetDistance(n);
267  fZ = source.GetZ(n);
268 }
269 
274 {
275  fPt = source.Pt();
276  fEta = source.Eta();
277  fPhi = source.Phi_0_2pi();
278  fR = 0;
279  fZ = 0;
280 }
281 
282 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonInfoSummary
283 
287 
292  fPt(0),
293  fEta(0),
294  fPhi(0)
295 {
296  Set(source);
297 }
298 
303 {
304  fPt = source.fD.Pt();
305  fEta = source.fD.Eta();
306  fPhi = source.fD.Phi_0_2pi();
307 }
308 
311 {
312  fPt = 0;
313  fEta = 0;
314  fPhi = 0;
315 }
316 
317 // Definitions of class AliAnalysisTaskDmesonJets::AliD0InfoSummary
318 
322 
327  AliDmesonInfoSummary(source),
328  fInvMass(source.fD.M())
329 {
330 }
331 
336 {
337  fInvMass = source.fD.M();
339 }
340 
343 {
345 
346  fInvMass = 0;
347 }
348 
349 // Definitions of class AliAnalysisTaskDmesonJets::AliDStarInfoSummary
350 
354 
359  AliDmesonInfoSummary(source),
360  f2ProngInvMass(source.fInvMass2Prong),
361  fDeltaInvMass(source.fD.M() - source.fInvMass2Prong)
362 {
363 }
364 
369 {
370  f2ProngInvMass = source.fInvMass2Prong;
371  fDeltaInvMass = source.fD.M() - source.fInvMass2Prong;
373 }
374 
377 {
379 
380  f2ProngInvMass = 0;
381  fDeltaInvMass = 0;
382 }
383 
384 // Definitions of class AliAnalysisTaskDmesonJets::AliJetDefinition
385 
389 
392  TObject(),
393  fJetType(AliJetContainer::kChargedJet),
394  fRadius(0),
395  fJetAlgo(AliJetContainer::antikt_algorithm),
396  fRecoScheme(AliJetContainer::pt_scheme),
397  fMinJetPt(0.),
398  fMaxJetPt(500.),
399  fMinJetPhi(0.),
400  fMaxJetPhi(0.),
401  fMinJetEta(-1.),
402  fMaxJetEta(1.),
403  fMinChargedPt(0.),
404  fMaxChargedPt(0.),
405  fMinNeutralPt(0.),
406  fMaxNeutralPt(0.)
407 {
408 }
409 
417  TObject(),
418  fJetType(type),
419  fRadius(r),
420  fJetAlgo(algo),
421  fRecoScheme(reco),
422  fMinJetPt(0.),
423  fMaxJetPt(500.),
424  fMinJetPhi(0.),
425  fMaxJetPhi(0.),
426  fMinJetEta(-1.),
427  fMaxJetEta(1.),
428  fMinChargedPt(0.),
429  fMaxChargedPt(0.),
430  fMinNeutralPt(0.),
431  fMaxNeutralPt(0.)
432 {
433 }
434 
439  TObject(),
440  fJetType(source.fJetType),
441  fRadius(source.fRadius),
442  fJetAlgo(source.fJetAlgo),
443  fRecoScheme(source.fRecoScheme),
444  fMinJetPt(source.fMinJetPt),
445  fMaxJetPt(source.fMaxJetPt),
446  fMinJetPhi(source.fMinJetPhi),
447  fMaxJetPhi(source.fMaxJetPhi),
448  fMinJetEta(source.fMinJetEta),
449  fMaxJetEta(source.fMaxJetEta),
450  fMinChargedPt(source.fMinChargedPt),
451  fMaxChargedPt(source.fMaxChargedPt),
452  fMinNeutralPt(source.fMinNeutralPt),
453  fMaxNeutralPt(source.fMaxNeutralPt)
454 {
455 }
456 
461 {
462  new (this) AliHFJetDefinition(source);
463  return *this;
464 }
465 
468 {
469  static TString name;
470 
471  name = AliJetContainer::GenerateJetName(fJetType, fJetAlgo, fRecoScheme, fRadius, 0, 0, "Jet");
472 
473  return name.Data();
474 }
475 
481 {
482  if (fMinJetEta < fMaxJetEta && (jet.Eta() < fMinJetEta || jet.Eta() > fMaxJetEta)) return kFALSE;
483  if (fMinJetPhi < fMaxJetPhi && (jet.Phi() < fMinJetPhi || jet.Phi() > fMaxJetPhi)) return kFALSE;
484  if (jet.Pt() > fMaxJetPt || jet.Pt() < fMinJetPt) return kFALSE;
485  if (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt) return kFALSE;
486  if (jet.fMaxNeutralPt < fMinNeutralPt || jet.fMaxNeutralPt > fMaxNeutralPt) return kFALSE;
487 
488  return kTRUE;
489 }
490 
496 {
497  const AliJetInfo* jet = dMesonJet.GetJet(n);
498  if (!jet) return kFALSE;
499  return IsJetInAcceptance((*jet));
500 }
501 
508 {
509  if (lhs.fJetType > rhs.fJetType) return false;
510  else if (lhs.fJetType < rhs.fJetType) return true;
511  else {
512  if (lhs.fRadius > rhs.fRadius) return false;
513  else if (lhs.fRadius < rhs.fRadius) return true;
514  else {
515  if (lhs.fJetAlgo > rhs.fJetAlgo) return false;
516  else if (lhs.fJetAlgo < rhs.fJetAlgo) return true;
517  else {
518  if (lhs.fRecoScheme < rhs.fRecoScheme) return true;
519  else return false;
520  }
521  }
522  }
523 }
524 
531 {
532  if (lhs.fJetType != rhs.fJetType) return false;
533  if (lhs.fRadius != rhs.fRadius) return false;
534  if (lhs.fJetAlgo != rhs.fJetAlgo) return false;
535  if (lhs.fRecoScheme != rhs.fRecoScheme) return false;
536  return true;
537 }
538 
539 // Definitions of class AliAnalysisTaskDmesonJets::AnalysisEngine
540 
544 
547  TObject(),
548  fCandidateType(kD0toKpi),
549  fCandidateName(),
550  fCandidatePDG(0),
551  fNDaughters(0),
552  fPDGdaughters(),
553  fBranchName(),
554  fMCMode(kNoMC),
555  fNMassBins(0),
556  fMinMass(0),
557  fMaxMass(0),
558  fRDHFCuts(0),
559  fRejectedOrigin(0),
560  fAcceptedDecay(0),
561  fInhibit(kFALSE),
562  fJetDefinitions(),
563  fPtBinWidth(0.5),
564  fMaxPt(100),
565  fDataSlotNumber(-1),
566  fTree(0),
567  fCurrentDmesonJetInfo(0),
568  fCurrentJetInfo(0),
569  fCandidateArray(0),
570  fMCContainer(0),
571  fTrackContainer(0),
572  fClusterContainer(0),
573  fAodEvent(0),
574  fFastJetWrapper(0),
575  fHistManager(0)
576 {
577 }
578 
587  TObject(),
588  fCandidateType(type),
589  fCandidateName(),
590  fCandidatePDG(0),
591  fNDaughters(0),
592  fPDGdaughters(),
593  fBranchName(),
594  fMCMode(MCmode),
595  fNMassBins(nMassBins),
596  fMinMass(0),
597  fMaxMass(0),
598  fRDHFCuts(cuts),
600  fAcceptedDecay(kAnyDecay),
601  fInhibit(kFALSE),
602  fJetDefinitions(),
603  fPtBinWidth(0.5),
604  fMaxPt(100),
605  fDataSlotNumber(-1),
606  fTree(0),
607  fCurrentDmesonJetInfo(0),
608  fCurrentJetInfo(0),
609  fCandidateArray(0),
610  fMCContainer(0),
611  fTrackContainer(0),
612  fClusterContainer(0),
613  fAodEvent(0),
614  fFastJetWrapper(0),
615  fHistManager(0)
616 {
617  SetCandidateProperties(range);
618 }
619 
624  TObject(source),
625  fCandidateType(source.fCandidateType),
626  fCandidateName(source.fCandidateName),
627  fCandidatePDG(source.fCandidatePDG),
628  fNDaughters(source.fNDaughters),
629  fPDGdaughters(source.fPDGdaughters),
630  fBranchName(source.fBranchName),
631  fMCMode(source.fMCMode),
632  fNMassBins(source.fNMassBins),
633  fMinMass(source.fMinMass),
634  fMaxMass(source.fMaxMass),
635  fRDHFCuts(),
636  fRejectedOrigin(source.fRejectedOrigin),
637  fAcceptedDecay(source.fAcceptedDecay),
638  fInhibit(source.fInhibit),
639  fJetDefinitions(source.fJetDefinitions),
640  fPtBinWidth(source.fPtBinWidth),
641  fMaxPt(source.fMaxPt),
642  fDataSlotNumber(-1),
643  fTree(0),
644  fCurrentDmesonJetInfo(0),
645  fCurrentJetInfo(0),
646  fCandidateArray(source.fCandidateArray),
647  fMCContainer(source.fMCContainer),
648  fTrackContainer(source.fTrackContainer),
649  fClusterContainer(source.fClusterContainer),
650  fAodEvent(source.fAodEvent),
652  fHistManager(source.fHistManager)
653 {
654  SetRDHFCuts(source.fRDHFCuts);
655 }
656 
657 // Destructor
659 {
660  delete fRDHFCuts;
661 }
662 
667 {
668  new (this) AnalysisEngine(source);
669  return *this;
670 }
671 
676 {
677  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
678  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
679  }
680 
681  return kFALSE;
682 }
683 
685 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const /*geom*/, Int_t /*runNumber*/)
686 {
687 }
688 
693 {
694  switch (fCandidateType) {
695  case kD0toKpi:
696  fCandidatePDG = 421;
697  fCandidateName = "D0";
698  fNDaughters = 2;
699  fPDGdaughters.Set(fNDaughters);
700  fPDGdaughters.Reset();
701  fPDGdaughters[0] = 211; // pi
702  fPDGdaughters[1] = 321; // K
703  fBranchName = "D0toKpi";
704  fAcceptedDecay = kDecayD0toKpi;
705  if (!fRDHFCuts) {
706  fRDHFCuts = new AliRDHFCutsD0toKpi();
707  fRDHFCuts->SetStandardCutsPP2010();
708  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
709  fRDHFCuts->SetTriggerClass("","");
710  }
711  break;
712  case kDstartoKpipi:
713  fCandidatePDG = 413;
714  fCandidateName = "DStar";
715  fNDaughters = 3;
716  fPDGdaughters.Set(fNDaughters);
717  fPDGdaughters.Reset();
718  fPDGdaughters[0] = 211; // pi soft
719  fPDGdaughters[1] = 211; // pi fromD0
720  fPDGdaughters[2] = 321; // K from D0
721  fBranchName = "Dstar";
722  fAcceptedDecay = kDecayDStartoKpipi;
723  if (!fRDHFCuts) {
724  fRDHFCuts = new AliRDHFCutsDStartoKpipi();
725  fRDHFCuts->SetStandardCutsPP2010();
726  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
727  fRDHFCuts->SetTriggerClass("","");
728  }
729  break;
730  default:
731  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
732  }
733 
734  CalculateMassLimits(range, fCandidatePDG, fNMassBins, fMinMass, fMaxMass);
735 }
736 
741 {
742  if (fRDHFCuts) delete fRDHFCuts;
743  fRDHFCuts = cuts;
744 }
745 
750 {
751  if (!cuts) return;
752  if (fRDHFCuts) delete fRDHFCuts;
753  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
754 }
755 
760 {
761  static TString name;
762 
763  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
764 
765  return name.Data();
766 }
767 
772 {
773  static TString name;
774 
775  name = fCandidateName;
776  switch (fMCMode) {
777  case kBackgroundOnly:
778  name += "_kBackgroundOnly";
779  break;
780  case kSignalOnly:
781  name += "_kSignalOnly";
782  break;
783  case kMCTruth:
784  name += "_MCTruth";
785  break;
786  default:
787  break;
788  }
789 
790  return name.Data();
791 }
792 
800 {
801  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
802 
803  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
804  fJetDefinitions.push_back(def);
805  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
806  "Total number of jet definitions is now %lu.",
807  def.GetName(), GetName(), fJetDefinitions.size());
808  // For detector level set maximum track pt to 100 GeV/c
809  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
810  }
811  else {
812  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
813  }
814 
815  return &(*it);
816 }
817 
829 {
830  AliHFJetDefinition def(type, r, algo, reco);
831 
832  return AddJetDefinition(def);
833 }
834 
840 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
841 {
842  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
843  while (it != fJetDefinitions.end() && (*it) != def) it++;
844  return it;
845 }
846 
853 {
854  if (lhs.fCandidateType > rhs.fCandidateType) return false;
855  else if (lhs.fCandidateType < rhs.fCandidateType) return true;
856  else {
857  if (lhs.fMCMode < rhs.fMCMode) return true;
858  else return false;
859  }
860 }
861 
868 {
869  if (lhs.fCandidateType != rhs.fCandidateType) return false;
870  if (lhs.fMCMode != rhs.fMCMode) return false;
871  return true;
872 }
873 
882 {
883  if (fCandidateType == kD0toKpi) { // D0 candidate
884  return ExtractD0Attributes(Dcand, DmesonJet, i);
885  }
886  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
887  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
888  }
889  else {
890  return kFALSE;
891  }
892 }
893 
902 {
903  Int_t MCtruthPdgCode = 0;
904 
905  Double_t invMassD = 0;
906 
907  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
908  // Checks also the origin, and if it matches the rejected origin mask, return false
909  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
910  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
911  DmesonJet.fMCLabel = mcLab;
912 
913  // Retrieve the generated particle (if exists) and its PDG code
914  if (mcLab >= 0) {
915  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
916 
917  if (aodMcPart) {
918  // Check origin and return false if it matches the rejected origin mask
919  if (fRejectedOrigin) {
920  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
921 
922  UInt_t rs = origin;
923  UShort_t p = 0;
924  while (rs >>= 1) { p++; }
925  TString hname = TString::Format("%s/fHistDmesonOrigin", GetName());
926  fHistManager->FillTH1(hname, p);
927 
928  if ((origin & fRejectedOrigin) == origin) return kFALSE;
929  }
930  MCtruthPdgCode = aodMcPart->PdgCode();
931  }
932  }
933  }
934 
935  //AliDebug(2,"Checking if D0 meson is selected");
936  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
937  if (isSelected == 1) { // selected as a D0
938  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
939 
940  if (fMCMode == kNoMC ||
941  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
942  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly)) {
943  // 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)
944  //AliDebug(2,"Selected as D0");
945  invMassD = Dcand->InvMassD0();
946  }
947  else { // conditions above not passed, so return FALSE
948  return kFALSE;
949  }
950  }
951  else if (isSelected == 2) { // selected as a D0bar
952  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
953 
954  if (fMCMode == kNoMC ||
955  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
956  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly)) {
957  // 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)
958  //AliDebug(2,"Selected as D0bar");
959  invMassD = Dcand->InvMassD0bar();
960  }
961  else { // conditions above not passed, so return FALSE
962  return kFALSE;
963  }
964  }
965  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
966  //AliDebug(2,"Selected as either D0 or D0bar");
967 
968  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
969  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
970  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kBackgroundOnly)) {
971  if (i > 0) return kFALSE;
972  //AliDebug(2, "MC truth is D0");
973  invMassD = Dcand->InvMassD0();
974  }
975  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
976  (MCtruthPdgCode == fCandidatePDG && fMCMode == kBackgroundOnly)) {
977  if (i > 0) return kFALSE;
978  //AliDebug(2, "MC truth is D0bar");
979  invMassD = Dcand->InvMassD0bar();
980  }
981  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
982 
983  // Only accept it if background-only OR background-and-signal was requested
984  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
985  // Select D0 or D0bar depending on the i-parameter
986  if (i == 0) {
987  //AliDebug(2, "Returning invariant mass with D0 hypothesis");
988  invMassD = Dcand->InvMassD0();
989  }
990  else if (i == 1) {
991  //AliDebug(2, "Returning invariant mass with D0bar hypothesis");
992  invMassD = Dcand->InvMassD0bar();
993  }
994  else { // i > 1
995  return kFALSE;
996  }
997  }
998  else { // signal-only was requested but this is not a true D0
999  return kFALSE;
1000  }
1001  }
1002  }
1003  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
1004  return kTRUE;
1005 }
1006 
1015 {
1016  if (i > 0) return kFALSE; // only one mass hypothesis for the D*
1017 
1018  Int_t MCtruthPdgCode = 0;
1019 
1020  Double_t invMassD = 0;
1021 
1022  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1023  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1024  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1025 
1026  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1027  //AliDebug(2, Form("MC label is %d", mcLab));
1028  DmesonJet.fMCLabel = mcLab;
1029  if (mcLab >= 0) {
1030  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1031 
1032  if (aodMcPart) {
1033  if (fRejectedOrigin) {
1034  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
1035 
1036  UInt_t rs = origin;
1037  UShort_t p = 0;
1038  while (rs >>= 1) { p++; }
1039  TString hname = TString::Format("%s/fHistDmesonOrigin", GetName());
1040  fHistManager->FillTH1(hname, p);
1041 
1042  if ((origin & fRejectedOrigin) == origin) return kFALSE;
1043  }
1044 
1045  MCtruthPdgCode = aodMcPart->PdgCode();
1046  //AliDebug(2, Form("MC truth pdg code is %d",MCtruthPdgCode));
1047  }
1048  }
1049  }
1050 
1051  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1052  if (fMCMode == kNoMC ||
1053  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1054  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1055  // 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)
1056  invMassD = DstarCand->InvMassDstarKpipi();
1057  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1058  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1059  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1060  return kTRUE;
1061  }
1062  else { // conditions above not passed, so return FALSE
1063  return kFALSE;
1064  }
1065 }
1066 
1074 {
1075  if (!part) return kDecayOther;
1076  if (!mcArray) return kDecayOther;
1077 
1079 
1080  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1081 
1082  if (part->GetNDaughters() == 2) {
1083 
1084  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1085  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1086 
1087  if (!d1 || !d2) {
1088  return decay;
1089  }
1090 
1091  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1092  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1093 
1094  if (absPdgPart == 421) { // D0 -> K pi
1095  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1096  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1097  decay = kDecayD0toKpi;
1098  }
1099  }
1100 
1101  if (absPdgPart == 413) { // D* -> D0 pi
1102  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1103  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1104  if (D0decay == kDecayD0toKpi) {
1105  decay = kDecayDStartoKpipi;
1106  }
1107  }
1108  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1109  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1110  if (D0decay == kDecayD0toKpi) {
1111  decay = kDecayDStartoKpipi;
1112  }
1113  }
1114  }
1115  }
1116 
1117  return decay;
1118 }
1119 
1127 {
1128  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1129 
1130  if (!part) return kUnknownQuark;
1131  if (!mcArray) return kUnknownQuark;
1132 
1133  Int_t mother = part->GetMother();
1134  while (mother >= 0) {
1135  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1136  if (mcGranma) {
1137  Int_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
1138 
1139  if (abspdgGranma == 1) return kFromDown;
1140  if (abspdgGranma == 2) return kFromUp;
1141  if (abspdgGranma == 3) return kFromStrange;
1142  if (abspdgGranma == 4) return kFromCharm;
1143  if (abspdgGranma == 5) return kFromBottom;
1144  if (abspdgGranma == 6) return kFromTop;
1145 
1146  mother = mcGranma->GetMother();
1147  }
1148  else {
1149  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1150  break;
1151  }
1152  }
1153 
1154  return kUnknownQuark;
1155 }
1156 
1159 {
1160  fDmesonJets.clear();
1161 
1162  for (auto& jetDef : fJetDefinitions) {
1163  jetDef.fJets.clear();
1164  }
1165 
1166  if (fMCMode == kMCTruth) {
1167  RunParticleLevelAnalysis();
1168  }
1169  else {
1170  RunDetectorLevelAnalysis();
1171  }
1172 }
1173 
1179 {
1181  fFastJetWrapper->SetR(jetDef.fRadius);
1184 
1185  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1186  fTrackContainer->SetDMesonCandidate(0);
1187  AddInputVectors(fTrackContainer, 100);
1188  }
1189 
1190  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1191  AddInputVectors(fClusterContainer, -100);
1192  }
1193 
1194  // run jet finder
1195  fFastJetWrapper->Run();
1196 
1197  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1198 
1199  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1200  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1201 
1202  Double_t maxChPt = 0;
1203  Double_t maxNePt = 0;
1204  Double_t totalNeutralPt = 0;
1205 
1206  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1207  if (constituents[ic].user_index() >= 100) {
1208  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1209  }
1210  else if (constituents[ic].user_index() <= -100) {
1211  totalNeutralPt += constituents[ic].pt();
1212  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1213  }
1214  }
1215 
1216  jetDef.fJets.push_back(
1217  AliJetInfo(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E(),
1218  constituents.size(), totalNeutralPt / jets_incl[ijet].pt(), maxChPt, maxNePt)
1219  );
1220  }
1221 }
1222 
1231 {
1232  if (jetDef.fJets.size() == 0) FindJets(jetDef);
1233 
1234  Double_t d_closest = 999;
1235  AliJetInfo* jet_closest = 0;
1236  const AliJetInfo* truth_jet = 0;
1237  try {
1238  truth_jet = &(dmeson.fJets.at(jetDef.GetName()));
1239  }
1240  catch(...) {
1241  return jet_distance_pair(0, 999);
1242  }
1243 
1244  for (auto& jet : jetDef.fJets) {
1245  Double_t d = truth_jet->GetDistance(jet);
1246 
1247  if (d > dMax) continue;
1248  if (d < d_closest) {
1249  d_closest = d;
1250  jet_closest = &jet;
1251  }
1252  }
1253 
1254  if (jet_closest && applyKinCuts) {
1255  if (!jetDef.IsJetInAcceptance(*jet_closest)) jet_closest = 0;
1256  }
1257 
1258  if (jet_closest) {
1259  AliDebug(2, Form("Found closest jet (pt=%.3f, eta=%.3f, phi=%.3f) to d meson (pt=%.3f, eta=%.3f, phi=%.3f) with d = %.3f",
1260  jet_closest->Pt(), jet_closest->Eta(), jet_closest->Phi_0_2pi(),
1261  dmeson.fD.Pt(), dmeson.fD.Eta(), dmeson.fD.Phi_0_2pi(),
1262  d_closest));
1263  }
1264 
1265  return jet_distance_pair(jet_closest, d_closest);
1266 }
1267 
1270 {
1271  const Int_t nD = fCandidateArray->GetEntriesFast();
1272 
1273  AliDmesonJetInfo DmesonJet;
1274 
1275  Int_t nAccCharm = 0;
1276  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1277  Int_t isSelected = 0;
1278 
1279  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1280  if (!charmCand) continue;
1281 
1282  Int_t nprongs = charmCand->GetNProngs();
1283 
1284  if (fCandidateType == kDstartoKpipi) {
1285  if (!charmCand->InheritsFrom("AliAODRecoCascadeHF")) {
1286  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::RunDetectorLevelAnalysis","Candidate type is D* but object type is wrong (should be AliAODRecoCascadeHF)");
1287  continue;
1288  }
1289  }
1290 
1291  // region of interest + cuts
1292  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1293 
1294  //candidate selected by cuts and PID
1295  isSelected = fRDHFCuts->IsSelected(charmCand, AliRDHFCuts::kAll, fAodEvent); //selected
1296 
1297  if (!isSelected) continue;
1298 
1299  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1300  DmesonJet.Reset();
1301  DmesonJet.fDmesonParticle = charmCand;
1302  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1303  for (std::vector<AliHFJetDefinition>::iterator itdef = fJetDefinitions.begin(); itdef != fJetDefinitions.end(); itdef++) {
1304  if (!FindJet(charmCand, DmesonJet, *itdef)) {
1305  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1306  (*itdef).GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1307  }
1308  }
1309  fDmesonJets[icharm] = DmesonJet;
1310  }
1311  }
1312  nAccCharm++;
1313  } // end of D cand loop
1314 
1315  TString hname;
1316 
1317  hname = TString::Format("%s/fHistNAcceptedDmesons", GetName());
1318  fHistManager->FillTH1(hname, nAccCharm);
1319 
1320  hname = TString::Format("%s/fHistNDmesons", GetName());
1321  fHistManager->FillTH1(hname, nD);
1322 }
1323 
1334 {
1335  TString hname;
1336 
1338  fFastJetWrapper->SetR(jetDef.fRadius);
1341 
1342  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1343 
1344  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1345  fTrackContainer->SetDMesonCandidate(Dcand);
1346  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1347  AddInputVectors(fTrackContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1348 
1349  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1350  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1351  const TObjArray& daughters = fTrackContainer->GetDaughterList();
1352  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1353  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1354  if (!fTrackContainer->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1355  }
1356  }
1357 
1358  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1359  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1360  AddInputVectors(fClusterContainer, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1361  }
1362 
1363  // run jet finder
1364  fFastJetWrapper->Run();
1365 
1366  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1367 
1368  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1369  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1370 
1371  Bool_t isDmesonJet = kFALSE;
1372 
1373  Double_t maxChPt = 0;
1374  Double_t maxNePt = 0;
1375  Double_t totalNeutralPt = 0;
1376 
1377  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1378  if (constituents[ic].user_index() == 0) {
1379  isDmesonJet = kTRUE;
1380  }
1381  else if (constituents[ic].user_index() >= 100) {
1382  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1383  }
1384  else if (constituents[ic].user_index() <= -100) {
1385  totalNeutralPt += constituents[ic].pt();
1386  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1387  }
1388  }
1389 
1390  if (isDmesonJet) {
1391  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1392  DmesonJet.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1393  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1394  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1395  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1396 
1397  return kTRUE;
1398  }
1399  }
1400 
1401  return kFALSE;
1402 }
1403 
1407 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist)
1408 {
1409  AliEmcalIterableMomentumContainer itcont = cont->all_momentum();
1410  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1411  UInt_t rejectionReason = 0;
1412  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1413  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1414  continue;
1415  }
1416  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1417  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1418  }
1419 }
1420 
1423 {
1424  TString hname;
1425 
1426  fMCContainer->SetSpecialPDG(fCandidatePDG);
1427  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1428  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1429 
1430  if (!fMCContainer->IsSpecialPDGFound()) return;
1431 
1432  for (auto &jetDef : fJetDefinitions) {
1433  switch (jetDef.fJetType) {
1435  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
1436  break;
1438  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
1439  break;
1441  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
1442  break;
1443  }
1444 
1446  fFastJetWrapper->SetR(jetDef.fRadius);
1449 
1450  hname = TString::Format("%s/fHistDmesonOrigin", GetName());
1451  fMCContainer->SetHistOrigin(static_cast<TH1*>(fHistManager->FindObject(hname)));
1452  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1453  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1454  fMCContainer->SetHistOrigin(0);
1455 
1456  fFastJetWrapper->Run();
1457 
1458  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1459 
1460  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1461  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1462 
1463  Bool_t isDmesonJet = kFALSE;
1464 
1465  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1466  Int_t iPart = constituents[ic].user_index() - 100;
1467  AliVParticle* part = fMCContainer->GetParticle(iPart);
1468  if (!part) {
1469  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1470  continue;
1471  }
1472  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1473  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1474  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1475  std::pair<int, AliDmesonJetInfo> element;
1476  element.first = iPart;
1477  element.second = AliDmesonJetInfo();
1478  dMesonJetIt = fDmesonJets.insert(element).first;
1479  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1480  (*dMesonJetIt).second.fDmesonParticle = part;
1481  }
1482 
1483  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1484  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1485  }
1486  }
1487  }
1488  }
1489 }
1490 
1495 {
1496  TString classname;
1497  switch (fCandidateType) {
1498  case kD0toKpi:
1499  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1500  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1501  break;
1502  case kDstartoKpipi:
1503  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1504  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1505  break;
1506  }
1507  TString treeName = TString::Format("%s_%s", taskName, GetName());
1508  fTree = new TTree(treeName, treeName);
1509  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
1510  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1511  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1512  fCurrentJetInfo[i] = new AliJetInfoSummary();
1513  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1514  }
1515 
1516  return fTree;
1517 }
1518 
1523 {
1524  TString hname;
1525 
1526  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
1527 
1528  for (auto &jetDef : fJetDefinitions) {
1529 
1530  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
1531 
1532  Double_t radius = jetDef.fRadius;
1533 
1534  TString title[30] = {""};
1535  Int_t nbins[30] = {0 };
1536  Double_t min [30] = {0.};
1537  Double_t max [30] = {0.};
1538  Int_t dim = 0 ;
1539 
1540  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1541  nbins[dim] = nPtBins;
1542  min[dim] = 0;
1543  max[dim] = fMaxPt;
1544  dim++;
1545 
1546  if ((enabledAxis & kPositionD) != 0) {
1547  title[dim] = "#eta_{D}";
1548  nbins[dim] = 50;
1549  min[dim] = -1;
1550  max[dim] = 1;
1551  dim++;
1552 
1553  title[dim] = "#phi_{D} (rad)";
1554  nbins[dim] = 150;
1555  min[dim] = 0;
1556  max[dim] = TMath::TwoPi();
1557  dim++;
1558  }
1559 
1560  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1561  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1562  nbins[dim] = fNMassBins;
1563  min[dim] = fMinMass;
1564  max[dim] = fMaxMass;
1565  dim++;
1566  }
1567 
1568  if (fCandidateType == kD0toKpi) {
1569  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1570  nbins[dim] = fNMassBins;
1571  min[dim] = fMinMass;
1572  max[dim] = fMaxMass;
1573  dim++;
1574  }
1575 
1576  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1577  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1578  nbins[dim] = fNMassBins;
1579  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1580  dim++;
1581  }
1582 
1583  if (fCandidateType == kDstartoKpipi) {
1584  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1585  nbins[dim] = fNMassBins*6;
1586  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1587 
1588  // subtract mass of D0
1589  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1590  min[dim] -= D0mass;
1591  max[dim] -= D0mass;
1592 
1593  dim++;
1594  }
1595 
1596  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1597  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1598  nbins[dim] = 100;
1599  min[dim] = 0;
1600  max[dim] = 25;
1601  dim++;
1602  }
1603 
1604  title[dim] = "#it{z}_{D}";
1605  nbins[dim] = 110;
1606  min[dim] = 0;
1607  max[dim] = 1.10;
1608  dim++;
1609 
1610  if ((enabledAxis & kDeltaR) != 0) {
1611  title[dim] = "#Delta R_{D-jet}";
1612  nbins[dim] = 100;
1613  min[dim] = 0;
1614  max[dim] = radius * 1.5;
1615  dim++;
1616  }
1617 
1618  if ((enabledAxis & kDeltaEta) != 0) {
1619  title[dim] = "#eta_{D} - #eta_{jet}";
1620  nbins[dim] = 100;
1621  min[dim] = -radius * 1.2;
1622  max[dim] = radius * 1.2;
1623  dim++;
1624  }
1625 
1626  if ((enabledAxis & kDeltaPhi) != 0) {
1627  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1628  nbins[dim] = 100;
1629  min[dim] = -radius * 1.2;
1630  max[dim] = radius * 1.2;
1631  dim++;
1632  }
1633 
1634  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1635  nbins[dim] = nPtBins;
1636  min[dim] = 0;
1637  max[dim] = fMaxPt;
1638  dim++;
1639 
1640  if ((enabledAxis & kPositionJet) != 0) {
1641  title[dim] = "#eta_{jet}";
1642  nbins[dim] = 50;
1643  min[dim] = -1;
1644  max[dim] = 1;
1645  dim++;
1646 
1647  title[dim] = "#phi_{jet} (rad)";
1648  nbins[dim] = 150;
1649  min[dim] = 0;
1650  max[dim] = TMath::TwoPi();
1651  dim++;
1652  }
1653 
1654  if ((enabledAxis & kJetConstituents) != 0) {
1655  title[dim] = "No. of constituents";
1656  nbins[dim] = 50;
1657  min[dim] = -0.5;
1658  max[dim] = 49.5;
1659  dim++;
1660  }
1661 
1662  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1663  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1664  for (Int_t j = 0; j < dim; j++) {
1665  h->GetAxis(j)->SetTitle(title[j]);
1666  }
1667  }
1668 }
1669 
1674 {
1675  TString hname;
1676 
1677  for (auto& dmeson_pair : fDmesonJets) {
1678  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
1679  Int_t accJets = 0;
1680  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1681  fCurrentJetInfo[ij]->Reset();
1682  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1683  if (!jet) continue;
1684  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1685  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1686  fHistManager->FillTH1(hname, jet->Pt());
1687  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1688  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1689  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1690  fHistManager->FillTH1(hname, jet->Eta());
1691  continue;
1692  }
1693  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
1694  accJets++;
1695  }
1696  if (accJets > 0) {
1697  fTree->Fill();
1698  }
1699  else {
1700  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1701  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1702  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1703  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1704  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1705  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1706  if (fCandidateType == kD0toKpi) {
1707  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1708  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1709  }
1710  else if (fCandidateType == kDstartoKpipi) {
1711  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1712  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1713 
1714  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1715  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1716  }
1717  }
1718  }
1719  return kTRUE;
1720 }
1721 
1727 {
1728  TString hname;
1729 
1730  for (auto& dmeson_pair : fDmesonJets) {
1731  Int_t accJets = 0;
1732  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1733  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1734  if (!jet) continue;
1735  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1736  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1737  fHistManager->FillTH1(hname, jet->Pt());
1738  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1739  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1740  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1741  fHistManager->FillTH1(hname, jet->Eta());
1742  continue;
1743  }
1744  accJets++;
1745  }
1746  if (accJets == 0) {
1747  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1748  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1749  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1750  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1751  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1752  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1753  if (fCandidateType == kD0toKpi) {
1754  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1755  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1756  }
1757  else if (fCandidateType == kDstartoKpipi) {
1758  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1759  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1760 
1761  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1762  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1763  }
1764  }
1765  }
1766  return kTRUE;
1767 }
1768 
1773 {
1774  TString hname;
1775 
1776  for (auto& dmeson_pair : fDmesonJets) {
1777  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
1778  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1779  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1780  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1781  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1782  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1783  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1784  }
1785  }
1786 
1787  for (auto &jetDef : fJetDefinitions) {
1788 
1789  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1790  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
1791 
1792  for (auto& dmeson_pair : fDmesonJets) {
1793  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
1794  if (!jet) continue;
1795  if (!jetDef.IsJetInAcceptance(*jet)) {
1796  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
1797  fHistManager->FillTH1(hname, jet->Pt());
1798  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
1799  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1800  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
1801  fHistManager->FillTH1(hname, jet->Eta());
1802  continue;
1803  }
1804  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
1805  }
1806  }
1807 
1808  return kTRUE;
1809 }
1810 
1816 Bool_t AliAnalysisTaskDmesonJets::AnalysisEngine::FillHnSparse(THnSparse* h, const AliDmesonJetInfo& DmesonJet, std::string n)
1817 {
1818  // Fill the THnSparse histogram.
1819 
1820  Double_t contents[30] = {0.};
1821 
1822  Double_t z = DmesonJet.GetZ(n);
1823  Double_t deltaPhi = 0;
1824  Double_t deltaEta = 0;
1825  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
1826 
1827  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
1828  if (it == DmesonJet.fJets.end()) return kFALSE;
1829 
1830  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
1831  TString title(h->GetAxis(i)->GetTitle());
1832  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
1833  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
1834  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
1835  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
1836  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
1837  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
1838  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
1839  else if (title=="#it{z}_{D}") contents[i] = z;
1840  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
1841  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
1842  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
1843  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
1844  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
1845  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
1846  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
1847  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
1848  }
1849 
1850  h->Fill(contents);
1851 
1852  return kTRUE;
1853 }
1854 
1855 // Definitions of class AliAnalysisTaskDmesonJets
1856 
1860 
1864  fAnalysisEngines(),
1865  fEnabledAxis(0),
1867  fHistManager(),
1868  fApplyKinematicCuts(kTRUE),
1869  fNOutputTrees(0),
1870  fAodEvent(0),
1871  fFastJetWrapper(0)
1872 {
1873  SetMakeGeneralHistograms(kTRUE);
1874 }
1875 
1879 AliAnalysisTaskDmesonJets::AliAnalysisTaskDmesonJets(const char* name, Int_t nOutputTrees) :
1880  AliAnalysisTaskEmcalLight(name, kTRUE),
1881  fAnalysisEngines(),
1882  fEnabledAxis(k2ProngInvMass),
1883  fOutputType(kTreeOutput),
1884  fHistManager(name),
1885  fApplyKinematicCuts(kTRUE),
1886  fNOutputTrees(nOutputTrees),
1887  fAodEvent(0),
1888  fFastJetWrapper(0)
1889 {
1890  SetMakeGeneralHistograms(kTRUE);
1891  for (Int_t i = 0; i < nOutputTrees; i++){
1892  DefineOutput(2+i, TTree::Class());
1893  }
1894 }
1895 
1898 {
1899  if (fFastJetWrapper) delete fFastJetWrapper;
1900 }
1901 
1909 {
1910  AliRDHFCuts* analysiscuts = 0;
1911  TFile* filecuts = TFile::Open(cutfname);
1912  if (!filecuts || filecuts->IsZombie()) {
1913  ::Warning("AddTaskDmesonJets", "Input file not found: will use std cuts.");
1914  filecuts = 0;
1915  }
1916 
1917  if (filecuts) {
1918  analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
1919  if (!analysiscuts) {
1920  ::Warning("AddTaskDmesonJetCorr", "Could not find analysis cuts '%s' in '%s'. Using std cuts.", cutsname.Data(), cutfname.Data());
1921  }
1922  }
1923 
1924  return analysiscuts;
1925 }
1926 
1936 {
1938  return AddAnalysisEngine(type, MCmode, jetDef, cutfname);
1939 }
1940 
1950 {
1951  AliRDHFCuts* cuts = 0;
1952 
1953  if (!cutfname.IsNull()) {
1954  TString cutsname;
1955 
1956  switch (type) {
1957  case kD0toKpi :
1958  cutsname = "D0toKpiCuts";
1959  break;
1960  case kDstartoKpipi :
1961  cutsname = "DStartoKpipiCuts";
1962  break;
1963  default:
1964  return 0;
1965  }
1966 
1967  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
1968  }
1969 
1970  AnalysisEngine eng(type, MCmode, cuts);
1971 
1972  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
1973 
1974  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
1975  eng.AddJetDefinition(jetDef);
1976  it = fAnalysisEngines.insert(it, eng);
1977  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(jetDef), fAnalysisEngines.size());
1978  }
1979  else {
1980  AnalysisEngine* found_eng = &(*it);
1981  ::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());
1982  found_eng->AddJetDefinition(jetDef);
1983 
1984  if (cuts && found_eng->fRDHFCuts != 0) {
1985  ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
1986  found_eng->SetRDHFCuts(cuts);
1987  }
1988  }
1989 
1990  return &(*it);
1991 }
1992 
1993 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
1994 {
1995  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
1996  while (it != fAnalysisEngines.end() && (*it) < eng) it++;
1997  return it;
1998 }
1999 
2002 {
2003  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
2004 
2006 
2007  // Define histograms
2008  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
2009 
2010  TString hname;
2011  TString htitle;
2012  TH1* h = 0;
2013  Int_t treeSlot = 0;
2014 
2015  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
2016  for (auto &param : fAnalysisEngines) {
2017  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
2018 
2019  fHistManager.CreateHistoGroup(param.GetName());
2020 
2021  param.fHistManager = &fHistManager;
2022 
2023  hname = TString::Format("%s/fHistNAcceptedDmesons", param.GetName());
2024  htitle = hname + ";Number of D accepted meson candidates;counts";
2025  h = fHistManager.CreateTH1(hname, htitle, 51, -0.5, 50.5);
2026 
2027  hname = TString::Format("%s/fHistNDmesons", param.GetName());
2028  htitle = hname + ";Number of D meson candidates;counts";
2029  h = fHistManager.CreateTH1(hname, htitle, 101, -0.5, 100.5);
2030 
2031  hname = TString::Format("%s/fHistNEvents", param.GetName());
2032  htitle = hname + ";Event status;counts";
2033  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
2034  h->GetXaxis()->SetBinLabel(1, "Accepted");
2035  h->GetXaxis()->SetBinLabel(2, "Rejected");
2036 
2037  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
2038  htitle = hname + ";Rejection reason;counts";
2039  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
2040 
2041  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
2042  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
2043  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2044 
2045  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
2046  htitle = hname + ";#it{#eta}_{D};counts";
2047  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2048 
2049  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
2050  htitle = hname + ";#it{#phi}_{D};counts";
2051  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2052 
2053  if (param.fCandidateType == kD0toKpi) {
2054  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
2055  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2056  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
2057  }
2058  else if (param.fCandidateType == kDstartoKpipi) {
2059  Double_t min = 0;
2060  Double_t max = 0;
2061 
2062  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
2063  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2064  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
2065  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
2066 
2067  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2068  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
2069  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2070  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
2071  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
2072  }
2073 
2074  if (param.fMCMode == kBackgroundOnly || param.fMCMode == kSignalOnly || param.fMCMode == kMCTruth) {
2075  hname = TString::Format("%s/fHistDmesonOrigin", param.GetName());
2076  htitle = hname + ";origin;counts";
2077  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2078  }
2079 
2080  for (std::vector<AliHFJetDefinition>::iterator itdef = param.fJetDefinitions.begin(); itdef != param.fJetDefinitions.end(); itdef++) {
2081  AliHFJetDefinition* jetDef = &(*itdef);
2082  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef->GetName());
2083 
2084  fHistManager.CreateHistoGroup(jetDef->GetName(), param.GetName());
2085 
2086  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef->GetName());
2087  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2088  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2089  SetRejectionReasonLabels(h->GetXaxis());
2090 
2091  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef->GetName());
2092  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2093  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2094  SetRejectionReasonLabels(h->GetXaxis());
2095 
2096  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef->GetName());
2097  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
2098  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2099  SetRejectionReasonLabels(h->GetXaxis());
2100 
2101  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef->GetName());
2102  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
2103  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
2104 
2105  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef->GetName());
2106  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
2107  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2108 
2109  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef->GetName());
2110  htitle = hname + ";#it{#eta}_{jet};counts";
2111  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2112 
2113  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef->GetName());
2114  htitle = hname + ";#it{#phi}_{jet};counts";
2115  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2116  }
2117  switch (fOutputType) {
2118  case kTreeOutput:
2119  param.BuildTree(GetName());
2120  if (treeSlot < fNOutputTrees) {
2121  param.AssignDataSlot(treeSlot+2);
2122  treeSlot++;
2124  }
2125  else {
2126  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
2127  }
2128  break;
2129  case kTHnOutput:
2130  param.BuildHnSparse(fEnabledAxis);
2131  break;
2132  case kNoOutput:
2133  break;
2134  }
2135  }
2136 
2138 
2139  PostData(1, fOutput);
2140 }
2141 
2145 {
2147 
2148  // Load the event
2149  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
2150 
2151  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
2152 
2153  fFastJetWrapper->SetAreaType(fastjet::active_area);
2155 
2156  if (!fAodEvent) {
2157  AliError(Form("This task need an AOD event! Task '%s' will be disabled!", GetName()));
2158  return;
2159  }
2160 
2161  for (auto &params : fAnalysisEngines) {
2162 
2163  params.fAodEvent = fAodEvent;
2164  params.fFastJetWrapper = fFastJetWrapper;
2165  params.Init(fGeom, fAodEvent->GetRunNumber());
2166 
2167  if (params.fMCMode != kMCTruth) {
2168  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
2169 
2170  if (params.fCandidateArray) {
2171  if (!params.fCandidateArray->GetClass()->InheritsFrom("AliAODRecoDecayHF2Prong")) {
2172  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2173  "%s: Objects of type %s in %s are not inherited from AliAODRecoDecayHF2Prong! Task will be disabled!",
2174  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName());
2175  params.fCandidateArray = 0;
2176  params.fInhibit = kTRUE;
2177  }
2178  }
2179  else {
2180  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2181  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
2182  params.fBranchName.Data(), params.GetName());
2183  params.fInhibit = kTRUE;
2184  }
2185  }
2186 
2187  if (params.fMCMode != kNoMC) {
2188  params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(0));
2189 
2190  if (!params.fMCContainer) params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(1));
2191 
2192  if (!params.fMCContainer) {
2193  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2194  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
2195  params.GetName());
2196  params.fInhibit = kTRUE;
2197  }
2198  }
2199 
2200  if (params.fMCMode != kMCTruth) {
2201  params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(0));
2202  if (!params.fTrackContainer) params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(1));
2203 
2204  params.fClusterContainer = GetClusterContainer(0);
2205 
2206  if (!params.fTrackContainer && !params.fClusterContainer) {
2207  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2208  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
2209  params.GetName());
2210  params.fInhibit = kTRUE;
2211  }
2212  }
2213  }
2214 }
2215 
2220 {
2221  if (!fAodEvent) return kFALSE;
2222 
2223  TString hname;
2224 
2225  // fix for temporary bug in ESDfilter
2226  // the AODs with null vertex pointer didn't pass the PhysSel
2227  // Now adding an entry in teh histogram so as to check that this is actually cutting anything out
2228  if (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001) {
2229  for (auto &eng : fAnalysisEngines) {
2230  if (eng.fInhibit) continue;
2231  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2232  fHistManager.FillTH1(hname, "ESDfilterBug");
2233  }
2234  return kFALSE;
2235  }
2236 
2237  for (auto &eng : fAnalysisEngines) {
2238 
2239  if (eng.fInhibit) continue;
2240 
2241  //Event selection
2242  hname = TString::Format("%s/fHistNEvents", eng.GetName());
2243  Bool_t iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
2244  if (!iseventselected) {
2245  fHistManager.FillTH1(hname, "Rejected");
2246  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2247  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
2248  TString label;
2249  do {
2250  label = GetHFEventRejectionReasonLabel(bitmap);
2251  if (label.IsNull()) break;
2252  fHistManager.FillTH1(hname, label);
2253  } while (true);
2254 
2255  continue;
2256  }
2257 
2258  fHistManager.FillTH1(hname, "Accepted");
2259 
2260  AliDebug(2, "Event selected");
2261 
2262  eng.RunAnalysis();
2263  }
2264  return kTRUE;
2265 }
2266 
2271 {
2272  TString hname;
2273  for (auto &param : fAnalysisEngines) {
2274 
2275  if (param.fInhibit) continue;
2276 
2277  if (fOutputType == kTreeOutput) {
2278  param.FillTree(fApplyKinematicCuts);
2280  }
2281  else if (fOutputType == kTHnOutput) {
2282  param.FillHnSparse(fApplyKinematicCuts);
2283  }
2284  }
2285  return kTRUE;
2286 }
2287 
2295 void AliAnalysisTaskDmesonJets::CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t& minMass, Double_t& maxMass)
2296 {
2297  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
2298 
2299  Double_t mass = part->Mass();
2300 
2301  // To make sure that the PDG mass value is not at the edge of a bin
2302  if (nbins % 2 == 0) {
2303  minMass = mass - range / 2 - range / nbins / 2;
2304  maxMass = mass + range / 2 - range / nbins / 2;
2305  }
2306  else {
2307  minMass = mass - range / 2;
2308  maxMass = mass + range / 2;
2309  }
2310 }
2311 
2318 {
2319  static TString label;
2320  label = "";
2321 
2322  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
2323  label = "NotSelTrigger";
2324  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
2325  return label.Data();
2326  }
2327  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
2328  label = "NoVertex";
2329  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
2330  return label.Data();
2331  }
2332  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
2333  label = "TooFewVtxContrib";
2334  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
2335  return label.Data();
2336  }
2337  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
2338  label = "ZVtxOutFid";
2339  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
2340  return label.Data();
2341  }
2342  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
2343  label = "Pileup";
2344  bitmap &= ~BIT(AliRDHFCuts::kPileup);
2345  return label.Data();
2346  }
2347  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
2348  label = "OutsideCentrality";
2349  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
2350  return label.Data();
2351  }
2352  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
2353  label = "PhysicsSelection";
2354  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
2355  return label.Data();
2356  }
2357  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
2358  label = "BadSPDVertex";
2359  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
2360  return label.Data();
2361  }
2362  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
2363  label = "ZVtxSPDOutFid";
2364  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
2365  return label.Data();
2366  }
2367  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
2368  label = "CentralityFlattening";
2369  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
2370  return label.Data();
2371  }
2372  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
2373  label = "BadTrackV0Correl";
2374  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
2375  return label.Data();
2376  }
2377 
2378  return label.Data();
2379 }
2380 
2387 {
2388  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
2389  PostData(eng.GetDataSlotNumber(), eng.GetTree());
2390  return eng.GetDataSlotNumber();
2391  }
2392  else {
2393  return -1;
2394  }
2395 }
Int_t pdg
void Print() const
Prints the content of this object in the standard output.
std::list< AnalysisEngine > fAnalysisEngines
Array of analysis parameters.
AliEMCALGeometry * fGeom
!emcal geometry
Analysis task for D meson jets.
virtual void Set(const AliDmesonJetInfo &source)
Invariant mass of the D0 meson candidate in GeV/c2.
virtual Int_t Run()
Bool_t IsJetInAcceptance(const AliJetInfo &jet) const
Double_t fSoftPionPt
! Transverse momentum of the soft pion of the D* candidate
const char * title
Definition: MakeQAPdf.C:26
AliRDHFCuts * fRDHFCuts
D meson candidates cuts.
Lightweight class that encapsulates D meson jets.
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
ERecoScheme_t fRecoScheme
Jet recombination scheme (pt scheme, E scheme, ...)
Double_t fMaxNeutralPt
Transverse momentum of the leading neutral particle (or cluster)
void RunDetectorLevelAnalysis()
Run a detector level analysis.
Double_t mass
jet_distance_pair FindJetMacthedToGeneratedDMeson(const AliDmesonJetInfo &dmeson, AliHFJetDefinition &jetDef, Double_t dMax, Bool_t applyKinCuts)
virtual void UserCreateOutputObjects()
Creates the output containers.
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.
EMCMode_t fMCMode
MC mode: No MC (data and MC detector level), background-only (MC), signal-only (MC), MC truth (particle level)
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:96
AnalysisEngine()
This is the default constructor, used for ROOT I/O purposes.
void Init(const AliEMCALGeometry *const geom, Int_t runNumber)
Initialize the analysis engine.
static FJRecoScheme ConvertToFJRecoScheme(ERecoScheme_t reco)
AliVParticle * fDmesonParticle
! pointer to the particle object
static TString GenerateJetName(EJetType_t jetType, EJetAlgo_t jetAlgo, ERecoScheme_t recoScheme, Double_t radius, AliParticleContainer *partCont, AliClusterContainer *clusCont, TString tag)
void SetR(Double_t r)
Definition: AliFJWrapper.h:101
Select tracks based on specific prescriptions of HF analysis.
ECandidateType_t fCandidateType
Candidate type.
Double_t GetDistance(std::string n, Double_t &deta, Double_t &dphi) const
Int_t fMCLabel
! MC label, i.e. index of the generator level D meson (only for detector level D meson candidates) ...
THashList * CreateHistoGroup(const char *groupname, const char *parent="/")
virtual void Set(const AliDmesonJetInfo &source)
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:33
EJetType_t fJetType
Jet type (charged, full, neutral)
std::pair< AliJetInfo *, Double_t > jet_distance_pair
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
virtual void Set(const AliDmesonJetInfo &source)
AliTLorentzVector fD
! 4-momentum of the D meson candidate
Double_t GetDistance(const AliJetInfo &jet, Double_t &deta, Double_t &dphi) const
TObject * FindObject(const char *name) const
THashList * GetListOfHistograms() const
Definition: THistManager.h:504
AliAnalysisTaskDmesonJets()
This is the default constructor, used for ROOT I/O purposes.
virtual void Reset()
Reset the current object.
AnalysisEngine & operator=(const AnalysisEngine &source)
Struct that encapsulates analysis parameters.
std::vector< AliHFJetDefinition >::iterator FindJetDefinition(const AliHFJetDefinition &eng)
Base task in the EMCAL framework (lighter version of AliAnalysisTaskEmcal)
Double_t Phi_0_2pi() const
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:95
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
AliClusterContainer * GetClusterContainer(Int_t i=0) const
AliAODTrack * GetBachelor() const
AliRDHFCuts * LoadDMesonCutsFromFile(TString cutfname, TString cutsname)
bool operator<(const AliAnalysisTaskDmesonJets::AliHFJetDefinition &lhs, const AliAnalysisTaskDmesonJets::AliHFJetDefinition &rhs)
std::list< AnalysisEngine >::iterator FindAnalysisEngine(const AnalysisEngine &eng)
std::vector< AliHFJetDefinition > fJetDefinitions
Jet definitions.
EJetAlgo_t fJetAlgo
Jet algorithm (kt, anti-kt,...)
static const char * GetHFEventRejectionReasonLabel(UInt_t &bitmap)
Lightweight class that encapsulates D*.
virtual void Clear(const Option_t *="")
Bool_t ExtractDstarAttributes(const AliAODRecoCascadeHF *DstarCand, AliDmesonJetInfo &DmesonJet, UInt_t i)
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
const AliJetInfo * GetJet(std::string n) const
virtual ~AliAnalysisTaskDmesonJets()
This is the standard destructor.
Bool_t ExtractD0Attributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
std::vector< AliJetInfo > fJets
! Inclusive jets reconstructed in the current event (includes D meson candidate daughters, if any)
Select MC particles based on specific prescriptions of HF analysis.
static EMesonDecayChannel_t CheckDecayChannel(const AliAODMCParticle *part, TClonesArray *mcArray)
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:99
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
void AddInputVectors(AliEmcalContainer *cont, Int_t offset, TH2 *rejectHist=0)
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Double_t minMass
void RunParticleLevelAnalysis()
Run a particle level analysis.
Bool_t fApplyKinematicCuts
Apply jet kinematic cuts.
const char * GetName() const
Generate a name for this jet definition.
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
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)
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:97
AliTLorentzVector fMomentum
4-momentum of the jet
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
AliHFJetDefinition * AddJetDefinition(EJetType_t type, Double_t r, EJetAlgo_t algo, ERecoScheme_t reco)
Lightweight class that encapsulates D0.
const Int_t nPtBins
Bool_t IsAnyJetInAcceptance(const AliDmesonJetInfo &dMesonJet) const
Container for jet within the EMCAL jet framework.
Int_t PostDataFromAnalysisEngine(const AnalysisEngine &eng)
static EMesonOrigin_t CheckOrigin(const AliAODMCParticle *part, TClonesArray *mcArray)
AliFJWrapper * fFastJetWrapper
! Fastjet wrapper
THistManager fHistManager
Histogram manager.
void Reset()
Reset all fields to their default values.
virtual void Set(const AliDmesonJetInfo &source, std::string n)