AliPhysics  63d3444 (63d3444)
 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 
89 {
90  fD.SetPtEtaPhiE(0,0,0,0);
91  fSoftPionPt = 0;
92  fInvMass2Prong = 0;
93  fDmesonParticle = 0;
94  fMCLabel = -1;
95  fReconstructed = kFALSE;
96  for (auto &jet : fJets) {
97  jet.second.fMomentum.SetPtEtaPhiE(0,0,0,0);
98  jet.second.fNConstituents = 0;
99  jet.second.fNEF = 0;
100  jet.second.fMaxChargedPt = 0;
101  jet.second.fMaxNeutralPt = 0;
102  }
103 }
104 
107 {
108  Printf("Printing D Meson Jet object.");
109  Printf("D Meson: pT = %.3f, eta = %.3f, phi = %.3f, inv. mass = %.3f", fD.Pt(), fD.Eta(), fD.Phi_0_2pi(), fD.M());
110  Printf("Soft pion pT: %.3f. 2-Prong Invariant mass = %.3f", fSoftPionPt, fInvMass2Prong);
111  for (auto &jet : fJets) {
112  Printf("Jet %s: pT = %.3f, eta = %.3f, phi = %.3f", jet.first.c_str(), jet.second.Pt(), jet.second.Eta(), jet.second.Phi_0_2pi());
113  Printf("Jet N Consituents = %d", jet.second.fNConstituents);
114  }
115 }
116 
121 {
122  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
123  if (it == fJets.end()) return 0;
124 
125  Double_t z = 0;
126 
127  if ((*it).second.Pt() > 0) {
128  TVector3 dvect = fD.Vect();
129  TVector3 jvect = (*it).second.fMomentum.Vect();
130 
131  Double_t jetMom = jvect * jvect;
132 
133  if (jetMom < 1e-6) {
134  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetZ", "Zero jet momentum!");
135  z = 0.999;
136  }
137  else {
138  z = (dvect * jvect) / jetMom;
139  }
140 
141  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
142  }
143 
144  return z;
145 }
146 
154 {
155  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
156  if (it == fJets.end()) return 0;
157 
158  return GetDistance((*it).second, deta, dphi);
159 }
160 
166 {
167  Double_t deta = 0;
168  Double_t dphi = 0;
169  return GetDistance(n, deta, dphi);
170 }
171 
179 {
180  dphi = TVector2::Phi_mpi_pi(fD.Phi() - jet.Phi());;
181  deta = fD.Eta() - jet.Eta();
182  return TMath::Sqrt(dphi*dphi + deta*deta);
183 }
184 
190 {
191  Double_t deta = 0;
192  Double_t dphi = 0;
193  return GetDistance(jet, deta, dphi);
194 }
195 
201 {
202  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
203  if (it == fJets.end()) {
204  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
205  n.c_str());
206  return 0;
207  }
208  return &((*it).second);
209 }
210 
216 {
217  std::map<std::string, AliJetInfo>::iterator it = fJets.find(n);
218  if (it == fJets.end()) {
219  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
220  n.c_str());
221  return 0;
222  }
223  return &((*it).second);
224 }
225 
226 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoSummary
227 
231 
237  fPt(0),
238  fEta(0),
239  fPhi(0),
240  fR(0),
241  fZ(0)
242 {
243  Set(source, n);
244 }
245 
248 {
249  fPt = 0;
250  fEta = 0;
251  fPhi = 0;
252  fR = 0;
253  fZ = 0;
254 }
255 
261 {
262  std::map<std::string, AliJetInfo>::const_iterator it = source.fJets.find(n);
263  if (it == source.fJets.end()) return;
264 
265  fPt = (*it).second.Pt();
266  fEta = (*it).second.Eta();
267  fPhi = (*it).second.Phi_0_2pi();
268  fR = source.GetDistance(n);
269  fZ = source.GetZ(n);
270 }
271 
276 {
277  fPt = source.Pt();
278  fEta = source.Eta();
279  fPhi = source.Phi_0_2pi();
280  fR = 0;
281  fZ = 0;
282 }
283 
284 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonInfoSummary
285 
289 
294  fPt(0),
295  fEta(0),
296  fPhi(0)
297 {
298  Set(source);
299 }
300 
305 {
306  fPt = source.fD.Pt();
307  fEta = source.fD.Eta();
308  fPhi = source.fD.Phi_0_2pi();
309 }
310 
313 {
314  fPt = 0;
315  fEta = 0;
316  fPhi = 0;
317 }
318 
319 // Definitions of class AliAnalysisTaskDmesonJets::AliD0InfoSummary
320 
324 
329  AliDmesonInfoSummary(source),
330  fInvMass(source.fD.M())
331 {
332 }
333 
338 {
339  fInvMass = source.fD.M();
341 }
342 
345 {
347 
348  fInvMass = 0;
349 }
350 
351 // Definitions of class AliAnalysisTaskDmesonJets::AliDStarInfoSummary
352 
356 
361  AliDmesonInfoSummary(source),
362  f2ProngInvMass(source.fInvMass2Prong),
363  fDeltaInvMass(source.fD.M() - source.fInvMass2Prong)
364 {
365 }
366 
371 {
372  f2ProngInvMass = source.fInvMass2Prong;
373  fDeltaInvMass = source.fD.M() - source.fInvMass2Prong;
375 }
376 
379 {
381 
382  f2ProngInvMass = 0;
383  fDeltaInvMass = 0;
384 }
385 
386 // Definitions of class AliAnalysisTaskDmesonJets::AliJetDefinition
387 
391 
394  TObject(),
395  fJetType(AliJetContainer::kChargedJet),
396  fRadius(0),
397  fJetAlgo(AliJetContainer::antikt_algorithm),
398  fRecoScheme(AliJetContainer::pt_scheme),
399  fMinJetPt(0.),
400  fMaxJetPt(500.),
401  fMinJetPhi(0.),
402  fMaxJetPhi(0.),
403  fMinJetEta(-1.),
404  fMaxJetEta(1.),
405  fMinChargedPt(0.),
406  fMaxChargedPt(0.),
407  fMinNeutralPt(0.),
408  fMaxNeutralPt(0.)
409 {
410 }
411 
419  TObject(),
420  fJetType(type),
421  fRadius(r),
422  fJetAlgo(algo),
423  fRecoScheme(reco),
424  fMinJetPt(0.),
425  fMaxJetPt(500.),
426  fMinJetPhi(0.),
427  fMaxJetPhi(0.),
428  fMinJetEta(-1.),
429  fMaxJetEta(1.),
430  fMinChargedPt(0.),
431  fMaxChargedPt(0.),
432  fMinNeutralPt(0.),
433  fMaxNeutralPt(0.)
434 {
435 }
436 
441  TObject(),
442  fJetType(source.fJetType),
443  fRadius(source.fRadius),
444  fJetAlgo(source.fJetAlgo),
445  fRecoScheme(source.fRecoScheme),
446  fMinJetPt(source.fMinJetPt),
447  fMaxJetPt(source.fMaxJetPt),
448  fMinJetPhi(source.fMinJetPhi),
449  fMaxJetPhi(source.fMaxJetPhi),
450  fMinJetEta(source.fMinJetEta),
451  fMaxJetEta(source.fMaxJetEta),
452  fMinChargedPt(source.fMinChargedPt),
453  fMaxChargedPt(source.fMaxChargedPt),
454  fMinNeutralPt(source.fMinNeutralPt),
455  fMaxNeutralPt(source.fMaxNeutralPt)
456 {
457 }
458 
463 {
464  new (this) AliHFJetDefinition(source);
465  return *this;
466 }
467 
470 {
471  static TString name;
472 
473  name = AliJetContainer::GenerateJetName(fJetType, fJetAlgo, fRecoScheme, fRadius, 0, 0, "Jet");
474 
475  return name.Data();
476 }
477 
483 {
484  if (fMinJetEta < fMaxJetEta && (jet.Eta() < fMinJetEta || jet.Eta() > fMaxJetEta)) return kFALSE;
485  if (fMinJetPhi < fMaxJetPhi && (jet.Phi() < fMinJetPhi || jet.Phi() > fMaxJetPhi)) return kFALSE;
486  if (jet.Pt() > fMaxJetPt || jet.Pt() < fMinJetPt) return kFALSE;
487  if (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt) return kFALSE;
488  if (jet.fMaxNeutralPt < fMinNeutralPt || jet.fMaxNeutralPt > fMaxNeutralPt) return kFALSE;
489 
490  return kTRUE;
491 }
492 
498 {
499  const AliJetInfo* jet = dMesonJet.GetJet(n);
500  if (!jet) return kFALSE;
501  return IsJetInAcceptance((*jet));
502 }
503 
510 {
511  if (lhs.fJetType > rhs.fJetType) return false;
512  else if (lhs.fJetType < rhs.fJetType) return true;
513  else {
514  if (lhs.fRadius > rhs.fRadius) return false;
515  else if (lhs.fRadius < rhs.fRadius) return true;
516  else {
517  if (lhs.fJetAlgo > rhs.fJetAlgo) return false;
518  else if (lhs.fJetAlgo < rhs.fJetAlgo) return true;
519  else {
520  if (lhs.fRecoScheme < rhs.fRecoScheme) return true;
521  else return false;
522  }
523  }
524  }
525 }
526 
533 {
534  if (lhs.fJetType != rhs.fJetType) return false;
535  if (lhs.fRadius != rhs.fRadius) return false;
536  if (lhs.fJetAlgo != rhs.fJetAlgo) return false;
537  if (lhs.fRecoScheme != rhs.fRecoScheme) return false;
538  return true;
539 }
540 
541 // Definitions of class AliAnalysisTaskDmesonJets::AnalysisEngine
542 
546 
549  TObject(),
550  fCandidateType(kD0toKpi),
551  fCandidateName(),
552  fCandidatePDG(0),
553  fNDaughters(0),
554  fPDGdaughters(),
555  fBranchName(),
556  fMCMode(kNoMC),
557  fNMassBins(0),
558  fMinMass(0),
559  fMaxMass(0),
560  fRDHFCuts(0),
561  fRejectedOrigin(0),
562  fAcceptedDecay(0),
563  fInhibit(kFALSE),
564  fJetDefinitions(),
565  fPtBinWidth(0.5),
566  fMaxPt(100),
567  fDataSlotNumber(-1),
568  fTree(0),
569  fCurrentDmesonJetInfo(0),
570  fCurrentJetInfo(0),
571  fCandidateArray(0),
572  fMCContainer(0),
573  fTrackContainer(0),
574  fClusterContainer(0),
575  fAodEvent(0),
576  fFastJetWrapper(0),
577  fHistManager(0)
578 {
579 }
580 
589  TObject(),
590  fCandidateType(type),
591  fCandidateName(),
592  fCandidatePDG(0),
593  fNDaughters(0),
594  fPDGdaughters(),
595  fBranchName(),
596  fMCMode(MCmode),
597  fNMassBins(nMassBins),
598  fMinMass(0),
599  fMaxMass(0),
600  fRDHFCuts(cuts),
602  fAcceptedDecay(kAnyDecay),
603  fInhibit(kFALSE),
604  fJetDefinitions(),
605  fPtBinWidth(0.5),
606  fMaxPt(100),
607  fDataSlotNumber(-1),
608  fTree(0),
609  fCurrentDmesonJetInfo(0),
610  fCurrentJetInfo(0),
611  fCandidateArray(0),
612  fMCContainer(0),
613  fTrackContainer(0),
614  fClusterContainer(0),
615  fAodEvent(0),
616  fFastJetWrapper(0),
617  fHistManager(0)
618 {
619  SetCandidateProperties(range);
620 }
621 
626  TObject(source),
627  fCandidateType(source.fCandidateType),
628  fCandidateName(source.fCandidateName),
629  fCandidatePDG(source.fCandidatePDG),
630  fNDaughters(source.fNDaughters),
631  fPDGdaughters(source.fPDGdaughters),
632  fBranchName(source.fBranchName),
633  fMCMode(source.fMCMode),
634  fNMassBins(source.fNMassBins),
635  fMinMass(source.fMinMass),
636  fMaxMass(source.fMaxMass),
637  fRDHFCuts(),
638  fRejectedOrigin(source.fRejectedOrigin),
639  fAcceptedDecay(source.fAcceptedDecay),
640  fInhibit(source.fInhibit),
641  fJetDefinitions(source.fJetDefinitions),
642  fPtBinWidth(source.fPtBinWidth),
643  fMaxPt(source.fMaxPt),
644  fDataSlotNumber(-1),
645  fTree(0),
646  fCurrentDmesonJetInfo(0),
647  fCurrentJetInfo(0),
648  fCandidateArray(source.fCandidateArray),
649  fMCContainer(source.fMCContainer),
650  fTrackContainer(source.fTrackContainer),
651  fClusterContainer(source.fClusterContainer),
652  fAodEvent(source.fAodEvent),
654  fHistManager(source.fHistManager)
655 {
656  SetRDHFCuts(source.fRDHFCuts);
657 }
658 
659 // Destructor
661 {
662  delete fRDHFCuts;
663 }
664 
669 {
670  new (this) AnalysisEngine(source);
671  return *this;
672 }
673 
678 {
679  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
680  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
681  }
682 
683  return kFALSE;
684 }
685 
687 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const /*geom*/, Int_t /*runNumber*/)
688 {
689 }
690 
695 {
696  switch (fCandidateType) {
697  case kD0toKpi:
698  fCandidatePDG = 421;
699  fCandidateName = "D0";
700  fNDaughters = 2;
701  fPDGdaughters.Set(fNDaughters);
702  fPDGdaughters.Reset();
703  fPDGdaughters[0] = 211; // pi
704  fPDGdaughters[1] = 321; // K
705  fBranchName = "D0toKpi";
706  fAcceptedDecay = kDecayD0toKpi;
707  if (!fRDHFCuts) {
708  fRDHFCuts = new AliRDHFCutsD0toKpi();
709  fRDHFCuts->SetStandardCutsPP2010();
710  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
711  fRDHFCuts->SetTriggerClass("","");
712  }
713  break;
714  case kD0toKpiLikeSign:
715  fCandidatePDG = 421;
716  fCandidateName = "2ProngLikeSign";
717  fNDaughters = 2;
718  fPDGdaughters.Set(fNDaughters);
719  fPDGdaughters.Reset();
720  fPDGdaughters[0] = 211; // pi
721  fPDGdaughters[1] = 321; // K
722  fBranchName = "LikeSign2Prong";
723  fAcceptedDecay = kAnyDecay;
724  if (!fRDHFCuts) {
725  fRDHFCuts = new AliRDHFCutsD0toKpi();
726  fRDHFCuts->SetStandardCutsPP2010();
727  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
728  fRDHFCuts->SetTriggerClass("","");
729  }
730  break;
731  case kDstartoKpipi:
732  fCandidatePDG = 413;
733  fCandidateName = "DStar";
734  fNDaughters = 3;
735  fPDGdaughters.Set(fNDaughters);
736  fPDGdaughters.Reset();
737  fPDGdaughters[0] = 211; // pi soft
738  fPDGdaughters[1] = 211; // pi fromD0
739  fPDGdaughters[2] = 321; // K from D0
740  fBranchName = "Dstar";
741  fAcceptedDecay = kDecayDStartoKpipi;
742  if (!fRDHFCuts) {
743  fRDHFCuts = new AliRDHFCutsDStartoKpipi();
744  fRDHFCuts->SetStandardCutsPP2010();
745  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
746  fRDHFCuts->SetTriggerClass("","");
747  }
748  break;
749  default:
750  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
751  }
752 
753  CalculateMassLimits(range, fCandidatePDG, fNMassBins, fMinMass, fMaxMass);
754 }
755 
760 {
761  if (fRDHFCuts) delete fRDHFCuts;
762  fRDHFCuts = cuts;
763 }
764 
769 {
770  if (!cuts) return;
771  if (fRDHFCuts) delete fRDHFCuts;
772  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
773 }
774 
779 {
780  static TString name;
781 
782  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
783 
784  return name.Data();
785 }
786 
791 {
792  static TString name;
793 
794  name = fCandidateName;
795  switch (fMCMode) {
796  case kBackgroundOnly:
797  name += "_kBackgroundOnly";
798  break;
799  case kSignalOnly:
800  name += "_kSignalOnly";
801  break;
802  case kMCTruth:
803  name += "_MCTruth";
804  break;
805  default:
806  break;
807  }
808 
809  return name.Data();
810 }
811 
819 {
820  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
821 
822  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
823  fJetDefinitions.push_back(def);
824  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
825  "Total number of jet definitions is now %lu.",
826  def.GetName(), GetName(), fJetDefinitions.size());
827  // For detector level set maximum track pt to 100 GeV/c
828  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
829  }
830  else {
831  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
832  }
833 
834  return &(*it);
835 }
836 
848 {
849  AliHFJetDefinition def(type, r, algo, reco);
850 
851  return AddJetDefinition(def);
852 }
853 
859 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
860 {
861  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
862  while (it != fJetDefinitions.end() && (*it) != def) it++;
863  return it;
864 }
865 
872 {
873  if (lhs.fCandidateType > rhs.fCandidateType) return false;
874  else if (lhs.fCandidateType < rhs.fCandidateType) return true;
875  else {
876  if (lhs.fMCMode < rhs.fMCMode) return true;
877  else return false;
878  }
879 }
880 
887 {
888  if (lhs.fCandidateType != rhs.fCandidateType) return false;
889  if (lhs.fMCMode != rhs.fMCMode) return false;
890  return true;
891 }
892 
901 {
902  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) { // D0 candidate
903  return ExtractD0Attributes(Dcand, DmesonJet, i);
904  }
905  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
906  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
907  }
908  else {
909  return kFALSE;
910  }
911 }
912 
921 {
922  //AliDebug(2,"Checking if D0 meson is selected");
923  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
924  if (isSelected == 0) return kFALSE;
925 
926  Int_t MCtruthPdgCode = 0;
927 
928  Double_t invMassD = 0;
929 
930  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
931  // Checks also the origin, and if it matches the rejected origin mask, return false
932  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
933  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
934  DmesonJet.fMCLabel = mcLab;
935 
936  // Retrieve the generated particle (if exists) and its PDG code
937  if (mcLab >= 0) {
938  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
939 
940  if (aodMcPart) {
941  // Check origin and return false if it matches the rejected origin mask
942  if (fRejectedOrigin) {
943  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
944 
945  UInt_t rs = origin;
946  UShort_t p = 0;
947  while (rs >>= 1) { p++; }
948  TString hname = TString::Format("%s/fHistDmesonOrigin", GetName());
949  fHistManager->FillTH1(hname, p);
950 
951  if ((origin & fRejectedOrigin) == origin) return kFALSE;
952  }
953  MCtruthPdgCode = aodMcPart->PdgCode();
954  }
955  }
956  }
957 
958  if (isSelected == 1) { // selected as a D0
959  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
960 
961  if (fMCMode == kNoMC ||
962  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
963  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly)) {
964  // 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)
965  //AliDebug(2,"Selected as D0");
966  invMassD = Dcand->InvMassD0();
967  }
968  else { // conditions above not passed, so return FALSE
969  return kFALSE;
970  }
971  }
972  else if (isSelected == 2) { // selected as a D0bar
973  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
974 
975  if (fMCMode == kNoMC ||
976  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
977  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly)) {
978  // 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)
979  //AliDebug(2,"Selected as D0bar");
980  invMassD = Dcand->InvMassD0bar();
981  }
982  else { // conditions above not passed, so return FALSE
983  return kFALSE;
984  }
985  }
986  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
987  //AliDebug(2,"Selected as either D0 or D0bar");
988 
989  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
990  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
991  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kBackgroundOnly)) {
992  if (i > 0) return kFALSE;
993  //AliDebug(2, "MC truth is D0");
994  invMassD = Dcand->InvMassD0();
995  }
996  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
997  (MCtruthPdgCode == fCandidatePDG && fMCMode == kBackgroundOnly)) {
998  if (i > 0) return kFALSE;
999  //AliDebug(2, "MC truth is D0bar");
1000  invMassD = Dcand->InvMassD0bar();
1001  }
1002  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
1003 
1004  // Only accept it if background-only OR background-and-signal was requested
1005  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
1006  // Select D0 or D0bar depending on the i-parameter
1007  if (i == 0) {
1008  //AliDebug(2, "Returning invariant mass with D0 hypothesis");
1009  invMassD = Dcand->InvMassD0();
1010  }
1011  else if (i == 1) {
1012  //AliDebug(2, "Returning invariant mass with D0bar hypothesis");
1013  invMassD = Dcand->InvMassD0bar();
1014  }
1015  else { // i > 1
1016  return kFALSE;
1017  }
1018  }
1019  else { // signal-only was requested but this is not a true D0
1020  return kFALSE;
1021  }
1022  }
1023  }
1024  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
1025  return kTRUE;
1026 }
1027 
1036 {
1037  //AliDebug(2,"Checking if D0 meson is selected");
1038  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoCascadeHF*>(DstarCand), AliRDHFCuts::kAll, fAodEvent);
1039  if (isSelected == 0) return kFALSE;
1040 
1041  if (i > 0) return kFALSE; // only one mass hypothesis for the D*
1042 
1043  Int_t MCtruthPdgCode = 0;
1044 
1045  Double_t invMassD = 0;
1046 
1047  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1048  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1049  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1050 
1051  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1052  //AliDebug(2, Form("MC label is %d", mcLab));
1053  DmesonJet.fMCLabel = mcLab;
1054  if (mcLab >= 0) {
1055  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1056 
1057  if (aodMcPart) {
1058  if (fRejectedOrigin) {
1059  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
1060 
1061  UInt_t rs = origin;
1062  UShort_t p = 0;
1063  while (rs >>= 1) { p++; }
1064  TString hname = TString::Format("%s/fHistDmesonOrigin", GetName());
1065  fHistManager->FillTH1(hname, p);
1066 
1067  if ((origin & fRejectedOrigin) == origin) return kFALSE;
1068  }
1069 
1070  MCtruthPdgCode = aodMcPart->PdgCode();
1071  //AliDebug(2, Form("MC truth pdg code is %d",MCtruthPdgCode));
1072  }
1073  }
1074  }
1075 
1076  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1077  if (fMCMode == kNoMC ||
1078  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1079  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1080  // 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)
1081  invMassD = DstarCand->InvMassDstarKpipi();
1082  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1083  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1084  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1085  return kTRUE;
1086  }
1087  else { // conditions above not passed, so return FALSE
1088  return kFALSE;
1089  }
1090 }
1091 
1099 {
1100  if (!part) return kUnknownDecay;
1101  if (!mcArray) return kUnknownDecay;
1102 
1104 
1105  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1106 
1107  if (part->GetNDaughters() == 2) {
1108 
1109  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1110  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1111 
1112  if (!d1 || !d2) {
1113  return decay;
1114  }
1115 
1116  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1117  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1118 
1119  if (absPdgPart == 421) { // D0 -> K pi
1120  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1121  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1122  decay = kDecayD0toKpi;
1123  }
1124  }
1125 
1126  if (absPdgPart == 413) { // D* -> D0 pi
1127  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1128  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1129  if (D0decay == kDecayD0toKpi) {
1130  decay = kDecayDStartoKpipi;
1131  }
1132  }
1133  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1134  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1135  if (D0decay == kDecayD0toKpi) {
1136  decay = kDecayDStartoKpipi;
1137  }
1138  }
1139  }
1140  }
1141 
1142  return decay;
1143 }
1144 
1152 {
1153  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1154 
1155  if (!part) return kUnknownQuark;
1156  if (!mcArray) return kUnknownQuark;
1157 
1158  Int_t mother = part->GetMother();
1159  while (mother >= 0) {
1160  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1161  if (mcGranma) {
1162  Int_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
1163 
1164  if (abspdgGranma == 1) return kFromDown;
1165  if (abspdgGranma == 2) return kFromUp;
1166  if (abspdgGranma == 3) return kFromStrange;
1167  if (abspdgGranma == 4) return kFromCharm;
1168  if (abspdgGranma == 5) return kFromBottom;
1169  if (abspdgGranma == 6) return kFromTop;
1170 
1171  mother = mcGranma->GetMother();
1172  }
1173  else {
1174  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1175  break;
1176  }
1177  }
1178 
1179  return kUnknownQuark;
1180 }
1181 
1184 {
1185  fDmesonJets.clear();
1186 
1187  for (auto& jetDef : fJetDefinitions) {
1188  jetDef.fJets.clear();
1189  }
1190 
1191  if (fMCMode == kMCTruth) {
1192  RunParticleLevelAnalysis();
1193  }
1194  else {
1195  RunDetectorLevelAnalysis();
1196  }
1197 }
1198 
1204 {
1206  fFastJetWrapper->SetR(jetDef.fRadius);
1209 
1210  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1211  fTrackContainer->SetDMesonCandidate(0);
1212  AddInputVectors(fTrackContainer, 100);
1213  }
1214 
1215  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1216  AddInputVectors(fClusterContainer, -100);
1217  }
1218 
1219  // run jet finder
1220  fFastJetWrapper->Run();
1221 
1222  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1223 
1224  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1225  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1226 
1227  Double_t maxChPt = 0;
1228  Double_t maxNePt = 0;
1229  Double_t totalNeutralPt = 0;
1230 
1231  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1232  if (constituents[ic].user_index() >= 100) {
1233  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1234  }
1235  else if (constituents[ic].user_index() <= -100) {
1236  totalNeutralPt += constituents[ic].pt();
1237  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1238  }
1239  }
1240 
1241  jetDef.fJets.push_back(
1242  AliJetInfo(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E(),
1243  constituents.size(), totalNeutralPt / jets_incl[ijet].pt(), maxChPt, maxNePt)
1244  );
1245  }
1246 }
1247 
1256 {
1257  if (jetDef.fJets.size() == 0) FindJets(jetDef);
1258 
1259  Double_t d_closest = 999;
1260  AliJetInfo* jet_closest = 0;
1261  const AliJetInfo* truth_jet = 0;
1262  try {
1263  truth_jet = &(dmeson.fJets.at(jetDef.GetName()));
1264  }
1265  catch(...) {
1266  return jet_distance_pair(0, 999);
1267  }
1268 
1269  for (auto& jet : jetDef.fJets) {
1270  Double_t d = truth_jet->GetDistance(jet);
1271 
1272  if (d > dMax) continue;
1273  if (d < d_closest) {
1274  d_closest = d;
1275  jet_closest = &jet;
1276  }
1277  }
1278 
1279  if (jet_closest && applyKinCuts) {
1280  if (!jetDef.IsJetInAcceptance(*jet_closest)) jet_closest = 0;
1281  }
1282 
1283  if (jet_closest) {
1284  AliDebug(2, Form("Found closest jet (pt=%.3f, eta=%.3f, phi=%.3f) to d meson (pt=%.3f, eta=%.3f, phi=%.3f) with d = %.3f",
1285  jet_closest->Pt(), jet_closest->Eta(), jet_closest->Phi_0_2pi(),
1286  dmeson.fD.Pt(), dmeson.fD.Eta(), dmeson.fD.Phi_0_2pi(),
1287  d_closest));
1288  }
1289 
1290  return jet_distance_pair(jet_closest, d_closest);
1291 }
1292 
1295 {
1296  const Int_t nD = fCandidateArray->GetEntriesFast();
1297 
1298  AliDmesonJetInfo DmesonJet;
1299 
1300  Int_t nAccCharm = 0;
1301  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1302  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1303  if (!charmCand) continue;
1304 
1305  // region of interest + cuts
1306  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1307 
1308  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1309  DmesonJet.Reset();
1310  DmesonJet.fDmesonParticle = charmCand;
1311  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1312  for (std::vector<AliHFJetDefinition>::iterator itdef = fJetDefinitions.begin(); itdef != fJetDefinitions.end(); itdef++) {
1313  if (!FindJet(charmCand, DmesonJet, *itdef)) {
1314  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1315  (*itdef).GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1316  }
1317  }
1318  fDmesonJets[icharm] = DmesonJet;
1319  }
1320  }
1321  nAccCharm++;
1322  } // end of D cand loop
1323 
1324  TString hname;
1325 
1326  hname = TString::Format("%s/fHistNAcceptedDmesons", GetName());
1327  fHistManager->FillTH1(hname, nAccCharm);
1328 
1329  hname = TString::Format("%s/fHistNDmesons", GetName());
1330  fHistManager->FillTH1(hname, nD);
1331 }
1332 
1343 {
1344  TString hname;
1345 
1347  fFastJetWrapper->SetR(jetDef.fRadius);
1350 
1351  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1352 
1353  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1354  fTrackContainer->SetDMesonCandidate(Dcand);
1355  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1356  AddInputVectors(fTrackContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1357 
1358  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1359  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1360  const TObjArray& daughters = fTrackContainer->GetDaughterList();
1361  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1362  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1363  if (!fTrackContainer->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1364  }
1365  }
1366 
1367  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1368  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1369  AddInputVectors(fClusterContainer, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1370  }
1371 
1372  // run jet finder
1373  fFastJetWrapper->Run();
1374 
1375  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1376 
1377  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1378  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1379 
1380  Bool_t isDmesonJet = kFALSE;
1381 
1382  Double_t maxChPt = 0;
1383  Double_t maxNePt = 0;
1384  Double_t totalNeutralPt = 0;
1385 
1386  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1387  if (constituents[ic].user_index() == 0) {
1388  isDmesonJet = kTRUE;
1389  }
1390  else if (constituents[ic].user_index() >= 100) {
1391  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1392  }
1393  else if (constituents[ic].user_index() <= -100) {
1394  totalNeutralPt += constituents[ic].pt();
1395  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1396  }
1397  }
1398 
1399  if (isDmesonJet) {
1400  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1401  DmesonJet.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1402  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1403  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1404  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1405 
1406  return kTRUE;
1407  }
1408  }
1409 
1410  return kFALSE;
1411 }
1412 
1416 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist)
1417 {
1418  AliEmcalIterableMomentumContainer itcont = cont->all_momentum();
1419  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1420  UInt_t rejectionReason = 0;
1421  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1422  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1423  continue;
1424  }
1425  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1426  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1427  }
1428 }
1429 
1432 {
1433  TString hname;
1434 
1435  fMCContainer->SetSpecialPDG(fCandidatePDG);
1436  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1437  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1438 
1439  if (!fMCContainer->IsSpecialPDGFound()) return;
1440 
1441  for (auto &jetDef : fJetDefinitions) {
1442  switch (jetDef.fJetType) {
1444  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
1445  break;
1447  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
1448  break;
1450  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
1451  break;
1452  }
1453 
1455  fFastJetWrapper->SetR(jetDef.fRadius);
1458 
1459  hname = TString::Format("%s/fHistDmesonOrigin", GetName());
1460  fMCContainer->SetHistOrigin(static_cast<TH1*>(fHistManager->FindObject(hname)));
1461  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1462  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1463  fMCContainer->SetHistOrigin(0);
1464 
1465  fFastJetWrapper->Run();
1466 
1467  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1468 
1469  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1470  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1471 
1472  Bool_t isDmesonJet = kFALSE;
1473 
1474  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1475  Int_t iPart = constituents[ic].user_index() - 100;
1476  AliVParticle* part = fMCContainer->GetParticle(iPart);
1477  if (!part) {
1478  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1479  continue;
1480  }
1481  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1482  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1483  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1484  std::pair<int, AliDmesonJetInfo> element;
1485  element.first = iPart;
1486  element.second = AliDmesonJetInfo();
1487  dMesonJetIt = fDmesonJets.insert(element).first;
1488  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1489  (*dMesonJetIt).second.fDmesonParticle = part;
1490  }
1491 
1492  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1493  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1494  }
1495  }
1496  }
1497  }
1498 }
1499 
1504 {
1505  TString classname;
1506  switch (fCandidateType) {
1507  case kD0toKpi:
1508  case kD0toKpiLikeSign:
1509  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1510  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1511  break;
1512  case kDstartoKpipi:
1513  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1514  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1515  break;
1516  }
1517  TString treeName = TString::Format("%s_%s", taskName, GetName());
1518  fTree = new TTree(treeName, treeName);
1519  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
1520  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1521  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1522  fCurrentJetInfo[i] = new AliJetInfoSummary();
1523  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1524  }
1525 
1526  return fTree;
1527 }
1528 
1533 {
1534  TString hname;
1535 
1536  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
1537 
1538  for (auto &jetDef : fJetDefinitions) {
1539 
1540  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
1541 
1542  Double_t radius = jetDef.fRadius;
1543 
1544  TString title[30] = {""};
1545  Int_t nbins[30] = {0 };
1546  Double_t min [30] = {0.};
1547  Double_t max [30] = {0.};
1548  Int_t dim = 0 ;
1549 
1550  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1551  nbins[dim] = nPtBins;
1552  min[dim] = 0;
1553  max[dim] = fMaxPt;
1554  dim++;
1555 
1556  if ((enabledAxis & kPositionD) != 0) {
1557  title[dim] = "#eta_{D}";
1558  nbins[dim] = 50;
1559  min[dim] = -1;
1560  max[dim] = 1;
1561  dim++;
1562 
1563  title[dim] = "#phi_{D} (rad)";
1564  nbins[dim] = 150;
1565  min[dim] = 0;
1566  max[dim] = TMath::TwoPi();
1567  dim++;
1568  }
1569 
1570  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1571  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1572  nbins[dim] = fNMassBins;
1573  min[dim] = fMinMass;
1574  max[dim] = fMaxMass;
1575  dim++;
1576  }
1577 
1578  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1579  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1580  nbins[dim] = fNMassBins;
1581  min[dim] = fMinMass;
1582  max[dim] = fMaxMass;
1583  dim++;
1584  }
1585 
1586  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1587  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1588  nbins[dim] = fNMassBins;
1589  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1590  dim++;
1591  }
1592 
1593  if (fCandidateType == kDstartoKpipi) {
1594  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1595  nbins[dim] = fNMassBins*6;
1596  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1597 
1598  // subtract mass of D0
1599  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1600  min[dim] -= D0mass;
1601  max[dim] -= D0mass;
1602 
1603  dim++;
1604  }
1605 
1606  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1607  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1608  nbins[dim] = 100;
1609  min[dim] = 0;
1610  max[dim] = 25;
1611  dim++;
1612  }
1613 
1614  title[dim] = "#it{z}_{D}";
1615  nbins[dim] = 110;
1616  min[dim] = 0;
1617  max[dim] = 1.10;
1618  dim++;
1619 
1620  if ((enabledAxis & kDeltaR) != 0) {
1621  title[dim] = "#Delta R_{D-jet}";
1622  nbins[dim] = 100;
1623  min[dim] = 0;
1624  max[dim] = radius * 1.5;
1625  dim++;
1626  }
1627 
1628  if ((enabledAxis & kDeltaEta) != 0) {
1629  title[dim] = "#eta_{D} - #eta_{jet}";
1630  nbins[dim] = 100;
1631  min[dim] = -radius * 1.2;
1632  max[dim] = radius * 1.2;
1633  dim++;
1634  }
1635 
1636  if ((enabledAxis & kDeltaPhi) != 0) {
1637  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1638  nbins[dim] = 100;
1639  min[dim] = -radius * 1.2;
1640  max[dim] = radius * 1.2;
1641  dim++;
1642  }
1643 
1644  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1645  nbins[dim] = nPtBins;
1646  min[dim] = 0;
1647  max[dim] = fMaxPt;
1648  dim++;
1649 
1650  if ((enabledAxis & kPositionJet) != 0) {
1651  title[dim] = "#eta_{jet}";
1652  nbins[dim] = 50;
1653  min[dim] = -1;
1654  max[dim] = 1;
1655  dim++;
1656 
1657  title[dim] = "#phi_{jet} (rad)";
1658  nbins[dim] = 150;
1659  min[dim] = 0;
1660  max[dim] = TMath::TwoPi();
1661  dim++;
1662  }
1663 
1664  if ((enabledAxis & kJetConstituents) != 0) {
1665  title[dim] = "No. of constituents";
1666  nbins[dim] = 50;
1667  min[dim] = -0.5;
1668  max[dim] = 49.5;
1669  dim++;
1670  }
1671 
1672  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1673  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1674  for (Int_t j = 0; j < dim; j++) {
1675  h->GetAxis(j)->SetTitle(title[j]);
1676  }
1677  }
1678 }
1679 
1684 {
1685  TString hname;
1686 
1687  for (auto& dmeson_pair : fDmesonJets) {
1688  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
1689  Int_t accJets = 0;
1690  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1691  fCurrentJetInfo[ij]->Reset();
1692  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1693  if (!jet) continue;
1694  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1695  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1696  fHistManager->FillTH1(hname, jet->Pt());
1697  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1698  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1699  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1700  fHistManager->FillTH1(hname, jet->Eta());
1701  continue;
1702  }
1703  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
1704  accJets++;
1705  }
1706  if (accJets > 0) {
1707  fTree->Fill();
1708  }
1709  else {
1710  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1711  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1712  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1713  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1714  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1715  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1716  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1717  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1718  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1719  }
1720  else if (fCandidateType == kDstartoKpipi) {
1721  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1722  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1723 
1724  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1725  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1726  }
1727  }
1728  }
1729  return kTRUE;
1730 }
1731 
1737 {
1738  TString hname;
1739 
1740  for (auto& dmeson_pair : fDmesonJets) {
1741  Int_t accJets = 0;
1742  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1743  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1744  if (!jet) continue;
1745  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1746  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1747  fHistManager->FillTH1(hname, jet->Pt());
1748  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1749  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1750  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1751  fHistManager->FillTH1(hname, jet->Eta());
1752  continue;
1753  }
1754  accJets++;
1755  }
1756  if (accJets == 0) {
1757  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1758  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1759  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1760  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1761  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1762  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1763  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1764  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1765  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1766  }
1767  else if (fCandidateType == kDstartoKpipi) {
1768  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1769  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1770 
1771  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1772  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1773  }
1774  }
1775  }
1776  return kTRUE;
1777 }
1778 
1783 {
1784  TString hname;
1785 
1786  for (auto& dmeson_pair : fDmesonJets) {
1787  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
1788  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1789  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1790  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1791  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1792  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1793  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1794  }
1795  }
1796 
1797  for (auto &jetDef : fJetDefinitions) {
1798 
1799  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1800  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
1801 
1802  for (auto& dmeson_pair : fDmesonJets) {
1803  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
1804  if (!jet) continue;
1805  if (!jetDef.IsJetInAcceptance(*jet)) {
1806  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
1807  fHistManager->FillTH1(hname, jet->Pt());
1808  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
1809  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1810  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
1811  fHistManager->FillTH1(hname, jet->Eta());
1812  continue;
1813  }
1814  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
1815  }
1816  }
1817 
1818  return kTRUE;
1819 }
1820 
1827 {
1828  // Fill the THnSparse histogram.
1829 
1830  Double_t contents[30] = {0.};
1831 
1832  Double_t z = DmesonJet.GetZ(n);
1833  Double_t deltaPhi = 0;
1834  Double_t deltaEta = 0;
1835  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
1836 
1837  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
1838  if (it == DmesonJet.fJets.end()) return kFALSE;
1839 
1840  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
1841  TString title(h->GetAxis(i)->GetTitle());
1842  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
1843  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
1844  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
1845  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
1846  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
1847  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
1848  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
1849  else if (title=="#it{z}_{D}") contents[i] = z;
1850  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
1851  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
1852  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
1853  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
1854  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
1855  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
1856  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
1857  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
1858  }
1859 
1860  h->Fill(contents);
1861 
1862  return kTRUE;
1863 }
1864 
1865 // Definitions of class AliAnalysisTaskDmesonJets
1866 
1870 
1874  fAnalysisEngines(),
1875  fEnabledAxis(0),
1877  fHistManager(),
1878  fApplyKinematicCuts(kTRUE),
1879  fNOutputTrees(0),
1880  fAodEvent(0),
1881  fFastJetWrapper(0)
1882 {
1883  SetMakeGeneralHistograms(kTRUE);
1884 }
1885 
1890  AliAnalysisTaskEmcalLight(name, kTRUE),
1891  fAnalysisEngines(),
1892  fEnabledAxis(k2ProngInvMass),
1893  fOutputType(kTreeOutput),
1894  fHistManager(name),
1895  fApplyKinematicCuts(kTRUE),
1896  fNOutputTrees(nOutputTrees),
1897  fAodEvent(0),
1898  fFastJetWrapper(0)
1899 {
1900  SetMakeGeneralHistograms(kTRUE);
1901  for (Int_t i = 0; i < nOutputTrees; i++){
1902  DefineOutput(2+i, TTree::Class());
1903  }
1904 }
1905 
1908 {
1909  if (fFastJetWrapper) delete fFastJetWrapper;
1910 }
1911 
1919 {
1920  AliRDHFCuts* analysiscuts = 0;
1921  TFile* filecuts = TFile::Open(cutfname);
1922  if (!filecuts || filecuts->IsZombie()) {
1923  ::Warning("AddTaskDmesonJets", "Input file not found: will use std cuts.");
1924  filecuts = 0;
1925  }
1926 
1927  if (filecuts) {
1928  analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
1929  if (!analysiscuts) {
1930  ::Warning("AddTaskDmesonJetCorr", "Could not find analysis cuts '%s' in '%s'. Using std cuts.", cutsname.Data(), cutfname.Data());
1931  }
1932  }
1933 
1934  return analysiscuts;
1935 }
1936 
1946 {
1948  return AddAnalysisEngine(type, MCmode, jetDef, cutfname);
1949 }
1950 
1960 {
1961  AliRDHFCuts* cuts = 0;
1962 
1963  if (!cutfname.IsNull()) {
1964  TString cutsname;
1965 
1966  switch (type) {
1967  case kD0toKpi:
1968  case kD0toKpiLikeSign:
1969  cutsname = "D0toKpiCuts";
1970  break;
1971  case kDstartoKpipi:
1972  cutsname = "DStartoKpipiCuts";
1973  break;
1974  default:
1975  return 0;
1976  }
1977 
1978  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
1979  }
1980 
1981  AnalysisEngine eng(type, MCmode, cuts);
1982 
1983  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
1984 
1985  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
1986  eng.AddJetDefinition(jetDef);
1987  it = fAnalysisEngines.insert(it, eng);
1988  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
1989  }
1990  else {
1991  AnalysisEngine* found_eng = &(*it);
1992  ::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());
1993  found_eng->AddJetDefinition(jetDef);
1994 
1995  if (cuts && found_eng->fRDHFCuts != 0) {
1996  ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
1997  found_eng->SetRDHFCuts(cuts);
1998  }
1999  }
2000 
2001  return &(*it);
2002 }
2003 
2004 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
2005 {
2006  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
2007  while (it != fAnalysisEngines.end() && (*it) < eng) it++;
2008  return it;
2009 }
2010 
2013 {
2014  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
2015 
2017 
2018  // Define histograms
2019  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
2020 
2021  TString hname;
2022  TString htitle;
2023  TH1* h = 0;
2024  Int_t treeSlot = 0;
2025 
2026  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
2027  for (auto &param : fAnalysisEngines) {
2028  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
2029 
2030  fHistManager.CreateHistoGroup(param.GetName());
2031 
2032  param.fHistManager = &fHistManager;
2033 
2034  hname = TString::Format("%s/fHistNAcceptedDmesons", param.GetName());
2035  htitle = hname + ";Number of D accepted meson candidates;counts";
2036  h = fHistManager.CreateTH1(hname, htitle, 51, -0.5, 50.5);
2037 
2038  hname = TString::Format("%s/fHistNDmesons", param.GetName());
2039  htitle = hname + ";Number of D meson candidates;counts";
2040  h = fHistManager.CreateTH1(hname, htitle, 101, -0.5, 100.5);
2041 
2042  hname = TString::Format("%s/fHistNEvents", param.GetName());
2043  htitle = hname + ";Event status;counts";
2044  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
2045  h->GetXaxis()->SetBinLabel(1, "Accepted");
2046  h->GetXaxis()->SetBinLabel(2, "Rejected");
2047 
2048  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
2049  htitle = hname + ";Rejection reason;counts";
2050  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
2051 
2052  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
2053  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
2054  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2055 
2056  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
2057  htitle = hname + ";#it{#eta}_{D};counts";
2058  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2059 
2060  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
2061  htitle = hname + ";#it{#phi}_{D};counts";
2062  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2063 
2064  if (param.fCandidateType == kD0toKpi || param.fCandidateType == kD0toKpiLikeSign) {
2065  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
2066  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2067  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
2068  }
2069  else if (param.fCandidateType == kDstartoKpipi) {
2070  Double_t min = 0;
2071  Double_t max = 0;
2072 
2073  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
2074  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2075  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
2076  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
2077 
2078  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2079  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
2080  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2081  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
2082  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
2083  }
2084 
2085  if (param.fMCMode == kBackgroundOnly || param.fMCMode == kSignalOnly || param.fMCMode == kMCTruth) {
2086  hname = TString::Format("%s/fHistDmesonOrigin", param.GetName());
2087  htitle = hname + ";origin;counts";
2088  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2089  }
2090 
2091  for (std::vector<AliHFJetDefinition>::iterator itdef = param.fJetDefinitions.begin(); itdef != param.fJetDefinitions.end(); itdef++) {
2092  AliHFJetDefinition* jetDef = &(*itdef);
2093  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef->GetName());
2094 
2095  fHistManager.CreateHistoGroup(jetDef->GetName(), param.GetName());
2096 
2097  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef->GetName());
2098  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2099  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2100  SetRejectionReasonLabels(h->GetXaxis());
2101 
2102  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef->GetName());
2103  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2104  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2105  SetRejectionReasonLabels(h->GetXaxis());
2106 
2107  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef->GetName());
2108  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
2109  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2110  SetRejectionReasonLabels(h->GetXaxis());
2111 
2112  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef->GetName());
2113  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
2114  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
2115 
2116  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef->GetName());
2117  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
2118  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2119 
2120  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef->GetName());
2121  htitle = hname + ";#it{#eta}_{jet};counts";
2122  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2123 
2124  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef->GetName());
2125  htitle = hname + ";#it{#phi}_{jet};counts";
2126  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2127  }
2128  switch (fOutputType) {
2129  case kTreeOutput:
2130  param.BuildTree(GetName());
2131  if (treeSlot < fNOutputTrees) {
2132  param.AssignDataSlot(treeSlot+2);
2133  treeSlot++;
2135  }
2136  else {
2137  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
2138  }
2139  break;
2140  case kTHnOutput:
2141  param.BuildHnSparse(fEnabledAxis);
2142  break;
2143  case kNoOutput:
2144  break;
2145  }
2146  }
2147 
2149 
2150  PostData(1, fOutput);
2151 }
2152 
2156 {
2158 
2159  // Load the event
2160  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
2161 
2162  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
2163 
2164  fFastJetWrapper->SetAreaType(fastjet::active_area);
2166 
2167  if (!fAodEvent) {
2168  AliError(Form("This task need an AOD event (Task '%s'). Expect troubles...", GetName()));
2169  //return;
2170  }
2171 
2172  for (auto &params : fAnalysisEngines) {
2173 
2174  params.fAodEvent = fAodEvent;
2175  params.fFastJetWrapper = fFastJetWrapper;
2176  if (fAodEvent) params.Init(fGeom, fAodEvent->GetRunNumber());
2177 
2178  if (params.fMCMode != kMCTruth && fAodEvent) {
2179  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
2180 
2181  if (params.fCandidateArray) {
2182  TString className;
2183  if (params.fCandidateType == kD0toKpi || params.fCandidateType == kD0toKpiLikeSign) {
2184  className = "AliAODRecoDecayHF2Prong";
2185  }
2186  else if (params.fCandidateType == kDstartoKpipi) {
2187  className = "AliAODRecoCascadeHF";
2188  }
2189  if (!params.fCandidateArray->GetClass()->InheritsFrom(className)) {
2190  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2191  "%s: Objects of type %s in %s are not inherited from %s! Task will be disabled!",
2192  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName(), className.Data());
2193  params.fCandidateArray = 0;
2194  params.fInhibit = kTRUE;
2195  }
2196  }
2197  else {
2198  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2199  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
2200  params.fBranchName.Data(), params.GetName());
2201  params.fInhibit = kTRUE;
2202  }
2203  }
2204 
2205  if (params.fMCMode != kNoMC) {
2206  params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(0));
2207 
2208  if (!params.fMCContainer) params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(1));
2209 
2210  if (!params.fMCContainer) {
2211  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2212  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
2213  params.GetName());
2214  params.fInhibit = kTRUE;
2215  }
2216  }
2217 
2218  if (params.fMCMode != kMCTruth) {
2219  params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(0));
2220  if (!params.fTrackContainer) params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(1));
2221 
2222  params.fClusterContainer = GetClusterContainer(0);
2223 
2224  if (!params.fTrackContainer && !params.fClusterContainer) {
2225  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2226  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
2227  params.GetName());
2228  params.fInhibit = kTRUE;
2229  }
2230  }
2231  }
2232 }
2233 
2238 {
2239  //if (!fAodEvent) return kFALSE;
2240 
2241  TString hname;
2242 
2243  // fix for temporary bug in ESDfilter
2244  // the AODs with null vertex pointer didn't pass the PhysSel
2245  // Now adding an entry in teh histogram so as to check that this is actually cutting anything out
2246  if (fAodEvent && (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001)) {
2247  for (auto &eng : fAnalysisEngines) {
2248  if (eng.fInhibit) continue;
2249  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2250  fHistManager.FillTH1(hname, "ESDfilterBug");
2251  }
2252  return kFALSE;
2253  }
2254 
2255  for (auto &eng : fAnalysisEngines) {
2256 
2257  if (eng.fInhibit) continue;
2258 
2259  //Event selection
2260  hname = TString::Format("%s/fHistNEvents", eng.GetName());
2261  if (fAodEvent) {
2262  Bool_t iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
2263  if (!iseventselected) {
2264  fHistManager.FillTH1(hname, "Rejected");
2265  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2266  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
2267  TString label;
2268  do {
2269  label = GetHFEventRejectionReasonLabel(bitmap);
2270  if (label.IsNull()) break;
2271  fHistManager.FillTH1(hname, label);
2272  } while (true);
2273 
2274  continue;
2275  }
2276  }
2277 
2278  fHistManager.FillTH1(hname, "Accepted");
2279 
2280  AliDebug(2, "Event selected");
2281 
2282  eng.RunAnalysis();
2283  }
2284  return kTRUE;
2285 }
2286 
2291 {
2292  TString hname;
2293  for (auto &param : fAnalysisEngines) {
2294 
2295  if (param.fInhibit) continue;
2296 
2297  if (fOutputType == kTreeOutput) {
2298  param.FillTree(fApplyKinematicCuts);
2300  }
2301  else if (fOutputType == kTHnOutput) {
2302  param.FillHnSparse(fApplyKinematicCuts);
2303  }
2304  }
2305  return kTRUE;
2306 }
2307 
2316 {
2317  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
2318 
2319  Double_t mass = part->Mass();
2320 
2321  // To make sure that the PDG mass value is not at the edge of a bin
2322  if (nbins % 2 == 0) {
2323  minMass = mass - range / 2 - range / nbins / 2;
2324  maxMass = mass + range / 2 - range / nbins / 2;
2325  }
2326  else {
2327  minMass = mass - range / 2;
2328  maxMass = mass + range / 2;
2329  }
2330 }
2331 
2338 {
2339  static TString label;
2340  label = "";
2341 
2342  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
2343  label = "NotSelTrigger";
2344  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
2345  return label.Data();
2346  }
2347  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
2348  label = "NoVertex";
2349  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
2350  return label.Data();
2351  }
2352  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
2353  label = "TooFewVtxContrib";
2354  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
2355  return label.Data();
2356  }
2357  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
2358  label = "ZVtxOutFid";
2359  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
2360  return label.Data();
2361  }
2362  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
2363  label = "Pileup";
2364  bitmap &= ~BIT(AliRDHFCuts::kPileup);
2365  return label.Data();
2366  }
2367  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
2368  label = "OutsideCentrality";
2369  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
2370  return label.Data();
2371  }
2372  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
2373  label = "PhysicsSelection";
2374  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
2375  return label.Data();
2376  }
2377  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
2378  label = "BadSPDVertex";
2379  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
2380  return label.Data();
2381  }
2382  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
2383  label = "ZVtxSPDOutFid";
2384  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
2385  return label.Data();
2386  }
2387  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
2388  label = "CentralityFlattening";
2389  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
2390  return label.Data();
2391  }
2392  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
2393  label = "BadTrackV0Correl";
2394  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
2395  return label.Data();
2396  }
2397 
2398  return label.Data();
2399 }
2400 
2407 {
2408  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
2409  PostData(eng.GetDataSlotNumber(), eng.GetTree());
2410  return eng.GetDataSlotNumber();
2411  }
2412  else {
2413  return -1;
2414  }
2415 }
2416 
2426 {
2427  // Get the pointer to the existing analysis manager via the static access method
2428  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
2429  if (!mgr) {
2430  ::Error("AddTaskDmesonJets", "No analysis manager to connect to.");
2431  return NULL;
2432  }
2433 
2434  // Check the analysis type using the event handlers connected to the analysis manager
2435  AliVEventHandler* handler = mgr->GetInputEventHandler();
2436  if (!handler) {
2437  ::Error("AddTaskEmcalJetSpectraQA", "This task requires an input event handler");
2438  return NULL;
2439  }
2440 
2441  EDataType_t dataType = kUnknownDataType;
2442 
2443  if (handler->InheritsFrom("AliESDInputHandler")) {
2444  dataType = kESD;
2445  }
2446  else if (handler->InheritsFrom("AliAODInputHandler")) {
2447  dataType = kAOD;
2448  }
2449 
2450  // Init the task and do settings
2451  if (ntracks == "usedefault") {
2452  if (dataType == kESD) {
2453  ntracks = "Tracks";
2454  }
2455  else if (dataType == kAOD) {
2456  ntracks = "tracks";
2457  }
2458  else {
2459  ntracks = "";
2460  }
2461  }
2462 
2463  if (nclusters == "usedefault") {
2464  if (dataType == kESD) {
2465  nclusters = "CaloClusters";
2466  }
2467  else if (dataType == kAOD) {
2468  nclusters = "caloClusters";
2469  }
2470  else {
2471  nclusters = "";
2472  }
2473  }
2474 
2475  if (nMCpart == "usedefault") {
2476  nMCpart = "mcparticles"; // Always needs AliAODMCParticle objects
2477  }
2478 
2479  TString name("AliAnalysisTaskDmesonJets");
2480  if (strcmp(suffix, "") != 0) {
2481  name += TString::Format("_%s", suffix.Data());
2482  }
2483 
2484  AliAnalysisTaskDmesonJets* jetTask = new AliAnalysisTaskDmesonJets(name, nMaxTrees);
2485  jetTask->SetVzRange(-10,10);
2486 
2487  if (!ntracks.IsNull()) {
2488  AliHFTrackContainer* trackCont = new AliHFTrackContainer(ntracks);
2489  jetTask->AdoptParticleContainer(trackCont);
2490  }
2491 
2492  if (!nMCpart.IsNull()) {
2493  AliMCParticleContainer* partCont = new AliHFAODMCParticleContainer(nMCpart);
2494  partCont->SetEtaLimits(-1.5, 1.5);
2495  jetTask->AdoptParticleContainer(partCont);
2496  }
2497 
2498  jetTask->AddClusterContainer(nclusters);
2499 
2500  // Final settings, pass to manager and set the containers
2501  mgr->AddTask(jetTask);
2502 
2503  // Create containers for input/output
2504  AliAnalysisDataContainer* cinput1 = mgr->GetCommonInputContainer();
2505  TString contname1(name);
2506  contname1 += "_histos";
2507  AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(contname1.Data(),
2508  TList::Class(), AliAnalysisManager::kOutputContainer,
2509  Form("%s", AliAnalysisManager::GetCommonFileName()));
2510 
2511  mgr->ConnectInput(jetTask, 0, cinput1);
2512  mgr->ConnectOutput(jetTask, 1, coutput1);
2513 
2514  for (Int_t i = 0; i < nMaxTrees; i++) {
2515  TString contname = TString::Format("%s_tree_%d", name.Data(), i);
2516  AliAnalysisDataContainer *coutput = mgr->CreateContainer(contname.Data(),
2517  TTree::Class(),AliAnalysisManager::kOutputContainer,
2518  Form("%s", AliAnalysisManager::GetCommonFileName()));
2519  mgr->ConnectOutput(jetTask, 2+i, coutput);
2520  }
2521  return jetTask;
2522 }
2523 
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:98
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:103
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:97
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:101
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:99
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)