AliPhysics  2797316 (2797316)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskDmesonJets.cxx
Go to the documentation of this file.
1 /*************************************************************************
2 * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15 
16 // Root
17 #include <TClonesArray.h>
18 #include <TDatabasePDG.h>
19 #include <TParticlePDG.h>
20 #include <TVector3.h>
21 #include <THnSparse.h>
22 #include <TParticle.h>
23 #include <TMath.h>
24 #include <THashList.h>
25 #include <TFile.h>
26 
27 // Aliroot general
28 #include "AliLog.h"
29 #include "AliEMCALGeometry.h"
30 
31 // Aliroot HF
32 #include "AliAODRecoCascadeHF.h"
34 #include "AliRDHFCutsD0toKpi.h"
37 #include "AliHFTrackContainer.h"
38 
39 // Aliroot EMCal jet framework
40 #include "AliEmcalJetTask.h"
41 #include "AliEmcalJet.h"
42 #include "AliJetContainer.h"
43 #include "AliParticleContainer.h"
44 #include "AliEmcalParticle.h"
45 #include "AliFJWrapper.h"
46 
48 
49 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfo
50 
54 
61 Double_t AliAnalysisTaskDmesonJets::AliJetInfo::GetDistance(const AliJetInfo& jet, Double_t& deta, Double_t& dphi) const
62 {
63  dphi = TVector2::Phi_mpi_pi(fMomentum.Phi() - jet.Phi());;
64  deta = fMomentum.Eta() - jet.Eta();
65  return TMath::Sqrt(dphi*dphi + deta*deta);
66 }
67 
73 {
74  Double_t deta = 0;
75  Double_t dphi = 0;
76  return GetDistance(jet, deta, dphi);
77 }
78 
79 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonJetInfo
80 
84 
87 {
88  fD.SetPtEtaPhiE(0,0,0,0);
89  fSoftPionPt = 0;
90  fInvMass2Prong = 0;
91  fDmesonParticle = 0;
92  fMCLabel = -1;
93  fReconstructed = kFALSE;
94  for (auto &jet : fJets) {
95  jet.second.fMomentum.SetPtEtaPhiE(0,0,0,0);
96  jet.second.fNConstituents = 0;
97  jet.second.fNEF = 0;
98  jet.second.fMaxChargedPt = 0;
99  jet.second.fMaxNeutralPt = 0;
100  }
101 }
102 
105 {
106  Printf("Printing D Meson Jet object.");
107  Printf("D Meson: pT = %.3f, eta = %.3f, phi = %.3f, inv. mass = %.3f", fD.Pt(), fD.Eta(), fD.Phi_0_2pi(), fD.M());
108  Printf("Soft pion pT: %.3f. 2-Prong Invariant mass = %.3f", fSoftPionPt, fInvMass2Prong);
109  for (auto &jet : fJets) {
110  Printf("Jet %s: pT = %.3f, eta = %.3f, phi = %.3f", jet.first.c_str(), jet.second.Pt(), jet.second.Eta(), jet.second.Phi_0_2pi());
111  Printf("Jet N Consituents = %d", jet.second.fNConstituents);
112  }
113 }
114 
119 {
120  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
121  if (it == fJets.end()) return 0;
122 
123  Double_t z = 0;
124 
125  if ((*it).second.Pt() > 0) {
126  TVector3 dvect = fD.Vect();
127  TVector3 jvect = (*it).second.fMomentum.Vect();
128 
129  Double_t jetMom = jvect * jvect;
130 
131  if (jetMom < 1e-6) {
132  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetZ", "Zero jet momentum!");
133  z = 0.999;
134  }
135  else {
136  z = (dvect * jvect) / jetMom;
137  }
138 
139  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin 0.9-1 rather than 1-1.1
140  }
141 
142  return z;
143 }
144 
151 Double_t AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetDistance(std::string n, Double_t& deta, Double_t& dphi) const
152 {
153  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
154  if (it == fJets.end()) return 0;
155 
156  return GetDistance((*it).second, deta, dphi);
157 }
158 
164 {
165  Double_t deta = 0;
166  Double_t dphi = 0;
167  return GetDistance(n, deta, dphi);
168 }
169 
176 Double_t AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetDistance(const AliJetInfo& jet, Double_t& deta, Double_t& dphi) const
177 {
178  dphi = TVector2::Phi_mpi_pi(fD.Phi() - jet.Phi());;
179  deta = fD.Eta() - jet.Eta();
180  return TMath::Sqrt(dphi*dphi + deta*deta);
181 }
182 
188 {
189  Double_t deta = 0;
190  Double_t dphi = 0;
191  return GetDistance(jet, deta, dphi);
192 }
193 
199 {
200  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
201  if (it == fJets.end()) {
202  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
203  n.c_str());
204  return 0;
205  }
206  return &((*it).second);
207 }
208 
214 {
215  std::map<std::string, AliJetInfo>::iterator it = fJets.find(n);
216  if (it == fJets.end()) {
217  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
218  n.c_str());
219  return 0;
220  }
221  return &((*it).second);
222 }
223 
224 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoSummary
225 
229 
235  fPt(0),
236  fEta(0),
237  fPhi(0),
238  fR(0),
239  fZ(0)
240 {
241  Set(source, n);
242 }
243 
246 {
247  fPt = 0;
248  fEta = 0;
249  fPhi = 0;
250  fR = 0;
251  fZ = 0;
252 }
253 
259 {
260  std::map<std::string, AliJetInfo>::const_iterator it = source.fJets.find(n);
261  if (it == source.fJets.end()) return;
262 
263  fPt = (*it).second.Pt();
264  fEta = (*it).second.Eta();
265  fPhi = (*it).second.Phi_0_2pi();
266  fR = source.GetDistance(n);
267  fZ = source.GetZ(n);
268 }
269 
274 {
275  fPt = source.Pt();
276  fEta = source.Eta();
277  fPhi = source.Phi_0_2pi();
278  fR = 0;
279  fZ = 0;
280 }
281 
282 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonInfoSummary
283 
287 
292  fPt(0),
293  fEta(0),
294  fPhi(0)
295 {
296  Set(source);
297 }
298 
303 {
304  fPt = source.fD.Pt();
305  fEta = source.fD.Eta();
306  fPhi = source.fD.Phi_0_2pi();
307 }
308 
311 {
312  fPt = 0;
313  fEta = 0;
314  fPhi = 0;
315 }
316 
317 // Definitions of class AliAnalysisTaskDmesonJets::AliD0InfoSummary
318 
322 
327  AliDmesonInfoSummary(source),
328  fInvMass(source.fD.M())
329 {
330 }
331 
336 {
337  fInvMass = source.fD.M();
339 }
340 
343 {
345 
346  fInvMass = 0;
347 }
348 
349 // Definitions of class AliAnalysisTaskDmesonJets::AliDStarInfoSummary
350 
354 
359  AliDmesonInfoSummary(source),
360  f2ProngInvMass(source.fInvMass2Prong),
361  fDeltaInvMass(source.fD.M() - source.fInvMass2Prong)
362 {
363 }
364 
369 {
370  f2ProngInvMass = source.fInvMass2Prong;
371  fDeltaInvMass = source.fD.M() - source.fInvMass2Prong;
373 }
374 
377 {
379 
380  f2ProngInvMass = 0;
381  fDeltaInvMass = 0;
382 }
383 
384 // Definitions of class AliAnalysisTaskDmesonJets::AliJetDefinition
385 
389 
392  TObject(),
393  fJetType(AliJetContainer::kChargedJet),
394  fRadius(0),
395  fJetAlgo(AliJetContainer::antikt_algorithm),
396  fRecoScheme(AliJetContainer::pt_scheme),
397  fAcceptance(AliEmcalJet::kUser),
398  fMinJetPt(0.),
399  fMaxJetPt(500.),
400  fMinJetPhi(0.),
401  fMaxJetPhi(0.),
402  fMinJetEta(0.),
403  fMaxJetEta(0.),
404  fMinChargedPt(0.),
405  fMaxChargedPt(0.),
406  fMinNeutralPt(0.),
407  fMaxNeutralPt(0.)
408 {
409 }
410 
418  TObject(),
419  fJetType(type),
420  fRadius(r),
421  fJetAlgo(algo),
422  fRecoScheme(reco),
423  fAcceptance(AliEmcalJet::kUser),
424  fMinJetPt(0.),
425  fMaxJetPt(500.),
426  fMinJetPhi(0.),
427  fMaxJetPhi(0.),
428  fMinJetEta(0.),
429  fMaxJetEta(0.),
430  fMinChargedPt(0.),
431  fMaxChargedPt(0.),
432  fMinNeutralPt(0.),
433  fMaxNeutralPt(0.)
434 {
435  // By default set detector fiducial acceptance
436  switch (type) {
440  break;
443  break;
444  }
445 }
446 
451  TObject(),
452  fJetType(source.fJetType),
453  fRadius(source.fRadius),
454  fJetAlgo(source.fJetAlgo),
455  fRecoScheme(source.fRecoScheme),
456  fAcceptance(source.fAcceptance),
457  fMinJetPt(source.fMinJetPt),
458  fMaxJetPt(source.fMaxJetPt),
459  fMinJetPhi(source.fMinJetPhi),
460  fMaxJetPhi(source.fMaxJetPhi),
461  fMinJetEta(source.fMinJetEta),
462  fMaxJetEta(source.fMaxJetEta),
463  fMinChargedPt(source.fMinChargedPt),
464  fMaxChargedPt(source.fMaxChargedPt),
465  fMinNeutralPt(source.fMinNeutralPt),
466  fMaxNeutralPt(source.fMaxNeutralPt)
467 {
468 }
469 
474 {
475  new (this) AliHFJetDefinition(source);
476  return *this;
477 }
478 
481 {
482  static TString name;
483 
484  name = AliJetContainer::GenerateJetName(fJetType, fJetAlgo, fRecoScheme, fRadius, 0, 0, "Jet");
485 
486  return name.Data();
487 }
488 
494 {
495  if (fMinJetEta < fMaxJetEta && (jet.Eta() < fMinJetEta || jet.Eta() > fMaxJetEta)) return kFALSE;
496  if (fMinJetPhi < fMaxJetPhi && (jet.Phi() < fMinJetPhi || jet.Phi() > fMaxJetPhi)) return kFALSE;
497  if (jet.Pt() > fMaxJetPt || jet.Pt() < fMinJetPt) return kFALSE;
498  if (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt) return kFALSE;
499  if (jet.fMaxNeutralPt < fMinNeutralPt || jet.fMaxNeutralPt > fMaxNeutralPt) return kFALSE;
500 
501  return kTRUE;
502 }
503 
509 {
510  const AliJetInfo* jet = dMesonJet.GetJet(n);
511  if (!jet) return kFALSE;
512  return IsJetInAcceptance((*jet));
513 }
514 
519 void AliAnalysisTaskDmesonJets::AliHFJetDefinition::SetDetectorJetEtaPhiRange(const AliEMCALGeometry* const geom, Int_t run)
520 {
521  Double_t r = 0;
522  switch (fAcceptance) {
524  r = fRadius;
525  // enforce fiducial acceptance
526  /* no break */
527  case AliEmcalJet::kTPC:
528  SetJetEtaRange(-0.9 + r, 0.9 - r);
529  SetJetPhiRange(0, 0); // No cut on phi
530  break;
531 
533  r = fRadius;
534  // enforce fiducial acceptance
535  /* no break */
536  case AliEmcalJet::kEMCAL:
537  if (geom) {
538  SetJetEtaRange(geom->GetArm1EtaMin() + r, geom->GetArm1EtaMax() - r);
539 
540  if(run>=177295 && run<=197470) {//small SM masked in 2012 and 2013
541  SetJetPhiRange(1.405 + r,3.135 - r);
542  }
543  else {
544  SetJetPhiRange(geom->GetArm1PhiMin() * TMath::DegToRad() + r, geom->GetEMCALPhiMax() * TMath::DegToRad() - r);
545  }
546  }
547  else {
548  AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
549  SetJetEtaRange(-0.7 + r, 0.7 - r);
550  SetJetPhiRange(1.405 + r, 3.135 - r);
551  }
552  break;
553 
555  r = fRadius;
556  // enforce fiducial acceptance
557  /* no break */
558  case AliEmcalJet::kDCAL:
559  if (geom) {
560  SetJetEtaRange(geom->GetArm1EtaMin() + r, geom->GetArm1EtaMax() - r);
561  SetJetPhiRange(geom->GetDCALPhiMin() * TMath::DegToRad() + r, geom->GetDCALPhiMax() * TMath::DegToRad() - r);
562  }
563  else {
564  AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for DCAL year 2015!!");
565  SetJetEtaRange(-0.7 + r, 0.7 - r);
566  SetJetPhiRange(4.538 + r, 5.727 - r);
567  }
568  break;
569 
570  case AliEmcalJet::kUser:
571  // Nothing to be done
572  break;
573  }
574 }
575 
582 {
583  if (lhs.fJetType > rhs.fJetType) return false;
584  else if (lhs.fJetType < rhs.fJetType) return true;
585  else {
586  if (lhs.fRadius > rhs.fRadius) return false;
587  else if (lhs.fRadius < rhs.fRadius) return true;
588  else {
589  if (lhs.fJetAlgo > rhs.fJetAlgo) return false;
590  else if (lhs.fJetAlgo < rhs.fJetAlgo) return true;
591  else {
592  if (lhs.fRecoScheme < rhs.fRecoScheme) return true;
593  else return false;
594  }
595  }
596  }
597 }
598 
605 {
606  if (lhs.fJetType != rhs.fJetType) return false;
607  if (lhs.fRadius != rhs.fRadius) return false;
608  if (lhs.fJetAlgo != rhs.fJetAlgo) return false;
609  if (lhs.fRecoScheme != rhs.fRecoScheme) return false;
610  return true;
611 }
612 
613 // Definitions of class AliAnalysisTaskDmesonJets::AnalysisEngine
614 
618 
621  TObject(),
622  fCandidateType(kD0toKpi),
623  fCandidateName(),
624  fCandidatePDG(0),
625  fNDaughters(0),
626  fPDGdaughters(),
627  fBranchName(),
628  fMCMode(kNoMC),
629  fNMassBins(0),
630  fMinMass(0),
631  fMaxMass(0),
632  fRDHFCuts(0),
633  fRejectedOrigin(0),
634  fAcceptedDecay(0),
635  fInhibit(kFALSE),
636  fJetDefinitions(),
637  fPtBinWidth(0.5),
638  fMaxPt(100),
639  fDataSlotNumber(-1),
640  fTree(0),
641  fCurrentDmesonJetInfo(0),
642  fCurrentJetInfo(0),
643  fCandidateArray(0),
644  fMCContainer(0),
645  fTrackContainer(0),
646  fClusterContainer(0),
647  fAodEvent(0),
648  fFastJetWrapper(0),
649  fHistManager(0)
650 {
651 }
652 
661  TObject(),
662  fCandidateType(type),
663  fCandidateName(),
664  fCandidatePDG(0),
665  fNDaughters(0),
666  fPDGdaughters(),
667  fBranchName(),
668  fMCMode(MCmode),
669  fNMassBins(nMassBins),
670  fMinMass(0),
671  fMaxMass(0),
672  fRDHFCuts(cuts),
674  fAcceptedDecay(kAnyDecay),
675  fInhibit(kFALSE),
676  fJetDefinitions(),
677  fPtBinWidth(0.5),
678  fMaxPt(100),
679  fDataSlotNumber(-1),
680  fTree(0),
681  fCurrentDmesonJetInfo(0),
682  fCurrentJetInfo(0),
683  fCandidateArray(0),
684  fMCContainer(0),
685  fTrackContainer(0),
686  fClusterContainer(0),
687  fAodEvent(0),
688  fFastJetWrapper(0),
689  fHistManager(0)
690 {
691  SetCandidateProperties(range);
692 }
693 
698  TObject(source),
699  fCandidateType(source.fCandidateType),
700  fCandidateName(source.fCandidateName),
701  fCandidatePDG(source.fCandidatePDG),
702  fNDaughters(source.fNDaughters),
703  fPDGdaughters(source.fPDGdaughters),
704  fBranchName(source.fBranchName),
705  fMCMode(source.fMCMode),
706  fNMassBins(source.fNMassBins),
707  fMinMass(source.fMinMass),
708  fMaxMass(source.fMaxMass),
709  fRDHFCuts(),
710  fRejectedOrigin(source.fRejectedOrigin),
711  fAcceptedDecay(source.fAcceptedDecay),
712  fInhibit(source.fInhibit),
713  fJetDefinitions(source.fJetDefinitions),
714  fPtBinWidth(source.fPtBinWidth),
715  fMaxPt(source.fMaxPt),
716  fDataSlotNumber(-1),
717  fTree(0),
718  fCurrentDmesonJetInfo(0),
719  fCurrentJetInfo(0),
720  fCandidateArray(source.fCandidateArray),
721  fMCContainer(source.fMCContainer),
722  fTrackContainer(source.fTrackContainer),
723  fClusterContainer(source.fClusterContainer),
724  fAodEvent(source.fAodEvent),
726  fHistManager(source.fHistManager)
727 {
728  SetRDHFCuts(source.fRDHFCuts);
729 }
730 
731 // Destructor
733 {
734  delete fRDHFCuts;
735 }
736 
741 {
742  new (this) AnalysisEngine(source);
743  return *this;
744 }
745 
750 {
751  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
752  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
753  }
754 
755  return kFALSE;
756 }
757 
759 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const geom, Int_t runNumber)
760 {
761  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
762  fJetDefinitions[i].SetDetectorJetEtaPhiRange(geom, runNumber);
763  }
764 }
765 
770 {
771  switch (fCandidateType) {
772  case kD0toKpi:
773  fCandidatePDG = 421;
774  fCandidateName = "D0";
775  fNDaughters = 2;
776  fPDGdaughters.Set(fNDaughters);
777  fPDGdaughters.Reset();
778  fPDGdaughters[0] = 211; // pi
779  fPDGdaughters[1] = 321; // K
780  fBranchName = "D0toKpi";
781  fAcceptedDecay = kDecayD0toKpi;
782  if (!fRDHFCuts) {
783  fRDHFCuts = new AliRDHFCutsD0toKpi();
784  fRDHFCuts->SetStandardCutsPP2010();
785  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
786  fRDHFCuts->SetTriggerClass("","");
787  }
788  break;
789  case kDstartoKpipi:
790  fCandidatePDG = 413;
791  fCandidateName = "DStar";
792  fNDaughters = 3;
793  fPDGdaughters.Set(fNDaughters);
794  fPDGdaughters.Reset();
795  fPDGdaughters[0] = 211; // pi soft
796  fPDGdaughters[1] = 211; // pi fromD0
797  fPDGdaughters[2] = 321; // K from D0
798  fBranchName = "Dstar";
799  fAcceptedDecay = kDecayDStartoKpipi;
800  if (!fRDHFCuts) {
801  fRDHFCuts = new AliRDHFCutsDStartoKpipi();
802  fRDHFCuts->SetStandardCutsPP2010();
803  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
804  fRDHFCuts->SetTriggerClass("","");
805  }
806  break;
807  default:
808  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
809  }
810 
811  CalculateMassLimits(range, fCandidatePDG, fNMassBins, fMinMass, fMaxMass);
812 }
813 
818 {
819  if (fRDHFCuts) delete fRDHFCuts;
820  fRDHFCuts = cuts;
821 }
822 
827 {
828  if (!cuts) return;
829  if (fRDHFCuts) delete fRDHFCuts;
830  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
831 }
832 
837 {
838  static TString name;
839 
840  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
841 
842  return name.Data();
843 }
844 
849 {
850  static TString name;
851 
852  name = fCandidateName;
853  switch (fMCMode) {
854  case kBackgroundOnly:
855  name += "_kBackgroundOnly";
856  break;
857  case kSignalOnly:
858  name += "_kSignalOnly";
859  break;
860  case kMCTruth:
861  name += "_MCTruth";
862  break;
863  default:
864  break;
865  }
866 
867  return name.Data();
868 }
869 
877 {
878  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
879 
880  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
881  fJetDefinitions.push_back(def);
882  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
883  "Total number of jet definitions is now %lu.",
884  def.GetName(), GetName(), fJetDefinitions.size());
885  // For detector level set maximum track pt to 100 GeV/c
886  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
887  }
888  else {
889  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
890  }
891 
892  return &(*it);
893 }
894 
906 {
907  AliHFJetDefinition def(type, r, algo, reco);
908 
909  return AddJetDefinition(def);
910 }
911 
917 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
918 {
919  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
920  while (it != fJetDefinitions.end() && (*it) != def) it++;
921  return it;
922 }
923 
930 {
931  if (lhs.fCandidateType > rhs.fCandidateType) return false;
932  else if (lhs.fCandidateType < rhs.fCandidateType) return true;
933  else {
934  if (lhs.fMCMode < rhs.fMCMode) return true;
935  else return false;
936  }
937 }
938 
945 {
946  if (lhs.fCandidateType != rhs.fCandidateType) return false;
947  if (lhs.fMCMode != rhs.fMCMode) return false;
948  return true;
949 }
950 
959 {
960  if (fCandidateType == kD0toKpi) { // D0 candidate
961  return ExtractD0Attributes(Dcand, DmesonJet, i);
962  }
963  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
964  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
965  }
966  else {
967  return kFALSE;
968  }
969 }
970 
979 {
980  Int_t MCtruthPdgCode = 0;
981 
982  Double_t invMassD = 0;
983 
984  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
985  // Checks also the origin, and if it matches the rejected origin mask, return false
986  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
987  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
988  DmesonJet.fMCLabel = mcLab;
989 
990  // Retrieve the generated particle (if exists) and its PDG code
991  if (mcLab >= 0) {
992  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
993 
994  if (aodMcPart) {
995  // Check origin and return false if it matches the rejected origin mask
996  if (fRejectedOrigin) {
997  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
998 
999  UInt_t rs = origin;
1000  UShort_t p = 0;
1001  while (rs >>= 1) { p++; }
1002  TString hname = TString::Format("%s/fHistDmesonOrigin", GetName());
1003  fHistManager->FillTH1(hname, p);
1004 
1005  if ((origin & fRejectedOrigin) == origin) return kFALSE;
1006  }
1007  MCtruthPdgCode = aodMcPart->PdgCode();
1008  }
1009  }
1010  }
1011 
1012  //AliDebug(2,"Checking if D0 meson is selected");
1013  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
1014  if (isSelected == 1) { // selected as a D0
1015  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
1016 
1017  if (fMCMode == kNoMC ||
1018  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1019  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly)) {
1020  // 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)
1021  //AliDebug(2,"Selected as D0");
1022  invMassD = Dcand->InvMassD0();
1023  }
1024  else { // conditions above not passed, so return FALSE
1025  return kFALSE;
1026  }
1027  }
1028  else if (isSelected == 2) { // selected as a D0bar
1029  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
1030 
1031  if (fMCMode == kNoMC ||
1032  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1033  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly)) {
1034  // 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)
1035  //AliDebug(2,"Selected as D0bar");
1036  invMassD = Dcand->InvMassD0bar();
1037  }
1038  else { // conditions above not passed, so return FALSE
1039  return kFALSE;
1040  }
1041  }
1042  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
1043  //AliDebug(2,"Selected as either D0 or D0bar");
1044 
1045  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
1046  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1047  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kBackgroundOnly)) {
1048  if (i > 0) return kFALSE;
1049  //AliDebug(2, "MC truth is D0");
1050  invMassD = Dcand->InvMassD0();
1051  }
1052  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1053  (MCtruthPdgCode == fCandidatePDG && fMCMode == kBackgroundOnly)) {
1054  if (i > 0) return kFALSE;
1055  //AliDebug(2, "MC truth is D0bar");
1056  invMassD = Dcand->InvMassD0bar();
1057  }
1058  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
1059 
1060  // Only accept it if background-only OR background-and-signal was requested
1061  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
1062  // Select D0 or D0bar depending on the i-parameter
1063  if (i == 0) {
1064  //AliDebug(2, "Returning invariant mass with D0 hypothesis");
1065  invMassD = Dcand->InvMassD0();
1066  }
1067  else if (i == 1) {
1068  //AliDebug(2, "Returning invariant mass with D0bar hypothesis");
1069  invMassD = Dcand->InvMassD0bar();
1070  }
1071  else { // i > 1
1072  return kFALSE;
1073  }
1074  }
1075  else { // signal-only was requested but this is not a true D0
1076  return kFALSE;
1077  }
1078  }
1079  }
1080  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
1081  return kTRUE;
1082 }
1083 
1092 {
1093  if (i > 0) return kFALSE; // only one mass hypothesis for the D*
1094 
1095  Int_t MCtruthPdgCode = 0;
1096 
1097  Double_t invMassD = 0;
1098 
1099  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1100  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1101  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1102 
1103  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1104  //AliDebug(2, Form("MC label is %d", mcLab));
1105  DmesonJet.fMCLabel = mcLab;
1106  if (mcLab >= 0) {
1107  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1108 
1109  if (aodMcPart) {
1110  if (fRejectedOrigin) {
1111  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
1112 
1113  UInt_t rs = origin;
1114  UShort_t p = 0;
1115  while (rs >>= 1) { p++; }
1116  TString hname = TString::Format("%s/fHistDmesonOrigin", GetName());
1117  fHistManager->FillTH1(hname, p);
1118 
1119  if ((origin & fRejectedOrigin) == origin) return kFALSE;
1120  }
1121 
1122  MCtruthPdgCode = aodMcPart->PdgCode();
1123  //AliDebug(2, Form("MC truth pdg code is %d",MCtruthPdgCode));
1124  }
1125  }
1126  }
1127 
1128  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1129  if (fMCMode == kNoMC ||
1130  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1131  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1132  // 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)
1133  invMassD = DstarCand->InvMassDstarKpipi();
1134  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1135  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1136  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1137  return kTRUE;
1138  }
1139  else { // conditions above not passed, so return FALSE
1140  return kFALSE;
1141  }
1142 }
1143 
1151 {
1152  if (!part) return kDecayOther;
1153  if (!mcArray) return kDecayOther;
1154 
1156 
1157  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1158 
1159  if (part->GetNDaughters() == 2) {
1160 
1161  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1162  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1163 
1164  if (!d1 || !d2) {
1165  return decay;
1166  }
1167 
1168  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1169  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1170 
1171  if (absPdgPart == 421) { // D0 -> K pi
1172  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1173  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1174  decay = kDecayD0toKpi;
1175  }
1176  }
1177 
1178  if (absPdgPart == 413) { // D* -> D0 pi
1179  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1180  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1181  if (D0decay == kDecayD0toKpi) {
1182  decay = kDecayDStartoKpipi;
1183  }
1184  }
1185  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1186  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1187  if (D0decay == kDecayD0toKpi) {
1188  decay = kDecayDStartoKpipi;
1189  }
1190  }
1191  }
1192  }
1193 
1194  return decay;
1195 }
1196 
1204 {
1205  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1206 
1207  if (!part) return kUnknownQuark;
1208  if (!mcArray) return kUnknownQuark;
1209 
1210  Int_t mother = part->GetMother();
1211  while (mother >= 0) {
1212  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1213  if (mcGranma) {
1214  Int_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
1215 
1216  if (abspdgGranma == 1) return kFromDown;
1217  if (abspdgGranma == 2) return kFromUp;
1218  if (abspdgGranma == 3) return kFromStrange;
1219  if (abspdgGranma == 4) return kFromCharm;
1220  if (abspdgGranma == 5) return kFromBottom;
1221  if (abspdgGranma == 6) return kFromTop;
1222 
1223  mother = mcGranma->GetMother();
1224  }
1225  else {
1226  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1227  break;
1228  }
1229  }
1230 
1231  return kUnknownQuark;
1232 }
1233 
1236 {
1237  fDmesonJets.clear();
1238 
1239  for (auto& jetDef : fJetDefinitions) {
1240  jetDef.fJets.clear();
1241  }
1242 
1243  if (fMCMode == kMCTruth) {
1244  RunParticleLevelAnalysis();
1245  }
1246  else {
1247  RunDetectorLevelAnalysis();
1248  }
1249 }
1250 
1256 {
1258  fFastJetWrapper->SetR(jetDef.fRadius);
1261 
1262  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1263  fTrackContainer->SetDMesonCandidate(0);
1264  AddInputVectors(fTrackContainer, 100);
1265  }
1266 
1267  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1268  AddInputVectors(fClusterContainer, -100);
1269  }
1270 
1271  // run jet finder
1272  fFastJetWrapper->Run();
1273 
1274  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1275 
1276  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1277  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1278 
1279  Double_t maxChPt = 0;
1280  Double_t maxNePt = 0;
1281  Double_t totalNeutralPt = 0;
1282 
1283  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1284  if (constituents[ic].user_index() >= 100) {
1285  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1286  }
1287  else if (constituents[ic].user_index() <= -100) {
1288  totalNeutralPt += constituents[ic].pt();
1289  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1290  }
1291  }
1292 
1293  jetDef.fJets.push_back(
1294  AliJetInfo(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E(),
1295  constituents.size(), totalNeutralPt / jets_incl[ijet].pt(), maxChPt, maxNePt)
1296  );
1297  }
1298 }
1299 
1308 {
1309  if (jetDef.fJets.size() == 0) FindJets(jetDef);
1310 
1311  Double_t d_closest = 999;
1312  AliJetInfo* jet_closest = 0;
1313  const AliJetInfo* truth_jet = 0;
1314  try {
1315  truth_jet = &(dmeson.fJets.at(jetDef.GetName()));
1316  }
1317  catch(...) {
1318  return jet_distance_pair(0, 999);
1319  }
1320 
1321  for (auto& jet : jetDef.fJets) {
1322  Double_t d = truth_jet->GetDistance(jet);
1323 
1324  if (d > dMax) continue;
1325  if (d < d_closest) {
1326  d_closest = d;
1327  jet_closest = &jet;
1328  }
1329  }
1330 
1331  if (jet_closest && applyKinCuts) {
1332  if (!jetDef.IsJetInAcceptance(*jet_closest)) jet_closest = 0;
1333  }
1334 
1335  if (jet_closest) {
1336  AliDebug(2, Form("Found closest jet (pt=%.3f, eta=%.3f, phi=%.3f) to d meson (pt=%.3f, eta=%.3f, phi=%.3f) with d = %.3f",
1337  jet_closest->Pt(), jet_closest->Eta(), jet_closest->Phi_0_2pi(),
1338  dmeson.fD.Pt(), dmeson.fD.Eta(), dmeson.fD.Phi_0_2pi(),
1339  d_closest));
1340  }
1341 
1342  return jet_distance_pair(jet_closest, d_closest);
1343 }
1344 
1347 {
1348  const Int_t nD = fCandidateArray->GetEntriesFast();
1349 
1350  AliDmesonJetInfo DmesonJet;
1351 
1352  Int_t nAccCharm = 0;
1353  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1354  Int_t isSelected = 0;
1355 
1356  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1357  if (!charmCand) continue;
1358 
1359  Int_t nprongs = charmCand->GetNProngs();
1360 
1361  if (fCandidateType == kDstartoKpipi) {
1362  if (!charmCand->InheritsFrom("AliAODRecoCascadeHF")) {
1363  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::RunDetectorLevelAnalysis","Candidate type is D* but object type is wrong (should be AliAODRecoCascadeHF)");
1364  continue;
1365  }
1366  }
1367 
1368  // region of interest + cuts
1369  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1370 
1371  //candidate selected by cuts and PID
1372  isSelected = fRDHFCuts->IsSelected(charmCand, AliRDHFCuts::kAll, fAodEvent); //selected
1373 
1374  if (!isSelected) continue;
1375 
1376  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1377  DmesonJet.Reset();
1378  DmesonJet.fDmesonParticle = charmCand;
1379  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1380  for (std::vector<AliHFJetDefinition>::iterator itdef = fJetDefinitions.begin(); itdef != fJetDefinitions.end(); itdef++) {
1381  if (!FindJet(charmCand, DmesonJet, *itdef)) {
1382  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1383  (*itdef).GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1384  }
1385  }
1386  fDmesonJets[icharm] = DmesonJet;
1387  }
1388  }
1389  nAccCharm++;
1390  } // end of D cand loop
1391 
1392  TString hname;
1393 
1394  hname = TString::Format("%s/fHistNAcceptedDmesons", GetName());
1395  fHistManager->FillTH1(hname, nAccCharm);
1396 
1397  hname = TString::Format("%s/fHistNDmesons", GetName());
1398  fHistManager->FillTH1(hname, nD);
1399 }
1400 
1411 {
1412  TString hname;
1413 
1415  fFastJetWrapper->SetR(jetDef.fRadius);
1418 
1419  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1420 
1421  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1422  fTrackContainer->SetDMesonCandidate(Dcand);
1423  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1424  AddInputVectors(fTrackContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1425 
1426  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1427  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1428  const TObjArray& daughters = fTrackContainer->GetDaughterList();
1429  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1430  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1431  if (!fTrackContainer->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1432  }
1433  }
1434 
1435  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1436  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1437  AddInputVectors(fClusterContainer, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1438  }
1439 
1440  // run jet finder
1441  fFastJetWrapper->Run();
1442 
1443  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1444 
1445  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1446  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1447 
1448  Bool_t isDmesonJet = kFALSE;
1449 
1450  Double_t maxChPt = 0;
1451  Double_t maxNePt = 0;
1452  Double_t totalNeutralPt = 0;
1453 
1454  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1455  if (constituents[ic].user_index() == 0) {
1456  isDmesonJet = kTRUE;
1457  }
1458  else if (constituents[ic].user_index() >= 100) {
1459  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1460  }
1461  else if (constituents[ic].user_index() <= -100) {
1462  totalNeutralPt += constituents[ic].pt();
1463  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1464  }
1465  }
1466 
1467  if (isDmesonJet) {
1468  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1469  DmesonJet.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1470  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1471  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1472  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1473 
1474  return kTRUE;
1475  }
1476  }
1477 
1478  return kFALSE;
1479 }
1480 
1484 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist)
1485 {
1486  AliEmcalIterableMomentumContainer itcont = cont->all_momentum();
1487  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1488  UInt_t rejectionReason = 0;
1489  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1490  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1491  continue;
1492  }
1493  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1494  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1495  }
1496 }
1497 
1500 {
1501  TString hname;
1502 
1503  fMCContainer->SetSpecialPDG(fCandidatePDG);
1504  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1505  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1506 
1507  if (!fMCContainer->IsSpecialPDGFound()) return;
1508 
1509  for (auto &jetDef : fJetDefinitions) {
1510  switch (jetDef.fJetType) {
1512  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
1513  break;
1515  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
1516  break;
1518  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
1519  break;
1520  }
1521 
1523  fFastJetWrapper->SetR(jetDef.fRadius);
1526 
1527  hname = TString::Format("%s/fHistDmesonOrigin", GetName());
1528  fMCContainer->SetHistOrigin(static_cast<TH1*>(fHistManager->FindObject(hname)));
1529  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1530  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1531  fMCContainer->SetHistOrigin(0);
1532 
1533  fFastJetWrapper->Run();
1534 
1535  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1536 
1537  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1538  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1539 
1540  Bool_t isDmesonJet = kFALSE;
1541 
1542  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1543  Int_t iPart = constituents[ic].user_index() - 100;
1544  AliVParticle* part = fMCContainer->GetParticle(iPart);
1545  if (!part) {
1546  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1547  continue;
1548  }
1549  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1550  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1551  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1552  std::pair<int, AliDmesonJetInfo> element;
1553  element.first = iPart;
1554  element.second = AliDmesonJetInfo();
1555  dMesonJetIt = fDmesonJets.insert(element).first;
1556  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1557  (*dMesonJetIt).second.fDmesonParticle = part;
1558  }
1559 
1560  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1561  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1562  }
1563  }
1564  }
1565  }
1566 }
1567 
1572 {
1573  TString classname;
1574  switch (fCandidateType) {
1575  case kD0toKpi:
1576  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1577  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1578  break;
1579  case kDstartoKpipi:
1580  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1581  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1582  break;
1583  }
1584  TString treeName = TString::Format("%s_%s", taskName, GetName());
1585  fTree = new TTree(treeName, treeName);
1586  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
1587  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1588  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1589  fCurrentJetInfo[i] = new AliJetInfoSummary();
1590  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1591  }
1592 
1593  return fTree;
1594 }
1595 
1600 {
1601  TString hname;
1602 
1603  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
1604 
1605  for (auto &jetDef : fJetDefinitions) {
1606 
1607  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
1608 
1609  Double_t radius = jetDef.fRadius;
1610 
1611  TString title[30] = {""};
1612  Int_t nbins[30] = {0 };
1613  Double_t min [30] = {0.};
1614  Double_t max [30] = {0.};
1615  Int_t dim = 0 ;
1616 
1617  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1618  nbins[dim] = nPtBins;
1619  min[dim] = 0;
1620  max[dim] = fMaxPt;
1621  dim++;
1622 
1623  if ((enabledAxis & kPositionD) != 0) {
1624  title[dim] = "#eta_{D}";
1625  nbins[dim] = 50;
1626  min[dim] = -1;
1627  max[dim] = 1;
1628  dim++;
1629 
1630  title[dim] = "#phi_{D} (rad)";
1631  nbins[dim] = 150;
1632  min[dim] = 0;
1633  max[dim] = TMath::TwoPi();
1634  dim++;
1635  }
1636 
1637  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1638  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1639  nbins[dim] = fNMassBins;
1640  min[dim] = fMinMass;
1641  max[dim] = fMaxMass;
1642  dim++;
1643  }
1644 
1645  if (fCandidateType == kD0toKpi) {
1646  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1647  nbins[dim] = fNMassBins;
1648  min[dim] = fMinMass;
1649  max[dim] = fMaxMass;
1650  dim++;
1651  }
1652 
1653  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1654  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1655  nbins[dim] = fNMassBins;
1656  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1657  dim++;
1658  }
1659 
1660  if (fCandidateType == kDstartoKpipi) {
1661  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1662  nbins[dim] = fNMassBins*6;
1663  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1664 
1665  // subtract mass of D0
1666  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1667  min[dim] -= D0mass;
1668  max[dim] -= D0mass;
1669 
1670  dim++;
1671  }
1672 
1673  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1674  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1675  nbins[dim] = 100;
1676  min[dim] = 0;
1677  max[dim] = 25;
1678  dim++;
1679  }
1680 
1681  title[dim] = "#it{z}_{D}";
1682  nbins[dim] = 110;
1683  min[dim] = 0;
1684  max[dim] = 1.10;
1685  dim++;
1686 
1687  if ((enabledAxis & kDeltaR) != 0) {
1688  title[dim] = "#Delta R_{D-jet}";
1689  nbins[dim] = 100;
1690  min[dim] = 0;
1691  max[dim] = radius * 1.5;
1692  dim++;
1693  }
1694 
1695  if ((enabledAxis & kDeltaEta) != 0) {
1696  title[dim] = "#eta_{D} - #eta_{jet}";
1697  nbins[dim] = 100;
1698  min[dim] = -radius * 1.2;
1699  max[dim] = radius * 1.2;
1700  dim++;
1701  }
1702 
1703  if ((enabledAxis & kDeltaPhi) != 0) {
1704  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1705  nbins[dim] = 100;
1706  min[dim] = -radius * 1.2;
1707  max[dim] = radius * 1.2;
1708  dim++;
1709  }
1710 
1711  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1712  nbins[dim] = nPtBins;
1713  min[dim] = 0;
1714  max[dim] = fMaxPt;
1715  dim++;
1716 
1717  if ((enabledAxis & kPositionJet) != 0) {
1718  title[dim] = "#eta_{jet}";
1719  nbins[dim] = 50;
1720  min[dim] = -1;
1721  max[dim] = 1;
1722  dim++;
1723 
1724  title[dim] = "#phi_{jet} (rad)";
1725  nbins[dim] = 150;
1726  min[dim] = 0;
1727  max[dim] = TMath::TwoPi();
1728  dim++;
1729  }
1730 
1731  if ((enabledAxis & kJetConstituents) != 0) {
1732  title[dim] = "No. of constituents";
1733  nbins[dim] = 50;
1734  min[dim] = -0.5;
1735  max[dim] = 49.5;
1736  dim++;
1737  }
1738 
1739  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1740  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1741  for (Int_t j = 0; j < dim; j++) {
1742  h->GetAxis(j)->SetTitle(title[j]);
1743  }
1744  }
1745 }
1746 
1751 {
1752  TString hname;
1753 
1754  for (auto& dmeson_pair : fDmesonJets) {
1755  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
1756  Int_t accJets = 0;
1757  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1758  fCurrentJetInfo[ij]->Reset();
1759  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1760  if (!jet) continue;
1761  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1762  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1763  fHistManager->FillTH1(hname, jet->Pt());
1764  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1765  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1766  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1767  fHistManager->FillTH1(hname, jet->Eta());
1768  continue;
1769  }
1770  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
1771  accJets++;
1772  }
1773  if (kTRUE || accJets > 0) { // always fill D meson tree, even if no jet was accepted
1774  fTree->Fill();
1775  }
1776  else { // this is never executed
1777  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1778  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1779  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1780  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1781  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1782  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1783  if (fCandidateType == kD0toKpi) {
1784  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1785  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1786  }
1787  else if (fCandidateType == kDstartoKpipi) {
1788  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1789  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1790 
1791  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1792  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1793  }
1794  }
1795  }
1796  return kTRUE;
1797 }
1798 
1804 {
1805  TString hname;
1806 
1807  for (auto& dmeson_pair : fDmesonJets) {
1808  Int_t accJets = 0;
1809  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1810  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1811  if (!jet) continue;
1812  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1813  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1814  fHistManager->FillTH1(hname, jet->Pt());
1815  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1816  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1817  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1818  fHistManager->FillTH1(hname, jet->Eta());
1819  continue;
1820  }
1821  accJets++;
1822  }
1823  if (accJets == 0) {
1824  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1825  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1826  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1827  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1828  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1829  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1830  if (fCandidateType == kD0toKpi) {
1831  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1832  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1833  }
1834  else if (fCandidateType == kDstartoKpipi) {
1835  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1836  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1837 
1838  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1839  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1840  }
1841  }
1842  }
1843  return kTRUE;
1844 }
1845 
1850 {
1851  TString hname;
1852 
1853  for (auto& dmeson_pair : fDmesonJets) {
1854  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
1855  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1856  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1857  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1858  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1859  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1860  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1861  }
1862  }
1863 
1864  for (auto &jetDef : fJetDefinitions) {
1865 
1866  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1867  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
1868 
1869  for (auto& dmeson_pair : fDmesonJets) {
1870  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
1871  if (!jet) continue;
1872  if (!jetDef.IsJetInAcceptance(*jet)) {
1873  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
1874  fHistManager->FillTH1(hname, jet->Pt());
1875  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
1876  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1877  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
1878  fHistManager->FillTH1(hname, jet->Eta());
1879  continue;
1880  }
1881  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
1882  }
1883  }
1884 
1885  return kTRUE;
1886 }
1887 
1893 Bool_t AliAnalysisTaskDmesonJets::AnalysisEngine::FillHnSparse(THnSparse* h, const AliDmesonJetInfo& DmesonJet, std::string n)
1894 {
1895  // Fill the THnSparse histogram.
1896 
1897  Double_t contents[30] = {0.};
1898 
1899  Double_t z = DmesonJet.GetZ(n);
1900  Double_t deltaPhi = 0;
1901  Double_t deltaEta = 0;
1902  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
1903 
1904  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
1905  if (it == DmesonJet.fJets.end()) return kFALSE;
1906 
1907  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
1908  TString title(h->GetAxis(i)->GetTitle());
1909  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
1910  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
1911  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
1912  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
1913  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
1914  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
1915  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
1916  else if (title=="#it{z}_{D}") contents[i] = z;
1917  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
1918  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
1919  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
1920  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
1921  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
1922  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
1923  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
1924  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
1925  }
1926 
1927  h->Fill(contents);
1928 
1929  return kTRUE;
1930 }
1931 
1932 // Definitions of class AliAnalysisTaskDmesonJets
1933 
1937 
1941  fAnalysisEngines(),
1942  fEnabledAxis(0),
1944  fHistManager(),
1945  fApplyKinematicCuts(kTRUE),
1946  fNOutputTrees(0),
1947  fAodEvent(0),
1948  fFastJetWrapper(0)
1949 {
1950  SetMakeGeneralHistograms(kTRUE);
1951 }
1952 
1956 AliAnalysisTaskDmesonJets::AliAnalysisTaskDmesonJets(const char* name, Int_t nOutputTrees) :
1957  AliAnalysisTaskEmcalLight(name, kTRUE),
1958  fAnalysisEngines(),
1959  fEnabledAxis(k2ProngInvMass),
1960  fOutputType(kTreeOutput),
1961  fHistManager(name),
1962  fApplyKinematicCuts(kTRUE),
1963  fNOutputTrees(nOutputTrees),
1964  fAodEvent(0),
1965  fFastJetWrapper(0)
1966 {
1967  SetMakeGeneralHistograms(kTRUE);
1968  for (Int_t i = 0; i < nOutputTrees; i++){
1969  DefineOutput(2+i, TTree::Class());
1970  }
1971 }
1972 
1975 {
1976  if (fFastJetWrapper) delete fFastJetWrapper;
1977 }
1978 
1986 {
1987  AliRDHFCuts* analysiscuts = 0;
1988  TFile* filecuts = TFile::Open(cutfname);
1989  if (!filecuts || filecuts->IsZombie()) {
1990  ::Warning("AddTaskDmesonJets", "Input file not found: will use std cuts.");
1991  filecuts = 0;
1992  }
1993 
1994  if (filecuts) {
1995  analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
1996  if (!analysiscuts) {
1997  ::Warning("AddTaskDmesonJetCorr", "Could not find analysis cuts '%s' in '%s'. Using std cuts.", cutsname.Data(), cutfname.Data());
1998  }
1999  }
2000 
2001  return analysiscuts;
2002 }
2003 
2013 {
2015  return AddAnalysisEngine(type, MCmode, jetDef, cutfname);
2016 }
2017 
2027 {
2028  AliRDHFCuts* cuts = 0;
2029 
2030  if (!cutfname.IsNull()) {
2031  TString cutsname;
2032 
2033  switch (type) {
2034  case kD0toKpi :
2035  cutsname = "D0toKpiCuts";
2036  break;
2037  case kDstartoKpipi :
2038  cutsname = "DStartoKpipiCuts";
2039  break;
2040  default:
2041  return 0;
2042  }
2043 
2044  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
2045  }
2046 
2047  AnalysisEngine eng(type, MCmode, cuts);
2048 
2049  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
2050 
2051  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
2052  eng.AddJetDefinition(jetDef);
2053  it = fAnalysisEngines.insert(it, eng);
2054  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(jetDef), fAnalysisEngines.size());
2055  }
2056  else {
2057  AnalysisEngine* found_eng = &(*it);
2058  ::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());
2059  found_eng->AddJetDefinition(jetDef);
2060 
2061  if (cuts && found_eng->fRDHFCuts != 0) {
2062  ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
2063  found_eng->SetRDHFCuts(cuts);
2064  }
2065  }
2066 
2067  return &(*it);
2068 }
2069 
2070 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
2071 {
2072  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
2073  while (it != fAnalysisEngines.end() && (*it) < eng) it++;
2074  return it;
2075 }
2076 
2079 {
2080  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
2081 
2083 
2084  // Define histograms
2085  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
2086 
2087  TString hname;
2088  TString htitle;
2089  TH1* h = 0;
2090  Int_t treeSlot = 0;
2091 
2092  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
2093  for (auto &param : fAnalysisEngines) {
2094  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
2095 
2096  fHistManager.CreateHistoGroup(param.GetName());
2097 
2098  param.fHistManager = &fHistManager;
2099 
2100  hname = TString::Format("%s/fHistNAcceptedDmesons", param.GetName());
2101  htitle = hname + ";Number of D accepted meson candidates;counts";
2102  h = fHistManager.CreateTH1(hname, htitle, 51, -0.5, 50.5);
2103 
2104  hname = TString::Format("%s/fHistNDmesons", param.GetName());
2105  htitle = hname + ";Number of D meson candidates;counts";
2106  h = fHistManager.CreateTH1(hname, htitle, 101, -0.5, 100.5);
2107 
2108  hname = TString::Format("%s/fHistNEvents", param.GetName());
2109  htitle = hname + ";Event status;counts";
2110  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
2111  h->GetXaxis()->SetBinLabel(1, "Accepted");
2112  h->GetXaxis()->SetBinLabel(2, "Rejected");
2113 
2114  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
2115  htitle = hname + ";Rejection reason;counts";
2116  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
2117 
2118  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
2119  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
2120  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2121 
2122  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
2123  htitle = hname + ";#it{#eta}_{D};counts";
2124  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2125 
2126  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
2127  htitle = hname + ";#it{#phi}_{D};counts";
2128  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2129 
2130  if (param.fCandidateType == kD0toKpi) {
2131  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
2132  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2133  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
2134  }
2135  else if (param.fCandidateType == kDstartoKpipi) {
2136  Double_t min = 0;
2137  Double_t max = 0;
2138 
2139  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
2140  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2141  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
2142  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
2143 
2144  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2145  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
2146  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2147  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
2148  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
2149  }
2150 
2151  if (param.fMCMode == kBackgroundOnly || param.fMCMode == kSignalOnly || param.fMCMode == kMCTruth) {
2152  hname = TString::Format("%s/fHistDmesonOrigin", param.GetName());
2153  htitle = hname + ";origin;counts";
2154  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2155  }
2156 
2157  for (std::vector<AliHFJetDefinition>::iterator itdef = param.fJetDefinitions.begin(); itdef != param.fJetDefinitions.end(); itdef++) {
2158  AliHFJetDefinition* jetDef = &(*itdef);
2159  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef->GetName());
2160 
2161  fHistManager.CreateHistoGroup(jetDef->GetName(), param.GetName());
2162 
2163  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef->GetName());
2164  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2165  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2166  SetRejectionReasonLabels(h->GetXaxis());
2167 
2168  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef->GetName());
2169  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2170  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2171  SetRejectionReasonLabels(h->GetXaxis());
2172 
2173  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef->GetName());
2174  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
2175  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2176  SetRejectionReasonLabels(h->GetXaxis());
2177 
2178  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef->GetName());
2179  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
2180  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
2181 
2182  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef->GetName());
2183  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
2184  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2185 
2186  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef->GetName());
2187  htitle = hname + ";#it{#eta}_{jet};counts";
2188  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2189 
2190  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef->GetName());
2191  htitle = hname + ";#it{#phi}_{jet};counts";
2192  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2193  }
2194  switch (fOutputType) {
2195  case kTreeOutput:
2196  param.BuildTree(GetName());
2197  if (treeSlot < fNOutputTrees) {
2198  param.AssignDataSlot(treeSlot+2);
2199  treeSlot++;
2201  }
2202  else {
2203  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
2204  }
2205  break;
2206  case kTHnOutput:
2207  param.BuildHnSparse(fEnabledAxis);
2208  break;
2209  case kNoOutput:
2210  break;
2211  }
2212  }
2213 
2215 
2216  PostData(1, fOutput);
2217 }
2218 
2222 {
2224 
2225  // Load the event
2226  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
2227 
2228  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
2229 
2230  fFastJetWrapper->SetAreaType(fastjet::active_area);
2232 
2233  if (!fAodEvent) {
2234  AliError(Form("This task need an AOD event! Task '%s' will be disabled!", GetName()));
2235  return;
2236  }
2237 
2238  for (auto &params : fAnalysisEngines) {
2239 
2240  params.fAodEvent = fAodEvent;
2241  params.fFastJetWrapper = fFastJetWrapper;
2242  params.Init(fGeom, fAodEvent->GetRunNumber());
2243 
2244  if (params.fMCMode != kMCTruth) {
2245  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
2246 
2247  if (params.fCandidateArray) {
2248  if (!params.fCandidateArray->GetClass()->InheritsFrom("AliAODRecoDecayHF2Prong")) {
2249  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2250  "%s: Objects of type %s in %s are not inherited from AliAODRecoDecayHF2Prong! Task will be disabled!",
2251  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName());
2252  params.fCandidateArray = 0;
2253  params.fInhibit = kTRUE;
2254  }
2255  }
2256  else {
2257  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2258  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
2259  params.fBranchName.Data(), params.GetName());
2260  params.fInhibit = kTRUE;
2261  }
2262  }
2263 
2264  if (params.fMCMode != kNoMC) {
2265  params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(0));
2266 
2267  if (!params.fMCContainer) params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(1));
2268 
2269  if (!params.fMCContainer) {
2270  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2271  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
2272  params.GetName());
2273  params.fInhibit = kTRUE;
2274  }
2275  }
2276 
2277  if (params.fMCMode != kMCTruth) {
2278  params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(0));
2279  if (!params.fTrackContainer) params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(1));
2280 
2281  params.fClusterContainer = GetClusterContainer(0);
2282 
2283  if (!params.fTrackContainer && !params.fClusterContainer) {
2284  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2285  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
2286  params.GetName());
2287  params.fInhibit = kTRUE;
2288  }
2289  }
2290  }
2291 }
2292 
2297 {
2298  if (!fAodEvent) return kFALSE;
2299 
2300  TString hname;
2301 
2302  // fix for temporary bug in ESDfilter
2303  // the AODs with null vertex pointer didn't pass the PhysSel
2304  if (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001) return kFALSE;
2305 
2306  for (auto &eng : fAnalysisEngines) {
2307 
2308  if (eng.fInhibit) continue;
2309 
2310  //Event selection
2311  hname = TString::Format("%s/fHistNEvents", eng.GetName());
2312  Bool_t iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
2313  if (!iseventselected) {
2314  fHistManager.FillTH1(hname, "Rejected");
2315  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2316  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
2317  TString label;
2318  do {
2319  label = GetHFEventRejectionReasonLabel(bitmap);
2320  if (label.IsNull()) break;
2321  fHistManager.FillTH1(hname, label);
2322  } while (true);
2323 
2324  continue;
2325  }
2326 
2327  fHistManager.FillTH1(hname, "Accepted");
2328 
2329  AliDebug(2, "Event selected");
2330 
2331  eng.RunAnalysis();
2332  }
2333  return kTRUE;
2334 }
2335 
2340 {
2341  TString hname;
2342  for (auto &param : fAnalysisEngines) {
2343 
2344  if (param.fInhibit) continue;
2345 
2346  if (fOutputType == kTreeOutput) {
2347  param.FillTree(fApplyKinematicCuts);
2349  }
2350  else if (fOutputType == kTHnOutput) {
2351  param.FillHnSparse(fApplyKinematicCuts);
2352  }
2353  }
2354  return kTRUE;
2355 }
2356 
2364 void AliAnalysisTaskDmesonJets::CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t& minMass, Double_t& maxMass)
2365 {
2366  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
2367 
2368  Double_t mass = part->Mass();
2369 
2370  // To make sure that the PDG mass value is not at the edge of a bin
2371  if (nbins % 2 == 0) {
2372  minMass = mass - range / 2 - range / nbins / 2;
2373  maxMass = mass + range / 2 - range / nbins / 2;
2374  }
2375  else {
2376  minMass = mass - range / 2;
2377  maxMass = mass + range / 2;
2378  }
2379 }
2380 
2387 {
2388  static TString label;
2389  label = "";
2390 
2391  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
2392  label = "NotSelTrigger";
2393  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
2394  return label.Data();
2395  }
2396  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
2397  label = "NoVertex";
2398  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
2399  return label.Data();
2400  }
2401  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
2402  label = "TooFewVtxContrib";
2403  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
2404  return label.Data();
2405  }
2406  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
2407  label = "ZVtxOutFid";
2408  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
2409  return label.Data();
2410  }
2411  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
2412  label = "Pileup";
2413  bitmap &= ~BIT(AliRDHFCuts::kPileup);
2414  return label.Data();
2415  }
2416  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
2417  label = "OutsideCentrality";
2418  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
2419  return label.Data();
2420  }
2421  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
2422  label = "PhysicsSelection";
2423  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
2424  return label.Data();
2425  }
2426  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
2427  label = "BadSPDVertex";
2428  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
2429  return label.Data();
2430  }
2431  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
2432  label = "ZVtxSPDOutFid";
2433  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
2434  return label.Data();
2435  }
2436  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
2437  label = "CentralityFlattening";
2438  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
2439  return label.Data();
2440  }
2441  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
2442  label = "BadTrackV0Correl";
2443  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
2444  return label.Data();
2445  }
2446 
2447  return label.Data();
2448 }
2449 
2456 {
2457  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
2458  PostData(eng.GetDataSlotNumber(), eng.GetTree());
2459  return eng.GetDataSlotNumber();
2460  }
2461  else {
2462  return -1;
2463  }
2464 }
Int_t pdg
void Print() const
Prints the content of this object in the standard output.
std::list< AnalysisEngine > fAnalysisEngines
Array of analysis parameters.
AliEMCALGeometry * fGeom
!emcal geometry
Analysis task for D meson jets.
virtual void Set(const AliDmesonJetInfo &source)
Invariant mass of the D0 meson candidate in GeV/c2.
virtual Int_t Run()
Bool_t IsJetInAcceptance(const AliJetInfo &jet) const
Double_t fSoftPionPt
! Transverse momentum of the soft pion of the D* candidate
const char * title
Definition: MakeQAPdf.C:26
AliRDHFCuts * fRDHFCuts
D meson candidates cuts.
Lightweight class that encapsulates D meson jets.
DCal acceptance.
Definition: AliEmcalJet.h:59
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.
Full acceptance, i.e. no acceptance cut applied – left to user.
Definition: AliEmcalJet.h:61
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="/")
TPC fiducial acceptance (each eta edge narrowed by jet R)
Definition: AliEmcalJet.h:56
virtual void Set(const AliDmesonJetInfo &source)
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:33
EJetType_t fJetType
Jet type (charged, full, neutral)
std::pair< AliJetInfo *, Double_t > jet_distance_pair
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
virtual void Set(const AliDmesonJetInfo &source)
AliTLorentzVector fD
! 4-momentum of the D meson candidate
Double_t GetDistance(const AliJetInfo &jet, Double_t &deta, Double_t &dphi) const
TObject * FindObject(const char *name) const
TPC acceptance.
Definition: AliEmcalJet.h:55
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
EMCal acceptance.
Definition: AliEmcalJet.h:57
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:95
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
AliClusterContainer * GetClusterContainer(Int_t i=0) const
AliAODTrack * GetBachelor() const
AliRDHFCuts * LoadDMesonCutsFromFile(TString cutfname, TString cutsname)
bool operator<(const AliAnalysisTaskDmesonJets::AliHFJetDefinition &lhs, const AliAnalysisTaskDmesonJets::AliHFJetDefinition &rhs)
std::list< AnalysisEngine >::iterator FindAnalysisEngine(const AnalysisEngine &eng)
std::vector< AliHFJetDefinition > fJetDefinitions
Jet definitions.
EJetAlgo_t fJetAlgo
Jet algorithm (kt, anti-kt,...)
static const char * GetHFEventRejectionReasonLabel(UInt_t &bitmap)
Lightweight class that encapsulates D*.
virtual void Clear(const Option_t *="")
Bool_t ExtractDstarAttributes(const AliAODRecoCascadeHF *DstarCand, AliDmesonJetInfo &DmesonJet, UInt_t i)
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
const AliJetInfo * GetJet(std::string n) const
virtual ~AliAnalysisTaskDmesonJets()
This is the standard destructor.
Bool_t ExtractD0Attributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
std::vector< AliJetInfo > fJets
! Inclusive jets reconstructed in the current event (includes D meson candidate daughters, if any)
Select MC particles based on specific prescriptions of HF analysis.
static EMesonDecayChannel_t CheckDecayChannel(const AliAODMCParticle *part, TClonesArray *mcArray)
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:99
Double_t fMaxChargedPt
Transverse momentum of the leading charged particle (or track)
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
static FJJetAlgo ConvertToFJAlgo(EJetAlgo_t algo)
decay
Definition: HFPtSpectrum.C:41
std::map< std::string, AliJetInfo > fJets
! list of jets
void AddInputVectors(AliEmcalContainer *cont, Int_t offset, TH2 *rejectHist=0)
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Double_t minMass
void RunParticleLevelAnalysis()
Run a particle level analysis.
Bool_t fApplyKinematicCuts
Apply jet kinematic cuts.
const char * GetName() const
Generate a name for this jet definition.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
void SetDetectorJetEtaPhiRange(const AliEMCALGeometry *const geom, Int_t run)
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="")
DCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:60
const Int_t nbins
Double_t maxMass
Int_t fNOutputTrees
Maximum number of output trees.
AliHFJetDefinition()
This is the default constructor, used for ROOT I/O purposes.
AliHFJetDefinition & operator=(const AliHFJetDefinition &source)
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:97
AliTLorentzVector fMomentum
4-momentum of the jet
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
AliHFJetDefinition * AddJetDefinition(EJetType_t type, Double_t r, EJetAlgo_t algo, ERecoScheme_t reco)
Lightweight class that encapsulates D0.
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:58
const Int_t nPtBins
Bool_t IsAnyJetInAcceptance(const AliDmesonJetInfo &dMesonJet) const
Container for jet within the EMCAL jet framework.
Int_t PostDataFromAnalysisEngine(const AnalysisEngine &eng)
static EMesonOrigin_t CheckOrigin(const AliAODMCParticle *part, TClonesArray *mcArray)
AliFJWrapper * fFastJetWrapper
! Fastjet wrapper
THistManager fHistManager
Histogram manager.
void Reset()
Reset all fields to their default values.
virtual void Set(const AliDmesonJetInfo &source, std::string n)