AliPhysics  e59a9ba (e59a9ba)
 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 
55 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonJetInfo
56 
60 
63 {
64  fD.SetPtEtaPhiE(0,0,0,0);
65  fSoftPionPt = 0;
66  fInvMass2Prong = 0;
67  fDmesonParticle = 0;
68  fMCLabel = -1;
69  fReconstructed = kFALSE;
70  for (auto &jet : fJets) {
71  jet.second.fMomentum.SetPtEtaPhiE(0,0,0,0);
72  jet.second.fNConstituents = 0;
73  jet.second.fNEF = 0;
74  jet.second.fMaxChargedPt = 0;
75  jet.second.fMaxNeutralPt = 0;
76  }
77 }
78 
81 {
82  Printf("Printing D Meson Jet object.");
83  Printf("D Meson: pT = %.3f, eta = %.3f, phi = %.3f, inv. mass = %.3f", fD.Pt(), fD.Eta(), fD.Phi_0_2pi(), fD.M());
84  Printf("Soft pion pT: %.3f. 2-Prong Invariant mass = %.3f", fSoftPionPt, fInvMass2Prong);
85  for (auto &jet : fJets) {
86  Printf("Jet %s: pT = %.3f, eta = %.3f, phi = %.3f", jet.first.c_str(), jet.second.Pt(), jet.second.Eta(), jet.second.Phi_0_2pi());
87  Printf("Jet N Consituents = %d", jet.second.fNConstituents);
88  }
89 }
90 
95 {
96  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
97  if (it == fJets.end()) return 0;
98 
99  Double_t z = 0;
100 
101  if ((*it).second.Pt() > 0) {
102  TVector3 dvect = fD.Vect();
103  TVector3 jvect = (*it).second.fMomentum.Vect();
104 
105  Double_t jetMom = jvect * jvect;
106 
107  if (jetMom < 1e-6) {
108  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetZ", "Zero jet momentum!");
109  z = 0.999;
110  }
111  else {
112  z = (dvect * jvect) / jetMom;
113  }
114 
115  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
116  }
117 
118  return z;
119 }
120 
126 Double_t AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetDistance(std::string n, Double_t& deta, Double_t& dphi) const
127 {
128  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
129  if (it == fJets.end()) return 0;
130 
131  dphi = TVector2::Phi_mpi_pi(fD.Phi() - (*it).second.Phi());;
132  deta = fD.Eta() - (*it).second.Eta();
133  return TMath::Sqrt(dphi*dphi + deta*deta);
134 }
135 
140 {
141  Double_t deta = 0;
142  Double_t dphi = 0;
143  return GetDistance(n, deta, dphi);
144 }
145 
151 {
152  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
153  if (it == fJets.end()) {
154  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
155  n.c_str());
156  return 0;
157  }
158  return &((*it).second);
159 }
160 
166 {
167  std::map<std::string, AliJetInfo>::iterator it = fJets.find(n);
168  if (it == fJets.end()) {
169  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
170  n.c_str());
171  return 0;
172  }
173  return &((*it).second);
174 }
175 
176 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoSummary
177 
181 
187  fPt(0),
188  fEta(0),
189  fPhi(0),
190  fR(0),
191  fZ(0)
192 {
193  Set(source, n);
194 }
195 
198 {
199  fPt = 0;
200  fEta = 0;
201  fPhi = 0;
202  fR = 0;
203  fZ = 0;
204 }
205 
211 {
212  std::map<std::string, AliJetInfo>::const_iterator it = source.fJets.find(n);
213  if (it == source.fJets.end()) return;
214 
215  fPt = (*it).second.Pt();
216  fEta = (*it).second.Eta();
217  fPhi = (*it).second.Phi_0_2pi();
218  fR = source.GetDistance(n);
219  fZ = source.GetZ(n);
220 }
221 
222 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonInfoSummary
223 
227 
232  fPt(0),
233  fEta(0),
234  fPhi(0)
235 {
236  Set(source);
237 }
238 
243 {
244  fPt = source.fD.Pt();
245  fEta = source.fD.Eta();
246  fPhi = source.fD.Phi_0_2pi();
247 }
248 
251 {
252  fPt = 0;
253  fEta = 0;
254  fPhi = 0;
255 }
256 
257 // Definitions of class AliAnalysisTaskDmesonJets::AliD0InfoSummary
258 
262 
267  AliDmesonInfoSummary(source),
268  fInvMass(source.fD.M())
269 {
270 }
271 
276 {
277  fInvMass = source.fD.M();
279 }
280 
283 {
285 
286  fInvMass = 0;
287 }
288 
289 // Definitions of class AliAnalysisTaskDmesonJets::AliDStarInfoSummary
290 
294 
299  AliDmesonInfoSummary(source),
300  f2ProngInvMass(source.fInvMass2Prong),
301  fDeltaInvMass(source.fD.M() - source.fInvMass2Prong)
302 {
303 }
304 
309 {
310  f2ProngInvMass = source.fInvMass2Prong;
311  fDeltaInvMass = source.fD.M() - source.fInvMass2Prong;
313 }
314 
317 {
319 
320  f2ProngInvMass = 0;
321  fDeltaInvMass = 0;
322 }
323 
324 // Definitions of class AliAnalysisTaskDmesonJets::AliJetDefinition
325 
329 
332  TObject(),
333  fJetType(AliJetContainer::kChargedJet),
334  fRadius(0),
335  fJetAlgo(AliJetContainer::antikt_algorithm),
336  fRecoScheme(AliJetContainer::pt_scheme),
337  fAcceptance(AliJetContainer::kUser),
338  fMinJetPt(0.),
339  fMaxJetPt(500.),
340  fMinJetPhi(0.),
341  fMaxJetPhi(0.),
342  fMinJetEta(0.),
343  fMaxJetEta(0.),
344  fMinChargedPt(0.),
345  fMaxChargedPt(0.),
346  fMinNeutralPt(0.),
347  fMaxNeutralPt(0.)
348 {
349 }
350 
358  TObject(),
359  fJetType(type),
360  fRadius(r),
361  fJetAlgo(algo),
362  fRecoScheme(reco),
363  fAcceptance(AliJetContainer::kUser),
364  fMinJetPt(0.),
365  fMaxJetPt(500.),
366  fMinJetPhi(0.),
367  fMaxJetPhi(0.),
368  fMinJetEta(0.),
369  fMaxJetEta(0.),
370  fMinChargedPt(0.),
371  fMaxChargedPt(0.),
372  fMinNeutralPt(0.),
373  fMaxNeutralPt(0.)
374 {
375  // By default set detector fiducial acceptance
376  switch (type) {
380  break;
383  break;
384  }
385 }
386 
391  TObject(),
392  fJetType(source.fJetType),
393  fRadius(source.fRadius),
394  fJetAlgo(source.fJetAlgo),
395  fRecoScheme(source.fRecoScheme),
396  fAcceptance(source.fAcceptance),
397  fMinJetPt(source.fMinJetPt),
398  fMaxJetPt(source.fMaxJetPt),
399  fMinJetPhi(source.fMinJetPhi),
400  fMaxJetPhi(source.fMaxJetPhi),
401  fMinJetEta(source.fMinJetEta),
402  fMaxJetEta(source.fMaxJetEta),
403  fMinChargedPt(source.fMinChargedPt),
404  fMaxChargedPt(source.fMaxChargedPt),
405  fMinNeutralPt(source.fMinNeutralPt),
406  fMaxNeutralPt(source.fMaxNeutralPt)
407 {
408 }
409 
414 {
415  new (this) AliHFJetDefinition(source);
416  return *this;
417 }
418 
421 {
422  static TString name;
423 
424  name = AliJetContainer::GenerateJetName(fJetType, fJetAlgo, fRecoScheme, fRadius, 0, 0, "Jet");
425 
426  return name.Data();
427 }
428 
434 {
435  if (fMinJetEta < fMaxJetEta && (jet.Eta() < fMinJetEta || jet.Eta() > fMaxJetEta)) return kFALSE;
436  if (fMinJetPhi < fMaxJetPhi && (jet.Phi() < fMinJetPhi || jet.Phi() > fMaxJetPhi)) return kFALSE;
437  if (jet.Pt() > fMaxJetPt || jet.Pt() < fMinJetPt) return kFALSE;
438  if (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt) return kFALSE;
439  if (jet.fMaxNeutralPt < fMinNeutralPt || jet.fMaxNeutralPt > fMaxNeutralPt) return kFALSE;
440 
441  return kTRUE;
442 }
443 
449 {
450  const AliJetInfo* jet = dMesonJet.GetJet(n);
451  if (!jet) return kFALSE;
452  return IsJetInAcceptance((*jet));
453 }
454 
459 void AliAnalysisTaskDmesonJets::AliHFJetDefinition::SetDetectorJetEtaPhiRange(const AliEMCALGeometry* const geom, Int_t run)
460 {
461  Double_t r = 0;
462  switch (fAcceptance) {
464  r = fRadius;
465  // enforce fiducial acceptance
466  /* no break */
468  SetJetEtaRange(-0.9 + r, 0.9 - r);
469  SetJetPhiRange(0, 0); // No cut on phi
470  break;
471 
473  r = fRadius;
474  // enforce fiducial acceptance
475  /* no break */
477  if (geom) {
478  SetJetEtaRange(geom->GetArm1EtaMin() + r, geom->GetArm1EtaMax() - r);
479 
480  if(run>=177295 && run<=197470) {//small SM masked in 2012 and 2013
481  SetJetPhiRange(1.405 + r,3.135 - r);
482  }
483  else {
484  SetJetPhiRange(geom->GetArm1PhiMin() * TMath::DegToRad() + r, geom->GetEMCALPhiMax() * TMath::DegToRad() - r);
485  }
486  }
487  else {
488  AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
489  SetJetEtaRange(-0.7 + r, 0.7 - r);
490  SetJetPhiRange(1.405 + r, 3.135 - r);
491  }
492  break;
493 
495  r = fRadius;
496  // enforce fiducial acceptance
497  /* no break */
499  if (geom) {
500  SetJetEtaRange(geom->GetArm1EtaMin() + r, geom->GetArm1EtaMax() - r);
501  SetJetPhiRange(geom->GetDCALPhiMin() * TMath::DegToRad() + r, geom->GetDCALPhiMax() * TMath::DegToRad() - r);
502  }
503  else {
504  AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for DCAL year 2015!!");
505  SetJetEtaRange(-0.7 + r, 0.7 - r);
506  SetJetPhiRange(4.538 + r, 5.727 - r);
507  }
508  break;
509 
511  // Nothing to be done
512  break;
513  }
514 }
515 
522 {
523  if (lhs.fJetType > rhs.fJetType) return false;
524  else if (lhs.fJetType < rhs.fJetType) return true;
525  else {
526  if (lhs.fRadius > rhs.fRadius) return false;
527  else if (lhs.fRadius < rhs.fRadius) return true;
528  else {
529  if (lhs.fJetAlgo > rhs.fJetAlgo) return false;
530  else if (lhs.fJetAlgo < rhs.fJetAlgo) return true;
531  else {
532  if (lhs.fRecoScheme < rhs.fRecoScheme) return true;
533  else return false;
534  }
535  }
536  }
537 }
538 
545 {
546  if (lhs.fJetType != rhs.fJetType) return false;
547  if (lhs.fRadius != rhs.fRadius) return false;
548  if (lhs.fJetAlgo != rhs.fJetAlgo) return false;
549  if (lhs.fRecoScheme != rhs.fRecoScheme) return false;
550  return true;
551 }
552 
553 // Definitions of class AliAnalysisTaskDmesonJets::AnalysisEngine
554 
558 
561  TObject(),
562  fCandidateType(kD0toKpi),
563  fCandidateName(),
564  fCandidatePDG(0),
565  fNDaughters(0),
566  fPDGdaughters(),
567  fBranchName(),
568  fMCMode(kNoMC),
569  fNMassBins(0),
570  fMinMass(0),
571  fMaxMass(0),
572  fRDHFCuts(0),
573  fRejectedOrigin(0),
574  fAcceptedDecay(0),
575  fInhibit(kFALSE),
576  fJetDefinitions(),
577  fPtBinWidth(0.5),
578  fMaxPt(100),
579  fDataSlotNumber(-1),
580  fTree(0),
581  fCurrentDmesonJetInfo(0),
582  fCurrentJetInfo(0),
583  fCandidateArray(0),
584  fMCContainer(0),
585  fTrackContainer(0),
586  fClusterContainer(0),
587  fAodEvent(0),
588  fFastJetWrapper(0),
589  fHistManager(0)
590 {
591 }
592 
601  TObject(),
602  fCandidateType(type),
603  fCandidateName(),
604  fCandidatePDG(0),
605  fNDaughters(0),
606  fPDGdaughters(),
607  fBranchName(),
608  fMCMode(MCmode),
609  fNMassBins(nMassBins),
610  fMinMass(0),
611  fMaxMass(0),
612  fRDHFCuts(cuts),
613  fRejectedOrigin(kUnknownQuark | kFromBottom),
614  fAcceptedDecay(kAnyDecay),
615  fInhibit(kFALSE),
616  fJetDefinitions(),
617  fPtBinWidth(0.5),
618  fMaxPt(100),
619  fDataSlotNumber(-1),
620  fTree(0),
621  fCurrentDmesonJetInfo(0),
622  fCurrentJetInfo(0),
623  fCandidateArray(0),
624  fMCContainer(0),
625  fTrackContainer(0),
626  fClusterContainer(0),
627  fAodEvent(0),
628  fFastJetWrapper(0),
629  fHistManager(0)
630 {
631  SetCandidateProperties(range);
632 }
633 
638  TObject(source),
639  fCandidateType(source.fCandidateType),
640  fCandidateName(source.fCandidateName),
641  fCandidatePDG(source.fCandidatePDG),
642  fNDaughters(source.fNDaughters),
643  fPDGdaughters(source.fPDGdaughters),
644  fBranchName(source.fBranchName),
645  fMCMode(source.fMCMode),
646  fNMassBins(source.fNMassBins),
647  fMinMass(source.fMinMass),
648  fMaxMass(source.fMaxMass),
649  fRDHFCuts(),
650  fRejectedOrigin(source.fRejectedOrigin),
651  fAcceptedDecay(source.fAcceptedDecay),
652  fInhibit(source.fInhibit),
653  fJetDefinitions(source.fJetDefinitions),
654  fPtBinWidth(source.fPtBinWidth),
655  fMaxPt(source.fMaxPt),
656  fDataSlotNumber(-1),
657  fTree(0),
658  fCurrentDmesonJetInfo(0),
659  fCurrentJetInfo(0),
660  fCandidateArray(source.fCandidateArray),
661  fMCContainer(source.fMCContainer),
662  fTrackContainer(source.fTrackContainer),
663  fClusterContainer(source.fClusterContainer),
664  fAodEvent(source.fAodEvent),
666  fHistManager(source.fHistManager)
667 {
668  SetRDHFCuts(source.fRDHFCuts);
669 }
670 
671 // Destructor
673 {
674  delete fRDHFCuts;
675 }
676 
681 {
682  new (this) AnalysisEngine(source);
683  return *this;
684 }
685 
690 {
691  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
692  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
693  }
694 
695  return kFALSE;
696 }
697 
699 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const geom, Int_t runNumber)
700 {
701  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
702  fJetDefinitions[i].SetDetectorJetEtaPhiRange(geom, runNumber);
703  }
704 }
705 
710 {
711  switch (fCandidateType) {
712  case kD0toKpi:
713  fCandidatePDG = 421;
714  fCandidateName = "D0";
715  fNDaughters = 2;
716  fPDGdaughters.Set(fNDaughters);
717  fPDGdaughters.Reset();
718  fPDGdaughters[0] = 211; // pi
719  fPDGdaughters[1] = 321; // K
720  fBranchName = "D0toKpi";
721  fAcceptedDecay = kD0toKpi;
722  if (!fRDHFCuts) {
723  fRDHFCuts = new AliRDHFCutsD0toKpi();
724  fRDHFCuts->SetStandardCutsPP2010();
725  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
726  fRDHFCuts->SetTriggerClass("","");
727  }
728  break;
729  case kDstartoKpipi:
730  fCandidatePDG = 413;
731  fCandidateName = "DStar";
732  fNDaughters = 3;
733  fPDGdaughters.Set(fNDaughters);
734  fPDGdaughters.Reset();
735  fPDGdaughters[0] = 211; // pi soft
736  fPDGdaughters[1] = 211; // pi fromD0
737  fPDGdaughters[2] = 321; // K from D0
738  fBranchName = "Dstar";
739  fAcceptedDecay = kDstartoKpipi;
740  if (!fRDHFCuts) {
741  fRDHFCuts = new AliRDHFCutsDStartoKpipi();
742  fRDHFCuts->SetStandardCutsPP2010();
743  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
744  fRDHFCuts->SetTriggerClass("","");
745  }
746  break;
747  default:
748  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
749  }
750 
751  CalculateMassLimits(range, fCandidatePDG, fNMassBins, fMinMass, fMaxMass);
752 }
753 
758 {
759  if (fRDHFCuts) delete fRDHFCuts;
760  fRDHFCuts = cuts;
761 }
762 
767 {
768  if (!cuts) return;
769  if (fRDHFCuts) delete fRDHFCuts;
770  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
771 }
772 
777 {
778  static TString name;
779 
780  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
781 
782  return name.Data();
783 }
784 
789 {
790  static TString name;
791 
792  name = fCandidateName;
793  switch (fMCMode) {
794  case kBackgroundOnly:
795  name += "_kBackgroundOnly";
796  break;
797  case kSignalOnly:
798  name += "_kSignalOnly";
799  break;
800  case kMCTruth:
801  name += "_MCTruth";
802  break;
803  default:
804  break;
805  }
806 
807  return name.Data();
808 }
809 
817 {
818  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
819 
820  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
821  fJetDefinitions.push_back(def);
822  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
823  "Total number of jet definitions is now %lu.",
824  def.GetName(), GetName(), fJetDefinitions.size());
825  // For detector level set maximum track pt to 100 GeV/c
826  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
827  }
828  else {
829  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
830  }
831 
832  return &(*it);
833 }
834 
846 {
847  AliHFJetDefinition def(type, r, algo, reco);
848 
849  return AddJetDefinition(def);
850 }
851 
857 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
858 {
859  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
860  while (it != fJetDefinitions.end() && (*it) != def) it++;
861  return it;
862 }
863 
870 {
871  if (lhs.fCandidateType > rhs.fCandidateType) return false;
872  else if (lhs.fCandidateType < rhs.fCandidateType) return true;
873  else {
874  if (lhs.fMCMode < rhs.fMCMode) return true;
875  else return false;
876  }
877 }
878 
885 {
886  if (lhs.fCandidateType != rhs.fCandidateType) return false;
887  if (lhs.fMCMode != rhs.fMCMode) return false;
888  return true;
889 }
890 
899 {
900  DmesonJet.fD.SetPtEtaPhiM(part->Pt(), part->Eta(), part->Phi(), part->M());
901  return kTRUE;
902 }
903 
912 {
913  if (fCandidateType == kD0toKpi) { // D0 candidate
914  return ExtractD0Attributes(Dcand, DmesonJet, i);
915  }
916  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
917  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
918  }
919  else {
920  return kFALSE;
921  }
922 }
923 
932 {
933  Int_t MCtruthPdgCode = 0;
934 
935  Double_t invMassD = 0;
936 
937  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
938  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
939  DmesonJet.fMCLabel = mcLab;
940  if (mcLab >= 0) {
941  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
942 
943  if (aodMcPart) {
944  if (fRejectedOrigin && fMCMode == kSignalOnly) {
945  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
946 
947  if ((origin & fRejectedOrigin) == origin) return kFALSE;
948  }
949  MCtruthPdgCode = aodMcPart->PdgCode();
950  }
951  }
952  }
953 
954  //AliDebug(2,"Checking if D0 meson is selected");
955  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
956  if (isSelected == 1) { // selected as a D0
957  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
958 
959  if (fMCMode == kNoMC ||
960  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
961  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly)) {
962  // 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)
963  //AliDebug(2,"Selected as D0");
964  invMassD = Dcand->InvMassD0();
965  }
966  else { // conditions above not passed, so return FALSE
967  return kFALSE;
968  }
969  }
970  else if (isSelected == 2) { // selected as a D0bar
971  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
972 
973  if (fMCMode == kNoMC ||
974  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
975  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly)) {
976  // 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)
977  //AliDebug(2,"Selected as D0bar");
978  invMassD = Dcand->InvMassD0bar();
979  }
980  else { // conditions above not passed, so return FALSE
981  return kFALSE;
982  }
983  }
984  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
985  //AliDebug(2,"Selected as either D0 or D0bar");
986 
987  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
988  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
989  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kBackgroundOnly)) {
990  if (i > 0) return kFALSE;
991  //AliDebug(2, "MC truth is D0");
992  invMassD = Dcand->InvMassD0();
993  }
994  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
995  (MCtruthPdgCode == fCandidatePDG && fMCMode == kBackgroundOnly)) {
996  if (i > 0) return kFALSE;
997  //AliDebug(2, "MC truth is D0bar");
998  invMassD = Dcand->InvMassD0bar();
999  }
1000  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
1001 
1002  // Only accept it if background-only OR background-and-signal was requested
1003  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
1004  // Select D0 or D0bar depending on the i-parameter
1005  if (i == 0) {
1006  //AliDebug(2, "Returning invariant mass with D0 hypothesis");
1007  invMassD = Dcand->InvMassD0();
1008  }
1009  else if (i == 1) {
1010  //AliDebug(2, "Returning invariant mass with D0bar hypothesis");
1011  invMassD = Dcand->InvMassD0bar();
1012  }
1013  else { // i > 1
1014  return kFALSE;
1015  }
1016  }
1017  else { // signal-only was requested but this is not a true D0
1018  return kFALSE;
1019  }
1020  }
1021  }
1022  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
1023  return kTRUE;
1024 }
1025 
1034 {
1035  if (i > 0) return kFALSE; // only one mass hypothesis for the D*
1036 
1037  Int_t MCtruthPdgCode = 0;
1038 
1039  Double_t invMassD = 0;
1040 
1041  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1042  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1043  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1044 
1045  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1046  //AliDebug(2, Form("MC label is %d", mcLab));
1047  DmesonJet.fMCLabel = mcLab;
1048  if (mcLab >= 0) {
1049  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1050 
1051  if (aodMcPart) {
1052  if (fRejectedOrigin && fMCMode == kSignalOnly) {
1053  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
1054 
1055  if ((origin & fRejectedOrigin) == origin) return kFALSE;
1056  }
1057 
1058  MCtruthPdgCode = aodMcPart->PdgCode();
1059  //AliDebug(2, Form("MC truth pdg code is %d",MCtruthPdgCode));
1060  }
1061  }
1062  }
1063 
1064  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1065  if (fMCMode == kNoMC ||
1066  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1067  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1068  // 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)
1069  invMassD = DstarCand->InvMassDstarKpipi();
1070  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1071  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1072  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1073  return kTRUE;
1074  }
1075  else { // conditions above not passed, so return FALSE
1076  return kFALSE;
1077  }
1078 }
1079 
1087 {
1088  if (!part) return kDecayOther;
1089  if (!mcArray) return kDecayOther;
1090 
1092 
1093  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1094 
1095  if (part->GetNDaughters() == 2) {
1096 
1097  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1098  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1099 
1100  if (!d1 || !d2) {
1101  return decay;
1102  }
1103 
1104  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1105  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1106 
1107  if (absPdgPart == 421) { // D0 -> K pi
1108 
1109  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1110  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1111  decay = kDecayD0toKpi;
1112  }
1113  }
1114 
1115  if (absPdgPart == 413) { // D* -> D0 pi
1116 
1117  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1118  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1119  if (D0decay == kDecayD0toKpi) {
1120  decay = kDecayDStartoKpipi;
1121  }
1122  }
1123 
1124  if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1125  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1126  if (D0decay == kDecayD0toKpi) {
1127  decay = kDecayDStartoKpipi;
1128  }
1129  }
1130  }
1131  }
1132 
1133  return decay;
1134 }
1135 
1143 {
1144  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1145 
1146  if (!part) return kUnknownQuark;
1147  if (!mcArray) return kUnknownQuark;
1148 
1149  Int_t pdgGranma = 0;
1150  Int_t mother = part->GetMother();
1151  Int_t istep = 0;
1152  Int_t abspdgGranma = 0;
1153  Bool_t isFromB = kFALSE;
1154  Bool_t isQuarkFound = kFALSE;
1155 
1156  while (mother >= 0) {
1157  istep++;
1158  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1159  if (mcGranma >= 0) {
1160  pdgGranma = mcGranma->GetPdgCode();
1161  abspdgGranma = TMath::Abs(pdgGranma);
1162  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1163  isFromB = kTRUE;
1164  }
1165 
1166  if (abspdgGranma == 4 || abspdgGranma == 5) isQuarkFound = kTRUE;
1167  mother = mcGranma->GetMother();
1168  }
1169  else {
1170  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1171  break;
1172  }
1173  }
1174 
1175  if (isQuarkFound) {
1176  if (isFromB) {
1177  return kFromBottom;
1178  }
1179  else {
1180  return kFromCharm;
1181  }
1182  }
1183  else {
1184  return kUnknownQuark;
1185  }
1186 }
1187 
1190 {
1191  fDmesonJets.clear();
1192 
1193  if (fMCMode == kMCTruth) {
1194  RunParticleLevelAnalysis();
1195  }
1196  else {
1197  RunDetectorLevelAnalysis();
1198  }
1199 }
1200 
1203 {
1204  const Int_t nD = fCandidateArray->GetEntriesFast();
1205 
1206  AliDmesonJetInfo DmesonJet;
1207 
1208  Int_t nAccCharm = 0;
1209  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1210  Int_t isSelected = 0;
1211 
1212  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1213  if (!charmCand) continue;
1214 
1215  Int_t nprongs = charmCand->GetNProngs();
1216 
1217  if (fCandidateType == kDstartoKpipi) {
1218  if (!charmCand->InheritsFrom("AliAODRecoCascadeHF")) {
1219  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::RunDetectorLevelAnalysis","Candidate type is D* but object type is wrong (should be AliAODRecoCascadeHF)");
1220  continue;
1221  }
1222  }
1223 
1224  // region of interest + cuts
1225  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1226 
1227  //candidate selected by cuts and PID
1228  isSelected = fRDHFCuts->IsSelected(charmCand, AliRDHFCuts::kAll, fAodEvent); //selected
1229 
1230  if (!isSelected) continue;
1231 
1232  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1233  DmesonJet.Reset();
1234  DmesonJet.fDmesonParticle = charmCand;
1235  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1236  for (std::vector<AliHFJetDefinition>::iterator itdef = fJetDefinitions.begin(); itdef != fJetDefinitions.end(); itdef++) {
1237  if (!FindJet(charmCand, DmesonJet, *itdef)) {
1238  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1239  (*itdef).GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1240  }
1241  }
1242  fDmesonJets[icharm] = DmesonJet;
1243  }
1244  }
1245  nAccCharm++;
1246  } // end of D cand loop
1247 
1248  TString hname;
1249 
1250  hname = TString::Format("%s/fHistNAcceptedDmesons", GetName());
1251  fHistManager->FillTH1(hname, nAccCharm);
1252 
1253  hname = TString::Format("%s/fHistNDmesons", GetName());
1254  fHistManager->FillTH1(hname, nD);
1255 }
1256 
1267 {
1268  TString hname;
1269 
1271  fFastJetWrapper->SetR(jetDef.fRadius);
1274 
1275  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1276 
1277  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1278  fTrackContainer->SetDMesonCandidate(Dcand);
1279  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1280  AddInputVectors(fTrackContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1281 
1282  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1283  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1284  const TObjArray& daughters = fTrackContainer->GetDaughterList();
1285  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1286  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1287  if (!fTrackContainer->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1288  }
1289  }
1290 
1291  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1292  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1293  AddInputVectors(fClusterContainer, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1294  }
1295 
1296  // run jet finder
1297  fFastJetWrapper->Run();
1298 
1299  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1300 
1301  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1302  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1303 
1304  Bool_t isDmesonJet = kFALSE;
1305 
1306  Double_t maxChPt = 0;
1307  Double_t maxNePt = 0;
1308  Double_t totalNeutralPt = 0;
1309 
1310  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1311  if (constituents[ic].user_index() == 0) {
1312  isDmesonJet = kTRUE;
1313  }
1314  else if (constituents[ic].user_index() >= 100) {
1315  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1316  }
1317  else if (constituents[ic].user_index() <= -100) {
1318  totalNeutralPt += constituents[ic].pt();
1319  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1320  }
1321  }
1322 
1323  if (isDmesonJet) {
1324  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1325  DmesonJet.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1326  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1327  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1328  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1329 
1330  return kTRUE;
1331  }
1332  }
1333 
1334  return kFALSE;
1335 }
1336 
1340 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist)
1341 {
1342  AliEmcalIterableMomentumContainer itcont = cont->all_momentum();
1343  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1344  UInt_t rejectionReason = 0;
1345  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1346  rejectHist->Fill(cont->GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1347  continue;
1348  }
1349  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1350  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1351  }
1352 }
1353 
1356 {
1357  TString hname;
1358 
1359  fMCContainer->SetSpecialPDG(fCandidatePDG);
1360  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1361  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1362 
1363  if (!fMCContainer->IsSpecialPDGFound()) return;
1364 
1365  for (auto &jetDef : fJetDefinitions) {
1366 
1368  fFastJetWrapper->SetR(jetDef.fRadius);
1371 
1372  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1373  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1374 
1375  fFastJetWrapper->Run();
1376 
1377  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1378 
1379  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1380  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1381 
1382  Bool_t isDmesonJet = kFALSE;
1383 
1384  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1385  Int_t iPart = constituents[ic].user_index() - 100;
1386  AliVParticle* part = fMCContainer->GetParticle(iPart);
1387  if (!part) {
1388  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1389  continue;
1390  }
1391  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1392  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1393  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1394  std::pair<int, AliDmesonJetInfo> element;
1395  element.first = iPart;
1396  element.second = AliDmesonJetInfo();
1397  dMesonJetIt = fDmesonJets.insert(element).first;
1398  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1399  (*dMesonJetIt).second.fDmesonParticle = part;
1400  }
1401 
1402  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1403  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1404  }
1405  }
1406  }
1407  }
1408 }
1409 
1414 {
1415  TString classname;
1416  switch (fCandidateType) {
1417  case kD0toKpi:
1418  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1419  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1420  break;
1421  case kDstartoKpipi:
1422  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1423  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1424  break;
1425  }
1426  TString treeName = TString::Format("%s_%s", taskName, GetName());
1427  fTree = new TTree(treeName, treeName);
1428  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
1429  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1430  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1431  fCurrentJetInfo[i] = new AliJetInfoSummary();
1432  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1433  }
1434 
1435  return fTree;
1436 }
1437 
1442 {
1443  TString hname;
1444 
1445  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
1446 
1447  for (auto &jetDef : fJetDefinitions) {
1448 
1449  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
1450 
1451  Double_t radius = jetDef.fRadius;
1452 
1453  TString title[30] = {""};
1454  Int_t nbins[30] = {0 };
1455  Double_t min [30] = {0.};
1456  Double_t max [30] = {0.};
1457  Int_t dim = 0 ;
1458 
1459  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1460  nbins[dim] = nPtBins;
1461  min[dim] = 0;
1462  max[dim] = fMaxPt;
1463  dim++;
1464 
1465  if ((enabledAxis & kPositionD) != 0) {
1466  title[dim] = "#eta_{D}";
1467  nbins[dim] = 50;
1468  min[dim] = -1;
1469  max[dim] = 1;
1470  dim++;
1471 
1472  title[dim] = "#phi_{D} (rad)";
1473  nbins[dim] = 150;
1474  min[dim] = 0;
1475  max[dim] = TMath::TwoPi();
1476  dim++;
1477  }
1478 
1479  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1480  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1481  nbins[dim] = fNMassBins;
1482  min[dim] = fMinMass;
1483  max[dim] = fMaxMass;
1484  dim++;
1485  }
1486 
1487  if (fCandidateType == kD0toKpi) {
1488  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1489  nbins[dim] = fNMassBins;
1490  min[dim] = fMinMass;
1491  max[dim] = fMaxMass;
1492  dim++;
1493  }
1494 
1495  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1496  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1497  nbins[dim] = fNMassBins;
1498  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1499  dim++;
1500  }
1501 
1502  if (fCandidateType == kDstartoKpipi) {
1503  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1504  nbins[dim] = fNMassBins*6;
1505  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1506 
1507  // subtract mass of D0
1508  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1509  min[dim] -= D0mass;
1510  max[dim] -= D0mass;
1511 
1512  dim++;
1513  }
1514 
1515  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1516  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1517  nbins[dim] = 100;
1518  min[dim] = 0;
1519  max[dim] = 25;
1520  dim++;
1521  }
1522 
1523  title[dim] = "#it{z}_{D}";
1524  nbins[dim] = 110;
1525  min[dim] = 0;
1526  max[dim] = 1.10;
1527  dim++;
1528 
1529  if ((enabledAxis & kDeltaR) != 0) {
1530  title[dim] = "#Delta R_{D-jet}";
1531  nbins[dim] = 100;
1532  min[dim] = 0;
1533  max[dim] = radius * 1.5;
1534  dim++;
1535  }
1536 
1537  if ((enabledAxis & kDeltaEta) != 0) {
1538  title[dim] = "#eta_{D} - #eta_{jet}";
1539  nbins[dim] = 100;
1540  min[dim] = -radius * 1.2;
1541  max[dim] = radius * 1.2;
1542  dim++;
1543  }
1544 
1545  if ((enabledAxis & kDeltaPhi) != 0) {
1546  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1547  nbins[dim] = 100;
1548  min[dim] = -radius * 1.2;
1549  max[dim] = radius * 1.2;
1550  dim++;
1551  }
1552 
1553  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1554  nbins[dim] = nPtBins;
1555  min[dim] = 0;
1556  max[dim] = fMaxPt;
1557  dim++;
1558 
1559  if ((enabledAxis & kPositionJet) != 0) {
1560  title[dim] = "#eta_{jet}";
1561  nbins[dim] = 50;
1562  min[dim] = -1;
1563  max[dim] = 1;
1564  dim++;
1565 
1566  title[dim] = "#phi_{jet} (rad)";
1567  nbins[dim] = 150;
1568  min[dim] = 0;
1569  max[dim] = TMath::TwoPi();
1570  dim++;
1571  }
1572 
1573  if ((enabledAxis & kJetConstituents) != 0) {
1574  title[dim] = "No. of constituents";
1575  nbins[dim] = 50;
1576  min[dim] = -0.5;
1577  max[dim] = 49.5;
1578  dim++;
1579  }
1580 
1581  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1582  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1583  for (Int_t j = 0; j < dim; j++) {
1584  h->GetAxis(j)->SetTitle(title[j]);
1585  }
1586  }
1587 }
1588 
1593 {
1594  TString hname;
1595 
1596  for (auto& dmeson_pair : fDmesonJets) {
1597  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
1598  Int_t accJets = 0;
1599  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1600  fCurrentJetInfo[ij]->Reset();
1601  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1602  if (!jet) continue;
1603  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1604  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1605  fHistManager->FillTH1(hname, jet->Pt());
1606  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1607  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1608  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1609  fHistManager->FillTH1(hname, jet->Eta());
1610  continue;
1611  }
1612  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
1613  accJets++;
1614  }
1615  if (kTRUE || accJets > 0) { // always fill D meson tree, even if no jet was accepted
1616  fTree->Fill();
1617  }
1618  else { // this is never executed
1619  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1620  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1621  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1622  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1623  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1624  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1625  if (fCandidateType == kD0toKpi) {
1626  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1627  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1628  }
1629  else if (fCandidateType == kDstartoKpipi) {
1630  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1631  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1632 
1633  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1634  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1635  }
1636  }
1637  }
1638  return kTRUE;
1639 }
1640 
1646 {
1647  TString hname;
1648 
1649  for (auto& dmeson_pair : fDmesonJets) {
1650  Int_t accJets = 0;
1651  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1652  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1653  if (!jet) continue;
1654  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1655  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1656  fHistManager->FillTH1(hname, jet->Pt());
1657  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1658  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1659  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1660  fHistManager->FillTH1(hname, jet->Eta());
1661  continue;
1662  }
1663  accJets++;
1664  }
1665  if (accJets == 0) {
1666  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1667  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1668  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1669  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1670  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1671  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1672  if (fCandidateType == kD0toKpi) {
1673  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1674  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1675  }
1676  else if (fCandidateType == kDstartoKpipi) {
1677  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1678  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1679 
1680  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1681  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1682  }
1683  }
1684  }
1685  return kTRUE;
1686 }
1687 
1692 {
1693  TString hname;
1694 
1695  for (auto& dmeson_pair : fDmesonJets) {
1696  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
1697  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1698  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1699  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1700  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1701  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1702  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1703  }
1704  }
1705 
1706  for (auto &jetDef : fJetDefinitions) {
1707 
1708  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1709  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
1710 
1711  for (auto& dmeson_pair : fDmesonJets) {
1712  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
1713  if (!jet) continue;
1714  if (!jetDef.IsJetInAcceptance(*jet)) {
1715  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
1716  fHistManager->FillTH1(hname, jet->Pt());
1717  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
1718  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1719  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
1720  fHistManager->FillTH1(hname, jet->Eta());
1721  continue;
1722  }
1723  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
1724  }
1725  }
1726 
1727  return kTRUE;
1728 }
1729 
1735 Bool_t AliAnalysisTaskDmesonJets::AnalysisEngine::FillHnSparse(THnSparse* h, const AliDmesonJetInfo& DmesonJet, std::string n)
1736 {
1737  // Fill the THnSparse histogram.
1738 
1739  Double_t contents[30] = {0.};
1740 
1741  Double_t z = DmesonJet.GetZ(n);
1742  Double_t deltaPhi = 0;
1743  Double_t deltaEta = 0;
1744  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
1745 
1746  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
1747  if (it == DmesonJet.fJets.end()) return kFALSE;
1748 
1749  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
1750  TString title(h->GetAxis(i)->GetTitle());
1751  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
1752  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
1753  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
1754  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
1755  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
1756  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
1757  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
1758  else if (title=="#it{z}_{D}") contents[i] = z;
1759  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
1760  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
1761  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
1762  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
1763  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
1764  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
1765  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
1766  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
1767  }
1768 
1769  h->Fill(contents);
1770 
1771  return kTRUE;
1772 }
1773 
1774 // Definitions of class AliAnalysisTaskDmesonJets
1775 
1779 
1783  fAnalysisEngines(),
1784  fEnabledAxis(0),
1786  fHistManager(),
1787  fApplyKinematicCuts(kTRUE),
1788  fNOutputTrees(0),
1789  fAodEvent(0),
1790  fFastJetWrapper(0)
1791 {
1792  SetMakeGeneralHistograms(kTRUE);
1793 }
1794 
1798 AliAnalysisTaskDmesonJets::AliAnalysisTaskDmesonJets(const char* name, Int_t nOutputTrees) :
1799  AliAnalysisTaskEmcalLight(name, kTRUE),
1800  fAnalysisEngines(),
1801  fEnabledAxis(k2ProngInvMass),
1802  fOutputType(kTreeOutput),
1803  fHistManager(name),
1804  fApplyKinematicCuts(kTRUE),
1805  fNOutputTrees(nOutputTrees),
1806  fAodEvent(0),
1807  fFastJetWrapper(0)
1808 {
1809  SetMakeGeneralHistograms(kTRUE);
1810  for (Int_t i = 0; i < nOutputTrees; i++){
1811  DefineOutput(2+i, TTree::Class());
1812  }
1813 }
1814 
1817 {
1818  if (fFastJetWrapper) delete fFastJetWrapper;
1819 }
1820 
1828 {
1829  AliRDHFCuts* analysiscuts = 0;
1830  TFile* filecuts = TFile::Open(cutfname);
1831  if (!filecuts || filecuts->IsZombie()) {
1832  ::Warning("AddTaskDmesonJets", "Input file not found: will use std cuts.");
1833  filecuts = 0;
1834  }
1835 
1836  if (filecuts) {
1837  analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
1838  if (!analysiscuts) {
1839  ::Warning("AddTaskDmesonJetCorr", "Could not find analysis cuts '%s' in '%s'. Using std cuts.", cutsname.Data(), cutfname.Data());
1840  }
1841  }
1842 
1843  return analysiscuts;
1844 }
1845 
1855 {
1857  return AddAnalysisEngine(type, MCmode, jetDef, cutfname);
1858 }
1859 
1869 {
1870  AliRDHFCuts* cuts = 0;
1871 
1872  if (!cutfname.IsNull()) {
1873  TString cutsname;
1874 
1875  switch (type) {
1876  case kD0toKpi :
1877  cutsname = "D0toKpiCuts";
1878  break;
1879  case kDstartoKpipi :
1880  cutsname = "DStartoKpipiCuts";
1881  break;
1882  default:
1883  return 0;
1884  }
1885 
1886  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
1887  }
1888 
1889  AnalysisEngine eng(type, MCmode, cuts);
1890 
1891  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
1892 
1893  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
1894  eng.AddJetDefinition(jetDef);
1895  it = fAnalysisEngines.insert(it, eng);
1896  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(jetDef), fAnalysisEngines.size());
1897  }
1898  else {
1899  AnalysisEngine* found_eng = &(*it);
1900  ::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());
1901  found_eng->AddJetDefinition(jetDef);
1902 
1903  if (cuts && found_eng->fRDHFCuts != 0) {
1904  ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
1905  found_eng->SetRDHFCuts(cuts);
1906  }
1907  }
1908 
1909  return &(*it);
1910 }
1911 
1912 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
1913 {
1914  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
1915  while (it != fAnalysisEngines.end() && (*it) < eng) it++;
1916  return it;
1917 }
1918 
1921 {
1922  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
1923 
1925 
1926  // Define histograms
1927  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
1928 
1929  TString hname;
1930  TString htitle;
1931  TH1* h = 0;
1932  Int_t treeSlot = 0;
1933 
1934  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
1935  for (auto &param : fAnalysisEngines) {
1936  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
1937 
1938  fHistManager.CreateHistoGroup(param.GetName());
1939 
1940  param.fHistManager = &fHistManager;
1941 
1942  hname = TString::Format("%s/fHistNAcceptedDmesons", param.GetName());
1943  htitle = hname + ";Number of D accepted meson candidates;counts";
1944  h = fHistManager.CreateTH1(hname, htitle, 51, -0.5, 50.5);
1945 
1946  hname = TString::Format("%s/fHistNDmesons", param.GetName());
1947  htitle = hname + ";Number of D meson candidates;counts";
1948  h = fHistManager.CreateTH1(hname, htitle, 101, -0.5, 100.5);
1949 
1950  hname = TString::Format("%s/fHistNEvents", param.GetName());
1951  htitle = hname + ";Event status;counts";
1952  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
1953  h->GetXaxis()->SetBinLabel(1, "Accepted");
1954  h->GetXaxis()->SetBinLabel(2, "Rejected");
1955 
1956  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
1957  htitle = hname + ";Rejection reason;counts";
1958  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
1959 
1960  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
1961  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
1962  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
1963 
1964  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
1965  htitle = hname + ";#it{#eta}_{D};counts";
1966  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
1967 
1968  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
1969  htitle = hname + ";#it{#phi}_{D};counts";
1970  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
1971 
1972  if (param.fCandidateType == kD0toKpi) {
1973  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
1974  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
1975  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
1976  }
1977  else if (param.fCandidateType == kDstartoKpipi) {
1978  Double_t min = 0;
1979  Double_t max = 0;
1980 
1981  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
1982  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
1983  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
1984  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
1985 
1986  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1987  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
1988  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
1989  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
1990  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
1991  }
1992 
1993  for (std::vector<AliHFJetDefinition>::iterator itdef = param.fJetDefinitions.begin(); itdef != param.fJetDefinitions.end(); itdef++) {
1994  AliHFJetDefinition* jetDef = &(*itdef);
1995  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef->GetName());
1996 
1997  fHistManager.CreateHistoGroup(jetDef->GetName(), param.GetName());
1998 
1999  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef->GetName());
2000  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2001  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2002  SetRejectionReasonLabels(h->GetXaxis());
2003 
2004  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef->GetName());
2005  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2006  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2007  SetRejectionReasonLabels(h->GetXaxis());
2008 
2009  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef->GetName());
2010  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
2011  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2012  SetRejectionReasonLabels(h->GetXaxis());
2013 
2014  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef->GetName());
2015  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
2016  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
2017 
2018  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef->GetName());
2019  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
2020  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2021 
2022  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef->GetName());
2023  htitle = hname + ";#it{#eta}_{jet};counts";
2024  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2025 
2026  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef->GetName());
2027  htitle = hname + ";#it{#phi}_{jet};counts";
2028  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2029  }
2030  switch (fOutputType) {
2031  case kTreeOutput:
2032  param.BuildTree(GetName());
2033  if (treeSlot < fNOutputTrees) {
2034  param.AssignDataSlot(treeSlot+2);
2035  treeSlot++;
2037  }
2038  else {
2039  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
2040  }
2041  break;
2042  case kTHnOutput:
2043  param.BuildHnSparse(fEnabledAxis);
2044  break;
2045  case kNoOutput:
2046  break;
2047  }
2048  }
2049 
2051 
2052  PostData(1, fOutput);
2053 }
2054 
2058 {
2060 
2061  // Load the event
2062  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
2063 
2064  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
2065 
2066  fFastJetWrapper->SetAreaType(fastjet::active_area);
2068 
2069  if (!fAodEvent) {
2070  AliError(Form("This task need an AOD event! Task '%s' will be disabled!", GetName()));
2071  return;
2072  }
2073 
2074  for (auto &params : fAnalysisEngines) {
2075 
2076  params.fAodEvent = fAodEvent;
2077  params.fFastJetWrapper = fFastJetWrapper;
2078  params.Init(fGeom, fAodEvent->GetRunNumber());
2079 
2080  if (params.fMCMode != kMCTruth) {
2081  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
2082 
2083  if (params.fCandidateArray) {
2084  if (!params.fCandidateArray->GetClass()->InheritsFrom("AliAODRecoDecayHF2Prong")) {
2085  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2086  "%s: Objects of type %s in %s are not inherited from AliAODRecoDecayHF2Prong! Task will be disabled!",
2087  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName());
2088  params.fCandidateArray = 0;
2089  params.fInhibit = kTRUE;
2090  }
2091  }
2092  else {
2093  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2094  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
2095  params.fBranchName.Data(), params.GetName());
2096  params.fInhibit = kTRUE;
2097  }
2098  }
2099 
2100  if (params.fMCMode != kNoMC) {
2101  params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(0));
2102 
2103  if (!params.fMCContainer) params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(1));
2104 
2105  if (!params.fMCContainer) {
2106  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2107  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
2108  params.GetName());
2109  params.fInhibit = kTRUE;
2110  }
2111  }
2112 
2113  if (params.fMCMode != kMCTruth) {
2114  params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(0));
2115  if (!params.fTrackContainer) params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(1));
2116 
2117  params.fClusterContainer = GetClusterContainer(0);
2118 
2119  if (!params.fTrackContainer && !params.fClusterContainer) {
2120  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2121  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
2122  params.GetName());
2123  params.fInhibit = kTRUE;
2124  }
2125  }
2126  }
2127 }
2128 
2133 {
2134  if (!fAodEvent) return kFALSE;
2135 
2136  TString hname;
2137 
2138  // fix for temporary bug in ESDfilter
2139  // the AODs with null vertex pointer didn't pass the PhysSel
2140  if (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001) return kFALSE;
2141 
2142  for (auto &eng : fAnalysisEngines) {
2143 
2144  if (eng.fInhibit) continue;
2145 
2146  //Event selection
2147  hname = TString::Format("%s/fHistNEvents", eng.GetName());
2148  Bool_t iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
2149  if (!iseventselected) {
2150  fHistManager.FillTH1(hname, "Rejected");
2151  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2152  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
2153  TString label;
2154  do {
2155  label = GetHFEventRejectionReasonLabel(bitmap);
2156  if (label.IsNull()) break;
2157  fHistManager.FillTH1(hname, label);
2158  } while (true);
2159 
2160  continue;
2161  }
2162 
2163  fHistManager.FillTH1(hname, "Accepted");
2164 
2165  AliDebug(2, "Event selected");
2166 
2167  eng.RunAnalysis();
2168  }
2169  return kTRUE;
2170 }
2171 
2176 {
2177  TString hname;
2178  for (auto &param : fAnalysisEngines) {
2179 
2180  if (param.fInhibit) continue;
2181 
2182  if (fOutputType == kTreeOutput) {
2183  param.FillTree(fApplyKinematicCuts);
2185  }
2186  else if (fOutputType == kTHnOutput) {
2187  param.FillHnSparse(fApplyKinematicCuts);
2188  }
2189  }
2190  return kTRUE;
2191 }
2192 
2200 void AliAnalysisTaskDmesonJets::CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t& minMass, Double_t& maxMass)
2201 {
2202  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
2203 
2204  Double_t mass = part->Mass();
2205 
2206  // To make sure that the PDG mass value is not at the edge of a bin
2207  if (nbins % 2 == 0) {
2208  minMass = mass - range / 2 - range / nbins / 2;
2209  maxMass = mass + range / 2 - range / nbins / 2;
2210  }
2211  else {
2212  minMass = mass - range / 2;
2213  maxMass = mass + range / 2;
2214  }
2215 }
2216 
2223 {
2224  static TString label;
2225  label = "";
2226 
2227  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
2228  label = "NotSelTrigger";
2229  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
2230  return label.Data();
2231  }
2232  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
2233  label = "NoVertex";
2234  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
2235  return label.Data();
2236  }
2237  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
2238  label = "TooFewVtxContrib";
2239  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
2240  return label.Data();
2241  }
2242  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
2243  label = "ZVtxOutFid";
2244  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
2245  return label.Data();
2246  }
2247  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
2248  label = "Pileup";
2249  bitmap &= ~BIT(AliRDHFCuts::kPileup);
2250  return label.Data();
2251  }
2252  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
2253  label = "OutsideCentrality";
2254  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
2255  return label.Data();
2256  }
2257  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
2258  label = "PhysicsSelection";
2259  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
2260  return label.Data();
2261  }
2262  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
2263  label = "BadSPDVertex";
2264  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
2265  return label.Data();
2266  }
2267  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
2268  label = "ZVtxSPDOutFid";
2269  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
2270  return label.Data();
2271  }
2272  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
2273  label = "CentralityFlattening";
2274  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
2275  return label.Data();
2276  }
2277  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
2278  label = "BadTrackV0Correl";
2279  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
2280  return label.Data();
2281  }
2282 
2283  return label.Data();
2284 }
2285 
2292 {
2293  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
2294  PostData(eng.GetDataSlotNumber(), eng.GetTree());
2295  return eng.GetDataSlotNumber();
2296  }
2297  else {
2298  return -1;
2299  }
2300 }
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.
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)
AliVParticle * fDmesonParticle
! pointer to the particle object
static TString GenerateJetName(EJetType_t jetType, EJetAlgo_t jetAlgo, ERecoScheme_t recoScheme, Double_t radius, AliParticleContainer *partCont, AliClusterContainer *clusCont, TString tag)
void SetR(Double_t r)
Definition: AliFJWrapper.h:101
Select tracks based on specific prescriptions of HF analysis.
ECandidateType_t fCandidateType
Candidate type.
Double_t GetDistance(std::string n, Double_t &deta, Double_t &dphi) const
Int_t fMCLabel
! MC label, i.e. index of the generator level D meson (only for detector level D meson candidates) ...
THashList * CreateHistoGroup(const char *groupname, const char *parent="/")
virtual void Set(const AliDmesonJetInfo &source)
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:33
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.
EJetAcceptanceType_t fAcceptance
Jet acceptance.
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)
Bool_t ExtractParticleLevelHFAttributes(const AliAODMCParticle *part, AliDmesonJetInfo &DmesonJet)
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
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 fReconstructed
! Whether this D meson was reconstructed (only for particle level D mesons)
Bool_t ExtractRecoDecayAttributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
Double_t InvMassD0() const
EOutputType_t fOutputType
Output type: none, TTree or THnSparse.
virtual void AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
bool operator==(const AliAnalysisTaskDmesonJets::AliHFJetDefinition &lhs, const AliAnalysisTaskDmesonJets::AliHFJetDefinition &rhs)
Double_t InvMassDstarKpipi() const
AnalysisEngine * AddAnalysisEngine(ECandidateType_t type, EMCMode_t bkgMode, EJetType_t jettype, Double_t jetradius, TString cutfname="")
const Int_t nbins
Double_t maxMass
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)
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)