AliPhysics  781d0c7 (781d0c7)
 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 
62 {
63  dphi = TVector2::Phi_mpi_pi(fMomentum.Phi() - jet.Phi());;
64  deta = fMomentum.Eta() - jet.Eta();
65  return TMath::Sqrt(dphi*dphi + deta*deta);
66 }
67 
73 {
74  Double_t deta = 0;
75  Double_t dphi = 0;
76  return GetDistance(jet, deta, dphi);
77 }
78 
79 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonJetInfo
80 
84 
87 {
88  fD.SetPtEtaPhiE(0,0,0,0);
89  fSoftPionPt = 0;
90  fInvMass2Prong = 0;
91  fDmesonParticle = 0;
92  fMCLabel = -1;
93  fReconstructed = kFALSE;
94  for (auto &jet : fJets) {
95  jet.second.fMomentum.SetPtEtaPhiE(0,0,0,0);
96  jet.second.fNConstituents = 0;
97  jet.second.fNEF = 0;
98  jet.second.fMaxChargedPt = 0;
99  jet.second.fMaxNeutralPt = 0;
100  }
101 }
102 
105 {
106  Printf("Printing D Meson Jet object.");
107  Printf("D Meson: pT = %.3f, eta = %.3f, phi = %.3f, inv. mass = %.3f", fD.Pt(), fD.Eta(), fD.Phi_0_2pi(), fD.M());
108  Printf("Soft pion pT: %.3f. 2-Prong Invariant mass = %.3f", fSoftPionPt, fInvMass2Prong);
109  for (auto &jet : fJets) {
110  Printf("Jet %s: pT = %.3f, eta = %.3f, phi = %.3f", jet.first.c_str(), jet.second.Pt(), jet.second.Eta(), jet.second.Phi_0_2pi());
111  Printf("Jet N Consituents = %d", jet.second.fNConstituents);
112  }
113 }
114 
119 {
120  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
121  if (it == fJets.end()) return 0;
122 
123  Double_t z = 0;
124 
125  if ((*it).second.Pt() > 0) {
126  TVector3 dvect = fD.Vect();
127  TVector3 jvect = (*it).second.fMomentum.Vect();
128 
129  Double_t jetMom = jvect * jvect;
130 
131  if (jetMom < 1e-6) {
132  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetZ", "Zero jet momentum!");
133  z = 0.999;
134  }
135  else {
136  z = (dvect * jvect) / jetMom;
137  }
138 
139  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin 0.9-1 rather than 1-1.1
140  }
141 
142  return z;
143 }
144 
152 {
153  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
154  if (it == fJets.end()) return 0;
155 
156  return GetDistance((*it).second, deta, dphi);
157 }
158 
164 {
165  Double_t deta = 0;
166  Double_t dphi = 0;
167  return GetDistance(n, deta, dphi);
168 }
169 
177 {
178  dphi = TVector2::Phi_mpi_pi(fD.Phi() - jet.Phi());;
179  deta = fD.Eta() - jet.Eta();
180  return TMath::Sqrt(dphi*dphi + deta*deta);
181 }
182 
188 {
189  Double_t deta = 0;
190  Double_t dphi = 0;
191  return GetDistance(jet, deta, dphi);
192 }
193 
199 {
200  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
201  if (it == fJets.end()) {
202  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
203  n.c_str());
204  return 0;
205  }
206  return &((*it).second);
207 }
208 
214 {
215  std::map<std::string, AliJetInfo>::iterator it = fJets.find(n);
216  if (it == fJets.end()) {
217  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
218  n.c_str());
219  return 0;
220  }
221  return &((*it).second);
222 }
223 
224 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoSummary
225 
229 
235  fPt(0),
236  fEta(0),
237  fPhi(0),
238  fR(0),
239  fZ(0)
240 {
241  Set(source, n);
242 }
243 
246 {
247  fPt = 0;
248  fEta = 0;
249  fPhi = 0;
250  fR = 0;
251  fZ = 0;
252 }
253 
259 {
260  std::map<std::string, AliJetInfo>::const_iterator it = source.fJets.find(n);
261  if (it == source.fJets.end()) return;
262 
263  fPt = (*it).second.Pt();
264  fEta = (*it).second.Eta();
265  fPhi = (*it).second.Phi_0_2pi();
266  fR = source.GetDistance(n);
267  fZ = source.GetZ(n);
268 }
269 
274 {
275  fPt = source.Pt();
276  fEta = source.Eta();
277  fPhi = source.Phi_0_2pi();
278  fR = 0;
279  fZ = 0;
280 }
281 
282 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonInfoSummary
283 
287 
292  fPt(0),
293  fEta(0),
294  fPhi(0)
295 {
296  Set(source);
297 }
298 
303 {
304  fPt = source.fD.Pt();
305  fEta = source.fD.Eta();
306  fPhi = source.fD.Phi_0_2pi();
307 }
308 
311 {
312  fPt = 0;
313  fEta = 0;
314  fPhi = 0;
315 }
316 
317 // Definitions of class AliAnalysisTaskDmesonJets::AliD0InfoSummary
318 
322 
327  AliDmesonInfoSummary(source),
328  fInvMass(source.fD.M())
329 {
330 }
331 
336 {
337  fInvMass = source.fD.M();
339 }
340 
343 {
345 
346  fInvMass = 0;
347 }
348 
349 // Definitions of class AliAnalysisTaskDmesonJets::AliDStarInfoSummary
350 
354 
359  AliDmesonInfoSummary(source),
360  f2ProngInvMass(source.fInvMass2Prong),
361  fDeltaInvMass(source.fD.M() - source.fInvMass2Prong)
362 {
363 }
364 
369 {
370  f2ProngInvMass = source.fInvMass2Prong;
371  fDeltaInvMass = source.fD.M() - source.fInvMass2Prong;
373 }
374 
377 {
379 
380  f2ProngInvMass = 0;
381  fDeltaInvMass = 0;
382 }
383 
384 // Definitions of class AliAnalysisTaskDmesonJets::AliJetDefinition
385 
389 
392  TObject(),
393  fJetType(AliJetContainer::kChargedJet),
394  fRadius(0),
395  fJetAlgo(AliJetContainer::antikt_algorithm),
396  fRecoScheme(AliJetContainer::pt_scheme),
397  fMinJetPt(0.),
398  fMaxJetPt(500.),
399  fMinJetPhi(0.),
400  fMaxJetPhi(0.),
401  fMinJetEta(-1.),
402  fMaxJetEta(1.),
403  fMinChargedPt(0.),
404  fMaxChargedPt(0.),
405  fMinNeutralPt(0.),
406  fMaxNeutralPt(0.)
407 {
408 }
409 
417  TObject(),
418  fJetType(type),
419  fRadius(r),
420  fJetAlgo(algo),
421  fRecoScheme(reco),
422  fMinJetPt(0.),
423  fMaxJetPt(500.),
424  fMinJetPhi(0.),
425  fMaxJetPhi(0.),
426  fMinJetEta(-1.),
427  fMaxJetEta(1.),
428  fMinChargedPt(0.),
429  fMaxChargedPt(0.),
430  fMinNeutralPt(0.),
431  fMaxNeutralPt(0.)
432 {
433 }
434 
439  TObject(),
440  fJetType(source.fJetType),
441  fRadius(source.fRadius),
442  fJetAlgo(source.fJetAlgo),
443  fRecoScheme(source.fRecoScheme),
444  fMinJetPt(source.fMinJetPt),
445  fMaxJetPt(source.fMaxJetPt),
446  fMinJetPhi(source.fMinJetPhi),
447  fMaxJetPhi(source.fMaxJetPhi),
448  fMinJetEta(source.fMinJetEta),
449  fMaxJetEta(source.fMaxJetEta),
450  fMinChargedPt(source.fMinChargedPt),
451  fMaxChargedPt(source.fMaxChargedPt),
452  fMinNeutralPt(source.fMinNeutralPt),
453  fMaxNeutralPt(source.fMaxNeutralPt)
454 {
455 }
456 
461 {
462  new (this) AliHFJetDefinition(source);
463  return *this;
464 }
465 
468 {
469  static TString name;
470 
471  name = AliJetContainer::GenerateJetName(fJetType, fJetAlgo, fRecoScheme, fRadius, 0, 0, "Jet");
472 
473  return name.Data();
474 }
475 
481 {
482  if (fMinJetEta < fMaxJetEta && (jet.Eta() < fMinJetEta || jet.Eta() > fMaxJetEta)) return kFALSE;
483  if (fMinJetPhi < fMaxJetPhi && (jet.Phi() < fMinJetPhi || jet.Phi() > fMaxJetPhi)) return kFALSE;
484  if (jet.Pt() > fMaxJetPt || jet.Pt() < fMinJetPt) return kFALSE;
485  if (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt) return kFALSE;
486  if (jet.fMaxNeutralPt < fMinNeutralPt || jet.fMaxNeutralPt > fMaxNeutralPt) return kFALSE;
487 
488  return kTRUE;
489 }
490 
496 {
497  const AliJetInfo* jet = dMesonJet.GetJet(n);
498  if (!jet) return kFALSE;
499  return IsJetInAcceptance((*jet));
500 }
501 
508 {
509  if (lhs.fJetType > rhs.fJetType) return false;
510  else if (lhs.fJetType < rhs.fJetType) return true;
511  else {
512  if (lhs.fRadius > rhs.fRadius) return false;
513  else if (lhs.fRadius < rhs.fRadius) return true;
514  else {
515  if (lhs.fJetAlgo > rhs.fJetAlgo) return false;
516  else if (lhs.fJetAlgo < rhs.fJetAlgo) return true;
517  else {
518  if (lhs.fRecoScheme < rhs.fRecoScheme) return true;
519  else return false;
520  }
521  }
522  }
523 }
524 
531 {
532  if (lhs.fJetType != rhs.fJetType) return false;
533  if (lhs.fRadius != rhs.fRadius) return false;
534  if (lhs.fJetAlgo != rhs.fJetAlgo) return false;
535  if (lhs.fRecoScheme != rhs.fRecoScheme) return false;
536  return true;
537 }
538 
539 // Definitions of class AliAnalysisTaskDmesonJets::AnalysisEngine
540 
544 
547  TObject(),
548  fCandidateType(kD0toKpi),
549  fCandidateName(),
550  fCandidatePDG(0),
551  fNDaughters(0),
552  fPDGdaughters(),
553  fBranchName(),
554  fMCMode(kNoMC),
555  fNMassBins(0),
556  fMinMass(0),
557  fMaxMass(0),
558  fRDHFCuts(0),
559  fRejectedOrigin(0),
560  fAcceptedDecay(0),
561  fInhibit(kFALSE),
562  fJetDefinitions(),
563  fPtBinWidth(0.5),
564  fMaxPt(100),
565  fDataSlotNumber(-1),
566  fTree(0),
567  fCurrentDmesonJetInfo(0),
568  fCurrentJetInfo(0),
569  fCandidateArray(0),
570  fMCContainer(0),
571  fTrackContainer(0),
572  fClusterContainer(0),
573  fAodEvent(0),
574  fFastJetWrapper(0),
575  fHistManager(0)
576 {
577 }
578 
587  TObject(),
588  fCandidateType(type),
589  fCandidateName(),
590  fCandidatePDG(0),
591  fNDaughters(0),
592  fPDGdaughters(),
593  fBranchName(),
594  fMCMode(MCmode),
595  fNMassBins(nMassBins),
596  fMinMass(0),
597  fMaxMass(0),
598  fRDHFCuts(cuts),
600  fAcceptedDecay(kAnyDecay),
601  fInhibit(kFALSE),
602  fJetDefinitions(),
603  fPtBinWidth(0.5),
604  fMaxPt(100),
605  fDataSlotNumber(-1),
606  fTree(0),
607  fCurrentDmesonJetInfo(0),
608  fCurrentJetInfo(0),
609  fCandidateArray(0),
610  fMCContainer(0),
611  fTrackContainer(0),
612  fClusterContainer(0),
613  fAodEvent(0),
614  fFastJetWrapper(0),
615  fHistManager(0)
616 {
617  SetCandidateProperties(range);
618 }
619 
624  TObject(source),
625  fCandidateType(source.fCandidateType),
626  fCandidateName(source.fCandidateName),
627  fCandidatePDG(source.fCandidatePDG),
628  fNDaughters(source.fNDaughters),
629  fPDGdaughters(source.fPDGdaughters),
630  fBranchName(source.fBranchName),
631  fMCMode(source.fMCMode),
632  fNMassBins(source.fNMassBins),
633  fMinMass(source.fMinMass),
634  fMaxMass(source.fMaxMass),
635  fRDHFCuts(),
636  fRejectedOrigin(source.fRejectedOrigin),
637  fAcceptedDecay(source.fAcceptedDecay),
638  fInhibit(source.fInhibit),
639  fJetDefinitions(source.fJetDefinitions),
640  fPtBinWidth(source.fPtBinWidth),
641  fMaxPt(source.fMaxPt),
642  fDataSlotNumber(-1),
643  fTree(0),
644  fCurrentDmesonJetInfo(0),
645  fCurrentJetInfo(0),
646  fCandidateArray(source.fCandidateArray),
647  fMCContainer(source.fMCContainer),
648  fTrackContainer(source.fTrackContainer),
649  fClusterContainer(source.fClusterContainer),
650  fAodEvent(source.fAodEvent),
652  fHistManager(source.fHistManager)
653 {
654  SetRDHFCuts(source.fRDHFCuts);
655 }
656 
657 // Destructor
659 {
660  delete fRDHFCuts;
661 }
662 
667 {
668  new (this) AnalysisEngine(source);
669  return *this;
670 }
671 
676 {
677  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
678  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
679  }
680 
681  return kFALSE;
682 }
683 
685 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const /*geom*/, Int_t /*runNumber*/)
686 {
687 }
688 
693 {
694  switch (fCandidateType) {
695  case kD0toKpi:
696  fCandidatePDG = 421;
697  fCandidateName = "D0";
698  fNDaughters = 2;
699  fPDGdaughters.Set(fNDaughters);
700  fPDGdaughters.Reset();
701  fPDGdaughters[0] = 211; // pi
702  fPDGdaughters[1] = 321; // K
703  fBranchName = "D0toKpi";
704  fAcceptedDecay = kDecayD0toKpi;
705  if (!fRDHFCuts) {
706  fRDHFCuts = new AliRDHFCutsD0toKpi();
707  fRDHFCuts->SetStandardCutsPP2010();
708  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
709  fRDHFCuts->SetTriggerClass("","");
710  }
711  break;
712  case kD0toKpiLikeSign:
713  fCandidatePDG = 421;
714  fCandidateName = "2ProngLikeSign";
715  fNDaughters = 2;
716  fPDGdaughters.Set(fNDaughters);
717  fPDGdaughters.Reset();
718  fPDGdaughters[0] = 211; // pi
719  fPDGdaughters[1] = 321; // K
720  fBranchName = "LikeSign2Prong";
721  fAcceptedDecay = kAnyDecay;
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 = kDecayDStartoKpipi;
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  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) { // D0 candidate
901  return ExtractD0Attributes(Dcand, DmesonJet, i);
902  }
903  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
904  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
905  }
906  else {
907  return kFALSE;
908  }
909 }
910 
919 {
920  //AliDebug(2,"Checking if D0 meson is selected");
921  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
922  if (isSelected == 0) return kFALSE;
923 
924  Int_t MCtruthPdgCode = 0;
925 
926  Double_t invMassD = 0;
927 
928  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
929  // Checks also the origin, and if it matches the rejected origin mask, return false
930  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
931  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
932  DmesonJet.fMCLabel = mcLab;
933 
934  // Retrieve the generated particle (if exists) and its PDG code
935  if (mcLab >= 0) {
936  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
937 
938  if (aodMcPart) {
939  // Check origin and return false if it matches the rejected origin mask
940  if (fRejectedOrigin) {
941  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
942 
943  UInt_t rs = origin;
944  UShort_t p = 0;
945  while (rs >>= 1) { p++; }
946  TString hname = TString::Format("%s/fHistDmesonOrigin", GetName());
947  fHistManager->FillTH1(hname, p);
948 
949  if ((origin & fRejectedOrigin) == origin) return kFALSE;
950  }
951  MCtruthPdgCode = aodMcPart->PdgCode();
952  }
953  }
954  }
955 
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  //AliDebug(2,"Checking if D0 meson is selected");
1036  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoCascadeHF*>(DstarCand), AliRDHFCuts::kAll, fAodEvent);
1037  if (isSelected == 0) return kFALSE;
1038 
1039  if (i > 0) return kFALSE; // only one mass hypothesis for the D*
1040 
1041  Int_t MCtruthPdgCode = 0;
1042 
1043  Double_t invMassD = 0;
1044 
1045  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1046  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1047  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1048 
1049  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1050  //AliDebug(2, Form("MC label is %d", mcLab));
1051  DmesonJet.fMCLabel = mcLab;
1052  if (mcLab >= 0) {
1053  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1054 
1055  if (aodMcPart) {
1056  if (fRejectedOrigin) {
1057  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
1058 
1059  UInt_t rs = origin;
1060  UShort_t p = 0;
1061  while (rs >>= 1) { p++; }
1062  TString hname = TString::Format("%s/fHistDmesonOrigin", GetName());
1063  fHistManager->FillTH1(hname, p);
1064 
1065  if ((origin & fRejectedOrigin) == origin) return kFALSE;
1066  }
1067 
1068  MCtruthPdgCode = aodMcPart->PdgCode();
1069  //AliDebug(2, Form("MC truth pdg code is %d",MCtruthPdgCode));
1070  }
1071  }
1072  }
1073 
1074  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1075  if (fMCMode == kNoMC ||
1076  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1077  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1078  // 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)
1079  invMassD = DstarCand->InvMassDstarKpipi();
1080  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1081  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1082  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1083  return kTRUE;
1084  }
1085  else { // conditions above not passed, so return FALSE
1086  return kFALSE;
1087  }
1088 }
1089 
1097 {
1098  if (!part) return kUnknownDecay;
1099  if (!mcArray) return kUnknownDecay;
1100 
1102 
1103  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1104 
1105  if (part->GetNDaughters() == 2) {
1106 
1107  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1108  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1109 
1110  if (!d1 || !d2) {
1111  return decay;
1112  }
1113 
1114  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1115  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1116 
1117  if (absPdgPart == 421) { // D0 -> K pi
1118  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1119  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1120  decay = kDecayD0toKpi;
1121  }
1122  }
1123 
1124  if (absPdgPart == 413) { // D* -> D0 pi
1125  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1126  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1127  if (D0decay == kDecayD0toKpi) {
1128  decay = kDecayDStartoKpipi;
1129  }
1130  }
1131  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1132  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1133  if (D0decay == kDecayD0toKpi) {
1134  decay = kDecayDStartoKpipi;
1135  }
1136  }
1137  }
1138  }
1139 
1140  return decay;
1141 }
1142 
1150 {
1151  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1152 
1153  if (!part) return kUnknownQuark;
1154  if (!mcArray) return kUnknownQuark;
1155 
1156  Int_t mother = part->GetMother();
1157  while (mother >= 0) {
1158  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1159  if (mcGranma) {
1160  Int_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
1161 
1162  if (abspdgGranma == 1) return kFromDown;
1163  if (abspdgGranma == 2) return kFromUp;
1164  if (abspdgGranma == 3) return kFromStrange;
1165  if (abspdgGranma == 4) return kFromCharm;
1166  if (abspdgGranma == 5) return kFromBottom;
1167  if (abspdgGranma == 6) return kFromTop;
1168 
1169  mother = mcGranma->GetMother();
1170  }
1171  else {
1172  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1173  break;
1174  }
1175  }
1176 
1177  return kUnknownQuark;
1178 }
1179 
1182 {
1183  fDmesonJets.clear();
1184 
1185  for (auto& jetDef : fJetDefinitions) {
1186  jetDef.fJets.clear();
1187  }
1188 
1189  if (fMCMode == kMCTruth) {
1190  RunParticleLevelAnalysis();
1191  }
1192  else {
1193  RunDetectorLevelAnalysis();
1194  }
1195 }
1196 
1202 {
1204  fFastJetWrapper->SetR(jetDef.fRadius);
1207 
1208  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1209  fTrackContainer->SetDMesonCandidate(0);
1210  AddInputVectors(fTrackContainer, 100);
1211  }
1212 
1213  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1214  AddInputVectors(fClusterContainer, -100);
1215  }
1216 
1217  // run jet finder
1218  fFastJetWrapper->Run();
1219 
1220  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1221 
1222  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1223  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1224 
1225  Double_t maxChPt = 0;
1226  Double_t maxNePt = 0;
1227  Double_t totalNeutralPt = 0;
1228 
1229  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1230  if (constituents[ic].user_index() >= 100) {
1231  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1232  }
1233  else if (constituents[ic].user_index() <= -100) {
1234  totalNeutralPt += constituents[ic].pt();
1235  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1236  }
1237  }
1238 
1239  jetDef.fJets.push_back(
1240  AliJetInfo(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E(),
1241  constituents.size(), totalNeutralPt / jets_incl[ijet].pt(), maxChPt, maxNePt)
1242  );
1243  }
1244 }
1245 
1254 {
1255  if (jetDef.fJets.size() == 0) FindJets(jetDef);
1256 
1257  Double_t d_closest = 999;
1258  AliJetInfo* jet_closest = 0;
1259  const AliJetInfo* truth_jet = 0;
1260  try {
1261  truth_jet = &(dmeson.fJets.at(jetDef.GetName()));
1262  }
1263  catch(...) {
1264  return jet_distance_pair(0, 999);
1265  }
1266 
1267  for (auto& jet : jetDef.fJets) {
1268  Double_t d = truth_jet->GetDistance(jet);
1269 
1270  if (d > dMax) continue;
1271  if (d < d_closest) {
1272  d_closest = d;
1273  jet_closest = &jet;
1274  }
1275  }
1276 
1277  if (jet_closest && applyKinCuts) {
1278  if (!jetDef.IsJetInAcceptance(*jet_closest)) jet_closest = 0;
1279  }
1280 
1281  if (jet_closest) {
1282  AliDebug(2, Form("Found closest jet (pt=%.3f, eta=%.3f, phi=%.3f) to d meson (pt=%.3f, eta=%.3f, phi=%.3f) with d = %.3f",
1283  jet_closest->Pt(), jet_closest->Eta(), jet_closest->Phi_0_2pi(),
1284  dmeson.fD.Pt(), dmeson.fD.Eta(), dmeson.fD.Phi_0_2pi(),
1285  d_closest));
1286  }
1287 
1288  return jet_distance_pair(jet_closest, d_closest);
1289 }
1290 
1293 {
1294  const Int_t nD = fCandidateArray->GetEntriesFast();
1295 
1296  AliDmesonJetInfo DmesonJet;
1297 
1298  Int_t nAccCharm = 0;
1299  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1300  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1301  if (!charmCand) continue;
1302 
1303  // region of interest + cuts
1304  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1305 
1306  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1307  DmesonJet.Reset();
1308  DmesonJet.fDmesonParticle = charmCand;
1309  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1310  for (std::vector<AliHFJetDefinition>::iterator itdef = fJetDefinitions.begin(); itdef != fJetDefinitions.end(); itdef++) {
1311  if (!FindJet(charmCand, DmesonJet, *itdef)) {
1312  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1313  (*itdef).GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1314  }
1315  }
1316  fDmesonJets[icharm] = DmesonJet;
1317  }
1318  }
1319  nAccCharm++;
1320  } // end of D cand loop
1321 
1322  TString hname;
1323 
1324  hname = TString::Format("%s/fHistNAcceptedDmesons", GetName());
1325  fHistManager->FillTH1(hname, nAccCharm);
1326 
1327  hname = TString::Format("%s/fHistNDmesons", GetName());
1328  fHistManager->FillTH1(hname, nD);
1329 }
1330 
1341 {
1342  TString hname;
1343 
1345  fFastJetWrapper->SetR(jetDef.fRadius);
1348 
1349  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1350 
1351  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1352  fTrackContainer->SetDMesonCandidate(Dcand);
1353  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1354  AddInputVectors(fTrackContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1355 
1356  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1357  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1358  const TObjArray& daughters = fTrackContainer->GetDaughterList();
1359  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1360  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1361  if (!fTrackContainer->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1362  }
1363  }
1364 
1365  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1366  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1367  AddInputVectors(fClusterContainer, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1368  }
1369 
1370  // run jet finder
1371  fFastJetWrapper->Run();
1372 
1373  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1374 
1375  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1376  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1377 
1378  Bool_t isDmesonJet = kFALSE;
1379 
1380  Double_t maxChPt = 0;
1381  Double_t maxNePt = 0;
1382  Double_t totalNeutralPt = 0;
1383 
1384  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1385  if (constituents[ic].user_index() == 0) {
1386  isDmesonJet = kTRUE;
1387  }
1388  else if (constituents[ic].user_index() >= 100) {
1389  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1390  }
1391  else if (constituents[ic].user_index() <= -100) {
1392  totalNeutralPt += constituents[ic].pt();
1393  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1394  }
1395  }
1396 
1397  if (isDmesonJet) {
1398  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1399  DmesonJet.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1400  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1401  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1402  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1403 
1404  return kTRUE;
1405  }
1406  }
1407 
1408  return kFALSE;
1409 }
1410 
1414 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist)
1415 {
1416  AliEmcalIterableMomentumContainer itcont = cont->all_momentum();
1417  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1418  UInt_t rejectionReason = 0;
1419  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1420  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1421  continue;
1422  }
1423  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1424  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1425  }
1426 }
1427 
1430 {
1431  TString hname;
1432 
1433  fMCContainer->SetSpecialPDG(fCandidatePDG);
1434  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1435  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1436 
1437  if (!fMCContainer->IsSpecialPDGFound()) return;
1438 
1439  for (auto &jetDef : fJetDefinitions) {
1440  switch (jetDef.fJetType) {
1442  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
1443  break;
1445  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
1446  break;
1448  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
1449  break;
1450  }
1451 
1453  fFastJetWrapper->SetR(jetDef.fRadius);
1456 
1457  hname = TString::Format("%s/fHistDmesonOrigin", GetName());
1458  fMCContainer->SetHistOrigin(static_cast<TH1*>(fHistManager->FindObject(hname)));
1459  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1460  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1461  fMCContainer->SetHistOrigin(0);
1462 
1463  fFastJetWrapper->Run();
1464 
1465  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1466 
1467  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1468  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1469 
1470  Bool_t isDmesonJet = kFALSE;
1471 
1472  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1473  Int_t iPart = constituents[ic].user_index() - 100;
1474  AliVParticle* part = fMCContainer->GetParticle(iPart);
1475  if (!part) {
1476  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1477  continue;
1478  }
1479  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1480  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1481  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1482  std::pair<int, AliDmesonJetInfo> element;
1483  element.first = iPart;
1484  element.second = AliDmesonJetInfo();
1485  dMesonJetIt = fDmesonJets.insert(element).first;
1486  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1487  (*dMesonJetIt).second.fDmesonParticle = part;
1488  }
1489 
1490  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1491  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1492  }
1493  }
1494  }
1495  }
1496 }
1497 
1502 {
1503  TString classname;
1504  switch (fCandidateType) {
1505  case kD0toKpi:
1506  case kD0toKpiLikeSign:
1507  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1508  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1509  break;
1510  case kDstartoKpipi:
1511  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1512  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1513  break;
1514  }
1515  TString treeName = TString::Format("%s_%s", taskName, GetName());
1516  fTree = new TTree(treeName, treeName);
1517  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
1518  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1519  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1520  fCurrentJetInfo[i] = new AliJetInfoSummary();
1521  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1522  }
1523 
1524  return fTree;
1525 }
1526 
1531 {
1532  TString hname;
1533 
1534  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
1535 
1536  for (auto &jetDef : fJetDefinitions) {
1537 
1538  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
1539 
1540  Double_t radius = jetDef.fRadius;
1541 
1542  TString title[30] = {""};
1543  Int_t nbins[30] = {0 };
1544  Double_t min [30] = {0.};
1545  Double_t max [30] = {0.};
1546  Int_t dim = 0 ;
1547 
1548  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1549  nbins[dim] = nPtBins;
1550  min[dim] = 0;
1551  max[dim] = fMaxPt;
1552  dim++;
1553 
1554  if ((enabledAxis & kPositionD) != 0) {
1555  title[dim] = "#eta_{D}";
1556  nbins[dim] = 50;
1557  min[dim] = -1;
1558  max[dim] = 1;
1559  dim++;
1560 
1561  title[dim] = "#phi_{D} (rad)";
1562  nbins[dim] = 150;
1563  min[dim] = 0;
1564  max[dim] = TMath::TwoPi();
1565  dim++;
1566  }
1567 
1568  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1569  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1570  nbins[dim] = fNMassBins;
1571  min[dim] = fMinMass;
1572  max[dim] = fMaxMass;
1573  dim++;
1574  }
1575 
1576  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1577  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1578  nbins[dim] = fNMassBins;
1579  min[dim] = fMinMass;
1580  max[dim] = fMaxMass;
1581  dim++;
1582  }
1583 
1584  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1585  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1586  nbins[dim] = fNMassBins;
1587  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1588  dim++;
1589  }
1590 
1591  if (fCandidateType == kDstartoKpipi) {
1592  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1593  nbins[dim] = fNMassBins*6;
1594  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1595 
1596  // subtract mass of D0
1597  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1598  min[dim] -= D0mass;
1599  max[dim] -= D0mass;
1600 
1601  dim++;
1602  }
1603 
1604  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1605  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1606  nbins[dim] = 100;
1607  min[dim] = 0;
1608  max[dim] = 25;
1609  dim++;
1610  }
1611 
1612  title[dim] = "#it{z}_{D}";
1613  nbins[dim] = 110;
1614  min[dim] = 0;
1615  max[dim] = 1.10;
1616  dim++;
1617 
1618  if ((enabledAxis & kDeltaR) != 0) {
1619  title[dim] = "#Delta R_{D-jet}";
1620  nbins[dim] = 100;
1621  min[dim] = 0;
1622  max[dim] = radius * 1.5;
1623  dim++;
1624  }
1625 
1626  if ((enabledAxis & kDeltaEta) != 0) {
1627  title[dim] = "#eta_{D} - #eta_{jet}";
1628  nbins[dim] = 100;
1629  min[dim] = -radius * 1.2;
1630  max[dim] = radius * 1.2;
1631  dim++;
1632  }
1633 
1634  if ((enabledAxis & kDeltaPhi) != 0) {
1635  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1636  nbins[dim] = 100;
1637  min[dim] = -radius * 1.2;
1638  max[dim] = radius * 1.2;
1639  dim++;
1640  }
1641 
1642  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1643  nbins[dim] = nPtBins;
1644  min[dim] = 0;
1645  max[dim] = fMaxPt;
1646  dim++;
1647 
1648  if ((enabledAxis & kPositionJet) != 0) {
1649  title[dim] = "#eta_{jet}";
1650  nbins[dim] = 50;
1651  min[dim] = -1;
1652  max[dim] = 1;
1653  dim++;
1654 
1655  title[dim] = "#phi_{jet} (rad)";
1656  nbins[dim] = 150;
1657  min[dim] = 0;
1658  max[dim] = TMath::TwoPi();
1659  dim++;
1660  }
1661 
1662  if ((enabledAxis & kJetConstituents) != 0) {
1663  title[dim] = "No. of constituents";
1664  nbins[dim] = 50;
1665  min[dim] = -0.5;
1666  max[dim] = 49.5;
1667  dim++;
1668  }
1669 
1670  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1671  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1672  for (Int_t j = 0; j < dim; j++) {
1673  h->GetAxis(j)->SetTitle(title[j]);
1674  }
1675  }
1676 }
1677 
1682 {
1683  TString hname;
1684 
1685  for (auto& dmeson_pair : fDmesonJets) {
1686  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
1687  Int_t accJets = 0;
1688  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1689  fCurrentJetInfo[ij]->Reset();
1690  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1691  if (!jet) continue;
1692  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1693  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1694  fHistManager->FillTH1(hname, jet->Pt());
1695  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1696  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1697  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1698  fHistManager->FillTH1(hname, jet->Eta());
1699  continue;
1700  }
1701  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
1702  accJets++;
1703  }
1704  if (accJets > 0) {
1705  fTree->Fill();
1706  }
1707  else {
1708  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1709  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1710  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1711  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1712  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1713  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1714  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1715  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1716  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1717  }
1718  else if (fCandidateType == kDstartoKpipi) {
1719  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1720  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1721 
1722  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1723  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1724  }
1725  }
1726  }
1727  return kTRUE;
1728 }
1729 
1735 {
1736  TString hname;
1737 
1738  for (auto& dmeson_pair : fDmesonJets) {
1739  Int_t accJets = 0;
1740  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1741  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1742  if (!jet) continue;
1743  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1744  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1745  fHistManager->FillTH1(hname, jet->Pt());
1746  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1747  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1748  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1749  fHistManager->FillTH1(hname, jet->Eta());
1750  continue;
1751  }
1752  accJets++;
1753  }
1754  if (accJets == 0) {
1755  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1756  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1757  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1758  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1759  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1760  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1761  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1762  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1763  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1764  }
1765  else if (fCandidateType == kDstartoKpipi) {
1766  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1767  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1768 
1769  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1770  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1771  }
1772  }
1773  }
1774  return kTRUE;
1775 }
1776 
1781 {
1782  TString hname;
1783 
1784  for (auto& dmeson_pair : fDmesonJets) {
1785  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
1786  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1787  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1788  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1789  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1790  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1791  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1792  }
1793  }
1794 
1795  for (auto &jetDef : fJetDefinitions) {
1796 
1797  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1798  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
1799 
1800  for (auto& dmeson_pair : fDmesonJets) {
1801  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
1802  if (!jet) continue;
1803  if (!jetDef.IsJetInAcceptance(*jet)) {
1804  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
1805  fHistManager->FillTH1(hname, jet->Pt());
1806  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
1807  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1808  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
1809  fHistManager->FillTH1(hname, jet->Eta());
1810  continue;
1811  }
1812  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
1813  }
1814  }
1815 
1816  return kTRUE;
1817 }
1818 
1825 {
1826  // Fill the THnSparse histogram.
1827 
1828  Double_t contents[30] = {0.};
1829 
1830  Double_t z = DmesonJet.GetZ(n);
1831  Double_t deltaPhi = 0;
1832  Double_t deltaEta = 0;
1833  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
1834 
1835  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
1836  if (it == DmesonJet.fJets.end()) return kFALSE;
1837 
1838  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
1839  TString title(h->GetAxis(i)->GetTitle());
1840  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
1841  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
1842  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
1843  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
1844  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
1845  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
1846  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
1847  else if (title=="#it{z}_{D}") contents[i] = z;
1848  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
1849  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
1850  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
1851  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
1852  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
1853  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
1854  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
1855  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
1856  }
1857 
1858  h->Fill(contents);
1859 
1860  return kTRUE;
1861 }
1862 
1863 // Definitions of class AliAnalysisTaskDmesonJets
1864 
1868 
1872  fAnalysisEngines(),
1873  fEnabledAxis(0),
1875  fHistManager(),
1876  fApplyKinematicCuts(kTRUE),
1877  fNOutputTrees(0),
1878  fAodEvent(0),
1879  fFastJetWrapper(0)
1880 {
1881  SetMakeGeneralHistograms(kTRUE);
1882 }
1883 
1888  AliAnalysisTaskEmcalLight(name, kTRUE),
1889  fAnalysisEngines(),
1890  fEnabledAxis(k2ProngInvMass),
1891  fOutputType(kTreeOutput),
1892  fHistManager(name),
1893  fApplyKinematicCuts(kTRUE),
1894  fNOutputTrees(nOutputTrees),
1895  fAodEvent(0),
1896  fFastJetWrapper(0)
1897 {
1898  SetMakeGeneralHistograms(kTRUE);
1899  for (Int_t i = 0; i < nOutputTrees; i++){
1900  DefineOutput(2+i, TTree::Class());
1901  }
1902 }
1903 
1906 {
1907  if (fFastJetWrapper) delete fFastJetWrapper;
1908 }
1909 
1917 {
1918  AliRDHFCuts* analysiscuts = 0;
1919  TFile* filecuts = TFile::Open(cutfname);
1920  if (!filecuts || filecuts->IsZombie()) {
1921  ::Warning("AddTaskDmesonJets", "Input file not found: will use std cuts.");
1922  filecuts = 0;
1923  }
1924 
1925  if (filecuts) {
1926  analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
1927  if (!analysiscuts) {
1928  ::Warning("AddTaskDmesonJetCorr", "Could not find analysis cuts '%s' in '%s'. Using std cuts.", cutsname.Data(), cutfname.Data());
1929  }
1930  }
1931 
1932  return analysiscuts;
1933 }
1934 
1944 {
1946  return AddAnalysisEngine(type, MCmode, jetDef, cutfname);
1947 }
1948 
1958 {
1959  AliRDHFCuts* cuts = 0;
1960 
1961  if (!cutfname.IsNull()) {
1962  TString cutsname;
1963 
1964  switch (type) {
1965  case kD0toKpi:
1966  case kD0toKpiLikeSign:
1967  cutsname = "D0toKpiCuts";
1968  break;
1969  case kDstartoKpipi:
1970  cutsname = "DStartoKpipiCuts";
1971  break;
1972  default:
1973  return 0;
1974  }
1975 
1976  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
1977  }
1978 
1979  AnalysisEngine eng(type, MCmode, cuts);
1980 
1981  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
1982 
1983  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
1984  eng.AddJetDefinition(jetDef);
1985  it = fAnalysisEngines.insert(it, eng);
1986  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
1987  }
1988  else {
1989  AnalysisEngine* found_eng = &(*it);
1990  ::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());
1991  found_eng->AddJetDefinition(jetDef);
1992 
1993  if (cuts && found_eng->fRDHFCuts != 0) {
1994  ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
1995  found_eng->SetRDHFCuts(cuts);
1996  }
1997  }
1998 
1999  return &(*it);
2000 }
2001 
2002 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
2003 {
2004  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
2005  while (it != fAnalysisEngines.end() && (*it) < eng) it++;
2006  return it;
2007 }
2008 
2011 {
2012  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
2013 
2015 
2016  // Define histograms
2017  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
2018 
2019  TString hname;
2020  TString htitle;
2021  TH1* h = 0;
2022  Int_t treeSlot = 0;
2023 
2024  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
2025  for (auto &param : fAnalysisEngines) {
2026  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
2027 
2028  fHistManager.CreateHistoGroup(param.GetName());
2029 
2030  param.fHistManager = &fHistManager;
2031 
2032  hname = TString::Format("%s/fHistNAcceptedDmesons", param.GetName());
2033  htitle = hname + ";Number of D accepted meson candidates;counts";
2034  h = fHistManager.CreateTH1(hname, htitle, 51, -0.5, 50.5);
2035 
2036  hname = TString::Format("%s/fHistNDmesons", param.GetName());
2037  htitle = hname + ";Number of D meson candidates;counts";
2038  h = fHistManager.CreateTH1(hname, htitle, 101, -0.5, 100.5);
2039 
2040  hname = TString::Format("%s/fHistNEvents", param.GetName());
2041  htitle = hname + ";Event status;counts";
2042  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
2043  h->GetXaxis()->SetBinLabel(1, "Accepted");
2044  h->GetXaxis()->SetBinLabel(2, "Rejected");
2045 
2046  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
2047  htitle = hname + ";Rejection reason;counts";
2048  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
2049 
2050  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
2051  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
2052  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2053 
2054  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
2055  htitle = hname + ";#it{#eta}_{D};counts";
2056  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2057 
2058  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
2059  htitle = hname + ";#it{#phi}_{D};counts";
2060  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2061 
2062  if (param.fCandidateType == kD0toKpi || param.fCandidateType == kD0toKpiLikeSign) {
2063  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
2064  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2065  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
2066  }
2067  else if (param.fCandidateType == kDstartoKpipi) {
2068  Double_t min = 0;
2069  Double_t max = 0;
2070 
2071  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
2072  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2073  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
2074  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
2075 
2076  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2077  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
2078  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2079  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
2080  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
2081  }
2082 
2083  if (param.fMCMode == kBackgroundOnly || param.fMCMode == kSignalOnly || param.fMCMode == kMCTruth) {
2084  hname = TString::Format("%s/fHistDmesonOrigin", param.GetName());
2085  htitle = hname + ";origin;counts";
2086  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2087  }
2088 
2089  for (std::vector<AliHFJetDefinition>::iterator itdef = param.fJetDefinitions.begin(); itdef != param.fJetDefinitions.end(); itdef++) {
2090  AliHFJetDefinition* jetDef = &(*itdef);
2091  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef->GetName());
2092 
2093  fHistManager.CreateHistoGroup(jetDef->GetName(), param.GetName());
2094 
2095  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef->GetName());
2096  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2097  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2098  SetRejectionReasonLabels(h->GetXaxis());
2099 
2100  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef->GetName());
2101  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2102  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2103  SetRejectionReasonLabels(h->GetXaxis());
2104 
2105  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef->GetName());
2106  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
2107  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2108  SetRejectionReasonLabels(h->GetXaxis());
2109 
2110  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef->GetName());
2111  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
2112  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
2113 
2114  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef->GetName());
2115  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
2116  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2117 
2118  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef->GetName());
2119  htitle = hname + ";#it{#eta}_{jet};counts";
2120  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2121 
2122  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef->GetName());
2123  htitle = hname + ";#it{#phi}_{jet};counts";
2124  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2125  }
2126  switch (fOutputType) {
2127  case kTreeOutput:
2128  param.BuildTree(GetName());
2129  if (treeSlot < fNOutputTrees) {
2130  param.AssignDataSlot(treeSlot+2);
2131  treeSlot++;
2133  }
2134  else {
2135  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
2136  }
2137  break;
2138  case kTHnOutput:
2139  param.BuildHnSparse(fEnabledAxis);
2140  break;
2141  case kNoOutput:
2142  break;
2143  }
2144  }
2145 
2147 
2148  PostData(1, fOutput);
2149 }
2150 
2154 {
2156 
2157  // Load the event
2158  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
2159 
2160  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
2161 
2162  fFastJetWrapper->SetAreaType(fastjet::active_area);
2164 
2165  if (!fAodEvent) {
2166  AliError(Form("This task need an AOD event! Task '%s' will be disabled!", GetName()));
2167  return;
2168  }
2169 
2170  for (auto &params : fAnalysisEngines) {
2171 
2172  params.fAodEvent = fAodEvent;
2173  params.fFastJetWrapper = fFastJetWrapper;
2174  params.Init(fGeom, fAodEvent->GetRunNumber());
2175 
2176  if (params.fMCMode != kMCTruth) {
2177  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
2178 
2179  if (params.fCandidateArray) {
2180  TString className;
2181  if (params.fCandidateType == kD0toKpi || params.fCandidateType == kD0toKpiLikeSign) {
2182  className = "AliAODRecoDecayHF2Prong";
2183  }
2184  else if (params.fCandidateType == kDstartoKpipi) {
2185  className = "AliAODRecoCascadeHF";
2186  }
2187  if (!params.fCandidateArray->GetClass()->InheritsFrom(className)) {
2188  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2189  "%s: Objects of type %s in %s are not inherited from %s! Task will be disabled!",
2190  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName(), className.Data());
2191  params.fCandidateArray = 0;
2192  params.fInhibit = kTRUE;
2193  }
2194  }
2195  else {
2196  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2197  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
2198  params.fBranchName.Data(), params.GetName());
2199  params.fInhibit = kTRUE;
2200  }
2201  }
2202 
2203  if (params.fMCMode != kNoMC) {
2204  params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(0));
2205 
2206  if (!params.fMCContainer) params.fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(1));
2207 
2208  if (!params.fMCContainer) {
2209  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2210  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
2211  params.GetName());
2212  params.fInhibit = kTRUE;
2213  }
2214  }
2215 
2216  if (params.fMCMode != kMCTruth) {
2217  params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(0));
2218  if (!params.fTrackContainer) params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(1));
2219 
2220  params.fClusterContainer = GetClusterContainer(0);
2221 
2222  if (!params.fTrackContainer && !params.fClusterContainer) {
2223  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2224  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
2225  params.GetName());
2226  params.fInhibit = kTRUE;
2227  }
2228  }
2229  }
2230 }
2231 
2236 {
2237  if (!fAodEvent) return kFALSE;
2238 
2239  TString hname;
2240 
2241  // fix for temporary bug in ESDfilter
2242  // the AODs with null vertex pointer didn't pass the PhysSel
2243  // Now adding an entry in teh histogram so as to check that this is actually cutting anything out
2244  if (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001) {
2245  for (auto &eng : fAnalysisEngines) {
2246  if (eng.fInhibit) continue;
2247  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2248  fHistManager.FillTH1(hname, "ESDfilterBug");
2249  }
2250  return kFALSE;
2251  }
2252 
2253  for (auto &eng : fAnalysisEngines) {
2254 
2255  if (eng.fInhibit) continue;
2256 
2257  //Event selection
2258  hname = TString::Format("%s/fHistNEvents", eng.GetName());
2259  Bool_t iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
2260  if (!iseventselected) {
2261  fHistManager.FillTH1(hname, "Rejected");
2262  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2263  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
2264  TString label;
2265  do {
2266  label = GetHFEventRejectionReasonLabel(bitmap);
2267  if (label.IsNull()) break;
2268  fHistManager.FillTH1(hname, label);
2269  } while (true);
2270 
2271  continue;
2272  }
2273 
2274  fHistManager.FillTH1(hname, "Accepted");
2275 
2276  AliDebug(2, "Event selected");
2277 
2278  eng.RunAnalysis();
2279  }
2280  return kTRUE;
2281 }
2282 
2287 {
2288  TString hname;
2289  for (auto &param : fAnalysisEngines) {
2290 
2291  if (param.fInhibit) continue;
2292 
2293  if (fOutputType == kTreeOutput) {
2294  param.FillTree(fApplyKinematicCuts);
2296  }
2297  else if (fOutputType == kTHnOutput) {
2298  param.FillHnSparse(fApplyKinematicCuts);
2299  }
2300  }
2301  return kTRUE;
2302 }
2303 
2312 {
2313  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
2314 
2315  Double_t mass = part->Mass();
2316 
2317  // To make sure that the PDG mass value is not at the edge of a bin
2318  if (nbins % 2 == 0) {
2319  minMass = mass - range / 2 - range / nbins / 2;
2320  maxMass = mass + range / 2 - range / nbins / 2;
2321  }
2322  else {
2323  minMass = mass - range / 2;
2324  maxMass = mass + range / 2;
2325  }
2326 }
2327 
2334 {
2335  static TString label;
2336  label = "";
2337 
2338  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
2339  label = "NotSelTrigger";
2340  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
2341  return label.Data();
2342  }
2343  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
2344  label = "NoVertex";
2345  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
2346  return label.Data();
2347  }
2348  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
2349  label = "TooFewVtxContrib";
2350  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
2351  return label.Data();
2352  }
2353  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
2354  label = "ZVtxOutFid";
2355  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
2356  return label.Data();
2357  }
2358  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
2359  label = "Pileup";
2360  bitmap &= ~BIT(AliRDHFCuts::kPileup);
2361  return label.Data();
2362  }
2363  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
2364  label = "OutsideCentrality";
2365  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
2366  return label.Data();
2367  }
2368  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
2369  label = "PhysicsSelection";
2370  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
2371  return label.Data();
2372  }
2373  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
2374  label = "BadSPDVertex";
2375  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
2376  return label.Data();
2377  }
2378  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
2379  label = "ZVtxSPDOutFid";
2380  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
2381  return label.Data();
2382  }
2383  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
2384  label = "CentralityFlattening";
2385  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
2386  return label.Data();
2387  }
2388  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
2389  label = "BadTrackV0Correl";
2390  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
2391  return label.Data();
2392  }
2393 
2394  return label.Data();
2395 }
2396 
2403 {
2404  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
2405  PostData(eng.GetDataSlotNumber(), eng.GetTree());
2406  return eng.GetDataSlotNumber();
2407  }
2408  else {
2409  return -1;
2410  }
2411 }
Int_t pdg
void Print() const
Prints the content of this object in the standard output.
std::list< AnalysisEngine > fAnalysisEngines
Array of analysis parameters.
AliEMCALGeometry * fGeom
!emcal geometry
Analysis task for D meson jets.
virtual void Set(const AliDmesonJetInfo &source)
Invariant mass of the D0 meson candidate in GeV/c2.
virtual Int_t Run()
double Double_t
Definition: External.C:58
Bool_t IsJetInAcceptance(const AliJetInfo &jet) const
Double_t fSoftPionPt
! Transverse momentum of the soft pion of the D* candidate
const char * title
Definition: MakeQAPdf.C:26
AliRDHFCuts * fRDHFCuts
D meson candidates cuts.
Lightweight class that encapsulates D meson jets.
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
ERecoScheme_t fRecoScheme
Jet recombination scheme (pt scheme, E scheme, ...)
Double_t fMaxNeutralPt
Transverse momentum of the leading neutral particle (or cluster)
void RunDetectorLevelAnalysis()
Run a detector level analysis.
Double_t mass
jet_distance_pair FindJetMacthedToGeneratedDMeson(const AliDmesonJetInfo &dmeson, AliHFJetDefinition &jetDef, Double_t dMax, Bool_t applyKinCuts)
virtual void UserCreateOutputObjects()
Creates the output containers.
Lightweight class that encapsulates D meson jets.
Double_t fInvMass2Prong
! 2-prong mass of the D* candidate (w/o the soft pion)
static void CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t &minMass, Double_t &maxMass)
AliAODEvent * fAodEvent
! AOD event
Bool_t FindJet(AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, AliHFJetDefinition &jetDef)
UInt_t fEnabledAxis
Use bit defined in EAxis_t to enable axis in the THnSparse.
EMCMode_t fMCMode
MC mode: No MC (data and MC detector level), background-only (MC), signal-only (MC), MC truth (particle level)
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h: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)
std::pair< AliJetInfo *, Double_t > jet_distance_pair
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
virtual void Set(const AliDmesonJetInfo &source)
AliTLorentzVector fD
! 4-momentum of the D meson candidate
Double_t GetDistance(const AliJetInfo &jet, Double_t &deta, Double_t &dphi) const
TObject * FindObject(const char *name) const
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Definition: THistManager.h:504
AliAnalysisTaskDmesonJets()
This is the default constructor, used for ROOT I/O purposes.
virtual void Reset()
Reset the current object.
AnalysisEngine & operator=(const AnalysisEngine &source)
Struct that encapsulates analysis parameters.
std::vector< AliHFJetDefinition >::iterator FindJetDefinition(const AliHFJetDefinition &eng)
Base task in the EMCAL framework (lighter version of AliAnalysisTaskEmcal)
Double_t Phi_0_2pi() const
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:95
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
AliClusterContainer * GetClusterContainer(Int_t i=0) const
AliAODTrack * GetBachelor() const
AliRDHFCuts * LoadDMesonCutsFromFile(TString cutfname, TString cutsname)
bool operator<(const AliAnalysisTaskDmesonJets::AliHFJetDefinition &lhs, const AliAnalysisTaskDmesonJets::AliHFJetDefinition &rhs)
std::list< AnalysisEngine >::iterator FindAnalysisEngine(const AnalysisEngine &eng)
std::vector< AliHFJetDefinition > fJetDefinitions
Jet definitions.
EJetAlgo_t fJetAlgo
Jet algorithm (kt, anti-kt,...)
static const char * GetHFEventRejectionReasonLabel(UInt_t &bitmap)
Lightweight class that encapsulates D*.
virtual void Clear(const Option_t *="")
Bool_t ExtractDstarAttributes(const AliAODRecoCascadeHF *DstarCand, AliDmesonJetInfo &DmesonJet, UInt_t i)
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
const AliJetInfo * GetJet(std::string n) const
virtual ~AliAnalysisTaskDmesonJets()
This is the standard destructor.
Bool_t ExtractD0Attributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
std::vector< AliJetInfo > fJets
! Inclusive jets reconstructed in the current event (includes D meson candidate daughters, if any)
Select MC particles based on specific prescriptions of HF analysis.
static EMesonDecayChannel_t CheckDecayChannel(const AliAODMCParticle *part, TClonesArray *mcArray)
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:99
Double_t fMaxChargedPt
Transverse momentum of the leading charged particle (or track)
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
static FJJetAlgo ConvertToFJAlgo(EJetAlgo_t algo)
decay
Definition: HFPtSpectrum.C:41
std::map< std::string, AliJetInfo > fJets
! list of jets
void AddInputVectors(AliEmcalContainer *cont, Int_t offset, TH2 *rejectHist=0)
Definition: External.C:220
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Double_t minMass
void RunParticleLevelAnalysis()
Run a particle level analysis.
Bool_t fApplyKinematicCuts
Apply jet kinematic cuts.
const char * GetName() const
Generate a name for this jet definition.
unsigned short UShort_t
Definition: External.C:28
Bool_t ExtractRecoDecayAttributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
Double_t InvMassD0() const
EOutputType_t fOutputType
Output type: none, TTree or THnSparse.
virtual void AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
bool operator==(const AliAnalysisTaskDmesonJets::AliHFJetDefinition &lhs, const AliAnalysisTaskDmesonJets::AliHFJetDefinition &rhs)
Double_t InvMassDstarKpipi() const
AnalysisEngine * AddAnalysisEngine(ECandidateType_t type, EMCMode_t bkgMode, EJetType_t jettype, Double_t jetradius, TString cutfname="")
const Int_t nbins
Double_t maxMass
bool Bool_t
Definition: External.C:53
Int_t fNOutputTrees
Maximum number of output trees.
AliHFJetDefinition()
This is the default constructor, used for ROOT I/O purposes.
AliHFJetDefinition & operator=(const AliHFJetDefinition &source)
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:97
AliTLorentzVector fMomentum
4-momentum of the jet
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
AliHFJetDefinition * AddJetDefinition(EJetType_t type, Double_t r, EJetAlgo_t algo, ERecoScheme_t reco)
Lightweight class that encapsulates D0.
const Int_t nPtBins
Bool_t IsAnyJetInAcceptance(const AliDmesonJetInfo &dMesonJet) const
Container for jet within the EMCAL jet framework.
Definition: External.C:196
Int_t PostDataFromAnalysisEngine(const AnalysisEngine &eng)
static EMesonOrigin_t CheckOrigin(const AliAODMCParticle *part, TClonesArray *mcArray)
AliFJWrapper * fFastJetWrapper
! Fastjet wrapper
THistManager fHistManager
Histogram manager.
void Reset()
Reset all fields to their default values.
virtual void Set(const AliDmesonJetInfo &source, std::string n)