AliPhysics  fde8a9f (fde8a9f)
 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 #include "AliAnalysisManager.h"
31 #include "AliVEventHandler.h"
32 
33 // Aliroot HF
34 #include "AliAODRecoCascadeHF.h"
36 #include "AliRDHFCutsD0toKpi.h"
39 #include "AliHFTrackContainer.h"
40 
41 // Aliroot EMCal jet framework
42 #include "AliEmcalJetTask.h"
43 #include "AliEmcalJet.h"
44 #include "AliJetContainer.h"
45 #include "AliParticleContainer.h"
46 #include "AliEmcalParticle.h"
47 #include "AliFJWrapper.h"
48 
50 
51 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfo
52 
56 
64 {
65  dphi = TVector2::Phi_mpi_pi(fMomentum.Phi() - jet.Phi());;
66  deta = fMomentum.Eta() - jet.Eta();
67  return TMath::Sqrt(dphi*dphi + deta*deta);
68 }
69 
75 {
76  Double_t deta = 0;
77  Double_t dphi = 0;
78  return GetDistance(jet, deta, dphi);
79 }
80 
81 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonJetInfo
82 
86 
89  fDmesonParticle(0),
90  fD(),
91  fSoftPionPt(0),
92  fInvMass2Prong(0),
93  fJets(),
94  fMCLabel(-1),
95  fReconstructed(kFALSE),
96  fFirstParton(0),
97  fFirstPartonType(0),
98  fLastParton(0),
99  fLastPartonType(0),
100  fSelectionType(0)
101 {
102 }
103 
108  fDmesonParticle(source.fDmesonParticle),
109  fD(source.fD),
110  fSoftPionPt(source.fSoftPionPt),
111  fInvMass2Prong(source.fInvMass2Prong),
112  fJets(source.fJets),
113  fMCLabel(source.fMCLabel),
114  fReconstructed(source.fReconstructed),
115  fFirstParton(source.fFirstParton),
116  fFirstPartonType(source.fFirstPartonType),
117  fLastParton(source.fLastParton),
118  fLastPartonType(source.fLastPartonType),
119  fSelectionType(source.fSelectionType)
120 {
121 }
122 
127 {
128  new (this) AliDmesonJetInfo(source);
129  return *this;
130 }
131 
134 {
135  fD.SetPtEtaPhiE(0,0,0,0);
136  fSoftPionPt = 0;
137  fInvMass2Prong = 0;
138  fDmesonParticle = 0;
139  fMCLabel = -1;
140  fReconstructed = kFALSE;
141  fFirstParton = 0;
142  fFirstPartonType = 0;
143  fLastParton = 0;
144  fLastPartonType = 0;
145  for (auto &jet : fJets) {
146  jet.second.fMomentum.SetPtEtaPhiE(0,0,0,0);
147  jet.second.fNConstituents = 0;
148  jet.second.fNEF = 0;
149  jet.second.fMaxChargedPt = 0;
150  jet.second.fMaxNeutralPt = 0;
151  }
152 }
153 
156 {
157  Printf("Printing D Meson Jet object.");
158  Printf("D Meson: pT = %.3f, eta = %.3f, phi = %.3f, inv. mass = %.3f", fD.Pt(), fD.Eta(), fD.Phi_0_2pi(), fD.M());
159  Printf("Soft pion pT: %.3f. 2-Prong Invariant mass = %.3f", fSoftPionPt, fInvMass2Prong);
160  for (auto &jet : fJets) {
161  Printf("Jet %s: pT = %.3f, eta = %.3f, phi = %.3f", jet.first.c_str(), jet.second.Pt(), jet.second.Eta(), jet.second.Phi_0_2pi());
162  Printf("Jet N Consituents = %d", jet.second.fNConstituents);
163  }
164 }
165 
170 {
171  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
172  if (it == fJets.end()) return 0;
173 
174  Double_t z = 0;
175 
176  if ((*it).second.Pt() > 0) {
177  TVector3 dvect = fD.Vect();
178  TVector3 jvect = (*it).second.fMomentum.Vect();
179 
180  Double_t jetMom = jvect * jvect;
181 
182  if (jetMom < 1e-6) {
183  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetZ", "Zero jet momentum!");
184  z = 0.999;
185  }
186  else {
187  z = (dvect * jvect) / jetMom;
188  }
189 
190  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
191  }
192 
193  return z;
194 }
195 
203 {
204  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
205  if (it == fJets.end()) return 0;
206 
207  return GetDistance((*it).second, deta, dphi);
208 }
209 
215 {
216  Double_t deta = 0;
217  Double_t dphi = 0;
218  return GetDistance(n, deta, dphi);
219 }
220 
228 {
229  dphi = TVector2::Phi_mpi_pi(fD.Phi() - jet.Phi());;
230  deta = fD.Eta() - jet.Eta();
231  return TMath::Sqrt(dphi*dphi + deta*deta);
232 }
233 
239 {
240  Double_t deta = 0;
241  Double_t dphi = 0;
242  return GetDistance(jet, deta, dphi);
243 }
244 
250 {
251  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
252  if (it == fJets.end()) {
253  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
254  n.c_str());
255  return 0;
256  }
257  return &((*it).second);
258 }
259 
265 {
266  std::map<std::string, AliJetInfo>::iterator it = fJets.find(n);
267  if (it == fJets.end()) {
268  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
269  n.c_str());
270  return 0;
271  }
272  return &((*it).second);
273 }
274 
275 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoSummary
276 
280 
286  fPt(0),
287  fEta(0),
288  fPhi(0),
289  fR(0),
290  fZ(0)
291 {
292  Set(source, n);
293 }
294 
297 {
298  fPt = 0;
299  fEta = 0;
300  fPhi = 0;
301  fR = 0;
302  fZ = 0;
303 }
304 
310 {
311  std::map<std::string, AliJetInfo>::const_iterator it = source.fJets.find(n);
312  if (it == source.fJets.end()) return;
313 
314  fPt = (*it).second.Pt();
315  fEta = (*it).second.Eta();
316  fPhi = (*it).second.Phi_0_2pi();
317  fR = source.GetDistance(n);
318  fZ = source.GetZ(n);
319 }
320 
325 {
326  fPt = source.Pt();
327  fEta = source.Eta();
328  fPhi = source.Phi_0_2pi();
329  fR = 0;
330  fZ = 0;
331 }
332 
333 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonInfoSummary
334 
338 
343  fPt(0),
344  fEta(0),
345  fPhi(0)
346 {
347  Set(source);
348 }
349 
354 {
355  fPt = source.fD.Pt();
356  fEta = source.fD.Eta();
357  fPhi = source.fD.Phi_0_2pi();
358 }
359 
362 {
363  fPt = 0;
364  fEta = 0;
365  fPhi = 0;
366 }
367 
368 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary
369 
373 
378  AliDmesonInfoSummary(source),
379  fFirstPartonType(0),
380  fFirstPartonPt(0),
381  fLastPartonType(0),
382  fLastPartonPt(0)
383 {
384  Set(source);
385 }
386 
391 {
393 
394  fFirstPartonType = source.fFirstPartonType;
395  fLastPartonType = source.fLastPartonType;
396 
397  if (source.fFirstParton) {
398  fFirstPartonPt = source.fFirstParton->Pt();
399  }
400  else {
401  fFirstPartonPt = 0.;
402  }
403 
404  if (source.fLastParton) {
405  fLastPartonPt = source.fLastParton->Pt();
406  }
407  else {
408  fLastPartonPt = 0.;
409  }
410 }
411 
414 {
416  fFirstPartonType = 0,
417  fFirstPartonPt = 0.;
418  fLastPartonType = 0,
419  fLastPartonPt = 0.;
420 }
421 
422 // Definitions of class AliAnalysisTaskDmesonJets::AliD0InfoSummary
423 
427 
432  AliDmesonInfoSummary(source),
433  fInvMass(source.fD.M()),
434  fSelectionType(0)
435 {
436 }
437 
442 {
443  fInvMass = source.fD.M();
444  fSelectionType = source.fSelectionType;
446 }
447 
450 {
452  fSelectionType = 0;
453  fInvMass = 0;
454 }
455 
456 // Definitions of class AliAnalysisTaskDmesonJets::AliDStarInfoSummary
457 
461 
466  AliDmesonInfoSummary(source),
467  f2ProngInvMass(source.fInvMass2Prong),
468  fDeltaInvMass(source.fD.M() - source.fInvMass2Prong)
469 {
470 }
471 
476 {
477  f2ProngInvMass = source.fInvMass2Prong;
478  fDeltaInvMass = source.fD.M() - source.fInvMass2Prong;
480 }
481 
484 {
486 
487  f2ProngInvMass = 0;
488  fDeltaInvMass = 0;
489 }
490 
491 // Definitions of class AliAnalysisTaskDmesonJets::AliJetDefinition
492 
496 
499  TObject(),
500  fJetType(AliJetContainer::kChargedJet),
501  fRadius(0),
502  fJetAlgo(AliJetContainer::antikt_algorithm),
503  fRecoScheme(AliJetContainer::pt_scheme),
504  fMinJetPt(0.),
505  fMaxJetPt(500.),
506  fMinJetPhi(0.),
507  fMaxJetPhi(0.),
508  fMinJetEta(-1.),
509  fMaxJetEta(1.),
510  fMinChargedPt(0.),
511  fMaxChargedPt(0.),
512  fMinNeutralPt(0.),
513  fMaxNeutralPt(0.)
514 {
515 }
516 
524  TObject(),
525  fJetType(type),
526  fRadius(r),
527  fJetAlgo(algo),
528  fRecoScheme(reco),
529  fMinJetPt(0.),
530  fMaxJetPt(500.),
531  fMinJetPhi(0.),
532  fMaxJetPhi(0.),
533  fMinJetEta(-1.),
534  fMaxJetEta(1.),
535  fMinChargedPt(0.),
536  fMaxChargedPt(0.),
537  fMinNeutralPt(0.),
538  fMaxNeutralPt(0.)
539 {
540 }
541 
546  TObject(),
547  fJetType(source.fJetType),
548  fRadius(source.fRadius),
549  fJetAlgo(source.fJetAlgo),
550  fRecoScheme(source.fRecoScheme),
551  fMinJetPt(source.fMinJetPt),
552  fMaxJetPt(source.fMaxJetPt),
553  fMinJetPhi(source.fMinJetPhi),
554  fMaxJetPhi(source.fMaxJetPhi),
555  fMinJetEta(source.fMinJetEta),
556  fMaxJetEta(source.fMaxJetEta),
557  fMinChargedPt(source.fMinChargedPt),
558  fMaxChargedPt(source.fMaxChargedPt),
559  fMinNeutralPt(source.fMinNeutralPt),
560  fMaxNeutralPt(source.fMaxNeutralPt)
561 {
562 }
563 
568 {
569  new (this) AliHFJetDefinition(source);
570  return *this;
571 }
572 
575 {
576  static TString name;
577 
578  name = AliJetContainer::GenerateJetName(fJetType, fJetAlgo, fRecoScheme, fRadius, 0, 0, "Jet");
579 
580  return name.Data();
581 }
582 
588 {
589  if (fMinJetEta < fMaxJetEta && (jet.Eta() < fMinJetEta || jet.Eta() > fMaxJetEta)) return kFALSE;
590  if (fMinJetPhi < fMaxJetPhi && (jet.Phi() < fMinJetPhi || jet.Phi() > fMaxJetPhi)) return kFALSE;
591  if (jet.Pt() > fMaxJetPt || jet.Pt() < fMinJetPt) return kFALSE;
592  if (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt) return kFALSE;
593  if (jet.fMaxNeutralPt < fMinNeutralPt || jet.fMaxNeutralPt > fMaxNeutralPt) return kFALSE;
594 
595  return kTRUE;
596 }
597 
603 {
604  const AliJetInfo* jet = dMesonJet.GetJet(n);
605  if (!jet) return kFALSE;
606  return IsJetInAcceptance((*jet));
607 }
608 
615 {
616  if (lhs.fJetType > rhs.fJetType) return false;
617  else if (lhs.fJetType < rhs.fJetType) return true;
618  else {
619  if (lhs.fRadius > rhs.fRadius) return false;
620  else if (lhs.fRadius < rhs.fRadius) return true;
621  else {
622  if (lhs.fJetAlgo > rhs.fJetAlgo) return false;
623  else if (lhs.fJetAlgo < rhs.fJetAlgo) return true;
624  else {
625  if (lhs.fRecoScheme < rhs.fRecoScheme) return true;
626  else return false;
627  }
628  }
629  }
630 }
631 
638 {
639  if (lhs.fJetType != rhs.fJetType) return false;
640  if (lhs.fRadius != rhs.fRadius) return false;
641  if (lhs.fJetAlgo != rhs.fJetAlgo) return false;
642  if (lhs.fRecoScheme != rhs.fRecoScheme) return false;
643  return true;
644 }
645 
646 // Definitions of class AliAnalysisTaskDmesonJets::AnalysisEngine
647 
651 
654  TObject(),
655  fFirstPartons(),
656  fLastPartons(),
657  fCandidateType(kD0toKpi),
658  fCandidateName(),
659  fCandidatePDG(0),
660  fNDaughters(0),
661  fPDGdaughters(),
662  fBranchName(),
663  fMCMode(kNoMC),
664  fNMassBins(0),
665  fMinMass(0),
666  fMaxMass(0),
667  fRDHFCuts(0),
668  fRejectedOrigin(0),
669  fAcceptedDecay(0),
670  fInhibit(kFALSE),
671  fJetDefinitions(),
672  fPtBinWidth(0.5),
673  fMaxPt(100),
674  fDataSlotNumber(-1),
675  fTree(0),
676  fCurrentDmesonJetInfo(0),
677  fCurrentJetInfo(0),
678  fCandidateArray(0),
679  fMCContainer(0),
680  fTrackContainer(0),
681  fClusterContainer(0),
682  fAodEvent(0),
683  fFastJetWrapper(0),
684  fHistManager(0)
685 {
686 }
687 
696  TObject(),
697  fFirstPartons(),
698  fLastPartons(),
699  fCandidateType(type),
700  fCandidateName(),
701  fCandidatePDG(0),
702  fNDaughters(0),
703  fPDGdaughters(),
704  fBranchName(),
705  fMCMode(MCmode),
706  fNMassBins(nMassBins),
707  fMinMass(0),
708  fMaxMass(0),
709  fRDHFCuts(cuts),
710  fRejectedOrigin(0),
711  fAcceptedDecay(kAnyDecay),
712  fInhibit(kFALSE),
713  fJetDefinitions(),
714  fPtBinWidth(0.5),
715  fMaxPt(100),
716  fDataSlotNumber(-1),
717  fTree(0),
718  fCurrentDmesonJetInfo(0),
719  fCurrentJetInfo(0),
720  fCandidateArray(0),
721  fMCContainer(0),
722  fTrackContainer(0),
723  fClusterContainer(0),
724  fAodEvent(0),
725  fFastJetWrapper(0),
726  fHistManager(0)
727 {
728  SetCandidateProperties(range);
729 }
730 
735  TObject(source),
736  fFirstPartons(source.fFirstPartons),
737  fLastPartons(source.fLastPartons),
738  fCandidateType(source.fCandidateType),
739  fCandidateName(source.fCandidateName),
740  fCandidatePDG(source.fCandidatePDG),
741  fNDaughters(source.fNDaughters),
742  fPDGdaughters(source.fPDGdaughters),
743  fBranchName(source.fBranchName),
744  fMCMode(source.fMCMode),
745  fNMassBins(source.fNMassBins),
746  fMinMass(source.fMinMass),
747  fMaxMass(source.fMaxMass),
748  fRDHFCuts(),
749  fRejectedOrigin(source.fRejectedOrigin),
750  fAcceptedDecay(source.fAcceptedDecay),
751  fInhibit(source.fInhibit),
752  fJetDefinitions(source.fJetDefinitions),
753  fPtBinWidth(source.fPtBinWidth),
754  fMaxPt(source.fMaxPt),
755  fDataSlotNumber(-1),
756  fTree(0),
757  fCurrentDmesonJetInfo(0),
758  fCurrentJetInfo(0),
759  fCandidateArray(source.fCandidateArray),
760  fMCContainer(source.fMCContainer),
761  fTrackContainer(source.fTrackContainer),
762  fClusterContainer(source.fClusterContainer),
763  fAodEvent(source.fAodEvent),
765  fHistManager(source.fHistManager)
766 {
767  SetRDHFCuts(source.fRDHFCuts);
768 }
769 
770 // Destructor
772 {
773  delete fRDHFCuts;
774 }
775 
780 {
781  new (this) AnalysisEngine(source);
782  return *this;
783 }
784 
789 {
790  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
791  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
792  }
793 
794  return kFALSE;
795 }
796 
798 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const /*geom*/, Int_t /*runNumber*/)
799 {
800 }
801 
806 {
807  switch (fCandidateType) {
808  case kD0toKpi:
809  fCandidatePDG = 421;
810  fCandidateName = "D0";
811  fNDaughters = 2;
812  fPDGdaughters.Set(fNDaughters);
813  fPDGdaughters.Reset();
814  fPDGdaughters[0] = 211; // pi
815  fPDGdaughters[1] = 321; // K
816  fBranchName = "D0toKpi";
817  fAcceptedDecay = kDecayD0toKpi;
818  if (!fRDHFCuts) {
819  fRDHFCuts = new AliRDHFCutsD0toKpi();
820  fRDHFCuts->SetStandardCutsPP2010();
821  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
822  fRDHFCuts->SetTriggerClass("","");
823  }
824  break;
825  case kD0toKpiLikeSign:
826  fCandidatePDG = 421;
827  fCandidateName = "2ProngLikeSign";
828  fNDaughters = 2;
829  fPDGdaughters.Set(fNDaughters);
830  fPDGdaughters.Reset();
831  fPDGdaughters[0] = 211; // pi
832  fPDGdaughters[1] = 321; // K
833  fBranchName = "LikeSign2Prong";
834  if (!fRDHFCuts) {
835  fRDHFCuts = new AliRDHFCutsD0toKpi();
836  fRDHFCuts->SetStandardCutsPP2010();
837  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
838  fRDHFCuts->SetTriggerClass("","");
839  }
840  break;
841  case kDstartoKpipi:
842  fCandidatePDG = 413;
843  fCandidateName = "DStar";
844  fNDaughters = 3;
845  fPDGdaughters.Set(fNDaughters);
846  fPDGdaughters.Reset();
847  fPDGdaughters[0] = 211; // pi soft
848  fPDGdaughters[1] = 211; // pi fromD0
849  fPDGdaughters[2] = 321; // K from D0
850  fBranchName = "Dstar";
851  fAcceptedDecay = kDecayDStartoKpipi;
852  if (!fRDHFCuts) {
853  fRDHFCuts = new AliRDHFCutsDStartoKpipi();
854  fRDHFCuts->SetStandardCutsPP2010();
855  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
856  fRDHFCuts->SetTriggerClass("","");
857  }
858  break;
859  default:
860  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
861  }
862 
863  CalculateMassLimits(range, fCandidatePDG, fNMassBins, fMinMass, fMaxMass);
864 }
865 
870 {
871  if (fRDHFCuts) delete fRDHFCuts;
872  fRDHFCuts = cuts;
873 }
874 
879 {
880  if (!cuts) return;
881  if (fRDHFCuts) delete fRDHFCuts;
882  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
883 }
884 
889 {
890  static TString name;
891 
892  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
893 
894  return name.Data();
895 }
896 
901 {
902  static TString name;
903 
904  name = fCandidateName;
905  switch (fMCMode) {
906  case kBackgroundOnly:
907  name += "_kBackgroundOnly";
908  break;
909  case kSignalOnly:
910  name += "_kSignalOnly";
911  break;
912  case kMCTruth:
913  name += "_MCTruth";
914  break;
915  default:
916  break;
917  }
918 
919  return name.Data();
920 }
921 
929 {
930  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
931 
932  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
933  fJetDefinitions.push_back(def);
934  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
935  "Total number of jet definitions is now %lu.",
936  def.GetName(), GetName(), fJetDefinitions.size());
937  // For detector level set maximum track pt to 100 GeV/c
938  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
939  }
940  else {
941  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
942  }
943 
944  return &(*it);
945 }
946 
958 {
959  AliHFJetDefinition def(type, r, algo, reco);
960 
961  return AddJetDefinition(def);
962 }
963 
969 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
970 {
971  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
972  while (it != fJetDefinitions.end() && (*it) != def) it++;
973  return it;
974 }
975 
982 {
983  if (lhs.fCandidateType > rhs.fCandidateType) return false;
984  else if (lhs.fCandidateType < rhs.fCandidateType) return true;
985  else {
986  if (lhs.fMCMode < rhs.fMCMode) return true;
987  else return false;
988  }
989 }
990 
997 {
998  if (lhs.fCandidateType != rhs.fCandidateType) return false;
999  if (lhs.fMCMode != rhs.fMCMode) return false;
1000  return true;
1001 }
1002 
1011 {
1012  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) { // D0 candidate
1013  return ExtractD0Attributes(Dcand, DmesonJet, i);
1014  }
1015  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
1016  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
1017  }
1018  else {
1019  return kFALSE;
1020  }
1021 }
1022 
1031 {
1032  AliDebug(10,"Checking if D0 meson is selected");
1033  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
1034  if (isSelected == 0) return kFALSE;
1035 
1036  Int_t MCtruthPdgCode = 0;
1037 
1038  Double_t invMassD = 0;
1039 
1040  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
1041  // Checks also the origin, and if it matches the rejected origin mask, return false
1042  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1043  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
1044  DmesonJet.fMCLabel = mcLab;
1045 
1046  // Retrieve the generated particle (if exists) and its PDG code
1047  if (mcLab >= 0) {
1048  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1049 
1050  if (aodMcPart) {
1051  // Check origin and return false if it matches the rejected origin mask
1052  if (fRejectedOrigin) {
1053  auto origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
1054  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1055  }
1056  MCtruthPdgCode = aodMcPart->PdgCode();
1057  }
1058  }
1059  }
1060 
1061  if (isSelected == 1) { // selected as a D0
1062  if (i != 0) return kFALSE; // only one mass hypothesis thanks to PID
1063 
1064  if (fMCMode == kNoMC ||
1065  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1066  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly)) {
1067  // 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)
1068  AliDebug(10,"Selected as D0");
1069  invMassD = Dcand->InvMassD0();
1070  }
1071  else { // conditions above not passed, so return FALSE
1072  return kFALSE;
1073  }
1074  }
1075  else if (isSelected == 2) { // selected as a D0bar
1076  if (i != 1) return kFALSE; // only one mass hypothesis thanks to PID
1077 
1078  if (fMCMode == kNoMC ||
1079  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1080  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly)) {
1081  // 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)
1082  AliDebug(10,"Selected as D0bar");
1083  invMassD = Dcand->InvMassD0bar();
1084  }
1085  else { // conditions above not passed, so return FALSE
1086  return kFALSE;
1087  }
1088  }
1089  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
1090  AliDebug(10,"Selected as either D0 or D0bar");
1091 
1092  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
1093  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
1094  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kBackgroundOnly)) {
1095  if (i != 0) return kFALSE;
1096  AliDebug(10, "MC truth is D0");
1097  invMassD = Dcand->InvMassD0();
1098  }
1099  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
1100  (MCtruthPdgCode == fCandidatePDG && fMCMode == kBackgroundOnly)) {
1101  if (i != 1) return kFALSE;
1102  AliDebug(10, "MC truth is D0bar");
1103  invMassD = Dcand->InvMassD0bar();
1104  }
1105  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
1106 
1107  // Only accept it if background-only OR background-and-signal was requested
1108  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
1109  // Select D0 or D0bar depending on the i-parameter
1110  if (i == 0) {
1111  AliDebug(10, "Returning invariant mass with D0 hypothesis");
1112  invMassD = Dcand->InvMassD0();
1113  }
1114  else if (i == 1) {
1115  AliDebug(10, "Returning invariant mass with D0bar hypothesis");
1116  invMassD = Dcand->InvMassD0bar();
1117  }
1118  else { // i > 1
1119  return kFALSE;
1120  }
1121  }
1122  else { // signal-only was requested but this is not a true D0
1123  return kFALSE;
1124  }
1125  }
1126  }
1127  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
1128  return kTRUE;
1129 }
1130 
1139 {
1140  AliDebug(10,"Checking if D* meson is selected");
1141  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoCascadeHF*>(DstarCand), AliRDHFCuts::kAll, fAodEvent);
1142  if (isSelected == 0) return kFALSE;
1143 
1144  if ((i == 1 && DstarCand->Charge()>0) || (i == 0 && DstarCand->Charge()<0) || i > 1) return kFALSE; // only one mass hypothesis for the D*
1145 
1146  Int_t MCtruthPdgCode = 0;
1147 
1148  Double_t invMassD = 0;
1149 
1150  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
1151  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
1152  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
1153 
1154  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
1155  AliDebug(10, Form("MC label is %d", mcLab));
1156  DmesonJet.fMCLabel = mcLab;
1157  if (mcLab >= 0) {
1158  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
1159 
1160  if (aodMcPart) {
1161  if (fRejectedOrigin) {
1162  auto origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
1163  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
1164  }
1165 
1166  MCtruthPdgCode = aodMcPart->PdgCode();
1167  AliDebug(10, Form("MC truth pdg code is %d",MCtruthPdgCode));
1168  }
1169  }
1170  }
1171 
1172  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1173  if (fMCMode == kNoMC ||
1174  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1175  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1176  // 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)
1177  invMassD = DstarCand->InvMassDstarKpipi();
1178  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1179  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1180  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1181  return kTRUE;
1182  }
1183  else { // conditions above not passed, so return FALSE
1184  return kFALSE;
1185  }
1186 }
1187 
1195 {
1196  if (!part) return kUnknownDecay;
1197  if (!mcArray) return kUnknownDecay;
1198 
1200 
1201  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1202 
1203  if (part->GetNDaughters() == 2) {
1204 
1205  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1206  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1207 
1208  if (!d1 || !d2) {
1209  return decay;
1210  }
1211 
1212  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1213  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1214 
1215  if (absPdgPart == 421) { // D0 -> K pi
1216  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1217  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1218  decay = kDecayD0toKpi;
1219  }
1220  }
1221 
1222  if (absPdgPart == 413) { // D* -> D0 pi
1223  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1224  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1225  if (D0decay == kDecayD0toKpi) {
1226  decay = kDecayDStartoKpipi;
1227  }
1228  }
1229  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1230  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1231  if (D0decay == kDecayD0toKpi) {
1232  decay = kDecayDStartoKpipi;
1233  }
1234  }
1235  }
1236  }
1237 
1238  return decay;
1239 }
1240 
1247 std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> AliAnalysisTaskDmesonJets::AnalysisEngine::CheckOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, Bool_t firstParton)
1248 {
1249  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1250 
1251  std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> result(kUnknownQuark, 0);
1252 
1253  if (!part) return result;
1254  if (!mcArray) return result;
1255 
1256  Int_t mother = part->GetMother();
1257  while (mother >= 0) {
1258  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1259  if (mcGranma) {
1260  Int_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
1261 
1262  if (abspdgGranma == 1) result = {kFromDown, mcGranma};
1263  if (abspdgGranma == 2) result = {kFromUp, mcGranma};
1264  if (abspdgGranma == 3) result = {kFromStrange, mcGranma};
1265  if (abspdgGranma == 4) result = {kFromCharm, mcGranma};
1266  if (abspdgGranma == 5) result = {kFromBottom, mcGranma};
1267  if (abspdgGranma == 6) result = {kFromTop, mcGranma};
1268  if (abspdgGranma == 9 || abspdgGranma == 21) result = {kFromGluon, mcGranma};
1269 
1270  // If looking for the very first parton in the hard scattering, it will continue the loop until it cannot find a mother particle
1271  if (result.first != kUnknownQuark && !firstParton) return result;
1272 
1273  mother = mcGranma->GetMother();
1274  }
1275  else {
1276  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1277  break;
1278  }
1279  }
1280 
1281  return result;
1282 }
1283 
1286 {
1287  for (auto& jetDef : fJetDefinitions) {
1288  jetDef.fJets.clear();
1289  }
1290 
1291  if (fMCMode == kMCTruth) {
1292  RunParticleLevelAnalysis();
1293  }
1294  else {
1295  RunDetectorLevelAnalysis();
1296  }
1297 }
1298 
1304 {
1306  fFastJetWrapper->SetR(jetDef.fRadius);
1309 
1310  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1311  fTrackContainer->SetDMesonCandidate(0);
1312  AddInputVectors(fTrackContainer, 100);
1313  }
1314 
1315  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1316  AddInputVectors(fClusterContainer, -100);
1317  }
1318 
1319  // run jet finder
1320  fFastJetWrapper->Run();
1321 
1322  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1323 
1324  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1325  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1326 
1327  Double_t maxChPt = 0;
1328  Double_t maxNePt = 0;
1329  Double_t totalNeutralPt = 0;
1330 
1331  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1332  if (constituents[ic].user_index() >= 100) {
1333  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1334  }
1335  else if (constituents[ic].user_index() <= -100) {
1336  totalNeutralPt += constituents[ic].pt();
1337  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1338  }
1339  }
1340 
1341  jetDef.fJets.push_back(
1342  AliJetInfo(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E(),
1343  constituents.size(), totalNeutralPt / jets_incl[ijet].pt(), maxChPt, maxNePt)
1344  );
1345  }
1346 }
1347 
1356 {
1357  if (jetDef.fJets.size() == 0) FindJets(jetDef);
1358 
1359  Double_t d_closest = 999;
1360  AliJetInfo* jet_closest = 0;
1361  const AliJetInfo* truth_jet = 0;
1362  try {
1363  truth_jet = &(dmeson.fJets.at(jetDef.GetName()));
1364  }
1365  catch(...) {
1366  return jet_distance_pair(0, 999);
1367  }
1368 
1369  for (auto& jet : jetDef.fJets) {
1370  Double_t d = truth_jet->GetDistance(jet);
1371 
1372  if (d > dMax) continue;
1373  if (d < d_closest) {
1374  d_closest = d;
1375  jet_closest = &jet;
1376  }
1377  }
1378 
1379  if (jet_closest && applyKinCuts) {
1380  if (!jetDef.IsJetInAcceptance(*jet_closest)) jet_closest = 0;
1381  }
1382 
1383  if (jet_closest) {
1384  AliDebug(2, Form("Found closest jet (pt=%.3f, eta=%.3f, phi=%.3f) to d meson (pt=%.3f, eta=%.3f, phi=%.3f) with d = %.3f",
1385  jet_closest->Pt(), jet_closest->Eta(), jet_closest->Phi_0_2pi(),
1386  dmeson.fD.Pt(), dmeson.fD.Eta(), dmeson.fD.Phi_0_2pi(),
1387  d_closest));
1388  }
1389 
1390  return jet_distance_pair(jet_closest, d_closest);
1391 }
1392 
1395 {
1396  const Int_t nD = fCandidateArray->GetEntriesFast();
1397 
1398  AliDmesonJetInfo DmesonJet;
1399 
1400  Int_t nAccCharm[3] = {0};
1401  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1402  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1403  if (!charmCand) continue;
1404 
1405  // region of interest + cuts
1406  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1407  Int_t nMassHypo = 0; // number of mass hypothesis accepted for this D meson
1408  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1409  DmesonJet.Reset();
1410  DmesonJet.fDmesonParticle = charmCand;
1411  DmesonJet.fSelectionType = im + 1;
1412  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1413  for (auto& def : fJetDefinitions) {
1414  if (!FindJet(charmCand, DmesonJet, def)) {
1415  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1416  def.GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1417  }
1418  }
1419  fDmesonJets[(icharm+1)*(1-(im*2))] = DmesonJet;
1420  nMassHypo++;
1421  nAccCharm[im]++;
1422  }
1423  }
1424  if (nMassHypo == 2) {
1425  nAccCharm[0]--;
1426  nAccCharm[1]--;
1427  nAccCharm[2] += 2;
1428  }
1429  if (nMassHypo == 2) { // both mass hypothesis accepted
1430  fDmesonJets[(icharm+1)].fSelectionType = 3;
1431  fDmesonJets[-(icharm+1)].fSelectionType = 3;
1432  }
1433  } // end of D cand loop
1434 
1435  TString hname;
1436 
1437  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1438  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1439  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1440  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1441 
1442  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1443  fHistManager->FillTH2(hname, fTrackContainer->GetNAcceptedTracks(), nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1444 
1445  hname = TString::Format("%s/fHistNDmesons", GetName());
1446  fHistManager->FillTH1(hname, nD);
1447 }
1448 
1459 {
1460  TString hname;
1461 
1463  fFastJetWrapper->SetR(jetDef.fRadius);
1466 
1467  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1468 
1469  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1470  fTrackContainer->SetDMesonCandidate(Dcand);
1471  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1472  AddInputVectors(fTrackContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1473 
1474  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1475  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1476  const TObjArray& daughters = fTrackContainer->GetDaughterList();
1477  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
1478  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
1479  if (!fTrackContainer->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
1480  }
1481  }
1482 
1483  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1484  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1485  AddInputVectors(fClusterContainer, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1486  }
1487 
1488  // run jet finder
1489  fFastJetWrapper->Run();
1490 
1491  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1492 
1493  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1494  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1495 
1496  Bool_t isDmesonJet = kFALSE;
1497 
1498  Double_t maxChPt = 0;
1499  Double_t maxNePt = 0;
1500  Double_t totalNeutralPt = 0;
1501 
1502  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1503  if (constituents[ic].user_index() == 0) {
1504  isDmesonJet = kTRUE;
1505  }
1506  else if (constituents[ic].user_index() >= 100) {
1507  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1508  }
1509  else if (constituents[ic].user_index() <= -100) {
1510  totalNeutralPt += constituents[ic].pt();
1511  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1512  }
1513  }
1514 
1515  if (isDmesonJet) {
1516  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1517  DmesonJet.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1518  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1519  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1520  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1521 
1522  return kTRUE;
1523  }
1524  }
1525 
1526  return kFALSE;
1527 }
1528 
1532 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist)
1533 {
1534  auto itcont = cont->all_momentum();
1535  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
1536  UInt_t rejectionReason = 0;
1537  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
1538  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
1539  continue;
1540  }
1541  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
1542  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
1543  }
1544 }
1545 
1548 {
1549  TString hname;
1550 
1551  fMCContainer->SetSpecialPDG(fCandidatePDG);
1552  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1553  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1554 
1555  if (!fMCContainer->IsSpecialPDGFound()) return;
1556 
1557  Int_t nAccCharm[3] = {0};
1558 
1559  for (auto &jetDef : fJetDefinitions) {
1560  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", GetName(), jetDef.GetName());
1561  TH1* histNDmesonsVsNconstituents = static_cast<TH1*>(fHistManager->FindObject(hname));
1562 
1563  switch (jetDef.fJetType) {
1565  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
1566  break;
1568  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
1569  break;
1571  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
1572  break;
1573  }
1574 
1576  fFastJetWrapper->SetR(jetDef.fRadius);
1579 
1580  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
1581  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1582 
1583  fFastJetWrapper->Run();
1584 
1585  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1586 
1587  for (auto jet : jets_incl) {
1588  Int_t nDmesonsInJet = 0;
1589 
1590  for (auto constituent : jet.constituents()) {
1591  Int_t iPart = constituent.user_index() - 100;
1592  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
1593  if (!part) {
1594  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1595  continue;
1596  }
1597  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
1598  nDmesonsInJet++;
1599  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
1600  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
1601  std::pair<int, AliDmesonJetInfo> element;
1602  element.first = iPart;
1603  dMesonJetIt = fDmesonJets.insert(element).first;
1604  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1605  (*dMesonJetIt).second.fDmesonParticle = part;
1606  (*dMesonJetIt).second.fSelectionType = part->PdgCode() > 0 ? 1 : 2;
1607 
1608  UShort_t p = 0;
1609  UInt_t rs = 0;
1610 
1611  auto firstParton = CheckOrigin(part, fMCContainer->GetArray(), kTRUE);
1612  p = 0;
1613  rs = firstParton.first;
1614  while (rs >>= 1) { p++; }
1615  (*dMesonJetIt).second.fFirstPartonType = p;
1616  (*dMesonJetIt).second.fFirstParton = firstParton.second;
1617 
1618  auto lastParton = CheckOrigin(part, fMCContainer->GetArray(), kFALSE);
1619  p = 0;
1620  rs = lastParton.first;
1621  while (rs >>= 1) { p++; }
1622  (*dMesonJetIt).second.fLastPartonType = p;
1623  (*dMesonJetIt).second.fLastParton = lastParton.second;
1624 
1625  if (part->PdgCode() > 0) {
1626  nAccCharm[0]++;
1627  }
1628  else {
1629  nAccCharm[1]++;
1630  }
1631  }
1632 
1633  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
1634  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = jet.constituents().size();
1635  } // if constituent is a D meson
1636  } // for each constituent
1637  if (nDmesonsInJet > 0) histNDmesonsVsNconstituents->Fill(jet.constituents().size(), nDmesonsInJet);
1638  } // for each jet
1639  } // for each jet definition
1640 
1641  if (fDmesonJets.size() != nAccCharm[0]+nAccCharm[1]) AliError(Form("I found %lu mesons (%d)?", fDmesonJets.size(), nAccCharm[0]+nAccCharm[1]));
1642  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
1643  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
1644  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
1645  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
1646 
1647  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
1648  fHistManager->FillTH2(hname, fMCContainer->GetNAcceptedParticles(), nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
1649 
1650  hname = TString::Format("%s/fHistNDmesons", GetName());
1651  fHistManager->FillTH1(hname, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]); // same as the number of accepted D mesons, since no selection is performed
1652 }
1653 
1658 {
1659  TString classname;
1660  if (fMCMode == kMCTruth) {
1661  classname = "AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary";
1662  fCurrentDmesonJetInfo = new AliDmesonMCInfoSummary();
1663  }
1664  else {
1665  switch (fCandidateType) {
1666  case kD0toKpi:
1667  case kD0toKpiLikeSign:
1668  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1669  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1670  break;
1671  case kDstartoKpipi:
1672  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1673  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1674  break;
1675  }
1676  }
1677  TString treeName = TString::Format("%s_%s", taskName, GetName());
1678  fTree = new TTree(treeName, treeName);
1679  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
1680  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1681  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1682  fCurrentJetInfo[i] = new AliJetInfoSummary();
1683  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1684  }
1685 
1686  return fTree;
1687 }
1688 
1693 {
1694  TString hname;
1695 
1696  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
1697 
1698  for (auto &jetDef : fJetDefinitions) {
1699 
1700  AliDebug(2,Form("Now working on '%s'", jetDef.GetName()));
1701 
1702  Double_t radius = jetDef.fRadius;
1703 
1704  TString title[30] = {""};
1705  Int_t nbins[30] = {0 };
1706  Double_t min [30] = {0.};
1707  Double_t max [30] = {0.};
1708  Int_t dim = 0 ;
1709 
1710  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1711  nbins[dim] = nPtBins;
1712  min[dim] = 0;
1713  max[dim] = fMaxPt;
1714  dim++;
1715 
1716  if ((enabledAxis & kPositionD) != 0) {
1717  title[dim] = "#eta_{D}";
1718  nbins[dim] = 50;
1719  min[dim] = -1;
1720  max[dim] = 1;
1721  dim++;
1722 
1723  title[dim] = "#phi_{D} (rad)";
1724  nbins[dim] = 150;
1725  min[dim] = 0;
1726  max[dim] = TMath::TwoPi();
1727  dim++;
1728  }
1729 
1730  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1731  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1732  nbins[dim] = fNMassBins;
1733  min[dim] = fMinMass;
1734  max[dim] = fMaxMass;
1735  dim++;
1736  }
1737 
1738  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1739  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1740  nbins[dim] = fNMassBins;
1741  min[dim] = fMinMass;
1742  max[dim] = fMaxMass;
1743  dim++;
1744  }
1745 
1746  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1747  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1748  nbins[dim] = fNMassBins;
1749  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1750  dim++;
1751  }
1752 
1753  if (fCandidateType == kDstartoKpipi) {
1754  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1755  nbins[dim] = fNMassBins*6;
1756  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1757 
1758  // subtract mass of D0
1759  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1760  min[dim] -= D0mass;
1761  max[dim] -= D0mass;
1762 
1763  dim++;
1764  }
1765 
1766  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1767  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1768  nbins[dim] = 100;
1769  min[dim] = 0;
1770  max[dim] = 25;
1771  dim++;
1772  }
1773 
1774  title[dim] = "#it{z}_{D}";
1775  nbins[dim] = 110;
1776  min[dim] = 0;
1777  max[dim] = 1.10;
1778  dim++;
1779 
1780  if ((enabledAxis & kDeltaR) != 0) {
1781  title[dim] = "#Delta R_{D-jet}";
1782  nbins[dim] = 100;
1783  min[dim] = 0;
1784  max[dim] = radius * 1.5;
1785  dim++;
1786  }
1787 
1788  if ((enabledAxis & kDeltaEta) != 0) {
1789  title[dim] = "#eta_{D} - #eta_{jet}";
1790  nbins[dim] = 100;
1791  min[dim] = -radius * 1.2;
1792  max[dim] = radius * 1.2;
1793  dim++;
1794  }
1795 
1796  if ((enabledAxis & kDeltaPhi) != 0) {
1797  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1798  nbins[dim] = 100;
1799  min[dim] = -radius * 1.2;
1800  max[dim] = radius * 1.2;
1801  dim++;
1802  }
1803 
1804  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1805  nbins[dim] = nPtBins;
1806  min[dim] = 0;
1807  max[dim] = fMaxPt;
1808  dim++;
1809 
1810  if ((enabledAxis & kPositionJet) != 0) {
1811  title[dim] = "#eta_{jet}";
1812  nbins[dim] = 50;
1813  min[dim] = -1;
1814  max[dim] = 1;
1815  dim++;
1816 
1817  title[dim] = "#phi_{jet} (rad)";
1818  nbins[dim] = 150;
1819  min[dim] = 0;
1820  max[dim] = TMath::TwoPi();
1821  dim++;
1822  }
1823 
1824  if ((enabledAxis & kJetConstituents) != 0) {
1825  title[dim] = "No. of constituents";
1826  nbins[dim] = 50;
1827  min[dim] = -0.5;
1828  max[dim] = 49.5;
1829  dim++;
1830  }
1831 
1832  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
1833  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1834  for (Int_t j = 0; j < dim; j++) {
1835  h->GetAxis(j)->SetTitle(title[j]);
1836  }
1837  }
1838 }
1839 
1844 {
1845  TString hname;
1846  fFirstPartons.clear();
1847  fLastPartons.clear();
1848  for (auto& dmeson_pair : fDmesonJets) {
1849  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
1850  Int_t accJets = 0;
1851  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1852  fCurrentJetInfo[ij]->Reset();
1853  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1854  if (!jet) continue;
1855  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1856  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1857  fHistManager->FillTH1(hname, jet->Pt());
1858  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1859  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1860  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1861  fHistManager->FillTH1(hname, jet->Eta());
1862  continue;
1863  }
1864  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions[ij].GetName());
1865  accJets++;
1866  }
1867  if (accJets > 0) {
1868  fFirstPartons[dmeson_pair.second.fFirstParton] = dmeson_pair.second.fFirstPartonType;
1869  fLastPartons[dmeson_pair.second.fLastParton] = dmeson_pair.second.fLastPartonType;
1870 
1871  fTree->Fill();
1872  }
1873  else {
1874  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1875  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1876  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1877  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1878  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1879  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1880  if (fMCMode != kMCTruth) {
1881  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1882  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1883  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1884  }
1885  else if (fCandidateType == kDstartoKpipi) {
1886  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1887  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1888 
1889  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1890  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1891  }
1892  }
1893  }
1894  }
1895 
1896  hname = TString::Format("%s/fHistFirstPartonPt", GetName());
1897  TH1* histFirstPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
1898  hname = TString::Format("%s/fHistFirstPartonEta", GetName());
1899  TH1* histFirstPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
1900  hname = TString::Format("%s/fHistFirstPartonPhi", GetName());
1901  TH1* histFirstPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
1902  hname = TString::Format("%s/fHistFirstPartonType", GetName());
1903  TH1* histFirstPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
1904 
1905  for (auto parton : fFirstPartons) {
1906  if (!parton.first) continue;
1907  histFirstPartonPt->Fill(parton.first->Pt());
1908  histFirstPartonEta->Fill(parton.first->Eta());
1909  histFirstPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
1910  histFirstPartonType->Fill(parton.second);
1911  }
1912 
1913  hname = TString::Format("%s/fHistLastPartonPt", GetName());
1914  TH1* histLastPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
1915  hname = TString::Format("%s/fHistLastPartonEta", GetName());
1916  TH1* histLastPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
1917  hname = TString::Format("%s/fHistLastPartonPhi", GetName());
1918  TH1* histLastPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
1919  hname = TString::Format("%s/fHistLastPartonType", GetName());
1920  TH1* histLastPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
1921 
1922  for (auto parton : fLastPartons) {
1923  if (!parton.first) continue;
1924  histLastPartonPt->Fill(parton.first->Pt());
1925  histLastPartonEta->Fill(parton.first->Eta());
1926  histLastPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
1927  histLastPartonType->Fill(parton.second);
1928  }
1929 
1930  return kTRUE;
1931 }
1932 
1938 {
1939  TString hname;
1940  fFirstPartons.clear();
1941  fLastPartons.clear();
1942  for (auto& dmeson_pair : fDmesonJets) {
1943  Int_t accJets = 0;
1944  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1945  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions[ij].GetName());
1946  if (!jet) continue;
1947  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1948  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1949  fHistManager->FillTH1(hname, jet->Pt());
1950  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1951  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1952  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1953  fHistManager->FillTH1(hname, jet->Eta());
1954  continue;
1955  }
1956  accJets++;
1957  }
1958  if (accJets > 0) {
1959  fFirstPartons[dmeson_pair.second.fFirstParton] = dmeson_pair.second.fFirstPartonType;
1960  fLastPartons[dmeson_pair.second.fLastParton] = dmeson_pair.second.fLastPartonType;
1961  }
1962  else {
1963  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1964  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1965  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1966  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1967  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1968  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1969  if (fMCMode != kMCTruth) {
1970  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) {
1971  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1972  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1973  }
1974  else if (fCandidateType == kDstartoKpipi) {
1975  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1976  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1977 
1978  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1979  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1980  }
1981  }
1982  }
1983  }
1984 
1985  hname = TString::Format("%s/fHistFirstPartonPt", GetName());
1986  TH1* histFirstPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
1987  hname = TString::Format("%s/fHistFirstPartonEta", GetName());
1988  TH1* histFirstPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
1989  hname = TString::Format("%s/fHistFirstPartonPhi", GetName());
1990  TH1* histFirstPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
1991  hname = TString::Format("%s/fHistFirstPartonType", GetName());
1992  TH1* histFirstPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
1993 
1994  for (auto parton : fFirstPartons) {
1995  if (!parton.first) continue;
1996  histFirstPartonPt->Fill(parton.first->Pt());
1997  histFirstPartonEta->Fill(parton.first->Eta());
1998  histFirstPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
1999  histFirstPartonType->Fill(parton.second);
2000  }
2001 
2002  hname = TString::Format("%s/fHistLastPartonPt", GetName());
2003  TH1* histLastPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
2004  hname = TString::Format("%s/fHistLastPartonEta", GetName());
2005  TH1* histLastPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
2006  hname = TString::Format("%s/fHistLastPartonPhi", GetName());
2007  TH1* histLastPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
2008  hname = TString::Format("%s/fHistLastPartonType", GetName());
2009  TH1* histLastPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
2010 
2011  for (auto parton : fLastPartons) {
2012  if (!parton.first) continue;
2013  histLastPartonPt->Fill(parton.first->Pt());
2014  histLastPartonEta->Fill(parton.first->Eta());
2015  histLastPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
2016  histLastPartonType->Fill(parton.second);
2017  }
2018 
2019  return kTRUE;
2020 }
2021 
2026 {
2027  TString hname;
2028 
2029  for (auto& dmeson_pair : fDmesonJets) {
2030  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
2031  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
2032  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
2033  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
2034  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
2035  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
2036  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
2037  }
2038  }
2039 
2040  for (auto &jetDef : fJetDefinitions) {
2041 
2042  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef.GetName());
2043  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
2044 
2045  for (auto& dmeson_pair : fDmesonJets) {
2046  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
2047  if (!jet) continue;
2048  if (!jetDef.IsJetInAcceptance(*jet)) {
2049  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef.GetName());
2050  fHistManager->FillTH1(hname, jet->Pt());
2051  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef.GetName());
2052  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
2053  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef.GetName());
2054  fHistManager->FillTH1(hname, jet->Eta());
2055  continue;
2056  }
2057  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
2058  }
2059  }
2060 
2061  return kTRUE;
2062 }
2063 
2070 {
2071  // Fill the THnSparse histogram.
2072 
2073  Double_t contents[30] = {0.};
2074 
2075  Double_t z = DmesonJet.GetZ(n);
2076  Double_t deltaPhi = 0;
2077  Double_t deltaEta = 0;
2078  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
2079 
2080  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
2081  if (it == DmesonJet.fJets.end()) return kFALSE;
2082 
2083  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
2084  TString title(h->GetAxis(i)->GetTitle());
2085  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
2086  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
2087  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
2088  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
2089  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
2090  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
2091  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
2092  else if (title=="#it{z}_{D}") contents[i] = z;
2093  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
2094  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
2095  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
2096  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
2097  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
2098  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
2099  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
2100  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
2101  }
2102 
2103  h->Fill(contents);
2104 
2105  return kTRUE;
2106 }
2107 
2108 // Definitions of class AliAnalysisTaskDmesonJets
2109 
2113 
2117  fAnalysisEngines(),
2118  fEnabledAxis(0),
2120  fHistManager(),
2121  fApplyKinematicCuts(kTRUE),
2122  fNOutputTrees(0),
2123  fAodEvent(0),
2124  fFastJetWrapper(0),
2125  fMCContainer(0)
2126 {
2127  SetMakeGeneralHistograms(kTRUE);
2128 }
2129 
2134  AliAnalysisTaskEmcalLight(name, kTRUE),
2135  fAnalysisEngines(),
2136  fEnabledAxis(k2ProngInvMass),
2137  fOutputType(kTreeOutput),
2138  fHistManager(name),
2139  fApplyKinematicCuts(kTRUE),
2140  fNOutputTrees(nOutputTrees),
2141  fAodEvent(0),
2142  fFastJetWrapper(0),
2143  fMCContainer(0)
2144 {
2145  SetMakeGeneralHistograms(kTRUE);
2146  for (Int_t i = 0; i < nOutputTrees; i++){
2147  DefineOutput(2+i, TTree::Class());
2148  }
2149 }
2150 
2153 {
2154  if (fFastJetWrapper) delete fFastJetWrapper;
2155 }
2156 
2164 {
2165  AliRDHFCuts* analysiscuts = 0;
2166  TFile* filecuts = TFile::Open(cutfname);
2167  if (!filecuts || filecuts->IsZombie()) {
2168  ::Warning("AddTaskDmesonJets", "Input file not found: will use std cuts.");
2169  filecuts = 0;
2170  }
2171 
2172  if (filecuts) {
2173  analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
2174  if (!analysiscuts) {
2175  ::Warning("AddTaskDmesonJetCorr", "Could not find analysis cuts '%s' in '%s'. Using std cuts.", cutsname.Data(), cutfname.Data());
2176  }
2177  }
2178 
2179  return analysiscuts;
2180 }
2181 
2191 {
2193  return AddAnalysisEngine(type, MCmode, jetDef, cutfname);
2194 }
2195 
2205 {
2206  AliRDHFCuts* cuts = 0;
2207 
2208  if (!cutfname.IsNull()) {
2209  TString cutsname;
2210 
2211  switch (type) {
2212  case kD0toKpi:
2213  case kD0toKpiLikeSign:
2214  cutsname = "D0toKpiCuts";
2215  break;
2216  case kDstartoKpipi:
2217  cutsname = "DStartoKpipiCuts";
2218  break;
2219  default:
2220  return 0;
2221  }
2222 
2223  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
2224  }
2225 
2226  AnalysisEngine eng(type, MCmode, cuts);
2227 
2228  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
2229 
2230  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
2231  eng.AddJetDefinition(jetDef);
2232  it = fAnalysisEngines.insert(it, eng);
2233  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
2234  }
2235  else {
2236  AnalysisEngine* found_eng = &(*it);
2237  ::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());
2238  found_eng->AddJetDefinition(jetDef);
2239 
2240  if (cuts && found_eng->fRDHFCuts != 0) {
2241  ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
2242  found_eng->SetRDHFCuts(cuts);
2243  }
2244  }
2245 
2246  return &(*it);
2247 }
2248 
2249 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
2250 {
2251  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
2252  while (it != fAnalysisEngines.end() && (*it) < eng) it++;
2253  return it;
2254 }
2255 
2258 {
2259  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
2260 
2262 
2263  // Define histograms
2264  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
2265 
2266  TString hname;
2267  TString htitle;
2268  TH1* h = 0;
2269  Int_t treeSlot = 0;
2270 
2271  hname = "fHistCharmPt";
2272  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2273  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2274 
2275  hname = "fHistCharmEta";
2276  htitle = hname + ";#eta_{charm};counts";
2277  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2278 
2279  hname = "fHistCharmPhi";
2280  htitle = hname + ";#phi_{charm};counts";
2281  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2282 
2283  hname = "fHistCharmPt_Eta05";
2284  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
2285  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2286 
2287  hname = "fHistBottomPt";
2288  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2289  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2290 
2291  hname = "fHistBottomEta";
2292  htitle = hname + ";#eta_{bottom};counts";
2293  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2294 
2295  hname = "fHistBottomPhi";
2296  htitle = hname + ";#phi_{bottom};counts";
2297  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2298 
2299  hname = "fHistBottomPt_Eta05";
2300  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2301  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2302 
2303  hname = "fHistHighestPartonPt";
2304  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
2305  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2306 
2307  hname = "fHistHighestPartonType";
2308  htitle = hname + ";type;counts";
2309  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2310 
2311  hname = "fHistNHeavyQuarks";
2312  htitle = hname + ";number of heavy-quarks;counts";
2313  fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
2314 
2315  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
2316  for (auto &param : fAnalysisEngines) {
2317  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
2318 
2319  param.fHistManager = &fHistManager;
2320 
2321  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", param.GetName());
2322  htitle = hname + ";#it{N}_{tracks};#it{N}_{D};events";
2323  h = fHistManager.CreateTH2(hname, htitle, 251, -0.5, 250.5, 21, -0.5, 20.5);
2324 
2325  hname = TString::Format("%s/fHistNTotAcceptedDmesons", param.GetName());
2326  htitle = hname + ";;#it{N}_{D}";
2327  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
2328 
2329  hname = TString::Format("%s/fHistNDmesons", param.GetName());
2330  htitle = hname + ";#it{N}_{D};events";
2331  h = fHistManager.CreateTH1(hname, htitle, 501, -0.5, 500.5);
2332 
2333  hname = TString::Format("%s/fHistNEvents", param.GetName());
2334  htitle = hname + ";Event status;counts";
2335  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
2336  h->GetXaxis()->SetBinLabel(1, "Accepted");
2337  h->GetXaxis()->SetBinLabel(2, "Rejected");
2338 
2339  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
2340  htitle = hname + ";Rejection reason;counts";
2341  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
2342 
2343  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
2344  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
2345  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2346 
2347  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
2348  htitle = hname + ";#it{#eta}_{D};counts";
2349  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2350 
2351  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
2352  htitle = hname + ";#it{#phi}_{D};counts";
2353  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2354 
2355  if (param.fMCMode != kMCTruth) {
2356  if (param.fCandidateType == kD0toKpi || param.fCandidateType == kD0toKpiLikeSign) {
2357  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
2358  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2359  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
2360  }
2361  else if (param.fCandidateType == kDstartoKpipi) {
2362  Double_t min = 0;
2363  Double_t max = 0;
2364 
2365  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
2366  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2367  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
2368  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
2369 
2370  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2371  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
2372  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
2373  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
2374  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
2375  }
2376  }
2377 
2378  if (param.fMCMode == kMCTruth) {
2379  hname = TString::Format("%s/fHistFirstPartonPt", param.GetName());
2380  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
2381  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2382 
2383  hname = TString::Format("%s/fHistFirstPartonEta", param.GetName());
2384  htitle = hname + ";#eta_{parton};counts";
2385  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2386 
2387  hname = TString::Format("%s/fHistFirstPartonPhi", param.GetName());
2388  htitle = hname + ";#phi_{parton};counts";
2389  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2390 
2391  hname = TString::Format("%s/fHistFirstPartonType", param.GetName());
2392  htitle = hname + ";type;counts";
2393  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2394 
2395  hname = TString::Format("%s/fHistLastPartonPt", param.GetName());
2396  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
2397  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
2398 
2399  hname = TString::Format("%s/fHistLastPartonEta", param.GetName());
2400  htitle = hname + ";#eta_{parton};counts";
2401  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
2402 
2403  hname = TString::Format("%s/fHistLastPartonPhi", param.GetName());
2404  htitle = hname + ";#phi_{parton};counts";
2405  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
2406 
2407  hname = TString::Format("%s/fHistLastPartonType", param.GetName());
2408  htitle = hname + ";type;counts";
2409  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
2410  }
2411 
2412  for (auto& jetDef : param.fJetDefinitions) {
2413  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef.GetName());
2414 
2415  if (param.fMCMode == kMCTruth) {
2416  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", param.GetName(), jetDef.GetName());
2417  htitle = hname + ";#it{N}_{constituents};#it{N}_{D};counts";
2418  h = fHistManager.CreateTH2(hname, htitle, 51, -0.5, 50.5, 10, 0.5, 10.5);
2419  }
2420 
2421  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef.GetName());
2422  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2423  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2424  SetRejectionReasonLabels(h->GetXaxis());
2425 
2426  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef.GetName());
2427  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
2428  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2429  SetRejectionReasonLabels(h->GetXaxis());
2430 
2431  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef.GetName());
2432  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
2433  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
2434  SetRejectionReasonLabels(h->GetXaxis());
2435 
2436  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef.GetName());
2437  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
2438  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
2439 
2440  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef.GetName());
2441  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
2442  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
2443 
2444  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef.GetName());
2445  htitle = hname + ";#it{#eta}_{jet};counts";
2446  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
2447 
2448  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef.GetName());
2449  htitle = hname + ";#it{#phi}_{jet};counts";
2450  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
2451  }
2452  switch (fOutputType) {
2453  case kTreeOutput:
2454  param.BuildTree(GetName());
2455  if (treeSlot < fNOutputTrees) {
2456  param.AssignDataSlot(treeSlot+2);
2457  treeSlot++;
2459  }
2460  else {
2461  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
2462  }
2463  break;
2464  case kTHnOutput:
2465  param.BuildHnSparse(fEnabledAxis);
2466  break;
2467  case kNoOutput:
2468  break;
2469  }
2470  }
2471 
2473 
2474  PostData(1, fOutput);
2475 }
2476 
2480 {
2482 
2483  // Load the event
2484  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
2485 
2486  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
2487 
2488  fFastJetWrapper->SetAreaType(fastjet::active_area);
2490 
2491  if (!fAodEvent) {
2492  AliError(Form("This task need an AOD event (Task '%s'). Expect troubles...", GetName()));
2493  //return;
2494  }
2495 
2498 
2499  for (auto &params : fAnalysisEngines) {
2500 
2501  params.fAodEvent = fAodEvent;
2502  params.fFastJetWrapper = fFastJetWrapper;
2503  if (fAodEvent) params.Init(fGeom, fAodEvent->GetRunNumber());
2504 
2505  if (params.fMCMode != kMCTruth && fAodEvent) {
2506  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
2507 
2508  if (params.fCandidateArray) {
2509  TString className;
2510  if (params.fCandidateType == kD0toKpi || params.fCandidateType == kD0toKpiLikeSign) {
2511  className = "AliAODRecoDecayHF2Prong";
2512  }
2513  else if (params.fCandidateType == kDstartoKpipi) {
2514  className = "AliAODRecoCascadeHF";
2515  }
2516  if (!params.fCandidateArray->GetClass()->InheritsFrom(className)) {
2517  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2518  "%s: Objects of type %s in %s are not inherited from %s! Task will be disabled!",
2519  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName(), className.Data());
2520  params.fCandidateArray = 0;
2521  params.fInhibit = kTRUE;
2522  }
2523  }
2524  else {
2525  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2526  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
2527  params.fBranchName.Data(), params.GetName());
2528  params.fInhibit = kTRUE;
2529  }
2530  }
2531 
2532  if (params.fMCMode != kNoMC) {
2533  if (!fMCContainer) {
2534  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2535  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
2536  params.GetName());
2537  params.fInhibit = kTRUE;
2538  }
2539  params.fMCContainer = fMCContainer;
2540  }
2541 
2542  if (params.fMCMode != kMCTruth) {
2543  params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(0));
2544  if (!params.fTrackContainer) params.fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(1));
2545 
2546  params.fClusterContainer = GetClusterContainer(0);
2547 
2548  if (!params.fTrackContainer && !params.fClusterContainer) {
2549  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2550  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
2551  params.GetName());
2552  params.fInhibit = kTRUE;
2553  }
2554  }
2555  }
2556 }
2557 
2562 {
2563  TString hname;
2564 
2565  // Fix for temporary bug in ESDfilter
2566  // The AODs with null vertex pointer didn't pass the PhysSel
2567  // Now adding an entry in the histogram so as to check that this is actually cutting anything out
2568  if (fAodEvent && (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001)) {
2569  for (auto &eng : fAnalysisEngines) {
2570  if (eng.fInhibit) continue;
2571  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2572  fHistManager.FillTH1(hname, "ESDfilterBug");
2573  }
2574  return kFALSE;
2575  }
2576 
2577  for (auto &eng : fAnalysisEngines) {
2578  eng.fDmesonJets.clear();
2579  if (eng.fInhibit) continue;
2580 
2581  //Event selection
2582  hname = TString::Format("%s/fHistNEvents", eng.GetName());
2583  if (fAodEvent) {
2584  Bool_t iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
2585  if (!iseventselected) {
2586  fHistManager.FillTH1(hname, "Rejected");
2587  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
2588  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
2589  TString label;
2590  do {
2591  label = GetHFEventRejectionReasonLabel(bitmap);
2592  if (label.IsNull()) break;
2593  fHistManager.FillTH1(hname, label);
2594  } while (true);
2595 
2596  continue;
2597  }
2598  }
2599 
2600  fHistManager.FillTH1(hname, "Accepted");
2601 
2602  AliDebug(2, "Event selected");
2603 
2604  eng.RunAnalysis();
2605  }
2606  return kTRUE;
2607 }
2608 
2613 {
2614  for (auto &param : fAnalysisEngines) {
2615  if (param.fInhibit) continue;
2616 
2617  if (fOutputType == kTreeOutput) {
2618  param.FillTree(fApplyKinematicCuts);
2620  }
2621  else if (fOutputType == kTHnOutput) {
2622  param.FillHnSparse(fApplyKinematicCuts);
2623  }
2624  }
2626  return kTRUE;
2627 }
2628 
2631 {
2632  auto itcont = fMCContainer->all_momentum();
2633  Int_t nHQ = 0;
2634  Double_t highestPartonPt = 0;
2635  Int_t absPdgHighParton = 0;
2636  for (auto part : itcont) {
2637  Int_t absPdgCode = TMath::Abs(part.second->GetPdgCode());
2638 
2639  // Skip all particles that are not either quarks or gluons
2640  if (absPdgCode > 9 && absPdgCode != 21) continue;
2641 
2642  // Look for highest momentum parton
2643  if (highestPartonPt < part.first.Pt()) {
2644  highestPartonPt = part.first.Pt();
2645  absPdgHighParton = absPdgCode;
2646  }
2647  /*
2648  // Look for the mother PDG code
2649  Int_t motherIndex = part.second->GetMother();
2650  AliAODMCParticle *mother = 0;
2651  Int_t motherPdg = 0;
2652  Double_t motherPt = 0;
2653  if (motherIndex >= 0) {
2654  mother = fMCContainer->GetMCParticle(motherIndex);
2655  if (motherIndex) {
2656  motherPdg = TMath::Abs(mother->GetPdgCode());
2657  motherPt = mother->Pt();
2658  }
2659  }
2660  */
2661  if (absPdgCode != 4 && absPdgCode != 5) continue;
2662  Bool_t notLastInPartonShower = kFALSE;
2663  for (Int_t idaugh = 0; idaugh < 2; idaugh++){
2664  Int_t daughterIndex = part.second->GetDaughter(idaugh);
2665  if (daughterIndex < 0) {
2666  AliDebug(10, Form("Could not find daughter of heavy quark (pdg=%d, pt=%.3f)!", absPdgCode, part.first.Pt()));
2667  continue;
2668  }
2669  AliAODMCParticle *daughter = fMCContainer->GetMCParticle(daughterIndex);
2670  if (!daughter) {
2671  AliDebug(10, Form("Could not find daughter %d of heavy quark (pdg=%d, pt=%.3f)!", daughterIndex, absPdgCode, part.first.Pt()));
2672  continue;
2673  }
2674  Int_t daughterAbsPdgCode = TMath::Abs(daughter->GetPdgCode());
2675  if (daughterAbsPdgCode <= 9 || daughterAbsPdgCode == 21) notLastInPartonShower = kTRUE; // this parton is not the last parton in the shower
2676  AliDebug(10, Form("Found daughter with PDG=%d, pt=%.3f", daughterAbsPdgCode, daughter->Pt()));
2677  }
2678  if (notLastInPartonShower) continue;
2679 
2680  if (absPdgCode == 4) {
2681  fHistManager.FillTH1("fHistCharmPt", part.first.Pt());
2682  fHistManager.FillTH1("fHistCharmEta", part.first.Eta());
2683  fHistManager.FillTH1("fHistCharmPhi", part.first.Phi_0_2pi());
2684  if (TMath::Abs(part.first.Eta()) < 0.5) fHistManager.FillTH1("fHistCharmPt_Eta05", part.first.Pt());
2685  }
2686  else if (absPdgCode == 5) {
2687  fHistManager.FillTH1("fHistBottomPt", part.first.Pt());
2688  fHistManager.FillTH1("fHistBottomEta", part.first.Eta());
2689  fHistManager.FillTH1("fHistBottomPhi", part.first.Phi_0_2pi());
2690  if (TMath::Abs(part.first.Eta()) < 0.5) fHistManager.FillTH1("fHistBottomPt_Eta05", part.first.Pt());
2691  }
2692  nHQ++;
2693  }
2694  fHistManager.FillTH1("fHistNHeavyQuarks", nHQ);
2695  fHistManager.FillTH1("fHistHighestPartonPt",highestPartonPt);
2696  Int_t partonType = 0;
2697  if (absPdgHighParton == 9 || absPdgHighParton == 21) {
2698  partonType = 7;
2699  }
2700  else {
2701  partonType = absPdgHighParton;
2702  }
2703  fHistManager.FillTH1("fHistHighestPartonType",partonType);
2704 }
2705 
2714 {
2715  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
2716 
2717  Double_t mass = part->Mass();
2718 
2719  // To make sure that the PDG mass value is not at the edge of a bin
2720  if (nbins % 2 == 0) {
2721  minMass = mass - range / 2 - range / nbins / 2;
2722  maxMass = mass + range / 2 - range / nbins / 2;
2723  }
2724  else {
2725  minMass = mass - range / 2;
2726  maxMass = mass + range / 2;
2727  }
2728 }
2729 
2736 {
2737  static TString label;
2738  label = "";
2739 
2740  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
2741  label = "NotSelTrigger";
2742  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
2743  return label.Data();
2744  }
2745  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
2746  label = "NoVertex";
2747  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
2748  return label.Data();
2749  }
2750  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
2751  label = "TooFewVtxContrib";
2752  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
2753  return label.Data();
2754  }
2755  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
2756  label = "ZVtxOutFid";
2757  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
2758  return label.Data();
2759  }
2760  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
2761  label = "Pileup";
2762  bitmap &= ~BIT(AliRDHFCuts::kPileup);
2763  return label.Data();
2764  }
2765  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
2766  label = "OutsideCentrality";
2767  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
2768  return label.Data();
2769  }
2770  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
2771  label = "PhysicsSelection";
2772  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
2773  return label.Data();
2774  }
2775  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
2776  label = "BadSPDVertex";
2777  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
2778  return label.Data();
2779  }
2780  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
2781  label = "ZVtxSPDOutFid";
2782  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
2783  return label.Data();
2784  }
2785  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
2786  label = "CentralityFlattening";
2787  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
2788  return label.Data();
2789  }
2790  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
2791  label = "BadTrackV0Correl";
2792  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
2793  return label.Data();
2794  }
2795 
2796  return label.Data();
2797 }
2798 
2805 {
2806  if (eng.GetDataSlotNumber() >= 0 && eng.GetTree()) {
2807  PostData(eng.GetDataSlotNumber(), eng.GetTree());
2808  return eng.GetDataSlotNumber();
2809  }
2810  else {
2811  return -1;
2812  }
2813 }
2814 
2824 {
2825  // Get the pointer to the existing analysis manager via the static access method
2826  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
2827  if (!mgr) {
2828  ::Error("AddTaskDmesonJets", "No analysis manager to connect to.");
2829  return NULL;
2830  }
2831 
2832  // Check the analysis type using the event handlers connected to the analysis manager
2833  AliVEventHandler* handler = mgr->GetInputEventHandler();
2834  if (!handler) {
2835  ::Error("AddTaskEmcalJetSpectraQA", "This task requires an input event handler");
2836  return NULL;
2837  }
2838 
2839  EDataType_t dataType = kUnknownDataType;
2840 
2841  if (handler->InheritsFrom("AliESDInputHandler")) {
2842  dataType = kESD;
2843  }
2844  else if (handler->InheritsFrom("AliAODInputHandler")) {
2845  dataType = kAOD;
2846  }
2847 
2848  // Init the task and do settings
2849  if (ntracks == "usedefault") {
2850  if (dataType == kESD) {
2851  ntracks = "Tracks";
2852  }
2853  else if (dataType == kAOD) {
2854  ntracks = "tracks";
2855  }
2856  else {
2857  ntracks = "";
2858  }
2859  }
2860 
2861  if (nclusters == "usedefault") {
2862  if (dataType == kESD) {
2863  nclusters = "CaloClusters";
2864  }
2865  else if (dataType == kAOD) {
2866  nclusters = "caloClusters";
2867  }
2868  else {
2869  nclusters = "";
2870  }
2871  }
2872 
2873  if (nMCpart == "usedefault") {
2874  nMCpart = "mcparticles"; // Always needs AliAODMCParticle objects
2875  }
2876 
2877  TString name("AliAnalysisTaskDmesonJets");
2878  if (strcmp(suffix, "") != 0) {
2879  name += TString::Format("_%s", suffix.Data());
2880  }
2881 
2882  AliAnalysisTaskDmesonJets* jetTask = new AliAnalysisTaskDmesonJets(name, nMaxTrees);
2883  jetTask->SetVzRange(-10,10);
2884 
2885  if (!ntracks.IsNull()) {
2886  AliHFTrackContainer* trackCont = new AliHFTrackContainer(ntracks);
2887  jetTask->AdoptParticleContainer(trackCont);
2888  }
2889 
2890  if (!nMCpart.IsNull()) {
2891  AliMCParticleContainer* partCont = new AliHFAODMCParticleContainer(nMCpart);
2892  partCont->SetEtaLimits(-1.5, 1.5);
2893  partCont->SetPtLimits(0, 1000);
2894  jetTask->AdoptParticleContainer(partCont);
2895  }
2896 
2897  jetTask->AddClusterContainer(nclusters);
2898 
2899  // Final settings, pass to manager and set the containers
2900  mgr->AddTask(jetTask);
2901 
2902  // Create containers for input/output
2903  AliAnalysisDataContainer* cinput1 = mgr->GetCommonInputContainer();
2904  TString contname1(name);
2905  contname1 += "_histos";
2906  AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(contname1.Data(),
2907  TList::Class(), AliAnalysisManager::kOutputContainer,
2908  Form("%s", AliAnalysisManager::GetCommonFileName()));
2909 
2910  mgr->ConnectInput(jetTask, 0, cinput1);
2911  mgr->ConnectOutput(jetTask, 1, coutput1);
2912 
2913  for (Int_t i = 0; i < nMaxTrees; i++) {
2914  TString contname = TString::Format("%s_tree_%d", name.Data(), i);
2915  AliAnalysisDataContainer *coutput = mgr->CreateContainer(contname.Data(),
2916  TTree::Class(),AliAnalysisManager::kOutputContainer,
2917  Form("%s", AliAnalysisManager::GetCommonFileName()));
2918  mgr->ConnectOutput(jetTask, 2+i, coutput);
2919  }
2920  return jetTask;
2921 }
2922 
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
EDataType_t
Switch for the data type.
virtual AliAODMCParticle * GetMCParticle(Int_t i=-1) const
const char * title
Definition: MakeQAPdf.C:26
AliRDHFCuts * fRDHFCuts
D meson candidates cuts.
virtual void Set(const AliDmesonJetInfo &source)
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
AliHFAODMCParticleContainer * fMCContainer
! MC particle container
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)
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
virtual void UserCreateOutputObjects()
Creates the output containers.
void FillPartonLevelHistograms()
Fill histograms with parton-level information.
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.
AliAODMCParticle * fFirstParton
! pointer to the first parton in the shower tree of the D meson (only for particle level D mesons) ...
EMCMode_t fMCMode
MC mode: No MC (data and MC detector level), background-only (MC), signal-only (MC), MC truth (particle level)
Short_t fLastPartonType
! type of the last parton in the shower tree (only for particle level D mesons)
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) ...
const Int_t nPtBins
virtual void Set(const AliDmesonJetInfo &source)
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:33
AliClusterContainer * AddClusterContainer(const char *n)
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
Short_t fFirstPartonType
! type of the first parton in the shower tree (only for particle level D mesons)
const AliMCParticleIterableMomentumContainer all_momentum() const
virtual ~AliAnalysisTaskDmesonJets()
This is the standard destructor.
Byte_t fSelectionType
! for D0: 0=not selected, 1=D0, 2=D0bar, 3=both
static std::pair< AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle * > CheckOrigin(const AliAODMCParticle *part, TClonesArray *mcArray, Bool_t firstParton=kFALSE)
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
AliDmesonJetInfo & operator=(const AliDmesonJetInfo &source)
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.
void SetVzRange(Double_t min, Double_t max)
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)
Int_t GetNAcceptedParticles() const
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)
Container for MC-true particles within the EMCAL framework.
Lightweight class that encapsulates D0.
Bool_t IsAnyJetInAcceptance(const AliDmesonJetInfo &dMesonJet) const
void AdoptParticleContainer(AliParticleContainer *cont)
void SetCharge(EChargeCut_t c)
Container for jet within the EMCAL jet framework.
Definition: External.C:196
AliAODMCParticle * fLastParton
! pointer to the last parton in the shower tree of the D meson (only for particle level D mesons) ...
Int_t PostDataFromAnalysisEngine(const AnalysisEngine &eng)
AliFJWrapper * fFastJetWrapper
! Fastjet wrapper
THistManager fHistManager
Histogram manager.
void Reset()
Reset all fields to their default values.
static AliAnalysisTaskDmesonJets * AddTaskDmesonJets(TString ntracks="usedefault", TString nclusters="usedefault", TString nMCpart="", Int_t nMaxTrees=2, TString suffix="")
virtual void Set(const AliDmesonJetInfo &source, std::string n)