AliPhysics  9fe175b (9fe175b)
 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),
673  fRejectedOrigin(kUnknownQuark | kFromBottom),
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  if ((origin & fRejectedOrigin) == origin) return kFALSE;
1000  }
1001  MCtruthPdgCode = aodMcPart->PdgCode();
1002  }
1003  }
1004  }
1005 
1006  //AliDebug(2,"Checking if D0 meson is selected");
1007  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
1008  if (isSelected == 1) { // selected as a D0
1009  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
1010 
1011  if (fMCMode == kNoMC ||
1012  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1013  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly)) {
1014  // 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)
1015  //AliDebug(2,"Selected as D0");
1016  invMassD = Dcand->InvMassD0();
1017  }
1018  else { // conditions above not passed, so return FALSE
1019  return kFALSE;
1020  }
1021  }
1022  else if (isSelected == 2) { // selected as a D0bar
1023  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
1024 
1025  if (fMCMode == kNoMC ||
1026  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1027  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly)) {
1028  // 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)
1029  //AliDebug(2,"Selected as D0bar");
1030  invMassD = Dcand->InvMassD0bar();
1031  }
1032  else { // conditions above not passed, so return FALSE
1033  return kFALSE;
1034  }
1035  }
1036  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
1037  //AliDebug(2,"Selected as either D0 or D0bar");
1038 
1039  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
1040  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1041  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kBackgroundOnly)) {
1042  if (i > 0) return kFALSE;
1043  //AliDebug(2, "MC truth is D0");
1044  invMassD = Dcand->InvMassD0();
1045  }
1046  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1047  (MCtruthPdgCode == fCandidatePDG && fMCMode == kBackgroundOnly)) {
1048  if (i > 0) return kFALSE;
1049  //AliDebug(2, "MC truth is D0bar");
1050  invMassD = Dcand->InvMassD0bar();
1051  }
1052  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
1053 
1054  // Only accept it if background-only OR background-and-signal was requested
1055  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
1056  // Select D0 or D0bar depending on the i-parameter
1057  if (i == 0) {
1058  //AliDebug(2, "Returning invariant mass with D0 hypothesis");
1059  invMassD = Dcand->InvMassD0();
1060  }
1061  else if (i == 1) {
1062  //AliDebug(2, "Returning invariant mass with D0bar hypothesis");
1063  invMassD = Dcand->InvMassD0bar();
1064  }
1065  else { // i > 1
1066  return kFALSE;
1067  }
1068  }
1069  else { // signal-only was requested but this is not a true D0
1070  return kFALSE;
1071  }
1072  }
1073  }
1074  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
1075  return kTRUE;
1076 }
1077 
1086 {
1087  if (i > 0) return kFALSE; // only one mass hypothesis for the D*
1088 
1089  Int_t MCtruthPdgCode = 0;
1090 
1091  Double_t invMassD = 0;
1092 
1093  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1094  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1095  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1096 
1097  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1098  //AliDebug(2, Form("MC label is %d", mcLab));
1099  DmesonJet.fMCLabel = mcLab;
1100  if (mcLab >= 0) {
1101  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1102 
1103  if (aodMcPart) {
1104  if (fRejectedOrigin) {
1105  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
1106 
1107  if ((origin & fRejectedOrigin) == origin) return kFALSE;
1108  }
1109 
1110  MCtruthPdgCode = aodMcPart->PdgCode();
1111  //AliDebug(2, Form("MC truth pdg code is %d",MCtruthPdgCode));
1112  }
1113  }
1114  }
1115 
1116  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1117  if (fMCMode == kNoMC ||
1118  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1119  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1120  // 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)
1121  invMassD = DstarCand->InvMassDstarKpipi();
1122  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1123  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1124  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1125  return kTRUE;
1126  }
1127  else { // conditions above not passed, so return FALSE
1128  return kFALSE;
1129  }
1130 }
1131 
1139 {
1140  if (!part) return kDecayOther;
1141  if (!mcArray) return kDecayOther;
1142 
1144 
1145  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1146 
1147  if (part->GetNDaughters() == 2) {
1148 
1149  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1150  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1151 
1152  if (!d1 || !d2) {
1153  return decay;
1154  }
1155 
1156  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1157  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1158 
1159  if (absPdgPart == 421) { // D0 -> K pi
1160 
1161  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1162  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1163  decay = kDecayD0toKpi;
1164  }
1165  }
1166 
1167  if (absPdgPart == 413) { // D* -> D0 pi
1168 
1169  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1170  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1171  if (D0decay == kDecayD0toKpi) {
1172  decay = kDecayDStartoKpipi;
1173  }
1174  }
1175 
1176  if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1177  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1178  if (D0decay == kDecayD0toKpi) {
1179  decay = kDecayDStartoKpipi;
1180  }
1181  }
1182  }
1183  }
1184 
1185  return decay;
1186 }
1187 
1195 {
1196  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1197 
1198  if (!part) return kUnknownQuark;
1199  if (!mcArray) return kUnknownQuark;
1200 
1201  Int_t pdgGranma = 0;
1202  Int_t mother = part->GetMother();
1203  Int_t istep = 0;
1204  Int_t abspdgGranma = 0;
1205  Bool_t isFromB = kFALSE;
1206  Bool_t isQuarkFound = kFALSE;
1207 
1208  while (mother >= 0) {
1209  istep++;
1210  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1211  if (mcGranma >= 0) {
1212  pdgGranma = mcGranma->GetPdgCode();
1213  abspdgGranma = TMath::Abs(pdgGranma);
1214  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1215  isFromB = kTRUE;
1216  }
1217 
1218  if (abspdgGranma == 4 || abspdgGranma == 5) isQuarkFound = kTRUE;
1219  mother = mcGranma->GetMother();
1220  }
1221  else {
1222  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1223  break;
1224  }
1225  }
1226 
1227  if (isQuarkFound) {
1228  if (isFromB) {
1229  return kFromBottom;
1230  }
1231  else {
1232  return kFromCharm;
1233  }
1234  }
1235  else {
1236  return kUnknownQuark;
1237  }
1238 }
1239 
1242 {
1243  fDmesonJets.clear();
1244 
1245  for (auto& jetDef : fJetDefinitions) {
1246  jetDef.fJets.clear();
1247  }
1248 
1249  if (fMCMode == kMCTruth) {
1250  RunParticleLevelAnalysis();
1251  }
1252  else {
1253  RunDetectorLevelAnalysis();
1254  }
1255 }
1256 
1262 {
1264  fFastJetWrapper->SetR(jetDef.fRadius);
1267 
1268  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1269  fTrackContainer->SetDMesonCandidate(0);
1270  AddInputVectors(fTrackContainer, 100);
1271  }
1272 
1273  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1274  AddInputVectors(fClusterContainer, -100);
1275  }
1276 
1277  // run jet finder
1278  fFastJetWrapper->Run();
1279 
1280  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1281 
1282  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1283  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1284 
1285  Double_t maxChPt = 0;
1286  Double_t maxNePt = 0;
1287  Double_t totalNeutralPt = 0;
1288 
1289  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1290  if (constituents[ic].user_index() >= 100) {
1291  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1292  }
1293  else if (constituents[ic].user_index() <= -100) {
1294  totalNeutralPt += constituents[ic].pt();
1295  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1296  }
1297  }
1298 
1299  jetDef.fJets.push_back(
1300  AliJetInfo(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E(),
1301  constituents.size(), totalNeutralPt / jets_incl[ijet].pt(), maxChPt, maxNePt)
1302  );
1303  }
1304 }
1305 
1314 {
1315  if (jetDef.fJets.size() == 0) FindJets(jetDef);
1316 
1317  Double_t d_closest = 999;
1318  AliJetInfo* jet_closest = 0;
1319  const AliJetInfo* truth_jet = 0;
1320  try {
1321  truth_jet = &(dmeson.fJets.at(jetDef.GetName()));
1322  }
1323  catch(...) {
1324  return jet_distance_pair(0, 999);
1325  }
1326 
1327  for (auto& jet : jetDef.fJets) {
1328  Double_t d = truth_jet->GetDistance(jet);
1329 
1330  if (d > dMax) continue;
1331  if (d < d_closest) {
1332  d_closest = d;
1333  jet_closest = &jet;
1334  }
1335  }
1336 
1337  if (jet_closest && applyKinCuts) {
1338  if (!jetDef.IsJetInAcceptance(*jet_closest)) jet_closest = 0;
1339  }
1340 
1341  if (jet_closest) {
1342  AliDebug(2, Form("Found closest jet (pt=%.3f, eta=%.3f, phi=%.3f) to d meson (pt=%.3f, eta=%.3f, phi=%.3f) with d = %.3f",
1343  jet_closest->Pt(), jet_closest->Eta(), jet_closest->Phi_0_2pi(),
1344  dmeson.fD.Pt(), dmeson.fD.Eta(), dmeson.fD.Phi_0_2pi(),
1345  d_closest));
1346  }
1347 
1348  return jet_distance_pair(jet_closest, d_closest);
1349 }
1350 
1353 {
1354  const Int_t nD = fCandidateArray->GetEntriesFast();
1355 
1356  AliDmesonJetInfo DmesonJet;
1357 
1358  Int_t nAccCharm = 0;
1359  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1360  Int_t isSelected = 0;
1361 
1362  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1363  if (!charmCand) continue;
1364 
1365  Int_t nprongs = charmCand->GetNProngs();
1366 
1367  if (fCandidateType == kDstartoKpipi) {
1368  if (!charmCand->InheritsFrom("AliAODRecoCascadeHF")) {
1369  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::RunDetectorLevelAnalysis","Candidate type is D* but object type is wrong (should be AliAODRecoCascadeHF)");
1370  continue;
1371  }
1372  }
1373 
1374  // region of interest + cuts
1375  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1376 
1377  //candidate selected by cuts and PID
1378  isSelected = fRDHFCuts->IsSelected(charmCand, AliRDHFCuts::kAll, fAodEvent); //selected
1379 
1380  if (!isSelected) continue;
1381 
1382  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1383  DmesonJet.Reset();
1384  DmesonJet.fDmesonParticle = charmCand;
1385  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1386  for (std::vector<AliHFJetDefinition>::iterator itdef = fJetDefinitions.begin(); itdef != fJetDefinitions.end(); itdef++) {
1387  if (!FindJet(charmCand, DmesonJet, *itdef)) {
1388  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1389  (*itdef).GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1390  }
1391  }
1392  fDmesonJets[icharm] = DmesonJet;
1393  }
1394  }
1395  nAccCharm++;
1396  } // end of D cand loop
1397 
1398  TString hname;
1399 
1400  hname = TString::Format("%s/fHistNAcceptedDmesons", GetName());
1401  fHistManager->FillTH1(hname, nAccCharm);
1402 
1403  hname = TString::Format("%s/fHistNDmesons", GetName());
1404  fHistManager->FillTH1(hname, nD);
1405 }
1406 
1417 {
1418  TString hname;
1419 
1421  fFastJetWrapper->SetR(jetDef.fRadius);
1424 
1425  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1426 
1427  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1428  fTrackContainer->SetDMesonCandidate(Dcand);
1429  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1430  AddInputVectors(fTrackContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1431 
1432  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1433  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1434  const TObjArray& daughters = fTrackContainer->GetDaughterList();
1435  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1436  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1437  if (!fTrackContainer->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1438  }
1439  }
1440 
1441  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1442  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1443  AddInputVectors(fClusterContainer, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1444  }
1445 
1446  // run jet finder
1447  fFastJetWrapper->Run();
1448 
1449  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1450 
1451  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1452  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1453 
1454  Bool_t isDmesonJet = kFALSE;
1455 
1456  Double_t maxChPt = 0;
1457  Double_t maxNePt = 0;
1458  Double_t totalNeutralPt = 0;
1459 
1460  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1461  if (constituents[ic].user_index() == 0) {
1462  isDmesonJet = kTRUE;
1463  }
1464  else if (constituents[ic].user_index() >= 100) {
1465  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1466  }
1467  else if (constituents[ic].user_index() <= -100) {
1468  totalNeutralPt += constituents[ic].pt();
1469  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1470  }
1471  }
1472 
1473  if (isDmesonJet) {
1474  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1475  DmesonJet.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1476  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1477  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1478  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1479 
1480  return kTRUE;
1481  }
1482  }
1483 
1484  return kFALSE;
1485 }
1486 
1490 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist)
1491 {
1492  AliEmcalIterableMomentumContainer itcont = cont->all_momentum();
1493  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1494  UInt_t rejectionReason = 0;
1495  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1496  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1497  continue;
1498  }
1499  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1500  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1501  }
1502 }
1503 
1506 {
1507  TString hname;
1508 
1509  fMCContainer->SetSpecialPDG(fCandidatePDG);
1510  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1511  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1512 
1513  if (!fMCContainer->IsSpecialPDGFound()) return;
1514 
1515  for (auto &jetDef : fJetDefinitions) {
1516 
1518  fFastJetWrapper->SetR(jetDef.fRadius);
1521 
1522  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1523  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1524 
1525  fFastJetWrapper->Run();
1526 
1527  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1528 
1529  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1530  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1531 
1532  Bool_t isDmesonJet = kFALSE;
1533 
1534  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1535  Int_t iPart = constituents[ic].user_index() - 100;
1536  AliVParticle* part = fMCContainer->GetParticle(iPart);
1537  if (!part) {
1538  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1539  continue;
1540  }
1541  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1542  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1543  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1544  std::pair<int, AliDmesonJetInfo> element;
1545  element.first = iPart;
1546  element.second = AliDmesonJetInfo();
1547  dMesonJetIt = fDmesonJets.insert(element).first;
1548  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1549  (*dMesonJetIt).second.fDmesonParticle = part;
1550  }
1551 
1552  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1553  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1554  }
1555  }
1556  }
1557  }
1558 }
1559 
1564 {
1565  TString classname;
1566  switch (fCandidateType) {
1567  case kD0toKpi:
1568  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1569  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1570  break;
1571  case kDstartoKpipi:
1572  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1573  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1574  break;
1575  }
1576  TString treeName = TString::Format("%s_%s", taskName, GetName());
1577  fTree = new TTree(treeName, treeName);
1578  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
1579  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1580  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1581  fCurrentJetInfo[i] = new AliJetInfoSummary();
1582  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1583  }
1584 
1585  return fTree;
1586 }
1587 
1592 {
1593  TString hname;
1594 
1595  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
1596 
1597  for (auto &jetDef : fJetDefinitions) {
1598 
1599  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
1600 
1601  Double_t radius = jetDef.fRadius;
1602 
1603  TString title[30] = {""};
1604  Int_t nbins[30] = {0 };
1605  Double_t min [30] = {0.};
1606  Double_t max [30] = {0.};
1607  Int_t dim = 0 ;
1608 
1609  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1610  nbins[dim] = nPtBins;
1611  min[dim] = 0;
1612  max[dim] = fMaxPt;
1613  dim++;
1614 
1615  if ((enabledAxis & kPositionD) != 0) {
1616  title[dim] = "#eta_{D}";
1617  nbins[dim] = 50;
1618  min[dim] = -1;
1619  max[dim] = 1;
1620  dim++;
1621 
1622  title[dim] = "#phi_{D} (rad)";
1623  nbins[dim] = 150;
1624  min[dim] = 0;
1625  max[dim] = TMath::TwoPi();
1626  dim++;
1627  }
1628 
1629  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1630  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1631  nbins[dim] = fNMassBins;
1632  min[dim] = fMinMass;
1633  max[dim] = fMaxMass;
1634  dim++;
1635  }
1636 
1637  if (fCandidateType == kD0toKpi) {
1638  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1639  nbins[dim] = fNMassBins;
1640  min[dim] = fMinMass;
1641  max[dim] = fMaxMass;
1642  dim++;
1643  }
1644 
1645  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1646  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1647  nbins[dim] = fNMassBins;
1648  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1649  dim++;
1650  }
1651 
1652  if (fCandidateType == kDstartoKpipi) {
1653  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1654  nbins[dim] = fNMassBins*6;
1655  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1656 
1657  // subtract mass of D0
1658  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1659  min[dim] -= D0mass;
1660  max[dim] -= D0mass;
1661 
1662  dim++;
1663  }
1664 
1665  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1666  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1667  nbins[dim] = 100;
1668  min[dim] = 0;
1669  max[dim] = 25;
1670  dim++;
1671  }
1672 
1673  title[dim] = "#it{z}_{D}";
1674  nbins[dim] = 110;
1675  min[dim] = 0;
1676  max[dim] = 1.10;
1677  dim++;
1678 
1679  if ((enabledAxis & kDeltaR) != 0) {
1680  title[dim] = "#Delta R_{D-jet}";
1681  nbins[dim] = 100;
1682  min[dim] = 0;
1683  max[dim] = radius * 1.5;
1684  dim++;
1685  }
1686 
1687  if ((enabledAxis & kDeltaEta) != 0) {
1688  title[dim] = "#eta_{D} - #eta_{jet}";
1689  nbins[dim] = 100;
1690  min[dim] = -radius * 1.2;
1691  max[dim] = radius * 1.2;
1692  dim++;
1693  }
1694 
1695  if ((enabledAxis & kDeltaPhi) != 0) {
1696  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1697  nbins[dim] = 100;
1698  min[dim] = -radius * 1.2;
1699  max[dim] = radius * 1.2;
1700  dim++;
1701  }
1702 
1703  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1704  nbins[dim] = nPtBins;
1705  min[dim] = 0;
1706  max[dim] = fMaxPt;
1707  dim++;
1708 
1709  if ((enabledAxis & kPositionJet) != 0) {
1710  title[dim] = "#eta_{jet}";
1711  nbins[dim] = 50;
1712  min[dim] = -1;
1713  max[dim] = 1;
1714  dim++;
1715 
1716  title[dim] = "#phi_{jet} (rad)";
1717  nbins[dim] = 150;
1718  min[dim] = 0;
1719  max[dim] = TMath::TwoPi();
1720  dim++;
1721  }
1722 
1723  if ((enabledAxis & kJetConstituents) != 0) {
1724  title[dim] = "No. of constituents";
1725  nbins[dim] = 50;
1726  min[dim] = -0.5;
1727  max[dim] = 49.5;
1728  dim++;
1729  }
1730 
1731  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1732  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1733  for (Int_t j = 0; j < dim; j++) {
1734  h->GetAxis(j)->SetTitle(title[j]);
1735  }
1736  }
1737 }
1738 
1743 {
1744  TString hname;
1745 
1746  for (auto& dmeson_pair : fDmesonJets) {
1747  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
1748  Int_t accJets = 0;
1749  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1750  fCurrentJetInfo[ij]->Reset();
1751  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1752  if (!jet) continue;
1753  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1754  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1755  fHistManager->FillTH1(hname, jet->Pt());
1756  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1757  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1758  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1759  fHistManager->FillTH1(hname, jet->Eta());
1760  continue;
1761  }
1762  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
1763  accJets++;
1764  }
1765  if (kTRUE || accJets > 0) { // always fill D meson tree, even if no jet was accepted
1766  fTree->Fill();
1767  }
1768  else { // this is never executed
1769  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1770  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1771  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1772  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1773  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1774  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1775  if (fCandidateType == kD0toKpi) {
1776  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1777  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1778  }
1779  else if (fCandidateType == kDstartoKpipi) {
1780  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1781  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1782 
1783  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1784  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1785  }
1786  }
1787  }
1788  return kTRUE;
1789 }
1790 
1796 {
1797  TString hname;
1798 
1799  for (auto& dmeson_pair : fDmesonJets) {
1800  Int_t accJets = 0;
1801  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1802  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1803  if (!jet) continue;
1804  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1805  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1806  fHistManager->FillTH1(hname, jet->Pt());
1807  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1808  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1809  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1810  fHistManager->FillTH1(hname, jet->Eta());
1811  continue;
1812  }
1813  accJets++;
1814  }
1815  if (accJets == 0) {
1816  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1817  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1818  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1819  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1820  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1821  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1822  if (fCandidateType == kD0toKpi) {
1823  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1824  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1825  }
1826  else if (fCandidateType == kDstartoKpipi) {
1827  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1828  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1829 
1830  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1831  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1832  }
1833  }
1834  }
1835  return kTRUE;
1836 }
1837 
1842 {
1843  TString hname;
1844 
1845  for (auto& dmeson_pair : fDmesonJets) {
1846  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
1847  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1848  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1849  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1850  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1851  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1852  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1853  }
1854  }
1855 
1856  for (auto &jetDef : fJetDefinitions) {
1857 
1858  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1859  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
1860 
1861  for (auto& dmeson_pair : fDmesonJets) {
1862  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
1863  if (!jet) continue;
1864  if (!jetDef.IsJetInAcceptance(*jet)) {
1865  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
1866  fHistManager->FillTH1(hname, jet->Pt());
1867  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
1868  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1869  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
1870  fHistManager->FillTH1(hname, jet->Eta());
1871  continue;
1872  }
1873  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
1874  }
1875  }
1876 
1877  return kTRUE;
1878 }
1879 
1885 Bool_t AliAnalysisTaskDmesonJets::AnalysisEngine::FillHnSparse(THnSparse* h, const AliDmesonJetInfo& DmesonJet, std::string n)
1886 {
1887  // Fill the THnSparse histogram.
1888 
1889  Double_t contents[30] = {0.};
1890 
1891  Double_t z = DmesonJet.GetZ(n);
1892  Double_t deltaPhi = 0;
1893  Double_t deltaEta = 0;
1894  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
1895 
1896  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
1897  if (it == DmesonJet.fJets.end()) return kFALSE;
1898 
1899  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
1900  TString title(h->GetAxis(i)->GetTitle());
1901  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
1902  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
1903  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
1904  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
1905  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
1906  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
1907  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
1908  else if (title=="#it{z}_{D}") contents[i] = z;
1909  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
1910  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
1911  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
1912  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
1913  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
1914  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
1915  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
1916  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
1917  }
1918 
1919  h->Fill(contents);
1920 
1921  return kTRUE;
1922 }
1923 
1924 // Definitions of class AliAnalysisTaskDmesonJets
1925 
1929 
1933  fAnalysisEngines(),
1934  fEnabledAxis(0),
1936  fHistManager(),
1937  fApplyKinematicCuts(kTRUE),
1938  fNOutputTrees(0),
1939  fAodEvent(0),
1940  fFastJetWrapper(0)
1941 {
1942  SetMakeGeneralHistograms(kTRUE);
1943 }
1944 
1948 AliAnalysisTaskDmesonJets::AliAnalysisTaskDmesonJets(const char* name, Int_t nOutputTrees) :
1949  AliAnalysisTaskEmcalLight(name, kTRUE),
1950  fAnalysisEngines(),
1951  fEnabledAxis(k2ProngInvMass),
1952  fOutputType(kTreeOutput),
1953  fHistManager(name),
1954  fApplyKinematicCuts(kTRUE),
1955  fNOutputTrees(nOutputTrees),
1956  fAodEvent(0),
1957  fFastJetWrapper(0)
1958 {
1959  SetMakeGeneralHistograms(kTRUE);
1960  for (Int_t i = 0; i < nOutputTrees; i++){
1961  DefineOutput(2+i, TTree::Class());
1962  }
1963 }
1964 
1967 {
1968  if (fFastJetWrapper) delete fFastJetWrapper;
1969 }
1970 
1978 {
1979  AliRDHFCuts* analysiscuts = 0;
1980  TFile* filecuts = TFile::Open(cutfname);
1981  if (!filecuts || filecuts->IsZombie()) {
1982  ::Warning("AddTaskDmesonJets", "Input file not found: will use std cuts.");
1983  filecuts = 0;
1984  }
1985 
1986  if (filecuts) {
1987  analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
1988  if (!analysiscuts) {
1989  ::Warning("AddTaskDmesonJetCorr", "Could not find analysis cuts '%s' in '%s'. Using std cuts.", cutsname.Data(), cutfname.Data());
1990  }
1991  }
1992 
1993  return analysiscuts;
1994 }
1995 
2005 {
2007  return AddAnalysisEngine(type, MCmode, jetDef, cutfname);
2008 }
2009 
2019 {
2020  AliRDHFCuts* cuts = 0;
2021 
2022  if (!cutfname.IsNull()) {
2023  TString cutsname;
2024 
2025  switch (type) {
2026  case kD0toKpi :
2027  cutsname = "D0toKpiCuts";
2028  break;
2029  case kDstartoKpipi :
2030  cutsname = "DStartoKpipiCuts";
2031  break;
2032  default:
2033  return 0;
2034  }
2035 
2036  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
2037  }
2038 
2039  AnalysisEngine eng(type, MCmode, cuts);
2040 
2041  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
2042 
2043  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
2044  eng.AddJetDefinition(jetDef);
2045  it = fAnalysisEngines.insert(it, eng);
2046  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(jetDef), fAnalysisEngines.size());
2047  }
2048  else {
2049  AnalysisEngine* found_eng = &(*it);
2050  ::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());
2051  found_eng->AddJetDefinition(jetDef);
2052 
2053  if (cuts && found_eng->fRDHFCuts != 0) {
2054  ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
2055  found_eng->SetRDHFCuts(cuts);
2056  }
2057  }
2058 
2059  return &(*it);
2060 }
2061 
2062 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
2063 {
2064  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
2065  while (it != fAnalysisEngines.end() && (*it) < eng) it++;
2066  return it;
2067 }
2068 
2071 {
2072  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
2073 
2075 
2076  // Define histograms
2077  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
2078 
2079  TString hname;
2080  TString htitle;
2081  TH1* h = 0;
2082  Int_t treeSlot = 0;
2083 
2084  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
2085  for (auto &param : fAnalysisEngines) {
2086  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
2087 
2088  fHistManager.CreateHistoGroup(param.GetName());
2089 
2090  param.fHistManager = &fHistManager;
2091 
2092  hname = TString::Format("%s/fHistNAcceptedDmesons", param.GetName());
2093  htitle = hname + ";Number of D accepted meson candidates;counts";
2094  h = fHistManager.CreateTH1(hname, htitle, 51, -0.5, 50.5);
2095 
2096  hname = TString::Format("%s/fHistNDmesons", param.GetName());
2097  htitle = hname + ";Number of D meson candidates;counts";
2098  h = fHistManager.CreateTH1(hname, htitle, 101, -0.5, 100.5);
2099 
2100  hname = TString::Format("%s/fHistNEvents", param.GetName());
2101  htitle = hname + ";Event status;counts";
2102  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
2103  h->GetXaxis()->SetBinLabel(1, "Accepted");
2104  h->GetXaxis()->SetBinLabel(2, "Rejected");
2105 
2106  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
2107  htitle = hname + ";Rejection reason;counts";
2108  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
2109 
2110  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
2111  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
2112  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2113 
2114  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
2115  htitle = hname + ";#it{#eta}_{D};counts";
2116  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2117 
2118  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
2119  htitle = hname + ";#it{#phi}_{D};counts";
2120  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2121 
2122  if (param.fCandidateType == kD0toKpi) {
2123  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
2124  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2125  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
2126  }
2127  else if (param.fCandidateType == kDstartoKpipi) {
2128  Double_t min = 0;
2129  Double_t max = 0;
2130 
2131  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
2132  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2133  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
2134  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
2135 
2136  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2137  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
2138  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2139  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
2140  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
2141  }
2142 
2143  for (std::vector<AliHFJetDefinition>::iterator itdef = param.fJetDefinitions.begin(); itdef != param.fJetDefinitions.end(); itdef++) {
2144  AliHFJetDefinition* jetDef = &(*itdef);
2145  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef->GetName());
2146 
2147  fHistManager.CreateHistoGroup(jetDef->GetName(), param.GetName());
2148 
2149  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef->GetName());
2150  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2151  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2152  SetRejectionReasonLabels(h->GetXaxis());
2153 
2154  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef->GetName());
2155  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2156  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2157  SetRejectionReasonLabels(h->GetXaxis());
2158 
2159  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef->GetName());
2160  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
2161  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2162  SetRejectionReasonLabels(h->GetXaxis());
2163 
2164  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef->GetName());
2165  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
2166  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
2167 
2168  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef->GetName());
2169  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
2170  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2171 
2172  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef->GetName());
2173  htitle = hname + ";#it{#eta}_{jet};counts";
2174  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2175 
2176  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef->GetName());
2177  htitle = hname + ";#it{#phi}_{jet};counts";
2178  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2179  }
2180  switch (fOutputType) {
2181  case kTreeOutput:
2182  param.BuildTree(GetName());
2183  if (treeSlot < fNOutputTrees) {
2184  param.AssignDataSlot(treeSlot+2);
2185  treeSlot++;
2187  }
2188  else {
2189  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
2190  }
2191  break;
2192  case kTHnOutput:
2193  param.BuildHnSparse(fEnabledAxis);
2194  break;
2195  case kNoOutput:
2196  break;
2197  }
2198  }
2199 
2201 
2202  PostData(1, fOutput);
2203 }
2204 
2208 {
2210 
2211  // Load the event
2212  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
2213 
2214  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
2215 
2216  fFastJetWrapper->SetAreaType(fastjet::active_area);
2218 
2219  if (!fAodEvent) {
2220  AliError(Form("This task need an AOD event! Task '%s' will be disabled!", GetName()));
2221  return;
2222  }
2223 
2224  for (auto &params : fAnalysisEngines) {
2225 
2226  params.fAodEvent = fAodEvent;
2227  params.fFastJetWrapper = fFastJetWrapper;
2228  params.Init(fGeom, fAodEvent->GetRunNumber());
2229 
2230  if (params.fMCMode != kMCTruth) {
2231  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
2232 
2233  if (params.fCandidateArray) {
2234  if (!params.fCandidateArray->GetClass()->InheritsFrom("AliAODRecoDecayHF2Prong")) {
2235  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2236  "%s: Objects of type %s in %s are not inherited from AliAODRecoDecayHF2Prong! Task will be disabled!",
2237  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName());
2238  params.fCandidateArray = 0;
2239  params.fInhibit = kTRUE;
2240  }
2241  }
2242  else {
2243  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2244  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
2245  params.fBranchName.Data(), params.GetName());
2246  params.fInhibit = kTRUE;
2247  }
2248  }
2249 
2250  if (params.fMCMode != kNoMC) {
2251  params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(0));
2252 
2253  if (!params.fMCContainer) params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(1));
2254 
2255  if (!params.fMCContainer) {
2256  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2257  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
2258  params.GetName());
2259  params.fInhibit = kTRUE;
2260  }
2261  }
2262 
2263  if (params.fMCMode != kMCTruth) {
2264  params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(0));
2265  if (!params.fTrackContainer) params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(1));
2266 
2267  params.fClusterContainer = GetClusterContainer(0);
2268 
2269  if (!params.fTrackContainer && !params.fClusterContainer) {
2270  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2271  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
2272  params.GetName());
2273  params.fInhibit = kTRUE;
2274  }
2275  }
2276  }
2277 }
2278 
2283 {
2284  if (!fAodEvent) return kFALSE;
2285 
2286  TString hname;
2287 
2288  // fix for temporary bug in ESDfilter
2289  // the AODs with null vertex pointer didn't pass the PhysSel
2290  if (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001) return kFALSE;
2291 
2292  for (auto &eng : fAnalysisEngines) {
2293 
2294  if (eng.fInhibit) continue;
2295 
2296  //Event selection
2297  hname = TString::Format("%s/fHistNEvents", eng.GetName());
2298  Bool_t iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
2299  if (!iseventselected) {
2300  fHistManager.FillTH1(hname, "Rejected");
2301  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2302  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
2303  TString label;
2304  do {
2305  label = GetHFEventRejectionReasonLabel(bitmap);
2306  if (label.IsNull()) break;
2307  fHistManager.FillTH1(hname, label);
2308  } while (true);
2309 
2310  continue;
2311  }
2312 
2313  fHistManager.FillTH1(hname, "Accepted");
2314 
2315  AliDebug(2, "Event selected");
2316 
2317  eng.RunAnalysis();
2318  }
2319  return kTRUE;
2320 }
2321 
2326 {
2327  TString hname;
2328  for (auto &param : fAnalysisEngines) {
2329 
2330  if (param.fInhibit) continue;
2331 
2332  if (fOutputType == kTreeOutput) {
2333  param.FillTree(fApplyKinematicCuts);
2335  }
2336  else if (fOutputType == kTHnOutput) {
2337  param.FillHnSparse(fApplyKinematicCuts);
2338  }
2339  }
2340  return kTRUE;
2341 }
2342 
2350 void AliAnalysisTaskDmesonJets::CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t& minMass, Double_t& maxMass)
2351 {
2352  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
2353 
2354  Double_t mass = part->Mass();
2355 
2356  // To make sure that the PDG mass value is not at the edge of a bin
2357  if (nbins % 2 == 0) {
2358  minMass = mass - range / 2 - range / nbins / 2;
2359  maxMass = mass + range / 2 - range / nbins / 2;
2360  }
2361  else {
2362  minMass = mass - range / 2;
2363  maxMass = mass + range / 2;
2364  }
2365 }
2366 
2373 {
2374  static TString label;
2375  label = "";
2376 
2377  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
2378  label = "NotSelTrigger";
2379  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
2380  return label.Data();
2381  }
2382  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
2383  label = "NoVertex";
2384  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
2385  return label.Data();
2386  }
2387  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
2388  label = "TooFewVtxContrib";
2389  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
2390  return label.Data();
2391  }
2392  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
2393  label = "ZVtxOutFid";
2394  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
2395  return label.Data();
2396  }
2397  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
2398  label = "Pileup";
2399  bitmap &= ~BIT(AliRDHFCuts::kPileup);
2400  return label.Data();
2401  }
2402  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
2403  label = "OutsideCentrality";
2404  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
2405  return label.Data();
2406  }
2407  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
2408  label = "PhysicsSelection";
2409  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
2410  return label.Data();
2411  }
2412  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
2413  label = "BadSPDVertex";
2414  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
2415  return label.Data();
2416  }
2417  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
2418  label = "ZVtxSPDOutFid";
2419  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
2420  return label.Data();
2421  }
2422  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
2423  label = "CentralityFlattening";
2424  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
2425  return label.Data();
2426  }
2427  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
2428  label = "BadTrackV0Correl";
2429  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
2430  return label.Data();
2431  }
2432 
2433  return label.Data();
2434 }
2435 
2442 {
2443  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
2444  PostData(eng.GetDataSlotNumber(), eng.GetTree());
2445  return eng.GetDataSlotNumber();
2446  }
2447  else {
2448  return -1;
2449  }
2450 }
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)