AliPhysics  c66ce6c (c66ce6c)
 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::AliDmesonJetInfo
50 
53 {
54  fD.SetPtEtaPhiE(0,0,0,0);
55  fSoftPionPt = 0;
56  fInvMass2Prong = 0;
57  for (auto &jet : fJets) {
58  jet.second.fMomentum.SetPtEtaPhiE(0,0,0,0);
59  jet.second.fNConstituents = 0;
60  jet.second.fNEF = 0;
61  jet.second.fMaxChargedPt = 0;
62  jet.second.fMaxNeutralPt = 0;
63  }
64 }
65 
68 {
69  Printf("Printing D Meson Jet object.");
70  Printf("D Meson: pT = %.3f, eta = %.3f, phi = %.3f, inv. mass = %.3f", fD.Pt(), fD.Eta(), fD.Phi_0_2pi(), fD.M());
71  Printf("Soft pion pT: %.3f. 2-Prong Invariant mass = %.3f", fSoftPionPt, fInvMass2Prong);
72  for (auto &jet : fJets) {
73  Printf("Jet %s: pT = %.3f, eta = %.3f, phi = %.3f", jet.first.c_str(), jet.second.Pt(), jet.second.Eta(), jet.second.Phi_0_2pi());
74  Printf("Jet N Consituents = %d", jet.second.fNConstituents);
75  }
76 }
77 
82 {
83  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
84  if (it == fJets.end()) return 0;
85 
86  Double_t z = 0;
87 
88  if ((*it).second.Pt() > 0) {
89  TVector3 dvect = fD.Vect();
90  TVector3 jvect = (*it).second.fMomentum.Vect();
91 
92  Double_t jetMom = jvect * jvect;
93 
94  if (jetMom < 1e-6) {
95  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetZ", "Zero jet momentum!");
96  z = 0.999;
97  }
98  else {
99  z = (dvect * jvect) / jetMom;
100  }
101 
102  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
103  }
104 
105  return z;
106 }
107 
113 Double_t AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetDistance(std::string n, Double_t& deta, Double_t& dphi) const
114 {
115  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
116  if (it == fJets.end()) return 0;
117 
118  dphi = TVector2::Phi_mpi_pi(fD.Phi() - (*it).second.Phi());;
119  deta = fD.Eta() - (*it).second.Eta();
120  return TMath::Sqrt(dphi*dphi + deta*deta);
121 }
122 
127 {
128  Double_t deta = 0;
129  Double_t dphi = 0;
130  return GetDistance(n, deta, dphi);
131 }
132 
134 {
135  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
136  if (it == fJets.end()) {
137  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!", n.c_str());
138  return 0;
139  }
140  return &((*it).second);
141 }
142 
144 {
145  std::map<std::string, AliJetInfo>::iterator it = fJets.find(n);
146  if (it == fJets.end()) {
147  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!", n.c_str());
148  return 0;
149  }
150  return &((*it).second);
151 }
152 
153 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoSummary
154 
158 
164  fPt(0),
165  fEta(0),
166  fPhi(0),
167  fR(0),
168  fZ(0)
169 {
170  Set(source, n);
171 }
172 
175 {
176  fPt = 0;
177  fEta = 0;
178  fPhi = 0;
179  fR = 0;
180  fZ = 0;
181 }
182 
188 {
189  std::map<std::string, AliJetInfo>::const_iterator it = source.fJets.find(n);
190  if (it == source.fJets.end()) return;
191 
192  fPt = (*it).second.Pt();
193  fEta = (*it).second.Eta();
194  fPhi = (*it).second.Phi_0_2pi();
195  fR = source.GetDistance(n);
196  fZ = source.GetZ(n);
197 }
198 
199 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonInfoSummary
200 
204 
209  fPt(0),
210  fEta(0),
211  fPhi(0)
212 {
213  Set(source);
214 }
215 
220 {
221  fPt = source.fD.Pt();
222  fEta = source.fD.Eta();
223  fPhi = source.fD.Phi_0_2pi();
224 }
225 
226 // Definitions of class AliAnalysisTaskDmesonJets::AliD0InfoSummary
227 
231 
236  AliDmesonInfoSummary(source),
237  fInvMass(source.fD.M())
238 {
239 }
240 
242 {
243  fInvMass = source.fD.M();
245 }
246 
247 // Definitions of class AliAnalysisTaskDmesonJets::AliDStarInfoSummary
248 
252 
257  AliDmesonInfoSummary(source),
258  f2ProngInvMass(source.fInvMass2Prong),
259  fDeltaInvMass(source.fD.M() - source.fInvMass2Prong)
260 {
261 }
262 
264 {
265  f2ProngInvMass = source.fInvMass2Prong;
266  fDeltaInvMass = source.fD.M() - source.fInvMass2Prong;
268 }
269 
270 // Definitions of class AliAnalysisTaskDmesonJets::AliJetDefinition
271 
275 
278  TObject(),
279  fJetType(AliJetContainer::kChargedJet),
280  fRadius(0),
281  fJetAlgo(AliJetContainer::antikt_algorithm),
282  fRecoScheme(AliJetContainer::pt_scheme),
283  fAcceptance(AliJetContainer::kUser),
284  fMinJetPt(0.),
285  fMinJetPhi(0.),
286  fMaxJetPhi(0.),
287  fMinJetEta(0.),
288  fMaxJetEta(0.),
289  fMinChargedPt(0.),
290  fMaxChargedPt(0.),
291  fMinNeutralPt(0.),
292  fMaxNeutralPt(0.)
293 {
294 }
295 
303  TObject(),
304  fJetType(type),
305  fRadius(r),
306  fJetAlgo(algo),
307  fRecoScheme(reco),
308  fAcceptance(AliJetContainer::kUser),
309  fMinJetPt(0.),
310  fMinJetPhi(0.),
311  fMaxJetPhi(0.),
312  fMinJetEta(0.),
313  fMaxJetEta(0.),
314  fMinChargedPt(0.),
315  fMaxChargedPt(0.),
316  fMinNeutralPt(0.),
317  fMaxNeutralPt(0.)
318 {
319  // By default set detector fiducial acceptance
320  switch (type) {
324  break;
327  break;
328  }
329 }
330 
335  TObject(),
336  fJetType(source.fJetType),
337  fRadius(source.fRadius),
338  fJetAlgo(source.fJetAlgo),
339  fRecoScheme(source.fRecoScheme),
340  fAcceptance(source.fAcceptance),
341  fMinJetPt(source.fMinJetPt),
342  fMinJetPhi(source.fMinJetPhi),
343  fMaxJetPhi(source.fMaxJetPhi),
344  fMinJetEta(source.fMinJetEta),
345  fMaxJetEta(source.fMaxJetEta),
346  fMinChargedPt(source.fMinChargedPt),
347  fMaxChargedPt(source.fMaxChargedPt),
348  fMinNeutralPt(source.fMinNeutralPt),
349  fMaxNeutralPt(source.fMaxNeutralPt)
350 {
351 }
352 
357 {
358  new (this) AliHFJetDefinition(source);
359  return *this;
360 }
361 
364 {
365  static TString name;
366 
367  name = AliJetContainer::GenerateJetName(fJetType, fJetAlgo, fRecoScheme, fRadius, 0, 0, "Jet");
368 
369  return name.Data();
370 }
371 
377 {
378  if (fMinJetEta < fMaxJetEta && (jet.Eta() < fMinJetEta || jet.Eta() > fMaxJetEta)) return kFALSE;
379  if (fMinJetPhi < fMaxJetPhi && (jet.Phi() < fMinJetPhi || jet.Phi() > fMaxJetPhi)) return kFALSE;
380  if (jet.Pt() < fMinJetPt) return kFALSE;
381  if (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt) return kFALSE;
382  if (jet.fMaxNeutralPt < fMinNeutralPt || jet.fMaxNeutralPt > fMaxNeutralPt) return kFALSE;
383 
384  return kTRUE;
385 }
386 
392 {
393  const AliJetInfo* jet = dMesonJet.GetJet(n);
394  if (!jet) return kFALSE;
395  return IsJetInAcceptance((*jet));
396 }
397 
402 void AliAnalysisTaskDmesonJets::AliHFJetDefinition::SetDetectorJetEtaPhiRange(const AliEMCALGeometry* const geom, Int_t run)
403 {
404  Double_t r = 0;
405  switch (fAcceptance) {
407  r = fRadius;
408  // enforce fiducial acceptance
409  /* no break */
411  SetJetEtaRange(-0.9 + r, 0.9 - r);
412  SetJetPhiRange(0, 0); // No cut on phi
413  break;
414 
416  r = fRadius;
417  // enforce fiducial acceptance
418  /* no break */
420  if (geom) {
421  SetJetEtaRange(geom->GetArm1EtaMin() + r, geom->GetArm1EtaMax() - r);
422 
423  if(run>=177295 && run<=197470) {//small SM masked in 2012 and 2013
424  SetJetPhiRange(1.405 + r,3.135 - r);
425  }
426  else {
427  SetJetPhiRange(geom->GetArm1PhiMin() * TMath::DegToRad() + r, geom->GetEMCALPhiMax() * TMath::DegToRad() - r);
428  }
429  }
430  else {
431  AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
432  SetJetEtaRange(-0.7 + r, 0.7 - r);
433  SetJetPhiRange(1.405 + r, 3.135 - r);
434  }
435  break;
436 
438  r = fRadius;
439  // enforce fiducial acceptance
440  /* no break */
442  if (geom) {
443  SetJetEtaRange(geom->GetArm1EtaMin() + r, geom->GetArm1EtaMax() - r);
444  SetJetPhiRange(geom->GetDCALPhiMin() * TMath::DegToRad() + r, geom->GetDCALPhiMax() * TMath::DegToRad() - r);
445  }
446  else {
447  AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for DCAL year 2015!!");
448  SetJetEtaRange(-0.7 + r, 0.7 - r);
449  SetJetPhiRange(4.538 + r, 5.727 - r);
450  }
451  break;
452 
454  // Nothing to be done
455  break;
456  }
457 }
458 
465 {
466  if (lhs.fJetType > rhs.fJetType) return false;
467  else if (lhs.fJetType < rhs.fJetType) return true;
468  else {
469  if (lhs.fRadius > rhs.fRadius) return false;
470  else if (lhs.fRadius < rhs.fRadius) return true;
471  else {
472  if (lhs.fJetAlgo > rhs.fJetAlgo) return false;
473  else if (lhs.fJetAlgo < rhs.fJetAlgo) return true;
474  else {
475  if (lhs.fRecoScheme < rhs.fRecoScheme) return true;
476  else return false;
477  }
478  }
479  }
480 }
481 
488 {
489  if (lhs.fJetType != rhs.fJetType) return false;
490  if (lhs.fRadius != rhs.fRadius) return false;
491  if (lhs.fJetAlgo != rhs.fJetAlgo) return false;
492  if (lhs.fRecoScheme != rhs.fRecoScheme) return false;
493  return true;
494 }
495 
496 // Definitions of class AliAnalysisTaskDmesonJets::AnalysisEngine
497 
501 
504  TObject(),
505  fCandidateType(kD0toKpi),
506  fCandidateName(),
507  fCandidatePDG(0),
508  fNDaughters(0),
509  fBranchName(),
510  fMCMode(kNoMC),
511  fNMassBins(0),
512  fMinMass(0),
513  fMaxMass(0),
514  fRDHFCuts(0),
515  fRejectedOrigin(0),
516  fAcceptedDecay(0),
517  fInhibit(kFALSE),
518  fJetDefinitions(),
519  fPtBinWidth(0.5),
520  fMaxPt(100),
521  fDataSlotNumber(-1),
522  fTree(0),
523  fCurrentDmesonJetInfo(0),
524  fCurrentJetInfo(0),
525  fCandidateArray(0),
526  fMCContainer(0),
527  fTrackContainer(0),
528  fClusterContainer(0),
529  fAodEvent(0),
530  fFastJetWrapper(0),
531  fHistManager(0)
532 {
533 }
534 
543  TObject(),
544  fCandidateType(type),
545  fCandidateName(),
546  fCandidatePDG(0),
547  fNDaughters(0),
548  fBranchName(),
549  fMCMode(MCmode),
550  fNMassBins(nMassBins),
551  fMinMass(0),
552  fMaxMass(0),
553  fRDHFCuts(cuts),
554  fRejectedOrigin(kUnknownQuark | kFromBottom),
555  fAcceptedDecay(kAnyDecay),
556  fInhibit(kFALSE),
557  fJetDefinitions(),
558  fPtBinWidth(0.5),
559  fMaxPt(100),
560  fDataSlotNumber(-1),
561  fTree(0),
562  fCurrentDmesonJetInfo(0),
563  fCurrentJetInfo(0),
564  fCandidateArray(0),
565  fMCContainer(0),
566  fTrackContainer(0),
567  fClusterContainer(0),
568  fAodEvent(0),
569  fFastJetWrapper(0),
570  fHistManager(0)
571 {
572  SetCandidateProperties(range);
573 }
574 
579  TObject(source),
580  fCandidateType(source.fCandidateType),
581  fCandidateName(source.fCandidateName),
582  fCandidatePDG(source.fCandidatePDG),
583  fNDaughters(source.fNDaughters),
584  fBranchName(source.fBranchName),
585  fMCMode(source.fMCMode),
586  fNMassBins(source.fNMassBins),
587  fMinMass(source.fMinMass),
588  fMaxMass(source.fMaxMass),
589  fRDHFCuts(),
590  fRejectedOrigin(source.fRejectedOrigin),
591  fAcceptedDecay(source.fAcceptedDecay),
592  fInhibit(source.fInhibit),
593  fJetDefinitions(source.fJetDefinitions),
594  fPtBinWidth(source.fPtBinWidth),
595  fMaxPt(source.fMaxPt),
596  fDataSlotNumber(-1),
597  fTree(0),
598  fCurrentDmesonJetInfo(0),
599  fCurrentJetInfo(0),
600  fCandidateArray(source.fCandidateArray),
601  fMCContainer(source.fMCContainer),
602  fTrackContainer(source.fTrackContainer),
603  fClusterContainer(source.fClusterContainer),
604  fAodEvent(source.fAodEvent),
606  fHistManager(source.fHistManager)
607 {
608  SetRDHFCuts(source.fRDHFCuts);
609 }
610 
611 // Destructor
613 {
614  delete fRDHFCuts;
615 }
616 
621 {
622  new (this) AnalysisEngine(source);
623  return *this;
624 }
625 
630 {
631  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
632  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
633  }
634 
635  return kFALSE;
636 }
637 
639 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const geom, Int_t runNumber)
640 {
641  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
642  fJetDefinitions[i].SetDetectorJetEtaPhiRange(geom, runNumber);
643  }
644 }
645 
650 {
651  switch (fCandidateType) {
652  case kD0toKpi:
653  fCandidatePDG = 421;
654  fCandidateName = "D0";
655  fNDaughters = 2;
656  fPDGdaughters.Set(fNDaughters);
657  fPDGdaughters.Reset();
658  fPDGdaughters[0] = 211; // pi
659  fPDGdaughters[1] = 321; // K
660  fBranchName = "D0toKpi";
661  fAcceptedDecay = kD0toKpi;
662  if (!fRDHFCuts) {
663  fRDHFCuts = new AliRDHFCutsD0toKpi();
664  fRDHFCuts->SetStandardCutsPP2010();
665  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
666  fRDHFCuts->SetTriggerClass("","");
667  }
668  break;
669  case kDstartoKpipi:
670  fCandidatePDG = 413;
671  fCandidateName = "DStar";
672  fNDaughters = 3;
673  fPDGdaughters.Set(fNDaughters);
674  fPDGdaughters.Reset();
675  fPDGdaughters[0] = 211; // pi soft
676  fPDGdaughters[1] = 211; // pi fromD0
677  fPDGdaughters[2] = 321; // K from D0
678  fBranchName = "Dstar";
679  fAcceptedDecay = kDstartoKpipi;
680  if (!fRDHFCuts) {
681  fRDHFCuts = new AliRDHFCutsDStartoKpipi();
682  fRDHFCuts->SetStandardCutsPP2010();
683  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
684  fRDHFCuts->SetTriggerClass("","");
685  }
686  break;
687  default:
688  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
689  }
690 
691  CalculateMassLimits(range, fCandidatePDG, fNMassBins, fMinMass, fMaxMass);
692 }
693 
698 {
699  if (fRDHFCuts) delete fRDHFCuts;
700  fRDHFCuts = cuts;
701 }
702 
707 {
708  if (!cuts) return;
709  if (fRDHFCuts) delete fRDHFCuts;
710  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
711 }
712 
717 {
718  static TString name;
719 
720  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
721 
722  return name.Data();
723 }
724 
729 {
730  static TString name;
731 
732  name = fCandidateName;
733  switch (fMCMode) {
734  case kBackgroundOnly:
735  name += "_kBackgroundOnly";
736  break;
737  case kSignalOnly:
738  name += "_kSignalOnly";
739  break;
740  case kMCTruth:
741  name += "_MCTruth";
742  break;
743  default:
744  break;
745  }
746 
747  return name.Data();
748 }
749 
757 {
758  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
759 
760  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
761  fJetDefinitions.push_back(def);
762  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
763  "Total number of jet definitions is now %lu.",
764  def.GetName(), GetName(), fJetDefinitions.size());
765  // For detector level set maximum track pt to 100 GeV/c
766  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
767  }
768  else {
769  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
770  }
771 
772  return &(*it);
773 }
774 
786 {
787  AliHFJetDefinition def(type, r, algo, reco);
788 
789  return AddJetDefinition(def);
790 }
791 
797 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
798 {
799  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
800  while (it != fJetDefinitions.end() && (*it) != def) it++;
801  return it;
802 }
803 
810 {
811  if (lhs.fCandidateType > rhs.fCandidateType) return false;
812  else if (lhs.fCandidateType < rhs.fCandidateType) return true;
813  else {
814  if (lhs.fMCMode < rhs.fMCMode) return true;
815  else return false;
816  }
817 }
818 
825 {
826  if (lhs.fCandidateType != rhs.fCandidateType) return false;
827  if (lhs.fMCMode != rhs.fMCMode) return false;
828  return true;
829 }
830 
839 {
840  DmesonJet.fD.SetPtEtaPhiM(part->Pt(), part->Eta(), part->Phi(), part->M());
841  return kTRUE;
842 }
843 
852 {
853  if (fCandidateType == kD0toKpi) { // D0 candidate
854  return ExtractD0Attributes(Dcand, DmesonJet, i);
855  }
856  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
857  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
858  }
859  else {
860  return kFALSE;
861  }
862 }
863 
872 {
873  Int_t MCtruthPdgCode = 0;
874 
875  Double_t invMassD = 0;
876 
877  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
878  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
879  if (mcLab >= 0) {
880  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
881 
882  if (aodMcPart) {
883  if (fRejectedOrigin && fMCMode == kSignalOnly) {
884  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
885 
886  if ((origin & fRejectedOrigin) == origin) return kFALSE;
887  }
888  MCtruthPdgCode = aodMcPart->PdgCode();
889  }
890  }
891  }
892 
893  //AliDebug(2,"Checking if D0 meson is selected");
894  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
895  if (isSelected == 1) { // selected as a D0
896  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
897 
898  if (fMCMode == kNoMC ||
899  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
900  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly)) {
901  // 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)
902  //AliDebug(2,"Selected as D0");
903  invMassD = Dcand->InvMassD0();
904  }
905  else { // conditions above not passed, so return FALSE
906  return kFALSE;
907  }
908  }
909  else if (isSelected == 2) { // selected as a D0bar
910  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
911 
912  if (fMCMode == kNoMC ||
913  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
914  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly)) {
915  // 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)
916  //AliDebug(2,"Selected as D0bar");
917  invMassD = Dcand->InvMassD0bar();
918  }
919  else { // conditions above not passed, so return FALSE
920  return kFALSE;
921  }
922  }
923  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
924  //AliDebug(2,"Selected as either D0 or D0bar");
925 
926  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
927  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
928  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kBackgroundOnly)) {
929  if (i > 0) return kFALSE;
930  //AliDebug(2, "MC truth is D0");
931  invMassD = Dcand->InvMassD0();
932  }
933  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
934  (MCtruthPdgCode == fCandidatePDG && fMCMode == kBackgroundOnly)) {
935  if (i > 0) return kFALSE;
936  //AliDebug(2, "MC truth is D0bar");
937  invMassD = Dcand->InvMassD0bar();
938  }
939  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
940 
941  // Only accept it if background-only OR background-and-signal was requested
942  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
943  // Select D0 or D0bar depending on the i-parameter
944  if (i == 0) {
945  //AliDebug(2, "Returning invariant mass with D0 hypothesis");
946  invMassD = Dcand->InvMassD0();
947  }
948  else if (i == 1) {
949  //AliDebug(2, "Returning invariant mass with D0bar hypothesis");
950  invMassD = Dcand->InvMassD0bar();
951  }
952  else { // i > 1
953  return kFALSE;
954  }
955  }
956  else { // signal-only was requested but this is not a true D0
957  return kFALSE;
958  }
959  }
960  }
961  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
962  return kTRUE;
963 }
964 
973 {
974  if (i > 0) return kFALSE; // only one mass hypothesis for the D*
975 
976  Int_t MCtruthPdgCode = 0;
977 
978  Double_t invMassD = 0;
979 
980  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
981  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
982  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
983 
984  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
985  //AliDebug(2, Form("MC label is %d", mcLab));
986  if (mcLab >= 0) {
987  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
988 
989  if (aodMcPart) {
990  if (fRejectedOrigin && fMCMode == kSignalOnly) {
991  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
992 
993  if ((origin & fRejectedOrigin) == origin) return kFALSE;
994  }
995 
996  MCtruthPdgCode = aodMcPart->PdgCode();
997  //AliDebug(2, Form("MC truth pdg code is %d",MCtruthPdgCode));
998  }
999  }
1000  }
1001 
1002  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1003  if (fMCMode == kNoMC ||
1004  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1005  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1006  // 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)
1007  invMassD = DstarCand->InvMassDstarKpipi();
1008  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1009  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1010  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1011  return kTRUE;
1012  }
1013  else { // conditions above not passed, so return FALSE
1014  return kFALSE;
1015  }
1016 }
1017 
1025 {
1026  if (!part) return kDecayOther;
1027  if (!mcArray) return kDecayOther;
1028 
1030 
1031  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1032 
1033  if (part->GetNDaughters() == 2) {
1034 
1035  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1036  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1037 
1038  if (!d1 || !d2) {
1039  return decay;
1040  }
1041 
1042  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1043  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1044 
1045  if (absPdgPart == 421) { // D0 -> K pi
1046 
1047  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1048  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1049  decay = kDecayD0toKpi;
1050  }
1051  }
1052 
1053  if (absPdgPart == 413) { // D* -> D0 pi
1054 
1055  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1056  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1057  if (D0decay == kDecayD0toKpi) {
1058  decay = kDecayDStartoKpipi;
1059  }
1060  }
1061 
1062  if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1063  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1064  if (D0decay == kDecayD0toKpi) {
1065  decay = kDecayDStartoKpipi;
1066  }
1067  }
1068  }
1069  }
1070 
1071  return decay;
1072 }
1073 
1081 {
1082  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1083 
1084  if (!part) return kUnknownQuark;
1085  if (!mcArray) return kUnknownQuark;
1086 
1087  Int_t pdgGranma = 0;
1088  Int_t mother = part->GetMother();
1089  Int_t istep = 0;
1090  Int_t abspdgGranma = 0;
1091  Bool_t isFromB = kFALSE;
1092  Bool_t isQuarkFound = kFALSE;
1093 
1094  while (mother >= 0) {
1095  istep++;
1096  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1097  if (mcGranma >= 0) {
1098  pdgGranma = mcGranma->GetPdgCode();
1099  abspdgGranma = TMath::Abs(pdgGranma);
1100  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1101  isFromB = kTRUE;
1102  }
1103 
1104  if (abspdgGranma == 4 || abspdgGranma == 5) isQuarkFound = kTRUE;
1105  mother = mcGranma->GetMother();
1106  }
1107  else {
1108  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1109  break;
1110  }
1111  }
1112 
1113  if (isQuarkFound) {
1114  if (isFromB) {
1115  return kFromBottom;
1116  }
1117  else {
1118  return kFromCharm;
1119  }
1120  }
1121  else {
1122  return kUnknownQuark;
1123  }
1124 }
1125 
1128 {
1129  fDmesonJets.clear();
1130 
1131  if (fMCMode == kMCTruth) {
1132  RunParticleLevelAnalysis();
1133  }
1134  else {
1135  RunDetectorLevelAnalysis();
1136  }
1137 }
1138 
1141 {
1142  const Int_t nD = fCandidateArray->GetEntriesFast();
1143 
1144  AliDmesonJetInfo DmesonJet;
1145 
1146  Int_t nAccCharm = 0;
1147  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1148  Int_t isSelected = 0;
1149 
1150  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1151  if (!charmCand) continue;
1152 
1153  Int_t nprongs = charmCand->GetNProngs();
1154 
1155  if (fCandidateType == kDstartoKpipi) {
1156  if (!charmCand->InheritsFrom("AliAODRecoCascadeHF")) {
1157  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::RunDetectorLevelAnalysis","Candidate type is D* but object type is wrong (should be AliAODRecoCascadeHF)");
1158  continue;
1159  }
1160  }
1161 
1162  // region of interest + cuts
1163  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1164 
1165  //candidate selected by cuts and PID
1166  isSelected = fRDHFCuts->IsSelected(charmCand, AliRDHFCuts::kAll, fAodEvent); //selected
1167 
1168  if (!isSelected) continue;
1169 
1170  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1171  DmesonJet.Reset();
1172  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1173  for (std::vector<AliHFJetDefinition>::iterator itdef = fJetDefinitions.begin(); itdef != fJetDefinitions.end(); itdef++) {
1174  if (!FindJet(charmCand, DmesonJet, *itdef)) {
1175  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1176  (*itdef).GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1177  }
1178  }
1179  fDmesonJets.push_back(DmesonJet);
1180  }
1181  }
1182  nAccCharm++;
1183  } // end of D cand loop
1184 
1185  TString hname;
1186 
1187  hname = TString::Format("%s/fHistNAcceptedDmesons", GetName());
1188  fHistManager->FillTH1(hname, nAccCharm);
1189 
1190  hname = TString::Format("%s/fHistNDmesons", GetName());
1191  fHistManager->FillTH1(hname, nD);
1192 }
1193 
1204 {
1205  TString hname;
1206 
1208  fFastJetWrapper->SetR(jetDef.fRadius);
1211 
1212  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1213 
1214  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1215  fTrackContainer->SetDMesonCandidate(Dcand);
1216  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1217  AddInputVectors(fTrackContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1218 
1219  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1220  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1221  const TObjArray daughterNotInJet = fTrackContainer->GetDaughterList();
1222  for (Int_t i = 0; i < daughterNotInJet.GetEntriesFast(); i++) {
1223  AliVParticle* daughter = static_cast<AliVParticle*>(daughterNotInJet.At(i));
1224  if (!daughter) continue;
1225  histDaughterNotInJet->Fill(daughter->Pt());
1226  }
1227  }
1228 
1229  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1230  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1231  AddInputVectors(fClusterContainer, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1232  }
1233 
1234  // run jet finder
1235  fFastJetWrapper->Run();
1236 
1237  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1238 
1239  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1240  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1241 
1242  Bool_t isDmesonJet = kFALSE;
1243 
1244  Double_t maxChPt = 0;
1245  Double_t maxNePt = 0;
1246  Double_t totalNeutralPt = 0;
1247 
1248  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1249  if (constituents[ic].user_index() == 0) {
1250  isDmesonJet = kTRUE;
1251  }
1252  else if (constituents[ic].user_index() >= 100) {
1253  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1254  }
1255  else if (constituents[ic].user_index() <= -100) {
1256  totalNeutralPt += constituents[ic].pt();
1257  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1258  }
1259  }
1260 
1261  if (isDmesonJet) {
1262  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1263  DmesonJet.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1264  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1265  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1266  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1267 
1268  return kTRUE;
1269  }
1270  }
1271 
1272  return kFALSE;
1273 }
1274 
1278 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist)
1279 {
1280  AliTLorentzVector part;
1281  for (Int_t i = 0; i < cont->GetNEntries(); i++) {
1282  cont->GetMomentum(part, i);
1283  UInt_t rejectionReason = 0;
1284  if (!cont->AcceptObject(i, rejectionReason)) {
1285  rejectHist->Fill(cont->GetRejectionReasonBitPosition(rejectionReason), part.Pt());
1286  continue;
1287  }
1288  Int_t uid = offset >= 0 ? i : -i;
1289  uid += offset;
1290  fFastJetWrapper->AddInputVector(part.Px(), part.Py(), part.Pz(), part.E(), uid);
1291  }
1292 }
1293 
1296 {
1297  TString hname;
1298 
1299  fMCContainer->SetSpecialPDG(fCandidatePDG);
1300  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1301  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1302 
1303  std::map<int, AliDmesonJetInfo> dMesonJets;
1304 
1305  for (auto &jetDef : fJetDefinitions) {
1306 
1308  fFastJetWrapper->SetR(jetDef.fRadius);
1311 
1312  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1313  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1314 
1315  fFastJetWrapper->Run();
1316 
1317  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1318 
1319  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1320  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1321 
1322  Bool_t isDmesonJet = kFALSE;
1323 
1324  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1325  Int_t iPart = constituents[ic].user_index() - 100;
1326  AliVParticle* part = fMCContainer->GetParticle(iPart);
1327  if (!part) {
1328  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1329  continue;
1330  }
1331  if (part->PdgCode() == fCandidatePDG) {
1332  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = dMesonJets.find(iPart);
1333  if (dMesonJetIt == dMesonJets.end()) {
1334  std::pair<int, AliDmesonJetInfo> element;
1335  element.first = iPart;
1336  element.second = AliDmesonJetInfo();
1337  dMesonJetIt = dMesonJets.insert(element).first;
1338  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1339  }
1340 
1341  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1342  break;
1343  }
1344  }
1345  }
1346  }
1347 
1348  for (std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = dMesonJets.begin();
1349  dMesonJetIt != dMesonJets.end();
1350  dMesonJetIt++) {
1351  fDmesonJets.push_back((*dMesonJetIt).second);
1352  }
1353 }
1354 
1359 {
1360  TString classname;
1361  switch (fCandidateType) {
1362  case kD0toKpi:
1363  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1364  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1365  break;
1366  case kDstartoKpipi:
1367  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1368  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1369  break;
1370  }
1371  TString treeName = TString::Format("AliAnalysisTaskDmesonJets_%s", GetName());
1372  fTree = new TTree(treeName, treeName);
1373  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo, 32000, 0);
1374  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1375  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1376  fCurrentJetInfo[i] = new AliJetInfoSummary();
1377  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1378  }
1379 
1380  return fTree;
1381 }
1382 
1387 {
1388  TString hname;
1389 
1390  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
1391 
1392  for (auto &jetDef : fJetDefinitions) {
1393 
1394  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
1395 
1396  Double_t radius = jetDef.fRadius;
1397 
1398  TString title[30] = {""};
1399  Int_t nbins[30] = {0 };
1400  Double_t min [30] = {0.};
1401  Double_t max [30] = {0.};
1402  Int_t dim = 0 ;
1403 
1404  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1405  nbins[dim] = nPtBins;
1406  min[dim] = 0;
1407  max[dim] = fMaxPt;
1408  dim++;
1409 
1410  if ((enabledAxis & kPositionD) != 0) {
1411  title[dim] = "#eta_{D}";
1412  nbins[dim] = 50;
1413  min[dim] = -1;
1414  max[dim] = 1;
1415  dim++;
1416 
1417  title[dim] = "#phi_{D} (rad)";
1418  nbins[dim] = 150;
1419  min[dim] = 0;
1420  max[dim] = TMath::TwoPi();
1421  dim++;
1422  }
1423 
1424  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1425  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1426  nbins[dim] = fNMassBins;
1427  min[dim] = fMinMass;
1428  max[dim] = fMaxMass;
1429  dim++;
1430  }
1431 
1432  if (fCandidateType == kD0toKpi) {
1433  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1434  nbins[dim] = fNMassBins;
1435  min[dim] = fMinMass;
1436  max[dim] = fMaxMass;
1437  dim++;
1438  }
1439 
1440  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1441  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1442  nbins[dim] = fNMassBins;
1443  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1444  dim++;
1445  }
1446 
1447  if (fCandidateType == kDstartoKpipi) {
1448  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1449  nbins[dim] = fNMassBins*6;
1450  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1451 
1452  // subtract mass of D0
1453  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1454  min[dim] -= D0mass;
1455  max[dim] -= D0mass;
1456 
1457  dim++;
1458  }
1459 
1460  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1461  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1462  nbins[dim] = 100;
1463  min[dim] = 0;
1464  max[dim] = 25;
1465  dim++;
1466  }
1467 
1468  title[dim] = "#it{z}_{D}";
1469  nbins[dim] = 110;
1470  min[dim] = 0;
1471  max[dim] = 1.10;
1472  dim++;
1473 
1474  if ((enabledAxis & kDeltaR) != 0) {
1475  title[dim] = "#Delta R_{D-jet}";
1476  nbins[dim] = 100;
1477  min[dim] = 0;
1478  max[dim] = radius * 1.5;
1479  dim++;
1480  }
1481 
1482  if ((enabledAxis & kDeltaEta) != 0) {
1483  title[dim] = "#eta_{D} - #eta_{jet}";
1484  nbins[dim] = 100;
1485  min[dim] = -radius * 1.2;
1486  max[dim] = radius * 1.2;
1487  dim++;
1488  }
1489 
1490  if ((enabledAxis & kDeltaPhi) != 0) {
1491  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1492  nbins[dim] = 100;
1493  min[dim] = -radius * 1.2;
1494  max[dim] = radius * 1.2;
1495  dim++;
1496  }
1497 
1498  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1499  nbins[dim] = nPtBins;
1500  min[dim] = 0;
1501  max[dim] = fMaxPt;
1502  dim++;
1503 
1504  if ((enabledAxis & kPositionJet) != 0) {
1505  title[dim] = "#eta_{jet}";
1506  nbins[dim] = 50;
1507  min[dim] = -1;
1508  max[dim] = 1;
1509  dim++;
1510 
1511  title[dim] = "#phi_{jet} (rad)";
1512  nbins[dim] = 150;
1513  min[dim] = 0;
1514  max[dim] = TMath::TwoPi();
1515  dim++;
1516  }
1517 
1518  if ((enabledAxis & kJetConstituents) != 0) {
1519  title[dim] = "No. of constituents";
1520  nbins[dim] = 50;
1521  min[dim] = -0.5;
1522  max[dim] = 49.5;
1523  dim++;
1524  }
1525 
1526  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1527  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1528  for (Int_t j = 0; j < dim; j++) {
1529  h->GetAxis(j)->SetTitle(title[j]);
1530  }
1531  }
1532 }
1533 
1538 {
1539  TString hname;
1540 
1541  for (Int_t id = 0; id < fDmesonJets.size(); id++) {
1542  fCurrentDmesonJetInfo->Set(fDmesonJets[id]);
1543  Int_t accJets = 0;
1544  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1545  fCurrentJetInfo[ij]->Reset();
1546  AliJetInfo* jet = fDmesonJets[id].GetJet(fJetDefinitions[ij].GetName());
1547  if (!jet) continue;
1548  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1549  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1550  fHistManager->FillTH1(hname, jet->Pt());
1551  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1552  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1553  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1554  fHistManager->FillTH1(hname, jet->Eta());
1555  continue;
1556  }
1557  fCurrentJetInfo[ij]->Set(fDmesonJets[id], fJetDefinitions[ij].GetName());
1558  accJets++;
1559  }
1560  if (accJets > 0) {
1561  fTree->Fill();
1562  }
1563  else {
1564  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1565  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Pt());
1566  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1567  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Phi_0_2pi());
1568  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1569  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Eta());
1570  if (fCandidateType == kD0toKpi) {
1571  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1572  fHistManager->FillTH1(hname, fDmesonJets[id].fD.M());
1573  }
1574  else if (fCandidateType == kDstartoKpipi) {
1575  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1576  fHistManager->FillTH1(hname, fDmesonJets[id].fInvMass2Prong);
1577 
1578  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1579  fHistManager->FillTH1(hname, fDmesonJets[id].fD.M() - fDmesonJets[id].fInvMass2Prong);
1580  }
1581  }
1582  }
1583  return kTRUE;
1584 }
1585 
1590 {
1591  TString hname;
1592 
1593  for (Int_t id = 0; id < fDmesonJets.size(); id++) {
1594  if (!IsAnyJetInAcceptance(fDmesonJets[id])) {
1595  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1596  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Pt());
1597  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1598  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Phi_0_2pi());
1599  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1600  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Eta());
1601  }
1602  }
1603 
1604  for (auto &jetDef : fJetDefinitions) {
1605 
1606  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1607  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
1608 
1609  for (Int_t id = 0; id < fDmesonJets.size(); id++) {
1610  const AliJetInfo* jet = fDmesonJets[id].GetJet(jetDef.GetName());
1611  if (!jet) continue;
1612  if (!jetDef.IsJetInAcceptance(*jet)) {
1613  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
1614  fHistManager->FillTH1(hname, jet->Pt());
1615  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
1616  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1617  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
1618  fHistManager->FillTH1(hname, jet->Eta());
1619  continue;
1620  }
1621  FillHnSparse(h, fDmesonJets[id], jetDef.GetName());
1622  }
1623  }
1624 
1625  return kTRUE;
1626 }
1627 
1633 Bool_t AliAnalysisTaskDmesonJets::AnalysisEngine::FillHnSparse(THnSparse* h, const AliDmesonJetInfo& DmesonJet, std::string n)
1634 {
1635  // Fill the THnSparse histogram.
1636 
1637  Double_t contents[30] = {0.};
1638 
1639  Double_t z = DmesonJet.GetZ(n);
1640  Double_t deltaPhi = 0;
1641  Double_t deltaEta = 0;
1642  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
1643 
1644  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
1645  if (it == DmesonJet.fJets.end()) return kFALSE;
1646 
1647  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
1648  TString title(h->GetAxis(i)->GetTitle());
1649  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
1650  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
1651  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
1652  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
1653  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
1654  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
1655  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
1656  else if (title=="#it{z}_{D}") contents[i] = z;
1657  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
1658  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
1659  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
1660  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
1661  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
1662  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
1663  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
1664  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
1665  }
1666 
1667  h->Fill(contents);
1668 
1669  return kTRUE;
1670 }
1671 
1672 // Definitions of class AliAnalysisTaskDmesonJets
1673 
1677 
1681  fAnalysisEngines(),
1682  fEnabledAxis(0),
1683  fTreeOutput(kFALSE),
1684  fHistManager(),
1685  fApplyKinematicCuts(kTRUE),
1686  fNOutputTrees(0),
1687  fAodEvent(0),
1688  fFastJetWrapper(0)
1689 {
1690  SetMakeGeneralHistograms(kTRUE);
1691 }
1692 
1696 AliAnalysisTaskDmesonJets::AliAnalysisTaskDmesonJets(const char* name, Int_t nOutputTrees) :
1697  AliAnalysisTaskEmcalLight(name, kTRUE),
1698  fAnalysisEngines(),
1699  fEnabledAxis(k2ProngInvMass),
1700  fTreeOutput(kFALSE),
1701  fHistManager(name),
1702  fApplyKinematicCuts(kTRUE),
1703  fNOutputTrees(nOutputTrees),
1704  fAodEvent(0),
1705  fFastJetWrapper(0)
1706 {
1707  SetMakeGeneralHistograms(kTRUE);
1708  for (Int_t i = 0; i < nOutputTrees; i++){
1709  DefineOutput(2+i, TTree::Class());
1710  }
1711 }
1712 
1715 {
1716  if (fFastJetWrapper) delete fFastJetWrapper;
1717 }
1718 
1726 {
1727  AliRDHFCuts* analysiscuts = 0;
1728  TFile* filecuts = TFile::Open(cutfname);
1729  if (!filecuts || filecuts->IsZombie()) {
1730  ::Warning("AddTaskDmesonJets", "Input file not found: will use std cuts.");
1731  filecuts = 0;
1732  }
1733 
1734  if (filecuts) {
1735  analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
1736  if (!analysiscuts) {
1737  ::Warning("AddTaskDmesonJetCorr", "Could not find analysis cuts '%s' in '%s'. Using std cuts.", cutsname.Data(), cutfname.Data());
1738  }
1739  }
1740 
1741  return analysiscuts;
1742 }
1743 
1753 {
1755  return AddAnalysisEngine(type, MCmode, jetDef, cutfname);
1756 }
1757 
1767 {
1768  AliRDHFCuts* cuts = 0;
1769 
1770  if (!cutfname.IsNull()) {
1771  TString cutsname;
1772 
1773  switch (type) {
1774  case kD0toKpi :
1775  cutsname = "D0toKpiCuts";
1776  break;
1777  case kDstartoKpipi :
1778  cutsname = "DStartoKpipiCuts";
1779  break;
1780  default:
1781  return 0;
1782  }
1783 
1784  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
1785  }
1786 
1787  AnalysisEngine eng(type, MCmode, cuts);
1788 
1789  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
1790 
1791  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
1792  eng.AddJetDefinition(jetDef);
1793  it = fAnalysisEngines.insert(it, eng);
1794  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(jetDef), fAnalysisEngines.size());
1795  }
1796  else {
1797  AnalysisEngine* found_eng = &(*it);
1798  ::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());
1799  found_eng->AddJetDefinition(jetDef);
1800 
1801  if (cuts && found_eng->fRDHFCuts != 0) {
1802  ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
1803  found_eng->SetRDHFCuts(cuts);
1804  }
1805  }
1806 
1807  return &(*it);
1808 }
1809 
1810 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
1811 {
1812  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
1813  while (it != fAnalysisEngines.end() && (*it) < eng) it++;
1814  return it;
1815 }
1816 
1819 {
1820  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
1821 
1823 
1824  // Define histograms
1825  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
1826 
1827  TString hname;
1828  TString htitle;
1829  TH1* h = 0;
1830  Int_t treeSlot = 0;
1831 
1832  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
1833  for (auto &param : fAnalysisEngines) {
1834  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
1835 
1836  fHistManager.CreateHistoGroup(param.GetName());
1837 
1838  param.fHistManager = &fHistManager;
1839 
1840  hname = TString::Format("%s/fHistNAcceptedDmesons", param.GetName());
1841  htitle = hname + ";Number of D accepted meson candidates;counts";
1842  h = fHistManager.CreateTH1(hname, htitle, 51, -0.5, 50.5);
1843 
1844  hname = TString::Format("%s/fHistNDmesons", param.GetName());
1845  htitle = hname + ";Number of D meson candidates;counts";
1846  h = fHistManager.CreateTH1(hname, htitle, 101, -0.5, 100.5);
1847 
1848  hname = TString::Format("%s/fHistNEvents", param.GetName());
1849  htitle = hname + ";Event status;counts";
1850  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
1851  h->GetXaxis()->SetBinLabel(1, "Accepted");
1852  h->GetXaxis()->SetBinLabel(2, "Rejected");
1853 
1854  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
1855  htitle = hname + ";Rejection reason;counts";
1856  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
1857 
1858  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
1859  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
1860  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
1861 
1862  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
1863  htitle = hname + ";#it{#eta}_{D};counts";
1864  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
1865 
1866  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
1867  htitle = hname + ";#it{#phi}_{D};counts";
1868  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
1869 
1870  if (param.fCandidateType == kD0toKpi) {
1871  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
1872  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
1873  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
1874  }
1875  else if (param.fCandidateType == kDstartoKpipi) {
1876  Double_t min = 0;
1877  Double_t max = 0;
1878 
1879  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
1880  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
1881  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
1882  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
1883 
1884  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1885  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
1886  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
1887  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
1888  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
1889  }
1890 
1891  for (std::vector<AliHFJetDefinition>::iterator itdef = param.fJetDefinitions.begin(); itdef != param.fJetDefinitions.end(); itdef++) {
1892  AliHFJetDefinition* jetDef = &(*itdef);
1893  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef->GetName());
1894 
1895  fHistManager.CreateHistoGroup(jetDef->GetName(), param.GetName());
1896 
1897  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef->GetName());
1898  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
1899  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
1900  SetRejectionReasonLabels(h->GetXaxis());
1901 
1902  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef->GetName());
1903  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
1904  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
1905  SetRejectionReasonLabels(h->GetXaxis());
1906 
1907  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef->GetName());
1908  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
1909  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
1910  SetRejectionReasonLabels(h->GetXaxis());
1911 
1912  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef->GetName());
1913  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
1914  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
1915 
1916  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef->GetName());
1917  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
1918  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
1919 
1920  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef->GetName());
1921  htitle = hname + ";#it{#eta}_{jet};counts";
1922  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
1923 
1924  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef->GetName());
1925  htitle = hname + ";#it{#phi}_{jet};counts";
1926  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
1927  }
1928  if (fTreeOutput) {
1929  param.BuildTree();
1930  if (treeSlot < fNOutputTrees) {
1931  param.AssignDataSlot(treeSlot+2);
1932  treeSlot++;
1934  }
1935  else {
1936  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
1937  }
1938  }
1939  else {
1940  param.BuildHnSparse(fEnabledAxis);
1941  }
1942  }
1943 
1945 
1946  PostData(1, fOutput);
1947 }
1948 
1952 {
1954 
1955  // Load the event
1956  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
1957 
1958  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
1959 
1960  fFastJetWrapper->SetAreaType(fastjet::active_area);
1962 
1963  if (!fAodEvent) {
1964  AliError(Form("This task need an AOD event! Task '%s' will be disabled!", GetName()));
1965  return;
1966  }
1967 
1968  for (auto &params : fAnalysisEngines) {
1969 
1970  params.fAodEvent = fAodEvent;
1971  params.fFastJetWrapper = fFastJetWrapper;
1972  params.Init(fGeom, fAodEvent->GetRunNumber());
1973 
1974  if (params.fMCMode != kMCTruth) {
1975  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
1976 
1977  if (params.fCandidateArray) {
1978  if (!params.fCandidateArray->GetClass()->InheritsFrom("AliAODRecoDecayHF2Prong")) {
1979  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
1980  "%s: Objects of type %s in %s are not inherited from AliAODRecoDecayHF2Prong! Task will be disabled!",
1981  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName());
1982  params.fCandidateArray = 0;
1983  params.fInhibit = kTRUE;
1984  }
1985  }
1986  else {
1987  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
1988  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
1989  params.fBranchName.Data(), params.GetName());
1990  params.fInhibit = kTRUE;
1991  }
1992  }
1993 
1994  if (params.fMCMode != kNoMC) {
1995  params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(0));
1996 
1997  if (!params.fMCContainer) params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(1));
1998 
1999  if (!params.fMCContainer) {
2000  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2001  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
2002  params.GetName());
2003  params.fInhibit = kTRUE;
2004  }
2005  }
2006 
2007  if (params.fMCMode != kMCTruth) {
2008  params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(0));
2009  if (!params.fTrackContainer) params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(1));
2010 
2011  params.fClusterContainer = GetClusterContainer(0);
2012 
2013  if (!params.fTrackContainer && !params.fClusterContainer) {
2014  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2015  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
2016  params.GetName());
2017  params.fInhibit = kTRUE;
2018  }
2019  }
2020  }
2021 }
2022 
2027 {
2028  if (!fAodEvent) return kFALSE;
2029 
2030  TString hname;
2031 
2032  // fix for temporary bug in ESDfilter
2033  // the AODs with null vertex pointer didn't pass the PhysSel
2034  if (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001) return kFALSE;
2035 
2036  for (auto &eng : fAnalysisEngines) {
2037 
2038  if (eng.fInhibit) continue;
2039 
2040  //Event selection
2041  hname = TString::Format("%s/fHistNEvents", eng.GetName());
2042  Bool_t iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
2043  if (!iseventselected) {
2044  fHistManager.FillTH1(hname, "Rejected");
2045  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2046  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
2047  TString label;
2048  do {
2049  label = GetHFEventRejectionReasonLabel(bitmap);
2050  if (label.IsNull()) break;
2051  fHistManager.FillTH1(hname, label);
2052  } while (true);
2053 
2054  continue;
2055  }
2056 
2057  fHistManager.FillTH1(hname, "Accepted");
2058 
2059  AliDebug(2, "Event selected");
2060 
2061  eng.RunAnalysis();
2062  }
2063  return kTRUE;
2064 }
2065 
2070 {
2071  TString hname;
2072  for (auto &param : fAnalysisEngines) {
2073 
2074  if (param.fInhibit) continue;
2075 
2076  if (fTreeOutput) {
2077  param.FillTree(fApplyKinematicCuts);
2078  }
2079  else {
2080  param.FillHnSparse(fApplyKinematicCuts);
2081  }
2082  }
2083  return kTRUE;
2084 }
2085 
2093 void AliAnalysisTaskDmesonJets::CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t& minMass, Double_t& maxMass)
2094 {
2095  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
2096 
2097  Double_t mass = part->Mass();
2098 
2099  // To make sure that the PDG mass value is not at the edge of a bin
2100  if (nbins % 2 == 0) {
2101  minMass = mass - range / 2 - range / nbins / 2;
2102  maxMass = mass + range / 2 - range / nbins / 2;
2103  }
2104  else {
2105  minMass = mass - range / 2;
2106  maxMass = mass + range / 2;
2107  }
2108 }
2109 
2116 {
2117  static TString label;
2118  label = "";
2119 
2120  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
2121  label = "NotSelTrigger";
2122  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
2123  return label.Data();
2124  }
2125  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
2126  label = "NoVertex";
2127  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
2128  return label.Data();
2129  }
2130  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
2131  label = "TooFewVtxContrib";
2132  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
2133  return label.Data();
2134  }
2135  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
2136  label = "ZVtxOutFid";
2137  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
2138  return label.Data();
2139  }
2140  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
2141  label = "Pileup";
2142  bitmap &= ~BIT(AliRDHFCuts::kPileup);
2143  return label.Data();
2144  }
2145  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
2146  label = "OutsideCentrality";
2147  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
2148  return label.Data();
2149  }
2150  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
2151  label = "PhysicsSelection";
2152  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
2153  return label.Data();
2154  }
2155  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
2156  label = "BadSPDVertex";
2157  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
2158  return label.Data();
2159  }
2160  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
2161  label = "ZVtxSPDOutFid";
2162  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
2163  return label.Data();
2164  }
2165  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
2166  label = "CentralityFlattening";
2167  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
2168  return label.Data();
2169  }
2170  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
2171  label = "BadTrackV0Correl";
2172  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
2173  return label.Data();
2174  }
2175 
2176  return label.Data();
2177 }
2178 
2185 {
2186  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
2187  PostData(eng.GetDataSlotNumber(), eng.GetTree());
2188  return eng.GetDataSlotNumber();
2189  }
2190  else {
2191  return -1;
2192  }
2193 }
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.
Bool_t fTreeOutput
If true, output will be posted in a TTree rather than a THnSparse.
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
virtual void UserCreateOutputObjects()
Creates the output containers.
Lightweight class that encapsulates D meson jets.
Double_t fInvMass2Prong
! 2-prong mass of the D* candidate (w/o the soft pion)
static void CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t &minMass, Double_t &maxMass)
AliAODEvent * fAodEvent
! AOD event
Bool_t FindJet(AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, AliHFJetDefinition &jetDef)
UInt_t fEnabledAxis
Use bit defined in EAxis_t to enable axis in the THnSparse.
EMCMode_t fMCMode
MC mode: No MC (data and MC detector level), background-only (MC), signal-only (MC), MC truth (particle level)
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:96
AnalysisEngine()
This is the default constructor, used for ROOT I/O purposes.
void Init(const AliEMCALGeometry *const geom, Int_t runNumber)
Initialize the analysis engine.
static FJRecoScheme ConvertToFJRecoScheme(ERecoScheme_t reco)
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.
static EMesonOrigin_t CheckOrigin(AliAODMCParticle *part, TClonesArray *mcArray)
Double_t GetDistance(std::string n, Double_t &deta, Double_t &dphi) const
THashList * CreateHistoGroup(const char *groupname, const char *parent="/")
virtual void Set(const AliDmesonJetInfo &source)
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:33
EJetType_t fJetType
Jet type (charged, full, neutral)
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
TObject * FindObject(const char *name) const
THashList * GetListOfHistograms() const
Definition: THistManager.h:504
AliAnalysisTaskDmesonJets()
This is the default constructor, used for ROOT I/O purposes.
virtual void Reset()
Reset the current object.
AnalysisEngine & operator=(const AnalysisEngine &source)
Struct that encapsulates analysis parameters.
std::vector< AliHFJetDefinition >::iterator FindJetDefinition(const AliHFJetDefinition &eng)
Base task in the EMCAL framework (lighter version of AliAnalysisTaskEmcal)
Double_t Phi_0_2pi() const
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:95
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
AliClusterContainer * GetClusterContainer(Int_t i=0) const
AliAODTrack * GetBachelor() const
AliRDHFCuts * LoadDMesonCutsFromFile(TString cutfname, TString cutsname)
bool operator<(const AliAnalysisTaskDmesonJets::AliHFJetDefinition &lhs, const AliAnalysisTaskDmesonJets::AliHFJetDefinition &rhs)
std::list< AnalysisEngine >::iterator FindAnalysisEngine(const AnalysisEngine &eng)
std::vector< AliHFJetDefinition > fJetDefinitions
Jet definitions.
static EMesonDecayChannel_t CheckDecayChannel(AliAODMCParticle *part, TClonesArray *mcArray)
EJetAcceptanceType_t fAcceptance
Jet acceptance.
EJetAlgo_t fJetAlgo
Jet algorithm (kt, anti-kt,...)
static const char * GetHFEventRejectionReasonLabel(UInt_t &bitmap)
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)
Bool_t ExtractParticleLevelHFAttributes(const AliAODMCParticle *part, AliDmesonJetInfo &DmesonJet)
Select MC particles based on specific prescriptions of HF analysis.
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
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.
void SetDetectorJetEtaPhiRange(const AliEMCALGeometry *const geom, Int_t run)
Bool_t ExtractRecoDecayAttributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
Double_t InvMassD0() const
virtual void AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
bool operator==(const AliAnalysisTaskDmesonJets::AliHFJetDefinition &lhs, const AliAnalysisTaskDmesonJets::AliHFJetDefinition &rhs)
Double_t InvMassDstarKpipi() const
AnalysisEngine * AddAnalysisEngine(ECandidateType_t type, EMCMode_t bkgMode, EJetType_t jettype, Double_t jetradius, TString cutfname="")
const Int_t nbins
Double_t maxMass
Int_t fNOutputTrees
Maximum number of output trees.
AliHFJetDefinition()
This is the default constructor, used for ROOT I/O purposes.
AliHFJetDefinition & operator=(const AliHFJetDefinition &source)
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:97
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)
void AddInputVectors(AliEmcalContainer *cont, Int_t offset, TH2 *rejectHist)
Lightweight class that encapsulates D0.
const Int_t nPtBins
Bool_t IsAnyJetInAcceptance(const AliDmesonJetInfo &dMesonJet) const
Container for jet within the EMCAL jet framework.
Int_t PostDataFromAnalysisEngine(const AnalysisEngine &eng)
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)