AliPhysics  a34469b (a34469b)
 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 #include "AliAnalysisManager.h"
31 #include "AliVEventHandler.h"
32 
33 // Aliroot HF
34 #include "AliAODRecoCascadeHF.h"
36 #include "AliRDHFCutsD0toKpi.h"
39 #include "AliHFTrackContainer.h"
40 
41 // Aliroot EMCal jet framework
42 #include "AliEmcalJetTask.h"
43 #include "AliEmcalJet.h"
44 #include "AliJetContainer.h"
45 #include "AliParticleContainer.h"
46 #include "AliEmcalParticle.h"
47 #include "AliFJWrapper.h"
48 
50 
51 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfo
52 
56 
64 {
65  dphi = TVector2::Phi_mpi_pi(fMomentum.Phi() - jet.Phi());;
66  deta = fMomentum.Eta() - jet.Eta();
67  return TMath::Sqrt(dphi*dphi + deta*deta);
68 }
69 
75 {
76  Double_t deta = 0;
77  Double_t dphi = 0;
78  return GetDistance(jet, deta, dphi);
79 }
80 
81 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonJetInfo
82 
86 
91  fDmesonParticle(source.fDmesonParticle),
92  fD(source.fD),
93  fSoftPionPt(source.fSoftPionPt),
94  fInvMass2Prong(source.fInvMass2Prong),
95  fJets(source.fJets),
96  fMCLabel(source.fMCLabel),
97  fReconstructed(source.fReconstructed)
98 {
99 }
100 
105 {
106  new (this) AliDmesonJetInfo(source);
107  return *this;
108 }
109 
112 {
113  fD.SetPtEtaPhiE(0,0,0,0);
114  fSoftPionPt = 0;
115  fInvMass2Prong = 0;
116  fDmesonParticle = 0;
117  fMCLabel = -1;
118  fReconstructed = kFALSE;
119  for (auto &jet : fJets) {
120  jet.second.fMomentum.SetPtEtaPhiE(0,0,0,0);
121  jet.second.fNConstituents = 0;
122  jet.second.fNEF = 0;
123  jet.second.fMaxChargedPt = 0;
124  jet.second.fMaxNeutralPt = 0;
125  }
126 }
127 
130 {
131  Printf("Printing D Meson Jet object.");
132  Printf("D Meson: pT = %.3f, eta = %.3f, phi = %.3f, inv. mass = %.3f", fD.Pt(), fD.Eta(), fD.Phi_0_2pi(), fD.M());
133  Printf("Soft pion pT: %.3f. 2-Prong Invariant mass = %.3f", fSoftPionPt, fInvMass2Prong);
134  for (auto &jet : fJets) {
135  Printf("Jet %s: pT = %.3f, eta = %.3f, phi = %.3f", jet.first.c_str(), jet.second.Pt(), jet.second.Eta(), jet.second.Phi_0_2pi());
136  Printf("Jet N Consituents = %d", jet.second.fNConstituents);
137  }
138 }
139 
144 {
145  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
146  if (it == fJets.end()) return 0;
147 
148  Double_t z = 0;
149 
150  if ((*it).second.Pt() > 0) {
151  TVector3 dvect = fD.Vect();
152  TVector3 jvect = (*it).second.fMomentum.Vect();
153 
154  Double_t jetMom = jvect * jvect;
155 
156  if (jetMom < 1e-6) {
157  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetZ", "Zero jet momentum!");
158  z = 0.999;
159  }
160  else {
161  z = (dvect * jvect) / jetMom;
162  }
163 
164  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
165  }
166 
167  return z;
168 }
169 
177 {
178  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
179  if (it == fJets.end()) return 0;
180 
181  return GetDistance((*it).second, deta, dphi);
182 }
183 
189 {
190  Double_t deta = 0;
191  Double_t dphi = 0;
192  return GetDistance(n, deta, dphi);
193 }
194 
202 {
203  dphi = TVector2::Phi_mpi_pi(fD.Phi() - jet.Phi());;
204  deta = fD.Eta() - jet.Eta();
205  return TMath::Sqrt(dphi*dphi + deta*deta);
206 }
207 
213 {
214  Double_t deta = 0;
215  Double_t dphi = 0;
216  return GetDistance(jet, deta, dphi);
217 }
218 
224 {
225  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
226  if (it == fJets.end()) {
227  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
228  n.c_str());
229  return 0;
230  }
231  return &((*it).second);
232 }
233 
239 {
240  std::map<std::string, AliJetInfo>::iterator it = fJets.find(n);
241  if (it == fJets.end()) {
242  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
243  n.c_str());
244  return 0;
245  }
246  return &((*it).second);
247 }
248 
249 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoSummary
250 
254 
260  fPt(0),
261  fEta(0),
262  fPhi(0),
263  fR(0),
264  fZ(0)
265 {
266  Set(source, n);
267 }
268 
271 {
272  fPt = 0;
273  fEta = 0;
274  fPhi = 0;
275  fR = 0;
276  fZ = 0;
277 }
278 
284 {
285  std::map<std::string, AliJetInfo>::const_iterator it = source.fJets.find(n);
286  if (it == source.fJets.end()) return;
287 
288  fPt = (*it).second.Pt();
289  fEta = (*it).second.Eta();
290  fPhi = (*it).second.Phi_0_2pi();
291  fR = source.GetDistance(n);
292  fZ = source.GetZ(n);
293 }
294 
299 {
300  fPt = source.Pt();
301  fEta = source.Eta();
302  fPhi = source.Phi_0_2pi();
303  fR = 0;
304  fZ = 0;
305 }
306 
307 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonInfoSummary
308 
312 
317  fPt(0),
318  fEta(0),
319  fPhi(0)
320 {
321  Set(source);
322 }
323 
328 {
329  fPt = source.fD.Pt();
330  fEta = source.fD.Eta();
331  fPhi = source.fD.Phi_0_2pi();
332 }
333 
336 {
337  fPt = 0;
338  fEta = 0;
339  fPhi = 0;
340 }
341 
342 // Definitions of class AliAnalysisTaskDmesonJets::AliD0InfoSummary
343 
347 
352  AliDmesonInfoSummary(source),
353  fInvMass(source.fD.M())
354 {
355 }
356 
361 {
362  fInvMass = source.fD.M();
364 }
365 
368 {
370 
371  fInvMass = 0;
372 }
373 
374 // Definitions of class AliAnalysisTaskDmesonJets::AliDStarInfoSummary
375 
379 
384  AliDmesonInfoSummary(source),
385  f2ProngInvMass(source.fInvMass2Prong),
386  fDeltaInvMass(source.fD.M() - source.fInvMass2Prong)
387 {
388 }
389 
394 {
395  f2ProngInvMass = source.fInvMass2Prong;
396  fDeltaInvMass = source.fD.M() - source.fInvMass2Prong;
398 }
399 
402 {
404 
405  f2ProngInvMass = 0;
406  fDeltaInvMass = 0;
407 }
408 
409 // Definitions of class AliAnalysisTaskDmesonJets::AliJetDefinition
410 
414 
417  TObject(),
418  fJetType(AliJetContainer::kChargedJet),
419  fRadius(0),
420  fJetAlgo(AliJetContainer::antikt_algorithm),
421  fRecoScheme(AliJetContainer::pt_scheme),
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 
442  TObject(),
443  fJetType(type),
444  fRadius(r),
445  fJetAlgo(algo),
446  fRecoScheme(reco),
447  fMinJetPt(0.),
448  fMaxJetPt(500.),
449  fMinJetPhi(0.),
450  fMaxJetPhi(0.),
451  fMinJetEta(-1.),
452  fMaxJetEta(1.),
453  fMinChargedPt(0.),
454  fMaxChargedPt(0.),
455  fMinNeutralPt(0.),
456  fMaxNeutralPt(0.)
457 {
458 }
459 
464  TObject(),
465  fJetType(source.fJetType),
466  fRadius(source.fRadius),
467  fJetAlgo(source.fJetAlgo),
468  fRecoScheme(source.fRecoScheme),
469  fMinJetPt(source.fMinJetPt),
470  fMaxJetPt(source.fMaxJetPt),
471  fMinJetPhi(source.fMinJetPhi),
472  fMaxJetPhi(source.fMaxJetPhi),
473  fMinJetEta(source.fMinJetEta),
474  fMaxJetEta(source.fMaxJetEta),
475  fMinChargedPt(source.fMinChargedPt),
476  fMaxChargedPt(source.fMaxChargedPt),
477  fMinNeutralPt(source.fMinNeutralPt),
478  fMaxNeutralPt(source.fMaxNeutralPt)
479 {
480 }
481 
486 {
487  new (this) AliHFJetDefinition(source);
488  return *this;
489 }
490 
493 {
494  static TString name;
495 
496  name = AliJetContainer::GenerateJetName(fJetType, fJetAlgo, fRecoScheme, fRadius, 0, 0, "Jet");
497 
498  return name.Data();
499 }
500 
506 {
507  if (fMinJetEta < fMaxJetEta && (jet.Eta() < fMinJetEta || jet.Eta() > fMaxJetEta)) return kFALSE;
508  if (fMinJetPhi < fMaxJetPhi && (jet.Phi() < fMinJetPhi || jet.Phi() > fMaxJetPhi)) return kFALSE;
509  if (jet.Pt() > fMaxJetPt || jet.Pt() < fMinJetPt) return kFALSE;
510  if (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt) return kFALSE;
511  if (jet.fMaxNeutralPt < fMinNeutralPt || jet.fMaxNeutralPt > fMaxNeutralPt) return kFALSE;
512 
513  return kTRUE;
514 }
515 
521 {
522  const AliJetInfo* jet = dMesonJet.GetJet(n);
523  if (!jet) return kFALSE;
524  return IsJetInAcceptance((*jet));
525 }
526 
533 {
534  if (lhs.fJetType > rhs.fJetType) return false;
535  else if (lhs.fJetType < rhs.fJetType) return true;
536  else {
537  if (lhs.fRadius > rhs.fRadius) return false;
538  else if (lhs.fRadius < rhs.fRadius) return true;
539  else {
540  if (lhs.fJetAlgo > rhs.fJetAlgo) return false;
541  else if (lhs.fJetAlgo < rhs.fJetAlgo) return true;
542  else {
543  if (lhs.fRecoScheme < rhs.fRecoScheme) return true;
544  else return false;
545  }
546  }
547  }
548 }
549 
556 {
557  if (lhs.fJetType != rhs.fJetType) return false;
558  if (lhs.fRadius != rhs.fRadius) return false;
559  if (lhs.fJetAlgo != rhs.fJetAlgo) return false;
560  if (lhs.fRecoScheme != rhs.fRecoScheme) return false;
561  return true;
562 }
563 
564 // Definitions of class AliAnalysisTaskDmesonJets::AnalysisEngine
565 
569 
572  TObject(),
573  fCandidateType(kD0toKpi),
574  fCandidateName(),
575  fCandidatePDG(0),
576  fNDaughters(0),
577  fPDGdaughters(),
578  fBranchName(),
579  fMCMode(kNoMC),
580  fNMassBins(0),
581  fMinMass(0),
582  fMaxMass(0),
583  fRDHFCuts(0),
584  fRejectedOrigin(0),
585  fAcceptedDecay(0),
586  fInhibit(kFALSE),
587  fJetDefinitions(),
588  fPtBinWidth(0.5),
589  fMaxPt(100),
590  fDataSlotNumber(-1),
591  fTree(0),
592  fCurrentDmesonJetInfo(0),
593  fCurrentJetInfo(0),
594  fCandidateArray(0),
595  fMCContainer(0),
596  fTrackContainer(0),
597  fClusterContainer(0),
598  fAodEvent(0),
599  fFastJetWrapper(0),
600  fHistManager(0)
601 {
602 }
603 
612  TObject(),
613  fCandidateType(type),
614  fCandidateName(),
615  fCandidatePDG(0),
616  fNDaughters(0),
617  fPDGdaughters(),
618  fBranchName(),
619  fMCMode(MCmode),
620  fNMassBins(nMassBins),
621  fMinMass(0),
622  fMaxMass(0),
623  fRDHFCuts(cuts),
625  fAcceptedDecay(kAnyDecay),
626  fInhibit(kFALSE),
627  fJetDefinitions(),
628  fPtBinWidth(0.5),
629  fMaxPt(100),
630  fDataSlotNumber(-1),
631  fTree(0),
632  fCurrentDmesonJetInfo(0),
633  fCurrentJetInfo(0),
634  fCandidateArray(0),
635  fMCContainer(0),
636  fTrackContainer(0),
637  fClusterContainer(0),
638  fAodEvent(0),
639  fFastJetWrapper(0),
640  fHistManager(0)
641 {
642  SetCandidateProperties(range);
643 }
644 
649  TObject(source),
650  fCandidateType(source.fCandidateType),
651  fCandidateName(source.fCandidateName),
652  fCandidatePDG(source.fCandidatePDG),
653  fNDaughters(source.fNDaughters),
654  fPDGdaughters(source.fPDGdaughters),
655  fBranchName(source.fBranchName),
656  fMCMode(source.fMCMode),
657  fNMassBins(source.fNMassBins),
658  fMinMass(source.fMinMass),
659  fMaxMass(source.fMaxMass),
660  fRDHFCuts(),
661  fRejectedOrigin(source.fRejectedOrigin),
662  fAcceptedDecay(source.fAcceptedDecay),
663  fInhibit(source.fInhibit),
664  fJetDefinitions(source.fJetDefinitions),
665  fPtBinWidth(source.fPtBinWidth),
666  fMaxPt(source.fMaxPt),
667  fDataSlotNumber(-1),
668  fTree(0),
669  fCurrentDmesonJetInfo(0),
670  fCurrentJetInfo(0),
671  fCandidateArray(source.fCandidateArray),
672  fMCContainer(source.fMCContainer),
673  fTrackContainer(source.fTrackContainer),
674  fClusterContainer(source.fClusterContainer),
675  fAodEvent(source.fAodEvent),
677  fHistManager(source.fHistManager)
678 {
679  SetRDHFCuts(source.fRDHFCuts);
680 }
681 
682 // Destructor
684 {
685  delete fRDHFCuts;
686 }
687 
692 {
693  new (this) AnalysisEngine(source);
694  return *this;
695 }
696 
701 {
702  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
703  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
704  }
705 
706  return kFALSE;
707 }
708 
710 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const /*geom*/, Int_t /*runNumber*/)
711 {
712 }
713 
718 {
719  switch (fCandidateType) {
720  case kD0toKpi:
721  fCandidatePDG = 421;
722  fCandidateName = "D0";
723  fNDaughters = 2;
724  fPDGdaughters.Set(fNDaughters);
725  fPDGdaughters.Reset();
726  fPDGdaughters[0] = 211; // pi
727  fPDGdaughters[1] = 321; // K
728  fBranchName = "D0toKpi";
729  fAcceptedDecay = kDecayD0toKpi;
730  if (!fRDHFCuts) {
731  fRDHFCuts = new AliRDHFCutsD0toKpi();
732  fRDHFCuts->SetStandardCutsPP2010();
733  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
734  fRDHFCuts->SetTriggerClass("","");
735  }
736  break;
737  case kD0toKpiLikeSign:
738  fCandidatePDG = 421;
739  fCandidateName = "2ProngLikeSign";
740  fNDaughters = 2;
741  fPDGdaughters.Set(fNDaughters);
742  fPDGdaughters.Reset();
743  fPDGdaughters[0] = 211; // pi
744  fPDGdaughters[1] = 321; // K
745  fBranchName = "LikeSign2Prong";
746  fAcceptedDecay = kAnyDecay;
747  if (!fRDHFCuts) {
748  fRDHFCuts = new AliRDHFCutsD0toKpi();
749  fRDHFCuts->SetStandardCutsPP2010();
750  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
751  fRDHFCuts->SetTriggerClass("","");
752  }
753  break;
754  case kDstartoKpipi:
755  fCandidatePDG = 413;
756  fCandidateName = "DStar";
757  fNDaughters = 3;
758  fPDGdaughters.Set(fNDaughters);
759  fPDGdaughters.Reset();
760  fPDGdaughters[0] = 211; // pi soft
761  fPDGdaughters[1] = 211; // pi fromD0
762  fPDGdaughters[2] = 321; // K from D0
763  fBranchName = "Dstar";
764  fAcceptedDecay = kDecayDStartoKpipi;
765  if (!fRDHFCuts) {
766  fRDHFCuts = new AliRDHFCutsDStartoKpipi();
767  fRDHFCuts->SetStandardCutsPP2010();
768  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
769  fRDHFCuts->SetTriggerClass("","");
770  }
771  break;
772  default:
773  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
774  }
775 
776  CalculateMassLimits(range, fCandidatePDG, fNMassBins, fMinMass, fMaxMass);
777 }
778 
783 {
784  if (fRDHFCuts) delete fRDHFCuts;
785  fRDHFCuts = cuts;
786 }
787 
792 {
793  if (!cuts) return;
794  if (fRDHFCuts) delete fRDHFCuts;
795  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
796 }
797 
802 {
803  static TString name;
804 
805  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
806 
807  return name.Data();
808 }
809 
814 {
815  static TString name;
816 
817  name = fCandidateName;
818  switch (fMCMode) {
819  case kBackgroundOnly:
820  name += "_kBackgroundOnly";
821  break;
822  case kSignalOnly:
823  name += "_kSignalOnly";
824  break;
825  case kMCTruth:
826  name += "_MCTruth";
827  break;
828  default:
829  break;
830  }
831 
832  return name.Data();
833 }
834 
842 {
843  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
844 
845  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
846  fJetDefinitions.push_back(def);
847  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
848  "Total number of jet definitions is now %lu.",
849  def.GetName(), GetName(), fJetDefinitions.size());
850  // For detector level set maximum track pt to 100 GeV/c
851  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
852  }
853  else {
854  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
855  }
856 
857  return &(*it);
858 }
859 
871 {
872  AliHFJetDefinition def(type, r, algo, reco);
873 
874  return AddJetDefinition(def);
875 }
876 
882 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
883 {
884  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
885  while (it != fJetDefinitions.end() && (*it) != def) it++;
886  return it;
887 }
888 
895 {
896  if (lhs.fCandidateType > rhs.fCandidateType) return false;
897  else if (lhs.fCandidateType < rhs.fCandidateType) return true;
898  else {
899  if (lhs.fMCMode < rhs.fMCMode) return true;
900  else return false;
901  }
902 }
903 
910 {
911  if (lhs.fCandidateType != rhs.fCandidateType) return false;
912  if (lhs.fMCMode != rhs.fMCMode) return false;
913  return true;
914 }
915 
924 {
925  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) { // D0 candidate
926  return ExtractD0Attributes(Dcand, DmesonJet, i);
927  }
928  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
929  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
930  }
931  else {
932  return kFALSE;
933  }
934 }
935 
944 {
945  //AliDebug(2,"Checking if D0 meson is selected");
946  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
947  if (isSelected == 0) return kFALSE;
948 
949  Int_t MCtruthPdgCode = 0;
950 
951  Double_t invMassD = 0;
952 
953  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
954  // Checks also the origin, and if it matches the rejected origin mask, return false
955  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
956  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
957  DmesonJet.fMCLabel = mcLab;
958 
959  // Retrieve the generated particle (if exists) and its PDG code
960  if (mcLab >= 0) {
961  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
962 
963  if (aodMcPart) {
964  // Check origin and return false if it matches the rejected origin mask
965  if (fRejectedOrigin) {
966  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
967 
968  UInt_t rs = origin;
969  UShort_t p = 0;
970  while (rs >>= 1) { p++; }
971  TString hname = TString::Format("%s/fHistDmesonOrigin", GetName());
972  fHistManager->FillTH1(hname, p);
973 
974  if ((origin & fRejectedOrigin) == origin) return kFALSE;
975  }
976  MCtruthPdgCode = aodMcPart->PdgCode();
977  }
978  }
979  }
980 
981  if (isSelected == 1) { // selected as a D0
982  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
983 
984  if (fMCMode == kNoMC ||
985  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
986  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly)) {
987  // 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)
988  //AliDebug(2,"Selected as D0");
989  invMassD = Dcand->InvMassD0();
990  }
991  else { // conditions above not passed, so return FALSE
992  return kFALSE;
993  }
994  }
995  else if (isSelected == 2) { // selected as a D0bar
996  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
997 
998  if (fMCMode == kNoMC ||
999  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1000  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly)) {
1001  // 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)
1002  //AliDebug(2,"Selected as D0bar");
1003  invMassD = Dcand->InvMassD0bar();
1004  }
1005  else { // conditions above not passed, so return FALSE
1006  return kFALSE;
1007  }
1008  }
1009  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
1010  //AliDebug(2,"Selected as either D0 or D0bar");
1011 
1012  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
1013  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1014  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kBackgroundOnly)) {
1015  if (i > 0) return kFALSE;
1016  //AliDebug(2, "MC truth is D0");
1017  invMassD = Dcand->InvMassD0();
1018  }
1019  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1020  (MCtruthPdgCode == fCandidatePDG && fMCMode == kBackgroundOnly)) {
1021  if (i > 0) return kFALSE;
1022  //AliDebug(2, "MC truth is D0bar");
1023  invMassD = Dcand->InvMassD0bar();
1024  }
1025  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
1026 
1027  // Only accept it if background-only OR background-and-signal was requested
1028  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
1029  // Select D0 or D0bar depending on the i-parameter
1030  if (i == 0) {
1031  //AliDebug(2, "Returning invariant mass with D0 hypothesis");
1032  invMassD = Dcand->InvMassD0();
1033  }
1034  else if (i == 1) {
1035  //AliDebug(2, "Returning invariant mass with D0bar hypothesis");
1036  invMassD = Dcand->InvMassD0bar();
1037  }
1038  else { // i > 1
1039  return kFALSE;
1040  }
1041  }
1042  else { // signal-only was requested but this is not a true D0
1043  return kFALSE;
1044  }
1045  }
1046  }
1047  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
1048  return kTRUE;
1049 }
1050 
1059 {
1060  //AliDebug(2,"Checking if D0 meson is selected");
1061  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoCascadeHF*>(DstarCand), AliRDHFCuts::kAll, fAodEvent);
1062  if (isSelected == 0) return kFALSE;
1063 
1064  if (i > 0) return kFALSE; // only one mass hypothesis for the D*
1065 
1066  Int_t MCtruthPdgCode = 0;
1067 
1068  Double_t invMassD = 0;
1069 
1070  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1071  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1072  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1073 
1074  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1075  //AliDebug(2, Form("MC label is %d", mcLab));
1076  DmesonJet.fMCLabel = mcLab;
1077  if (mcLab >= 0) {
1078  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1079 
1080  if (aodMcPart) {
1081  if (fRejectedOrigin) {
1082  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
1083 
1084  UInt_t rs = origin;
1085  UShort_t p = 0;
1086  while (rs >>= 1) { p++; }
1087  TString hname = TString::Format("%s/fHistDmesonOrigin", GetName());
1088  fHistManager->FillTH1(hname, p);
1089 
1090  if ((origin & fRejectedOrigin) == origin) return kFALSE;
1091  }
1092 
1093  MCtruthPdgCode = aodMcPart->PdgCode();
1094  //AliDebug(2, Form("MC truth pdg code is %d",MCtruthPdgCode));
1095  }
1096  }
1097  }
1098 
1099  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1100  if (fMCMode == kNoMC ||
1101  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1102  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1103  // 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)
1104  invMassD = DstarCand->InvMassDstarKpipi();
1105  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1106  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1107  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1108  return kTRUE;
1109  }
1110  else { // conditions above not passed, so return FALSE
1111  return kFALSE;
1112  }
1113 }
1114 
1122 {
1123  if (!part) return kUnknownDecay;
1124  if (!mcArray) return kUnknownDecay;
1125 
1127 
1128  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1129 
1130  if (part->GetNDaughters() == 2) {
1131 
1132  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1133  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1134 
1135  if (!d1 || !d2) {
1136  return decay;
1137  }
1138 
1139  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1140  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1141 
1142  if (absPdgPart == 421) { // D0 -> K pi
1143  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1144  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1145  decay = kDecayD0toKpi;
1146  }
1147  }
1148 
1149  if (absPdgPart == 413) { // D* -> D0 pi
1150  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1151  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1152  if (D0decay == kDecayD0toKpi) {
1153  decay = kDecayDStartoKpipi;
1154  }
1155  }
1156  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1157  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1158  if (D0decay == kDecayD0toKpi) {
1159  decay = kDecayDStartoKpipi;
1160  }
1161  }
1162  }
1163  }
1164 
1165  return decay;
1166 }
1167 
1175 {
1176  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1177 
1178  if (!part) return kUnknownQuark;
1179  if (!mcArray) return kUnknownQuark;
1180 
1181  Int_t mother = part->GetMother();
1182  while (mother >= 0) {
1183  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1184  if (mcGranma) {
1185  Int_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
1186 
1187  if (abspdgGranma == 1) return kFromDown;
1188  if (abspdgGranma == 2) return kFromUp;
1189  if (abspdgGranma == 3) return kFromStrange;
1190  if (abspdgGranma == 4) return kFromCharm;
1191  if (abspdgGranma == 5) return kFromBottom;
1192  if (abspdgGranma == 6) return kFromTop;
1193 
1194  mother = mcGranma->GetMother();
1195  }
1196  else {
1197  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1198  break;
1199  }
1200  }
1201 
1202  return kUnknownQuark;
1203 }
1204 
1207 {
1208  fDmesonJets.clear();
1209 
1210  for (auto& jetDef : fJetDefinitions) {
1211  jetDef.fJets.clear();
1212  }
1213 
1214  if (fMCMode == kMCTruth) {
1215  RunParticleLevelAnalysis();
1216  }
1217  else {
1218  RunDetectorLevelAnalysis();
1219  }
1220 }
1221 
1227 {
1229  fFastJetWrapper->SetR(jetDef.fRadius);
1232 
1233  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1234  fTrackContainer->SetDMesonCandidate(0);
1235  AddInputVectors(fTrackContainer, 100);
1236  }
1237 
1238  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1239  AddInputVectors(fClusterContainer, -100);
1240  }
1241 
1242  // run jet finder
1243  fFastJetWrapper->Run();
1244 
1245  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1246 
1247  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1248  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1249 
1250  Double_t maxChPt = 0;
1251  Double_t maxNePt = 0;
1252  Double_t totalNeutralPt = 0;
1253 
1254  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1255  if (constituents[ic].user_index() >= 100) {
1256  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1257  }
1258  else if (constituents[ic].user_index() <= -100) {
1259  totalNeutralPt += constituents[ic].pt();
1260  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1261  }
1262  }
1263 
1264  jetDef.fJets.push_back(
1265  AliJetInfo(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E(),
1266  constituents.size(), totalNeutralPt / jets_incl[ijet].pt(), maxChPt, maxNePt)
1267  );
1268  }
1269 }
1270 
1279 {
1280  if (jetDef.fJets.size() == 0) FindJets(jetDef);
1281 
1282  Double_t d_closest = 999;
1283  AliJetInfo* jet_closest = 0;
1284  const AliJetInfo* truth_jet = 0;
1285  try {
1286  truth_jet = &(dmeson.fJets.at(jetDef.GetName()));
1287  }
1288  catch(...) {
1289  return jet_distance_pair(0, 999);
1290  }
1291 
1292  for (auto& jet : jetDef.fJets) {
1293  Double_t d = truth_jet->GetDistance(jet);
1294 
1295  if (d > dMax) continue;
1296  if (d < d_closest) {
1297  d_closest = d;
1298  jet_closest = &jet;
1299  }
1300  }
1301 
1302  if (jet_closest && applyKinCuts) {
1303  if (!jetDef.IsJetInAcceptance(*jet_closest)) jet_closest = 0;
1304  }
1305 
1306  if (jet_closest) {
1307  AliDebug(2, Form("Found closest jet (pt=%.3f, eta=%.3f, phi=%.3f) to d meson (pt=%.3f, eta=%.3f, phi=%.3f) with d = %.3f",
1308  jet_closest->Pt(), jet_closest->Eta(), jet_closest->Phi_0_2pi(),
1309  dmeson.fD.Pt(), dmeson.fD.Eta(), dmeson.fD.Phi_0_2pi(),
1310  d_closest));
1311  }
1312 
1313  return jet_distance_pair(jet_closest, d_closest);
1314 }
1315 
1318 {
1319  const Int_t nD = fCandidateArray->GetEntriesFast();
1320 
1321  AliDmesonJetInfo DmesonJet;
1322 
1323  Int_t nAccCharm = 0;
1324  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1325  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1326  if (!charmCand) continue;
1327 
1328  // region of interest + cuts
1329  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1330 
1331  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1332  DmesonJet.Reset();
1333  DmesonJet.fDmesonParticle = charmCand;
1334  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1335  for (std::vector<AliHFJetDefinition>::iterator itdef = fJetDefinitions.begin(); itdef != fJetDefinitions.end(); itdef++) {
1336  if (!FindJet(charmCand, DmesonJet, *itdef)) {
1337  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1338  (*itdef).GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1339  }
1340  }
1341  fDmesonJets[icharm] = DmesonJet;
1342  }
1343  }
1344  nAccCharm++;
1345  } // end of D cand loop
1346 
1347  TString hname;
1348 
1349  hname = TString::Format("%s/fHistNAcceptedDmesons", GetName());
1350  fHistManager->FillTH1(hname, nAccCharm);
1351 
1352  hname = TString::Format("%s/fHistNDmesons", GetName());
1353  fHistManager->FillTH1(hname, nD);
1354 }
1355 
1366 {
1367  TString hname;
1368 
1370  fFastJetWrapper->SetR(jetDef.fRadius);
1373 
1374  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1375 
1376  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1377  fTrackContainer->SetDMesonCandidate(Dcand);
1378  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1379  AddInputVectors(fTrackContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1380 
1381  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1382  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1383  const TObjArray& daughters = fTrackContainer->GetDaughterList();
1384  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1385  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1386  if (!fTrackContainer->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1387  }
1388  }
1389 
1390  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1391  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1392  AddInputVectors(fClusterContainer, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1393  }
1394 
1395  // run jet finder
1396  fFastJetWrapper->Run();
1397 
1398  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1399 
1400  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1401  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1402 
1403  Bool_t isDmesonJet = kFALSE;
1404 
1405  Double_t maxChPt = 0;
1406  Double_t maxNePt = 0;
1407  Double_t totalNeutralPt = 0;
1408 
1409  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1410  if (constituents[ic].user_index() == 0) {
1411  isDmesonJet = kTRUE;
1412  }
1413  else if (constituents[ic].user_index() >= 100) {
1414  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1415  }
1416  else if (constituents[ic].user_index() <= -100) {
1417  totalNeutralPt += constituents[ic].pt();
1418  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1419  }
1420  }
1421 
1422  if (isDmesonJet) {
1423  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1424  DmesonJet.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1425  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1426  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1427  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1428 
1429  return kTRUE;
1430  }
1431  }
1432 
1433  return kFALSE;
1434 }
1435 
1439 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist)
1440 {
1441  AliEmcalIterableMomentumContainer itcont = cont->all_momentum();
1442  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1443  UInt_t rejectionReason = 0;
1444  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1445  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1446  continue;
1447  }
1448  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1449  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1450  }
1451 }
1452 
1455 {
1456  TString hname;
1457 
1458  fMCContainer->SetSpecialPDG(fCandidatePDG);
1459  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1460  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1461 
1462  if (!fMCContainer->IsSpecialPDGFound()) return;
1463 
1464  for (auto &jetDef : fJetDefinitions) {
1465  switch (jetDef.fJetType) {
1467  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
1468  break;
1470  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
1471  break;
1473  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
1474  break;
1475  }
1476 
1478  fFastJetWrapper->SetR(jetDef.fRadius);
1481 
1482  hname = TString::Format("%s/fHistDmesonOrigin", GetName());
1483  fMCContainer->SetHistOrigin(static_cast<TH1*>(fHistManager->FindObject(hname)));
1484  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1485  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1486  fMCContainer->SetHistOrigin(0);
1487 
1488  fFastJetWrapper->Run();
1489 
1490  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1491 
1492  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1493  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1494 
1495  Bool_t isDmesonJet = kFALSE;
1496 
1497  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1498  Int_t iPart = constituents[ic].user_index() - 100;
1499  AliVParticle* part = fMCContainer->GetParticle(iPart);
1500  if (!part) {
1501  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1502  continue;
1503  }
1504  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1505  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1506  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1507  std::pair<int, AliDmesonJetInfo> element;
1508  element.first = iPart;
1509  element.second = AliDmesonJetInfo();
1510  dMesonJetIt = fDmesonJets.insert(element).first;
1511  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1512  (*dMesonJetIt).second.fDmesonParticle = part;
1513  }
1514 
1515  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1516  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1517  }
1518  }
1519  }
1520  }
1521 }
1522 
1527 {
1528  TString classname;
1529  switch (fCandidateType) {
1530  case kD0toKpi:
1531  case kD0toKpiLikeSign:
1532  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1533  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1534  break;
1535  case kDstartoKpipi:
1536  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1537  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1538  break;
1539  }
1540  TString treeName = TString::Format("%s_%s", taskName, GetName());
1541  fTree = new TTree(treeName, treeName);
1542  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
1543  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1544  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1545  fCurrentJetInfo[i] = new AliJetInfoSummary();
1546  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1547  }
1548 
1549  return fTree;
1550 }
1551 
1556 {
1557  TString hname;
1558 
1559  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
1560 
1561  for (auto &jetDef : fJetDefinitions) {
1562 
1563  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
1564 
1565  Double_t radius = jetDef.fRadius;
1566 
1567  TString title[30] = {""};
1568  Int_t nbins[30] = {0 };
1569  Double_t min [30] = {0.};
1570  Double_t max [30] = {0.};
1571  Int_t dim = 0 ;
1572 
1573  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1574  nbins[dim] = nPtBins;
1575  min[dim] = 0;
1576  max[dim] = fMaxPt;
1577  dim++;
1578 
1579  if ((enabledAxis & kPositionD) != 0) {
1580  title[dim] = "#eta_{D}";
1581  nbins[dim] = 50;
1582  min[dim] = -1;
1583  max[dim] = 1;
1584  dim++;
1585 
1586  title[dim] = "#phi_{D} (rad)";
1587  nbins[dim] = 150;
1588  min[dim] = 0;
1589  max[dim] = TMath::TwoPi();
1590  dim++;
1591  }
1592 
1593  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1594  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1595  nbins[dim] = fNMassBins;
1596  min[dim] = fMinMass;
1597  max[dim] = fMaxMass;
1598  dim++;
1599  }
1600 
1601  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1602  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1603  nbins[dim] = fNMassBins;
1604  min[dim] = fMinMass;
1605  max[dim] = fMaxMass;
1606  dim++;
1607  }
1608 
1609  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1610  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1611  nbins[dim] = fNMassBins;
1612  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1613  dim++;
1614  }
1615 
1616  if (fCandidateType == kDstartoKpipi) {
1617  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1618  nbins[dim] = fNMassBins*6;
1619  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1620 
1621  // subtract mass of D0
1622  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1623  min[dim] -= D0mass;
1624  max[dim] -= D0mass;
1625 
1626  dim++;
1627  }
1628 
1629  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1630  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1631  nbins[dim] = 100;
1632  min[dim] = 0;
1633  max[dim] = 25;
1634  dim++;
1635  }
1636 
1637  title[dim] = "#it{z}_{D}";
1638  nbins[dim] = 110;
1639  min[dim] = 0;
1640  max[dim] = 1.10;
1641  dim++;
1642 
1643  if ((enabledAxis & kDeltaR) != 0) {
1644  title[dim] = "#Delta R_{D-jet}";
1645  nbins[dim] = 100;
1646  min[dim] = 0;
1647  max[dim] = radius * 1.5;
1648  dim++;
1649  }
1650 
1651  if ((enabledAxis & kDeltaEta) != 0) {
1652  title[dim] = "#eta_{D} - #eta_{jet}";
1653  nbins[dim] = 100;
1654  min[dim] = -radius * 1.2;
1655  max[dim] = radius * 1.2;
1656  dim++;
1657  }
1658 
1659  if ((enabledAxis & kDeltaPhi) != 0) {
1660  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1661  nbins[dim] = 100;
1662  min[dim] = -radius * 1.2;
1663  max[dim] = radius * 1.2;
1664  dim++;
1665  }
1666 
1667  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1668  nbins[dim] = nPtBins;
1669  min[dim] = 0;
1670  max[dim] = fMaxPt;
1671  dim++;
1672 
1673  if ((enabledAxis & kPositionJet) != 0) {
1674  title[dim] = "#eta_{jet}";
1675  nbins[dim] = 50;
1676  min[dim] = -1;
1677  max[dim] = 1;
1678  dim++;
1679 
1680  title[dim] = "#phi_{jet} (rad)";
1681  nbins[dim] = 150;
1682  min[dim] = 0;
1683  max[dim] = TMath::TwoPi();
1684  dim++;
1685  }
1686 
1687  if ((enabledAxis & kJetConstituents) != 0) {
1688  title[dim] = "No. of constituents";
1689  nbins[dim] = 50;
1690  min[dim] = -0.5;
1691  max[dim] = 49.5;
1692  dim++;
1693  }
1694 
1695  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1696  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1697  for (Int_t j = 0; j < dim; j++) {
1698  h->GetAxis(j)->SetTitle(title[j]);
1699  }
1700  }
1701 }
1702 
1707 {
1708  TString hname;
1709 
1710  for (auto& dmeson_pair : fDmesonJets) {
1711  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
1712  Int_t accJets = 0;
1713  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1714  fCurrentJetInfo[ij]->Reset();
1715  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1716  if (!jet) continue;
1717  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1718  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1719  fHistManager->FillTH1(hname, jet->Pt());
1720  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1721  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1722  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1723  fHistManager->FillTH1(hname, jet->Eta());
1724  continue;
1725  }
1726  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
1727  accJets++;
1728  }
1729  if (accJets > 0) {
1730  fTree->Fill();
1731  }
1732  else {
1733  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1734  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1735  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1736  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1737  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1738  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1739  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1740  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1741  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1742  }
1743  else if (fCandidateType == kDstartoKpipi) {
1744  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1745  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1746 
1747  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1748  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1749  }
1750  }
1751  }
1752  return kTRUE;
1753 }
1754 
1760 {
1761  TString hname;
1762 
1763  for (auto& dmeson_pair : fDmesonJets) {
1764  Int_t accJets = 0;
1765  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1766  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1767  if (!jet) continue;
1768  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1769  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1770  fHistManager->FillTH1(hname, jet->Pt());
1771  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1772  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1773  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1774  fHistManager->FillTH1(hname, jet->Eta());
1775  continue;
1776  }
1777  accJets++;
1778  }
1779  if (accJets == 0) {
1780  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1781  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1782  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1783  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1784  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1785  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1786  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1787  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1788  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1789  }
1790  else if (fCandidateType == kDstartoKpipi) {
1791  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1792  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1793 
1794  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1795  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1796  }
1797  }
1798  }
1799  return kTRUE;
1800 }
1801 
1806 {
1807  TString hname;
1808 
1809  for (auto& dmeson_pair : fDmesonJets) {
1810  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
1811  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1812  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1813  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1814  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1815  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1816  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1817  }
1818  }
1819 
1820  for (auto &jetDef : fJetDefinitions) {
1821 
1822  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1823  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
1824 
1825  for (auto& dmeson_pair : fDmesonJets) {
1826  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
1827  if (!jet) continue;
1828  if (!jetDef.IsJetInAcceptance(*jet)) {
1829  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
1830  fHistManager->FillTH1(hname, jet->Pt());
1831  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
1832  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1833  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
1834  fHistManager->FillTH1(hname, jet->Eta());
1835  continue;
1836  }
1837  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
1838  }
1839  }
1840 
1841  return kTRUE;
1842 }
1843 
1850 {
1851  // Fill the THnSparse histogram.
1852 
1853  Double_t contents[30] = {0.};
1854 
1855  Double_t z = DmesonJet.GetZ(n);
1856  Double_t deltaPhi = 0;
1857  Double_t deltaEta = 0;
1858  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
1859 
1860  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
1861  if (it == DmesonJet.fJets.end()) return kFALSE;
1862 
1863  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
1864  TString title(h->GetAxis(i)->GetTitle());
1865  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
1866  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
1867  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
1868  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
1869  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
1870  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
1871  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
1872  else if (title=="#it{z}_{D}") contents[i] = z;
1873  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
1874  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
1875  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
1876  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
1877  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
1878  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
1879  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
1880  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
1881  }
1882 
1883  h->Fill(contents);
1884 
1885  return kTRUE;
1886 }
1887 
1888 // Definitions of class AliAnalysisTaskDmesonJets
1889 
1893 
1897  fAnalysisEngines(),
1898  fEnabledAxis(0),
1900  fHistManager(),
1901  fApplyKinematicCuts(kTRUE),
1902  fNOutputTrees(0),
1903  fAodEvent(0),
1904  fFastJetWrapper(0)
1905 {
1906  SetMakeGeneralHistograms(kTRUE);
1907 }
1908 
1913  AliAnalysisTaskEmcalLight(name, kTRUE),
1914  fAnalysisEngines(),
1915  fEnabledAxis(k2ProngInvMass),
1916  fOutputType(kTreeOutput),
1917  fHistManager(name),
1918  fApplyKinematicCuts(kTRUE),
1919  fNOutputTrees(nOutputTrees),
1920  fAodEvent(0),
1921  fFastJetWrapper(0)
1922 {
1923  SetMakeGeneralHistograms(kTRUE);
1924  for (Int_t i = 0; i < nOutputTrees; i++){
1925  DefineOutput(2+i, TTree::Class());
1926  }
1927 }
1928 
1931 {
1932  if (fFastJetWrapper) delete fFastJetWrapper;
1933 }
1934 
1942 {
1943  AliRDHFCuts* analysiscuts = 0;
1944  TFile* filecuts = TFile::Open(cutfname);
1945  if (!filecuts || filecuts->IsZombie()) {
1946  ::Warning("AddTaskDmesonJets", "Input file not found: will use std cuts.");
1947  filecuts = 0;
1948  }
1949 
1950  if (filecuts) {
1951  analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
1952  if (!analysiscuts) {
1953  ::Warning("AddTaskDmesonJetCorr", "Could not find analysis cuts '%s' in '%s'. Using std cuts.", cutsname.Data(), cutfname.Data());
1954  }
1955  }
1956 
1957  return analysiscuts;
1958 }
1959 
1969 {
1971  return AddAnalysisEngine(type, MCmode, jetDef, cutfname);
1972 }
1973 
1983 {
1984  AliRDHFCuts* cuts = 0;
1985 
1986  if (!cutfname.IsNull()) {
1987  TString cutsname;
1988 
1989  switch (type) {
1990  case kD0toKpi:
1991  case kD0toKpiLikeSign:
1992  cutsname = "D0toKpiCuts";
1993  break;
1994  case kDstartoKpipi:
1995  cutsname = "DStartoKpipiCuts";
1996  break;
1997  default:
1998  return 0;
1999  }
2000 
2001  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
2002  }
2003 
2004  AnalysisEngine eng(type, MCmode, cuts);
2005 
2006  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
2007 
2008  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
2009  eng.AddJetDefinition(jetDef);
2010  it = fAnalysisEngines.insert(it, eng);
2011  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
2012  }
2013  else {
2014  AnalysisEngine* found_eng = &(*it);
2015  ::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());
2016  found_eng->AddJetDefinition(jetDef);
2017 
2018  if (cuts && found_eng->fRDHFCuts != 0) {
2019  ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
2020  found_eng->SetRDHFCuts(cuts);
2021  }
2022  }
2023 
2024  return &(*it);
2025 }
2026 
2027 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
2028 {
2029  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
2030  while (it != fAnalysisEngines.end() && (*it) < eng) it++;
2031  return it;
2032 }
2033 
2036 {
2037  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
2038 
2040 
2041  // Define histograms
2042  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
2043 
2044  TString hname;
2045  TString htitle;
2046  TH1* h = 0;
2047  Int_t treeSlot = 0;
2048 
2049  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
2050  for (auto &param : fAnalysisEngines) {
2051  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
2052 
2053  fHistManager.CreateHistoGroup(param.GetName());
2054 
2055  param.fHistManager = &fHistManager;
2056 
2057  hname = TString::Format("%s/fHistNAcceptedDmesons", param.GetName());
2058  htitle = hname + ";Number of D accepted meson candidates;counts";
2059  h = fHistManager.CreateTH1(hname, htitle, 51, -0.5, 50.5);
2060 
2061  hname = TString::Format("%s/fHistNDmesons", param.GetName());
2062  htitle = hname + ";Number of D meson candidates;counts";
2063  h = fHistManager.CreateTH1(hname, htitle, 101, -0.5, 100.5);
2064 
2065  hname = TString::Format("%s/fHistNEvents", param.GetName());
2066  htitle = hname + ";Event status;counts";
2067  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
2068  h->GetXaxis()->SetBinLabel(1, "Accepted");
2069  h->GetXaxis()->SetBinLabel(2, "Rejected");
2070 
2071  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
2072  htitle = hname + ";Rejection reason;counts";
2073  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
2074 
2075  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
2076  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
2077  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2078 
2079  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
2080  htitle = hname + ";#it{#eta}_{D};counts";
2081  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2082 
2083  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
2084  htitle = hname + ";#it{#phi}_{D};counts";
2085  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2086 
2087  if (param.fCandidateType == kD0toKpi || param.fCandidateType == kD0toKpiLikeSign) {
2088  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
2089  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2090  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
2091  }
2092  else if (param.fCandidateType == kDstartoKpipi) {
2093  Double_t min = 0;
2094  Double_t max = 0;
2095 
2096  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
2097  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2098  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
2099  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
2100 
2101  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2102  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
2103  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2104  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
2105  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
2106  }
2107 
2108  if (param.fMCMode == kBackgroundOnly || param.fMCMode == kSignalOnly || param.fMCMode == kMCTruth) {
2109  hname = TString::Format("%s/fHistDmesonOrigin", param.GetName());
2110  htitle = hname + ";origin;counts";
2111  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2112  }
2113 
2114  for (std::vector<AliHFJetDefinition>::iterator itdef = param.fJetDefinitions.begin(); itdef != param.fJetDefinitions.end(); itdef++) {
2115  AliHFJetDefinition* jetDef = &(*itdef);
2116  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef->GetName());
2117 
2118  fHistManager.CreateHistoGroup(jetDef->GetName(), param.GetName());
2119 
2120  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef->GetName());
2121  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2122  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2123  SetRejectionReasonLabels(h->GetXaxis());
2124 
2125  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef->GetName());
2126  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2127  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2128  SetRejectionReasonLabels(h->GetXaxis());
2129 
2130  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef->GetName());
2131  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
2132  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2133  SetRejectionReasonLabels(h->GetXaxis());
2134 
2135  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef->GetName());
2136  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
2137  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
2138 
2139  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef->GetName());
2140  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
2141  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2142 
2143  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef->GetName());
2144  htitle = hname + ";#it{#eta}_{jet};counts";
2145  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2146 
2147  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef->GetName());
2148  htitle = hname + ";#it{#phi}_{jet};counts";
2149  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2150  }
2151  switch (fOutputType) {
2152  case kTreeOutput:
2153  param.BuildTree(GetName());
2154  if (treeSlot < fNOutputTrees) {
2155  param.AssignDataSlot(treeSlot+2);
2156  treeSlot++;
2158  }
2159  else {
2160  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
2161  }
2162  break;
2163  case kTHnOutput:
2164  param.BuildHnSparse(fEnabledAxis);
2165  break;
2166  case kNoOutput:
2167  break;
2168  }
2169  }
2170 
2172 
2173  PostData(1, fOutput);
2174 }
2175 
2179 {
2181 
2182  // Load the event
2183  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
2184 
2185  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
2186 
2187  fFastJetWrapper->SetAreaType(fastjet::active_area);
2189 
2190  if (!fAodEvent) {
2191  AliError(Form("This task need an AOD event (Task '%s'). Expect troubles...", GetName()));
2192  //return;
2193  }
2194 
2195  for (auto &params : fAnalysisEngines) {
2196 
2197  params.fAodEvent = fAodEvent;
2198  params.fFastJetWrapper = fFastJetWrapper;
2199  if (fAodEvent) params.Init(fGeom, fAodEvent->GetRunNumber());
2200 
2201  if (params.fMCMode != kMCTruth && fAodEvent) {
2202  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
2203 
2204  if (params.fCandidateArray) {
2205  TString className;
2206  if (params.fCandidateType == kD0toKpi || params.fCandidateType == kD0toKpiLikeSign) {
2207  className = "AliAODRecoDecayHF2Prong";
2208  }
2209  else if (params.fCandidateType == kDstartoKpipi) {
2210  className = "AliAODRecoCascadeHF";
2211  }
2212  if (!params.fCandidateArray->GetClass()->InheritsFrom(className)) {
2213  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2214  "%s: Objects of type %s in %s are not inherited from %s! Task will be disabled!",
2215  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName(), className.Data());
2216  params.fCandidateArray = 0;
2217  params.fInhibit = kTRUE;
2218  }
2219  }
2220  else {
2221  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2222  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
2223  params.fBranchName.Data(), params.GetName());
2224  params.fInhibit = kTRUE;
2225  }
2226  }
2227 
2228  if (params.fMCMode != kNoMC) {
2229  params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(0));
2230 
2231  if (!params.fMCContainer) params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(1));
2232 
2233  if (!params.fMCContainer) {
2234  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2235  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
2236  params.GetName());
2237  params.fInhibit = kTRUE;
2238  }
2239  }
2240 
2241  if (params.fMCMode != kMCTruth) {
2242  params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(0));
2243  if (!params.fTrackContainer) params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(1));
2244 
2245  params.fClusterContainer = GetClusterContainer(0);
2246 
2247  if (!params.fTrackContainer && !params.fClusterContainer) {
2248  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2249  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
2250  params.GetName());
2251  params.fInhibit = kTRUE;
2252  }
2253  }
2254  }
2255 }
2256 
2261 {
2262  //if (!fAodEvent) return kFALSE;
2263 
2264  TString hname;
2265 
2266  // fix for temporary bug in ESDfilter
2267  // the AODs with null vertex pointer didn't pass the PhysSel
2268  // Now adding an entry in teh histogram so as to check that this is actually cutting anything out
2269  if (fAodEvent && (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001)) {
2270  for (auto &eng : fAnalysisEngines) {
2271  if (eng.fInhibit) continue;
2272  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2273  fHistManager.FillTH1(hname, "ESDfilterBug");
2274  }
2275  return kFALSE;
2276  }
2277 
2278  for (auto &eng : fAnalysisEngines) {
2279 
2280  if (eng.fInhibit) continue;
2281 
2282  //Event selection
2283  hname = TString::Format("%s/fHistNEvents", eng.GetName());
2284  if (fAodEvent) {
2285  Bool_t iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
2286  if (!iseventselected) {
2287  fHistManager.FillTH1(hname, "Rejected");
2288  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2289  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
2290  TString label;
2291  do {
2292  label = GetHFEventRejectionReasonLabel(bitmap);
2293  if (label.IsNull()) break;
2294  fHistManager.FillTH1(hname, label);
2295  } while (true);
2296 
2297  continue;
2298  }
2299  }
2300 
2301  fHistManager.FillTH1(hname, "Accepted");
2302 
2303  AliDebug(2, "Event selected");
2304 
2305  eng.RunAnalysis();
2306  }
2307  return kTRUE;
2308 }
2309 
2314 {
2315  TString hname;
2316  for (auto &param : fAnalysisEngines) {
2317 
2318  if (param.fInhibit) continue;
2319 
2320  if (fOutputType == kTreeOutput) {
2321  param.FillTree(fApplyKinematicCuts);
2323  }
2324  else if (fOutputType == kTHnOutput) {
2325  param.FillHnSparse(fApplyKinematicCuts);
2326  }
2327  }
2328  return kTRUE;
2329 }
2330 
2339 {
2340  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
2341 
2342  Double_t mass = part->Mass();
2343 
2344  // To make sure that the PDG mass value is not at the edge of a bin
2345  if (nbins % 2 == 0) {
2346  minMass = mass - range / 2 - range / nbins / 2;
2347  maxMass = mass + range / 2 - range / nbins / 2;
2348  }
2349  else {
2350  minMass = mass - range / 2;
2351  maxMass = mass + range / 2;
2352  }
2353 }
2354 
2361 {
2362  static TString label;
2363  label = "";
2364 
2365  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
2366  label = "NotSelTrigger";
2367  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
2368  return label.Data();
2369  }
2370  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
2371  label = "NoVertex";
2372  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
2373  return label.Data();
2374  }
2375  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
2376  label = "TooFewVtxContrib";
2377  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
2378  return label.Data();
2379  }
2380  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
2381  label = "ZVtxOutFid";
2382  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
2383  return label.Data();
2384  }
2385  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
2386  label = "Pileup";
2387  bitmap &= ~BIT(AliRDHFCuts::kPileup);
2388  return label.Data();
2389  }
2390  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
2391  label = "OutsideCentrality";
2392  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
2393  return label.Data();
2394  }
2395  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
2396  label = "PhysicsSelection";
2397  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
2398  return label.Data();
2399  }
2400  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
2401  label = "BadSPDVertex";
2402  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
2403  return label.Data();
2404  }
2405  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
2406  label = "ZVtxSPDOutFid";
2407  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
2408  return label.Data();
2409  }
2410  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
2411  label = "CentralityFlattening";
2412  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
2413  return label.Data();
2414  }
2415  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
2416  label = "BadTrackV0Correl";
2417  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
2418  return label.Data();
2419  }
2420 
2421  return label.Data();
2422 }
2423 
2430 {
2431  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
2432  PostData(eng.GetDataSlotNumber(), eng.GetTree());
2433  return eng.GetDataSlotNumber();
2434  }
2435  else {
2436  return -1;
2437  }
2438 }
2439 
2449 {
2450  // Get the pointer to the existing analysis manager via the static access method
2451  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
2452  if (!mgr) {
2453  ::Error("AddTaskDmesonJets", "No analysis manager to connect to.");
2454  return NULL;
2455  }
2456 
2457  // Check the analysis type using the event handlers connected to the analysis manager
2458  AliVEventHandler* handler = mgr->GetInputEventHandler();
2459  if (!handler) {
2460  ::Error("AddTaskEmcalJetSpectraQA", "This task requires an input event handler");
2461  return NULL;
2462  }
2463 
2464  EDataType_t dataType = kUnknownDataType;
2465 
2466  if (handler->InheritsFrom("AliESDInputHandler")) {
2467  dataType = kESD;
2468  }
2469  else if (handler->InheritsFrom("AliAODInputHandler")) {
2470  dataType = kAOD;
2471  }
2472 
2473  // Init the task and do settings
2474  if (ntracks == "usedefault") {
2475  if (dataType == kESD) {
2476  ntracks = "Tracks";
2477  }
2478  else if (dataType == kAOD) {
2479  ntracks = "tracks";
2480  }
2481  else {
2482  ntracks = "";
2483  }
2484  }
2485 
2486  if (nclusters == "usedefault") {
2487  if (dataType == kESD) {
2488  nclusters = "CaloClusters";
2489  }
2490  else if (dataType == kAOD) {
2491  nclusters = "caloClusters";
2492  }
2493  else {
2494  nclusters = "";
2495  }
2496  }
2497 
2498  if (nMCpart == "usedefault") {
2499  nMCpart = "mcparticles"; // Always needs AliAODMCParticle objects
2500  }
2501 
2502  TString name("AliAnalysisTaskDmesonJets");
2503  if (strcmp(suffix, "") != 0) {
2504  name += TString::Format("_%s", suffix.Data());
2505  }
2506 
2507  AliAnalysisTaskDmesonJets* jetTask = new AliAnalysisTaskDmesonJets(name, nMaxTrees);
2508  jetTask->SetVzRange(-10,10);
2509 
2510  if (!ntracks.IsNull()) {
2511  AliHFTrackContainer* trackCont = new AliHFTrackContainer(ntracks);
2512  jetTask->AdoptParticleContainer(trackCont);
2513  }
2514 
2515  if (!nMCpart.IsNull()) {
2516  AliMCParticleContainer* partCont = new AliHFAODMCParticleContainer(nMCpart);
2517  partCont->SetEtaLimits(-1.5, 1.5);
2518  jetTask->AdoptParticleContainer(partCont);
2519  }
2520 
2521  jetTask->AddClusterContainer(nclusters);
2522 
2523  // Final settings, pass to manager and set the containers
2524  mgr->AddTask(jetTask);
2525 
2526  // Create containers for input/output
2527  AliAnalysisDataContainer* cinput1 = mgr->GetCommonInputContainer();
2528  TString contname1(name);
2529  contname1 += "_histos";
2530  AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(contname1.Data(),
2531  TList::Class(), AliAnalysisManager::kOutputContainer,
2532  Form("%s", AliAnalysisManager::GetCommonFileName()));
2533 
2534  mgr->ConnectInput(jetTask, 0, cinput1);
2535  mgr->ConnectOutput(jetTask, 1, coutput1);
2536 
2537  for (Int_t i = 0; i < nMaxTrees; i++) {
2538  TString contname = TString::Format("%s_tree_%d", name.Data(), i);
2539  AliAnalysisDataContainer *coutput = mgr->CreateContainer(contname.Data(),
2540  TTree::Class(),AliAnalysisManager::kOutputContainer,
2541  Form("%s", AliAnalysisManager::GetCommonFileName()));
2542  mgr->ConnectOutput(jetTask, 2+i, coutput);
2543  }
2544  return jetTask;
2545 }
2546 
Int_t pdg
void Print() const
Prints the content of this object in the standard output.
std::list< AnalysisEngine > fAnalysisEngines
Array of analysis parameters.
AliEMCALGeometry * fGeom
!emcal geometry
Analysis task for D meson jets.
virtual void Set(const AliDmesonJetInfo &source)
Invariant mass of the D0 meson candidate in GeV/c2.
virtual Int_t Run()
double Double_t
Definition: External.C:58
Bool_t IsJetInAcceptance(const AliJetInfo &jet) const
Double_t fSoftPionPt
! Transverse momentum of the soft pion of the D* candidate
EDataType_t
Switch for the data type.
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
AliClusterContainer * AddClusterContainer(const char *n)
EJetType_t fJetType
Jet type (charged, full, neutral)
std::pair< AliJetInfo *, Double_t > jet_distance_pair
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
virtual void Set(const AliDmesonJetInfo &source)
AliTLorentzVector fD
! 4-momentum of the D meson candidate
Double_t GetDistance(const AliJetInfo &jet, Double_t &deta, Double_t &dphi) const
TObject * FindObject(const char *name) const
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Definition: THistManager.h:504
AliAnalysisTaskDmesonJets()
This is the default constructor, used for ROOT I/O purposes.
virtual void Reset()
Reset the current object.
AnalysisEngine & operator=(const AnalysisEngine &source)
Struct that encapsulates analysis parameters.
std::vector< AliHFJetDefinition >::iterator FindJetDefinition(const AliHFJetDefinition &eng)
Base task in the EMCAL framework (lighter version of AliAnalysisTaskEmcal)
Double_t Phi_0_2pi() const
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:95
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
AliClusterContainer * GetClusterContainer(Int_t i=0) const
AliAODTrack * GetBachelor() const
AliRDHFCuts * LoadDMesonCutsFromFile(TString cutfname, TString cutsname)
bool operator<(const AliAnalysisTaskDmesonJets::AliHFJetDefinition &lhs, const AliAnalysisTaskDmesonJets::AliHFJetDefinition &rhs)
std::list< AnalysisEngine >::iterator FindAnalysisEngine(const AnalysisEngine &eng)
std::vector< AliHFJetDefinition > fJetDefinitions
Jet definitions.
EJetAlgo_t fJetAlgo
Jet algorithm (kt, anti-kt,...)
static const char * GetHFEventRejectionReasonLabel(UInt_t &bitmap)
Lightweight class that encapsulates D*.
virtual void Clear(const Option_t *="")
Bool_t ExtractDstarAttributes(const AliAODRecoCascadeHF *DstarCand, AliDmesonJetInfo &DmesonJet, UInt_t i)
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
const AliJetInfo * GetJet(std::string n) const
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
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
void AddInputVectors(AliEmcalContainer *cont, Int_t offset, TH2 *rejectHist=0)
Definition: External.C:220
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Double_t minMass
void RunParticleLevelAnalysis()
Run a particle level analysis.
Bool_t fApplyKinematicCuts
Apply jet kinematic cuts.
void SetVzRange(Double_t min, Double_t max)
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)
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:97
AliTLorentzVector fMomentum
4-momentum of the jet
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
AliHFJetDefinition * AddJetDefinition(EJetType_t type, Double_t r, EJetAlgo_t algo, ERecoScheme_t reco)
Container for MC-true particles within the EMCAL framework.
Lightweight class that encapsulates D0.
const Int_t nPtBins
Bool_t IsAnyJetInAcceptance(const AliDmesonJetInfo &dMesonJet) const
void AdoptParticleContainer(AliParticleContainer *cont)
Container for jet within the EMCAL jet framework.
Definition: External.C:196
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.
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)