AliPhysics  8bb951a (8bb951a)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskDmesonJets.cxx
Go to the documentation of this file.
1 /*************************************************************************
2 * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15 
16 // Root
17 #include <TClonesArray.h>
18 #include <TDatabasePDG.h>
19 #include <TParticlePDG.h>
20 #include <TVector3.h>
21 #include <THnSparse.h>
22 #include <TParticle.h>
23 #include <TMath.h>
24 #include <THashList.h>
25 #include <TFile.h>
26 
27 // Aliroot general
28 #include "AliLog.h"
29 #include "AliEMCALGeometry.h"
30 
31 // Aliroot HF
32 #include "AliAODRecoCascadeHF.h"
34 #include "AliRDHFCutsD0toKpi.h"
37 #include "AliHFTrackContainer.h"
38 
39 // Aliroot EMCal jet framework
40 #include "AliEmcalJetTask.h"
41 #include "AliEmcalJet.h"
42 #include "AliJetContainer.h"
43 #include "AliParticleContainer.h"
44 #include "AliEmcalParticle.h"
45 #include "AliFJWrapper.h"
46 
48 
49 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonJetInfo
50 
53 {
54  fD.SetPtEtaPhiE(0,0,0,0);
55  fSoftPionPt = 0;
56  fInvMass2Prong = 0;
57  for (std::map<std::string, AliJetInfo>::iterator it = fJets.begin(); it != fJets.end(); it++) {
58  (*it).second.fMomentum.SetPtEtaPhiE(0,0,0,0);
59  (*it).second.fNConstituents = 0;
60  (*it).second.fNEF = 0;
61  (*it).second.fMaxChargedPt = 0;
62  (*it).second.fMaxNeutralPt = 0;
63  }
64 }
65 
68 {
69  Printf("Printing D Meson Jet object.");
70  Printf("D Meson: pT = %.3f, eta = %.3f, phi = %.3f, inv. mass = %.3f", fD.Pt(), fD.Eta(), fD.Phi_0_2pi(), fD.M());
71  Printf("Soft pion pT: %.3f. 2-Prong Invariant mass = %.3f", fSoftPionPt, fInvMass2Prong);
72  for (std::map<std::string, AliJetInfo>::const_iterator it = fJets.begin(); it != fJets.end(); it++) {
73  Printf("Jet %s: pT = %.3f, eta = %.3f, phi = %.3f", (*it).first.c_str(), (*it).second.Pt(), (*it).second.Eta(), (*it).second.Phi_0_2pi());
74  Printf("Jet N Consituents = %d", (*it).second.fNConstituents);
75  }
76 }
77 
82 {
83  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
84  if (it == fJets.end()) return 0;
85 
86  Double_t z = 0;
87 
88  if ((*it).second.Pt() > 0) {
89  TVector3 dvect = fD.Vect();
90  TVector3 jvect = (*it).second.fMomentum.Vect();
91 
92  Double_t jetMom = jvect * jvect;
93 
94  if (jetMom < 1e-6) {
95  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetZ", "Zero jet momentum!");
96  z = 0.999;
97  }
98  else {
99  z = (dvect * jvect) / jetMom;
100  }
101 
102  if (z == 1 || (z > 1 && z - 1 < 1e-3)) z = 0.999; // so that it will contribute to the bin 0.9-1 rather than 1-1.1
103  }
104 
105  return z;
106 }
107 
113 Double_t AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetDistance(std::string n, Double_t& deta, Double_t& dphi) const
114 {
115  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
116  if (it == fJets.end()) return 0;
117 
118  dphi = TVector2::Phi_mpi_pi(fD.Phi() - (*it).second.Phi());;
119  deta = fD.Eta() - (*it).second.Eta();
120  return TMath::Sqrt(dphi*dphi + deta*deta);
121 }
122 
127 {
128  Double_t deta = 0;
129  Double_t dphi = 0;
130  return GetDistance(n, deta, dphi);
131 }
132 
134 {
135  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
136  if (it == fJets.end()) {
137  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!", n.c_str());
138  return 0;
139  }
140  return &((*it).second);
141 }
142 
144 {
145  std::map<std::string, AliJetInfo>::iterator it = fJets.find(n);
146  if (it == fJets.end()) {
147  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!", n.c_str());
148  return 0;
149  }
150  return &((*it).second);
151 }
152 
153 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoSummary
154 
158 
164  fPt(0),
165  fEta(0),
166  fPhi(0),
167  fR(0),
168  fZ(0)
169 {
170  Set(source, n);
171 }
172 
175 {
176  fPt = 0;
177  fEta = 0;
178  fPhi = 0;
179  fR = 0;
180  fZ = 0;
181 }
182 
188 {
189  std::map<std::string, AliJetInfo>::const_iterator it = source.fJets.find(n);
190  if (it == source.fJets.end()) return;
191 
192  fPt = (*it).second.Pt();
193  fEta = (*it).second.Eta();
194  fPhi = (*it).second.Phi_0_2pi();
195  fR = source.GetDistance(n);
196  fZ = source.GetZ(n);
197 }
198 
199 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonInfoSummary
200 
204 
209  fPt(0),
210  fEta(0),
211  fPhi(0)
212 {
213  Set(source);
214 }
215 
220 {
221  fPt = source.fD.Pt();
222  fEta = source.fD.Eta();
223  fPhi = source.fD.Phi_0_2pi();
224 }
225 
226 // Definitions of class AliAnalysisTaskDmesonJets::AliD0InfoSummary
227 
231 
236  AliDmesonInfoSummary(source),
237  fInvMass(source.fD.M())
238 {
239 }
240 
242 {
243  fInvMass = source.fD.M();
245 }
246 
247 // Definitions of class AliAnalysisTaskDmesonJets::AliDStarInfoSummary
248 
252 
257  AliDmesonInfoSummary(source),
258  f2ProngInvMass(source.fInvMass2Prong),
259  fDeltaInvMass(source.fD.M() - source.fInvMass2Prong)
260 {
261 }
262 
264 {
265  f2ProngInvMass = source.fInvMass2Prong;
266  fDeltaInvMass = source.fD.M() - source.fInvMass2Prong;
268 }
269 
270 // Definitions of class AliAnalysisTaskDmesonJets::AliJetDefinition
271 
275 
278  TObject(),
279  fJetType(AliJetContainer::kChargedJet),
280  fRadius(0),
281  fJetAlgo(AliJetContainer::antikt_algorithm),
282  fRecoScheme(AliJetContainer::pt_scheme),
283  fAcceptance(AliJetContainer::kUser),
284  fMinJetPt(0.),
285  fMinJetPhi(0.),
286  fMaxJetPhi(0.),
287  fMinJetEta(0.),
288  fMaxJetEta(0.),
289  fMinChargedPt(0.),
290  fMaxChargedPt(0.),
291  fMinNeutralPt(0.),
292  fMaxNeutralPt(0.)
293 {
294 }
295 
303  TObject(),
304  fJetType(type),
305  fRadius(r),
306  fJetAlgo(algo),
307  fRecoScheme(reco),
308  fAcceptance(AliJetContainer::kUser),
309  fMinJetPt(0.),
310  fMinJetPhi(0.),
311  fMaxJetPhi(0.),
312  fMinJetEta(0.),
313  fMaxJetEta(0.),
314  fMinChargedPt(0.),
315  fMaxChargedPt(0.),
316  fMinNeutralPt(0.),
317  fMaxNeutralPt(0.)
318 {
319  // By default set detector fiducial acceptance
320  switch (type) {
324  break;
327  break;
328  }
329 }
330 
335  TObject(),
336  fJetType(source.fJetType),
337  fRadius(source.fRadius),
338  fJetAlgo(source.fJetAlgo),
339  fRecoScheme(source.fRecoScheme),
340  fAcceptance(source.fAcceptance),
341  fMinJetPt(source.fMinJetPt),
342  fMinJetPhi(source.fMinJetPhi),
343  fMaxJetPhi(source.fMaxJetPhi),
344  fMinJetEta(source.fMinJetEta),
345  fMaxJetEta(source.fMaxJetEta),
346  fMinChargedPt(source.fMinChargedPt),
347  fMaxChargedPt(source.fMaxChargedPt),
348  fMinNeutralPt(source.fMinNeutralPt),
349  fMaxNeutralPt(source.fMaxNeutralPt)
350 {
351 }
352 
357 {
358  new (this) AliHFJetDefinition(source);
359  return *this;
360 }
361 
364 {
365  static TString name;
366 
367  name = AliJetContainer::GenerateJetName(fJetType, fJetAlgo, fRecoScheme, fRadius, 0, 0, "Jet");
368 
369  return name.Data();
370 }
371 
377 {
378  if (fMinJetEta < fMaxJetEta && (jet.Eta() < fMinJetEta || jet.Eta() > fMaxJetEta)) return kFALSE;
379  if (fMinJetPhi < fMaxJetPhi && (jet.Phi() < fMinJetPhi || jet.Phi() > fMaxJetPhi)) return kFALSE;
380  if (jet.Pt() < fMinJetPt) return kFALSE;
381  if (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt) return kFALSE;
382  if (jet.fMaxNeutralPt < fMinNeutralPt || jet.fMaxNeutralPt > fMaxNeutralPt) return kFALSE;
383 
384  return kTRUE;
385 }
386 
392 {
393  const AliJetInfo* jet = dMesonJet.GetJet(n);
394  if (!jet) return kFALSE;
395  return IsJetInAcceptance((*jet));
396 }
397 
402 void AliAnalysisTaskDmesonJets::AliHFJetDefinition::SetDetectorJetEtaPhiRange(const AliEMCALGeometry* const geom, Int_t run)
403 {
404  Double_t r = 0;
405  switch (fAcceptance) {
407  r = fRadius;
408  // enforce fiducial acceptance
409  /* no break */
411  SetJetEtaRange(-0.9 + r, 0.9 - r);
412  SetJetPhiRange(0, 0); // No cut on phi
413  break;
414 
416  r = fRadius;
417  // enforce fiducial acceptance
418  /* no break */
420  if (geom) {
421  SetJetEtaRange(geom->GetArm1EtaMin() + r, geom->GetArm1EtaMax() - r);
422 
423  if(run>=177295 && run<=197470) {//small SM masked in 2012 and 2013
424  SetJetPhiRange(1.405 + r,3.135 - r);
425  }
426  else {
427  SetJetPhiRange(geom->GetArm1PhiMin() * TMath::DegToRad() + r, geom->GetEMCALPhiMax() * TMath::DegToRad() - r);
428  }
429  }
430  else {
431  AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for EMCAL year 2011!!");
432  SetJetEtaRange(-0.7 + r, 0.7 - r);
433  SetJetPhiRange(1.405 + r, 3.135 - r);
434  }
435  break;
436 
438  r = fRadius;
439  // enforce fiducial acceptance
440  /* no break */
442  if (geom) {
443  SetJetEtaRange(geom->GetArm1EtaMin() + r, geom->GetArm1EtaMax() - r);
444  SetJetPhiRange(geom->GetDCALPhiMin() * TMath::DegToRad() + r, geom->GetDCALPhiMax() * TMath::DegToRad() - r);
445  }
446  else {
447  AliWarning("Could not get instance of AliEMCALGeometry. Using manual settings for DCAL year 2015!!");
448  SetJetEtaRange(-0.7 + r, 0.7 - r);
449  SetJetPhiRange(4.538 + r, 5.727 - r);
450  }
451  break;
452 
454  // Nothing to be done
455  break;
456  }
457 }
458 
465 {
466  if (lhs.fJetType > rhs.fJetType) return false;
467  else if (lhs.fJetType < rhs.fJetType) return true;
468  else {
469  if (lhs.fRadius > rhs.fRadius) return false;
470  else if (lhs.fRadius < rhs.fRadius) return true;
471  else {
472  if (lhs.fJetAlgo > rhs.fJetAlgo) return false;
473  else if (lhs.fJetAlgo < rhs.fJetAlgo) return true;
474  else {
475  if (lhs.fRecoScheme < rhs.fRecoScheme) return true;
476  else return false;
477  }
478  }
479  }
480 }
481 
488 {
489  if (lhs.fJetType != rhs.fJetType) return false;
490  if (lhs.fRadius != rhs.fRadius) return false;
491  if (lhs.fJetAlgo != rhs.fJetAlgo) return false;
492  if (lhs.fRecoScheme != rhs.fRecoScheme) return false;
493  return true;
494 }
495 
496 // Definitions of class AliAnalysisTaskDmesonJets::AnalysisEngine
497 
501 
504  TObject(),
505  fCandidateType(kD0toKpi),
506  fCandidateName(),
507  fCandidatePDG(0),
508  fNDaughters(0),
509  fBranchName(),
510  fMCMode(kNoMC),
511  fNMassBins(0),
512  fMinMass(0),
513  fMaxMass(0),
514  fRDHFCuts(0),
515  fRejectedOrigin(0),
516  fAcceptedDecay(0),
517  fInhibit(kFALSE),
518  fJetDefinitions(),
519  fTree(0),
520  fCurrentDmesonJetInfo(0),
521  fCurrentJetInfo(0),
522  fCandidateArray(0),
523  fMCContainer(0),
524  fTrackContainer(0),
525  fClusterContainer(0),
526  fAodEvent(0),
527  fFastJetWrapper(0),
528  fHistManager(0)
529 {
530 }
531 
540  TObject(),
541  fCandidateType(type),
542  fCandidateName(),
543  fCandidatePDG(0),
544  fNDaughters(0),
545  fBranchName(),
546  fMCMode(MCmode),
547  fNMassBins(nMassBins),
548  fMinMass(0),
549  fMaxMass(0),
550  fRDHFCuts(cuts),
551  fRejectedOrigin(kUnknownQuark | kFromBottom),
552  fAcceptedDecay(kAnyDecay),
553  fInhibit(kFALSE),
554  fJetDefinitions(),
555  fTree(0),
556  fCurrentDmesonJetInfo(0),
557  fCurrentJetInfo(0),
558  fCandidateArray(0),
559  fMCContainer(0),
560  fTrackContainer(0),
561  fClusterContainer(0),
562  fAodEvent(0),
563  fFastJetWrapper(0),
564  fHistManager(0)
565 {
566  SetCandidateProperties(range);
567 }
568 
573  TObject(source),
574  fCandidateType(source.fCandidateType),
575  fCandidateName(source.fCandidateName),
576  fCandidatePDG(source.fCandidatePDG),
577  fNDaughters(source.fNDaughters),
578  fBranchName(source.fBranchName),
579  fMCMode(source.fMCMode),
580  fNMassBins(source.fNMassBins),
581  fMinMass(source.fMinMass),
582  fMaxMass(source.fMaxMass),
583  fRDHFCuts(),
584  fRejectedOrigin(source.fRejectedOrigin),
585  fAcceptedDecay(source.fAcceptedDecay),
586  fInhibit(source.fInhibit),
587  fJetDefinitions(source.fJetDefinitions),
588  fTree(0),
589  fCurrentDmesonJetInfo(0),
590  fCurrentJetInfo(0),
591  fCandidateArray(source.fCandidateArray),
592  fMCContainer(source.fMCContainer),
593  fTrackContainer(source.fTrackContainer),
594  fClusterContainer(source.fClusterContainer),
595  fAodEvent(source.fAodEvent),
597  fHistManager(source.fHistManager)
598 {
599  SetRDHFCuts(source.fRDHFCuts);
600 }
601 
602 // Destructor
604 {
605  if (fRDHFCuts) delete fRDHFCuts;
606 }
607 
612 {
613  new (this) AnalysisEngine(source);
614  return *this;
615 }
616 
621 {
622  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
623  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
624  }
625 
626  return kFALSE;
627 }
628 
630 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const geom, Int_t runNumber)
631 {
632  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
633  fJetDefinitions[i].SetDetectorJetEtaPhiRange(geom, runNumber);
634  }
635 }
636 
641 {
642  switch (fCandidateType) {
643  case kD0toKpi:
644  fCandidatePDG = 421;
645  fCandidateName = "D0";
646  fNDaughters = 2;
647  fPDGdaughters.Set(fNDaughters);
648  fPDGdaughters.Reset();
649  fPDGdaughters[0] = 211; // pi
650  fPDGdaughters[1] = 321; // K
651  fBranchName = "D0toKpi";
652  fAcceptedDecay = kD0toKpi;
653  if (!fRDHFCuts) {
654  fRDHFCuts = new AliRDHFCutsD0toKpi();
655  fRDHFCuts->SetStandardCutsPP2010();
656  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
657  fRDHFCuts->SetTriggerClass("","");
658  }
659  break;
660  case kDstartoKpipi:
661  fCandidatePDG = 413;
662  fCandidateName = "DStar";
663  fNDaughters = 3;
664  fPDGdaughters.Set(fNDaughters);
665  fPDGdaughters.Reset();
666  fPDGdaughters[0] = 211; // pi soft
667  fPDGdaughters[1] = 211; // pi fromD0
668  fPDGdaughters[2] = 321; // K from D0
669  fBranchName = "Dstar";
670  fAcceptedDecay = kDstartoKpipi;
671  if (!fRDHFCuts) {
672  fRDHFCuts = new AliRDHFCutsDStartoKpipi();
673  fRDHFCuts->SetStandardCutsPP2010();
674  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
675  fRDHFCuts->SetTriggerClass("","");
676  }
677  break;
678  default:
679  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
680  }
681 
682  CalculateMassLimits(range, fCandidatePDG, fNMassBins, fMinMass, fMaxMass);
683 }
684 
689 {
690  if (fRDHFCuts) delete fRDHFCuts;
691  fRDHFCuts = cuts;
692 }
693 
698 {
699  if (!cuts) return;
700  if (fRDHFCuts) delete fRDHFCuts;
701  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
702 }
703 
708 {
709  static TString name;
710 
711  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
712 
713  return name.Data();
714 }
715 
720 {
721  static TString name;
722 
723  name = fCandidateName;
724  switch (fMCMode) {
725  case kBackgroundOnly:
726  name += "_kBackgroundOnly";
727  break;
728  case kSignalOnly:
729  name += "_kSignalOnly";
730  break;
731  case kMCTruth:
732  name += "_MCTruth";
733  break;
734  default:
735  break;
736  }
737 
738  return name.Data();
739 }
740 
748 {
749  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
750 
751  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
752  fJetDefinitions.push_back(def);
753  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
754  "Total number of jet definitions is now %lu.",
755  def.GetName(), GetName(), fJetDefinitions.size());
756  // For detector level set maximum track pt to 100 GeV/c
757  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
758  }
759  else {
760  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
761  }
762 
763  return &(*it);
764 }
765 
777 {
778  AliHFJetDefinition def(type, r, algo, reco);
779 
780  return AddJetDefinition(def);
781 }
782 
788 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
789 {
790  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
791  while (it != fJetDefinitions.end() && (*it) != def) it++;
792  return it;
793 }
794 
801 {
802  if (lhs.fCandidateType > rhs.fCandidateType) return false;
803  else if (lhs.fCandidateType < rhs.fCandidateType) return true;
804  else {
805  if (lhs.fMCMode < rhs.fMCMode) return true;
806  else return false;
807  }
808 }
809 
816 {
817  if (lhs.fCandidateType != rhs.fCandidateType) return false;
818  if (lhs.fMCMode != rhs.fMCMode) return false;
819  return true;
820 }
821 
830 {
831  DmesonJet.fD.SetPtEtaPhiM(part->Pt(), part->Eta(), part->Phi(), part->M());
832  return kTRUE;
833 }
834 
843 {
844  if (fCandidateType == kD0toKpi) { // D0 candidate
845  return ExtractD0Attributes(Dcand, DmesonJet, i);
846  }
847  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
848  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
849  }
850  else {
851  return kFALSE;
852  }
853 }
854 
863 {
864  Int_t MCtruthPdgCode = 0;
865 
866  Double_t invMassD = 0;
867 
868  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
869  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
870  if (mcLab >= 0) {
871  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
872 
873  if (aodMcPart) {
874  if (fRejectedOrigin && fMCMode == kSignalOnly) {
875  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
876 
877  if ((origin & fRejectedOrigin) == origin) return kFALSE;
878  }
879  MCtruthPdgCode = aodMcPart->PdgCode();
880  }
881  }
882  }
883 
884  //AliDebug(2,"Checking if D0 meson is selected");
885  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
886  if (isSelected == 1) { // selected as a D0
887  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
888 
889  if (fMCMode == kNoMC ||
890  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
891  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly)) {
892  // 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)
893  //AliDebug(2,"Selected as D0");
894  invMassD = Dcand->InvMassD0();
895  }
896  else { // conditions above not passed, so return FALSE
897  return kFALSE;
898  }
899  }
900  else if (isSelected == 2) { // selected as a D0bar
901  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
902 
903  if (fMCMode == kNoMC ||
904  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
905  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly)) {
906  // 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)
907  //AliDebug(2,"Selected as D0bar");
908  invMassD = Dcand->InvMassD0bar();
909  }
910  else { // conditions above not passed, so return FALSE
911  return kFALSE;
912  }
913  }
914  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
915  //AliDebug(2,"Selected as either D0 or D0bar");
916 
917  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
918  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
919  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kBackgroundOnly)) {
920  if (i > 0) return kFALSE;
921  //AliDebug(2, "MC truth is D0");
922  invMassD = Dcand->InvMassD0();
923  }
924  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
925  (MCtruthPdgCode == fCandidatePDG && fMCMode == kBackgroundOnly)) {
926  if (i > 0) return kFALSE;
927  //AliDebug(2, "MC truth is D0bar");
928  invMassD = Dcand->InvMassD0bar();
929  }
930  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
931 
932  // Only accept it if background-only OR background-and-signal was requested
933  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
934  // Select D0 or D0bar depending on the i-parameter
935  if (i == 0) {
936  //AliDebug(2, "Returning invariant mass with D0 hypothesis");
937  invMassD = Dcand->InvMassD0();
938  }
939  else if (i == 1) {
940  //AliDebug(2, "Returning invariant mass with D0bar hypothesis");
941  invMassD = Dcand->InvMassD0bar();
942  }
943  else { // i > 1
944  return kFALSE;
945  }
946  }
947  else { // signal-only was requested but this is not a true D0
948  return kFALSE;
949  }
950  }
951  }
952  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
953  return kTRUE;
954 }
955 
964 {
965  if (i > 0) return kFALSE; // only one mass hypothesis for the D*
966 
967  Int_t MCtruthPdgCode = 0;
968 
969  Double_t invMassD = 0;
970 
971  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
972  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
973  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
974 
975  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
976  //AliDebug(2, Form("MC label is %d", mcLab));
977  if (mcLab >= 0) {
978  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
979 
980  if (aodMcPart) {
981  if (fRejectedOrigin && fMCMode == kSignalOnly) {
982  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
983 
984  if ((origin & fRejectedOrigin) == origin) return kFALSE;
985  }
986 
987  MCtruthPdgCode = aodMcPart->PdgCode();
988  //AliDebug(2, Form("MC truth pdg code is %d",MCtruthPdgCode));
989  }
990  }
991  }
992 
993  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
994  if (fMCMode == kNoMC ||
995  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
996  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
997  // 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)
998  invMassD = DstarCand->InvMassDstarKpipi();
999  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1000  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1001  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1002  return kTRUE;
1003  }
1004  else { // conditions above not passed, so return FALSE
1005  return kFALSE;
1006  }
1007 }
1008 
1016 {
1017  if (!part) return kDecayOther;
1018  if (!mcArray) return kDecayOther;
1019 
1021 
1022  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1023 
1024  if (part->GetNDaughters() == 2) {
1025 
1026  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1027  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1028 
1029  if (!d1 || !d2) {
1030  return decay;
1031  }
1032 
1033  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1034  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1035 
1036  if (absPdgPart == 421) { // D0 -> K pi
1037 
1038  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1039  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1040  decay = kDecayD0toKpi;
1041  }
1042  }
1043 
1044  if (absPdgPart == 413) { // D* -> D0 pi
1045 
1046  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1047  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1048  if (D0decay == kDecayD0toKpi) {
1049  decay = kDecayDStartoKpipi;
1050  }
1051  }
1052 
1053  if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1054  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1055  if (D0decay == kDecayD0toKpi) {
1056  decay = kDecayDStartoKpipi;
1057  }
1058  }
1059  }
1060  }
1061 
1062  return decay;
1063 }
1064 
1072 {
1073  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1074 
1075  if (!part) return kUnknownQuark;
1076  if (!mcArray) return kUnknownQuark;
1077 
1078  Int_t pdgGranma = 0;
1079  Int_t mother = part->GetMother();
1080  Int_t istep = 0;
1081  Int_t abspdgGranma = 0;
1082  Bool_t isFromB = kFALSE;
1083  Bool_t isQuarkFound = kFALSE;
1084 
1085  while (mother >= 0) {
1086  istep++;
1087  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1088  if (mcGranma >= 0) {
1089  pdgGranma = mcGranma->GetPdgCode();
1090  abspdgGranma = TMath::Abs(pdgGranma);
1091  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1092  isFromB = kTRUE;
1093  }
1094 
1095  if (abspdgGranma == 4 || abspdgGranma == 5) isQuarkFound = kTRUE;
1096  mother = mcGranma->GetMother();
1097  }
1098  else {
1099  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1100  break;
1101  }
1102  }
1103 
1104  if (isQuarkFound) {
1105  if (isFromB) {
1106  return kFromBottom;
1107  }
1108  else {
1109  return kFromCharm;
1110  }
1111  }
1112  else {
1113  return kUnknownQuark;
1114  }
1115 }
1116 
1119 {
1120  fDmesonJets.clear();
1121 
1122  if (fMCMode == kMCTruth) {
1123  RunParticleLevelAnalysis();
1124  }
1125  else {
1126  RunDetectorLevelAnalysis();
1127  }
1128 }
1129 
1132 {
1133  const Int_t nD = fCandidateArray->GetEntriesFast();
1134 
1135  AliDmesonJetInfo DmesonJet;
1136 
1137  Int_t nAccCharm = 0;
1138  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1139  Int_t isSelected = 0;
1140 
1141  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1142  if (!charmCand) continue;
1143 
1144  Int_t nprongs = charmCand->GetNProngs();
1145 
1146  if (fCandidateType == kDstartoKpipi) {
1147  if (!charmCand->InheritsFrom("AliAODRecoCascadeHF")) {
1148  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::RunDetectorLevelAnalysis","Candidate type is D* but object type is wrong (should be AliAODRecoCascadeHF)");
1149  continue;
1150  }
1151  }
1152 
1153  // region of interest + cuts
1154  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1155 
1156  //candidate selected by cuts and PID
1157  isSelected = fRDHFCuts->IsSelected(charmCand, AliRDHFCuts::kAll, fAodEvent); //selected
1158 
1159  if (!isSelected) continue;
1160 
1161  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1162  DmesonJet.Reset();
1163  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1164  for (std::vector<AliHFJetDefinition>::iterator itdef = fJetDefinitions.begin(); itdef != fJetDefinitions.end(); itdef++) {
1165  if (!FindJet(charmCand, DmesonJet, *itdef)) {
1166  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1167  (*itdef).GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1168  }
1169  }
1170  fDmesonJets.push_back(DmesonJet);
1171  }
1172  }
1173  nAccCharm++;
1174  } // end of D cand loop
1175 
1176  TString hname;
1177 
1178  hname = TString::Format("%s/fHistNAcceptedDmesons", GetName());
1179  fHistManager->FillTH1(hname, nAccCharm);
1180 
1181  hname = TString::Format("%s/fHistNDmesons", GetName());
1182  fHistManager->FillTH1(hname, nD);
1183 }
1184 
1195 {
1196  TString hname;
1197 
1199  fFastJetWrapper->SetR(jetDef.fRadius);
1202 
1203  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1204 
1205  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1206  fTrackContainer->SetDMesonCandidate(Dcand);
1207  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1208  AddInputVectors(fTrackContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1209 
1210  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1211  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1212  const TObjArray daughterNotInJet = fTrackContainer->GetDaughterList();
1213  for (Int_t i = 0; i < daughterNotInJet.GetEntriesFast(); i++) {
1214  AliVParticle* daughter = static_cast<AliVParticle*>(daughterNotInJet.At(i));
1215  if (!daughter) continue;
1216  histDaughterNotInJet->Fill(daughter->Pt());
1217  }
1218  }
1219 
1220  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1221  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1222  AddInputVectors(fClusterContainer, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1223  }
1224 
1225  // run jet finder
1226  fFastJetWrapper->Run();
1227 
1228  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1229 
1230  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1231  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1232 
1233  Bool_t isDmesonJet = kFALSE;
1234 
1235  Double_t maxChPt = 0;
1236  Double_t maxNePt = 0;
1237  Double_t totalNeutralPt = 0;
1238 
1239  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1240  if (constituents[ic].user_index() == 0) {
1241  isDmesonJet = kTRUE;
1242  }
1243  else if (constituents[ic].user_index() >= 100) {
1244  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1245  }
1246  else if (constituents[ic].user_index() <= -100) {
1247  totalNeutralPt += constituents[ic].pt();
1248  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1249  }
1250  }
1251 
1252  if (isDmesonJet) {
1253  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1254  DmesonJet.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1255  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1256  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1257  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1258 
1259  return kTRUE;
1260  }
1261  }
1262 
1263  return kFALSE;
1264 }
1265 
1270 {
1271  AliTLorentzVector part;
1272  for (Int_t i = 0; i < cont->GetNEntries(); i++) {
1273  cont->GetMomentum(part, i);
1274  UInt_t rejectionReason = 0;
1275  if (!cont->AcceptObject(i, rejectionReason)) {
1276  rejectHist->Fill(cont->GetRejectionReasonBitPosition(rejectionReason), part.Pt());
1277  continue;
1278  }
1279  Int_t uid = offset >= 0 ? i : -i;
1280  uid += offset;
1281  fFastJetWrapper->AddInputVector(part.Px(), part.Py(), part.Pz(), part.E(), uid);
1282  }
1283 }
1284 
1287 {
1288  TString hname;
1289 
1290  fMCContainer->SetSpecialPDG(fCandidatePDG);
1291  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1292  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1293 
1294  std::map<int, AliDmesonJetInfo> dMesonJets;
1295 
1296  for (std::vector<AliHFJetDefinition>::iterator itdef = fJetDefinitions.begin(); itdef != fJetDefinitions.end(); itdef++) {
1297  AliHFJetDefinition* jetDef = &(*itdef);
1298 
1300  fFastJetWrapper->SetR(jetDef->fRadius);
1303 
1304  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef->GetName());
1305  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1306 
1307  fFastJetWrapper->Run();
1308 
1309  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1310 
1311  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1312  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1313 
1314  Bool_t isDmesonJet = kFALSE;
1315 
1316  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1317  Int_t iPart = constituents[ic].user_index() - 100;
1318  AliVParticle* part = fMCContainer->GetParticle(iPart);
1319  if (!part) {
1320  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1321  continue;
1322  }
1323  if (part->PdgCode() == fCandidatePDG) {
1324  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = dMesonJets.find(iPart);
1325  if (dMesonJetIt == dMesonJets.end()) {
1326  std::pair<int, AliDmesonJetInfo> element;
1327  element.first = iPart;
1328  element.second = AliDmesonJetInfo();
1329  dMesonJetIt = dMesonJets.insert(element).first;
1330  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1331  }
1332 
1333  (*dMesonJetIt).second.fJets[jetDef->GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1334  break;
1335  }
1336  }
1337  }
1338  }
1339 
1340  for (std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = dMesonJets.begin();
1341  dMesonJetIt != dMesonJets.end();
1342  dMesonJetIt++) {
1343  fDmesonJets.push_back((*dMesonJetIt).second);
1344  }
1345 }
1346 
1351 {
1352  TString classname;
1353  switch (fCandidateType) {
1354  case kD0toKpi:
1355  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1356  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1357  break;
1358  case kDstartoKpipi:
1359  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1360  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1361  break;
1362  }
1363  fTree = new TTree(GetName(), GetName());
1364  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo, 32000, 0);
1365  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1366  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1367  fCurrentJetInfo[i] = new AliJetInfoSummary();
1368  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1369  }
1370 
1371  return fTree;
1372 }
1373 
1377 void AliAnalysisTaskDmesonJets::AnalysisEngine::BuildHnSparse(UInt_t enabledAxis, Int_t nBins, Double_t minBinPt, Double_t maxBinPt)
1378 {
1379  TString hname;
1380 
1381  for (std::vector<AliHFJetDefinition>::const_iterator it = fJetDefinitions.begin(); it != fJetDefinitions.end(); it++) {
1382  const AliHFJetDefinition* jetDef = &(*it);
1383 
1384  AliDebug(2,Form("Now working on '%s'", jetDef->GetName()));
1385 
1386  Double_t radius = jetDef->fRadius;
1387 
1388  TString title[30] = {""};
1389  Int_t nbins[30] = {0 };
1390  Double_t min [30] = {0.};
1391  Double_t max [30] = {0.};
1392  Int_t dim = 0 ;
1393 
1394  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1395  nbins[dim] = nBins;
1396  min[dim] = 0;
1397  max[dim] = 100;
1398  dim++;
1399 
1400  if ((enabledAxis & kPositionD) != 0) {
1401  title[dim] = "#eta_{D}";
1402  nbins[dim] = 50;
1403  min[dim] = -1;
1404  max[dim] = 1;
1405  dim++;
1406 
1407  title[dim] = "#phi_{D} (rad)";
1408  nbins[dim] = 150;
1409  min[dim] = 0;
1410  max[dim] = TMath::TwoPi();
1411  dim++;
1412  }
1413 
1414  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1415  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1416  nbins[dim] = fNMassBins;
1417  min[dim] = fMinMass;
1418  max[dim] = fMaxMass;
1419  dim++;
1420  }
1421 
1422  if (fCandidateType == kD0toKpi) {
1423  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1424  nbins[dim] = fNMassBins;
1425  min[dim] = fMinMass;
1426  max[dim] = fMaxMass;
1427  dim++;
1428  }
1429 
1430  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1431  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1432  nbins[dim] = fNMassBins;
1433  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1434  dim++;
1435  }
1436 
1437  if (fCandidateType == kDstartoKpipi) {
1438  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1439  nbins[dim] = fNMassBins*6;
1440  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1441 
1442  // subtract mass of D0
1443  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1444  min[dim] -= D0mass;
1445  max[dim] -= D0mass;
1446 
1447  dim++;
1448  }
1449 
1450  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1451  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1452  nbins[dim] = 100;
1453  min[dim] = 0;
1454  max[dim] = 25;
1455  dim++;
1456  }
1457 
1458  title[dim] = "#it{z}_{D}";
1459  nbins[dim] = 110;
1460  min[dim] = 0;
1461  max[dim] = 1.10;
1462  dim++;
1463 
1464  if ((enabledAxis & kDeltaR) != 0) {
1465  title[dim] = "#Delta R_{D-jet}";
1466  nbins[dim] = 100;
1467  min[dim] = 0;
1468  max[dim] = radius * 1.5;
1469  dim++;
1470  }
1471 
1472  if ((enabledAxis & kDeltaEta) != 0) {
1473  title[dim] = "#eta_{D} - #eta_{jet}";
1474  nbins[dim] = 100;
1475  min[dim] = -radius * 1.2;
1476  max[dim] = radius * 1.2;
1477  dim++;
1478  }
1479 
1480  if ((enabledAxis & kDeltaPhi) != 0) {
1481  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1482  nbins[dim] = 100;
1483  min[dim] = -radius * 1.2;
1484  max[dim] = radius * 1.2;
1485  dim++;
1486  }
1487 
1488  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1489  nbins[dim] = nBins;
1490  min[dim] = minBinPt;
1491  max[dim] = maxBinPt;
1492  dim++;
1493 
1494  if ((enabledAxis & kPositionJet) != 0) {
1495  title[dim] = "#eta_{jet}";
1496  nbins[dim] = 50;
1497  min[dim] = -1;
1498  max[dim] = 1;
1499  dim++;
1500 
1501  title[dim] = "#phi_{jet} (rad)";
1502  nbins[dim] = 150;
1503  min[dim] = 0;
1504  max[dim] = TMath::TwoPi();
1505  dim++;
1506  }
1507 
1508  if ((enabledAxis & kJetConstituents) != 0) {
1509  title[dim] = "No. of constituents";
1510  nbins[dim] = 50;
1511  min[dim] = -0.5;
1512  max[dim] = 49.5;
1513  dim++;
1514  }
1515 
1516  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef->GetName());
1517  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1518  for (Int_t j = 0; j < dim; j++) {
1519  h->GetAxis(j)->SetTitle(title[j]);
1520  }
1521  }
1522 }
1523 
1528 {
1529  TString hname;
1530 
1531  for (Int_t id = 0; id < fDmesonJets.size(); id++) {
1532  fCurrentDmesonJetInfo->Set(fDmesonJets[id]);
1533  Int_t accJets = 0;
1534  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1535  fCurrentJetInfo[ij]->Reset();
1536  AliJetInfo* jet = fDmesonJets[id].GetJet(fJetDefinitions[ij].GetName());
1537  if (!jet) continue;
1538  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1539  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1540  fHistManager->FillTH1(hname, jet->Pt());
1541  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1542  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1543  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1544  fHistManager->FillTH1(hname, jet->Eta());
1545  continue;
1546  }
1547  fCurrentJetInfo[ij]->Set(fDmesonJets[id], fJetDefinitions[ij].GetName());
1548  accJets++;
1549  }
1550  if (accJets > 0) {
1551  fTree->Fill();
1552  }
1553  else {
1554  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1555  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Pt());
1556  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1557  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Phi_0_2pi());
1558  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1559  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Eta());
1560  if (fCandidateType == kD0toKpi) {
1561  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1562  fHistManager->FillTH1(hname, fDmesonJets[id].fD.M());
1563  }
1564  else if (fCandidateType == kDstartoKpipi) {
1565  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1566  fHistManager->FillTH1(hname, fDmesonJets[id].fInvMass2Prong);
1567 
1568  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1569  fHistManager->FillTH1(hname, fDmesonJets[id].fD.M() - fDmesonJets[id].fInvMass2Prong);
1570  }
1571  }
1572  }
1573  return kTRUE;
1574 }
1575 
1580 {
1581  TString hname;
1582 
1583  for (Int_t id = 0; id < fDmesonJets.size(); id++) {
1584  if (!IsAnyJetInAcceptance(fDmesonJets[id])) {
1585  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1586  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Pt());
1587  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1588  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Phi_0_2pi());
1589  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1590  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Eta());
1591  }
1592  }
1593 
1594  for (std::vector<AliHFJetDefinition>::iterator itdef = fJetDefinitions.begin(); itdef != fJetDefinitions.end(); itdef++) {
1595 
1596  AliHFJetDefinition* jetDef = &(*itdef);
1597 
1598  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef->GetName());
1599  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
1600 
1601  for (Int_t id = 0; id < fDmesonJets.size(); id++) {
1602  const AliJetInfo* jet = fDmesonJets[id].GetJet(jetDef->GetName());
1603  if (!jet) continue;
1604  if (!jetDef->IsJetInAcceptance(*jet)) {
1605  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef->GetName());
1606  fHistManager->FillTH1(hname, jet->Pt());
1607  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef->GetName());
1608  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1609  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef->GetName());
1610  fHistManager->FillTH1(hname, jet->Eta());
1611  continue;
1612  }
1613  FillHnSparse(h, fDmesonJets[id], (*itdef).GetName());
1614  }
1615  }
1616 
1617  return kTRUE;
1618 }
1619 
1625 Bool_t AliAnalysisTaskDmesonJets::AnalysisEngine::FillHnSparse(THnSparse* h, const AliDmesonJetInfo& DmesonJet, std::string n)
1626 {
1627  // Fill the THnSparse histogram.
1628 
1629  Double_t contents[30] = {0.};
1630 
1631  Double_t z = DmesonJet.GetZ(n);
1632  Double_t deltaPhi = 0;
1633  Double_t deltaEta = 0;
1634  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
1635 
1636  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
1637  if (it == DmesonJet.fJets.end()) return kFALSE;
1638 
1639  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
1640  TString title(h->GetAxis(i)->GetTitle());
1641  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
1642  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
1643  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
1644  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
1645  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
1646  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
1647  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
1648  else if (title=="#it{z}_{D}") contents[i] = z;
1649  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
1650  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
1651  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
1652  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
1653  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
1654  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
1655  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
1656  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
1657  }
1658 
1659  h->Fill(contents);
1660 
1661  return kTRUE;
1662 }
1663 
1664 // Definitions of class AliAnalysisTaskDmesonJets
1665 
1669 
1673  fAnalysisEngines(),
1674  fEnabledAxis(0),
1675  fTreeOutput(kFALSE),
1676  fHistManager(),
1677  fApplyKinematicCuts(kTRUE),
1678  fAodEvent(0),
1679  fFastJetWrapper(0)
1680 {
1681  SetMakeGeneralHistograms(kTRUE);
1682 }
1683 
1688  AliAnalysisTaskEmcal(name, kTRUE),
1689  fAnalysisEngines(),
1690  fEnabledAxis(k2ProngInvMass),
1691  fTreeOutput(kFALSE),
1692  fHistManager(name),
1693  fApplyKinematicCuts(kTRUE),
1694  fAodEvent(0),
1695  fFastJetWrapper(0)
1696 {
1697  SetMakeGeneralHistograms(kTRUE);
1698 }
1699 
1702 {
1703  if (fFastJetWrapper) delete fFastJetWrapper;
1704 }
1705 
1713 {
1714  AliRDHFCuts* analysiscuts = 0;
1715  TFile* filecuts = TFile::Open(cutfname);
1716  if (!filecuts || filecuts->IsZombie()) {
1717  ::Warning("AddTaskDmesonJets", "Input file not found: will use std cuts.");
1718  filecuts = 0;
1719  }
1720 
1721  if (filecuts) {
1722  analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
1723  if (!analysiscuts) {
1724  ::Warning("AddTaskDmesonJetCorr", "Could not find analysis cuts '%s' in '%s'. Using std cuts.", cutsname.Data(), cutfname.Data());
1725  }
1726  }
1727 
1728  return analysiscuts;
1729 }
1730 
1740 {
1742  return AddAnalysisEngine(type, MCmode, jetDef, cutfname);
1743 }
1744 
1754 {
1755  AliRDHFCuts* cuts = 0;
1756 
1757  if (!cutfname.IsNull()) {
1758  TString cutsname;
1759 
1760  switch (type) {
1761  case kD0toKpi :
1762  cutsname = "D0toKpiCuts";
1763  break;
1764  case kDstartoKpipi :
1765  cutsname = "DStartoKpipiCuts";
1766  break;
1767  default:
1768  return 0;
1769  }
1770 
1771  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
1772  }
1773 
1774  AnalysisEngine eng(type, MCmode, cuts);
1775 
1776  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
1777 
1778  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
1779  eng.AddJetDefinition(jetDef);
1780  it = fAnalysisEngines.insert(it, eng);
1781  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(jetDef), fAnalysisEngines.size());
1782  }
1783  else {
1784  AnalysisEngine* found_eng = &(*it);
1785  ::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());
1786  found_eng->AddJetDefinition(jetDef);
1787 
1788  if (cuts && found_eng->fRDHFCuts != 0) {
1789  ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
1790  found_eng->SetRDHFCuts(cuts);
1791  }
1792  }
1793 
1794  return &(*it);
1795 }
1796 
1797 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
1798 {
1799  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
1800  while (it != fAnalysisEngines.end() && (*it) < eng) it++;
1801  return it;
1802 }
1803 
1806 {
1807  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
1808 
1810 
1811  // Define histograms
1812  // the TList fOutput is already defined in AliAnalysisTaskEmcal::UserCreateOutputObjects()
1813 
1814  TString hname;
1815  TString htitle;
1816  TH1* h = 0;
1817 
1818  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
1819  for (std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin(); it != fAnalysisEngines.end(); it++) {
1820  AnalysisEngine* param = &(*it);
1821  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param->GetName(), param->fJetDefinitions.size());
1822 
1823  fHistManager.CreateHistoGroup(param->GetName());
1824 
1825  param->fHistManager = &fHistManager;
1826 
1827  hname = TString::Format("%s/fHistNAcceptedDmesons", param->GetName());
1828  htitle = hname + ";Number of D accepted meson candidates;counts";
1829  h = fHistManager.CreateTH1(hname, htitle, 51, -0.5, 50.5);
1830 
1831  hname = TString::Format("%s/fHistNDmesons", param->GetName());
1832  htitle = hname + ";Number of D meson candidates;counts";
1833  h = fHistManager.CreateTH1(hname, htitle, 101, -0.5, 100.5);
1834 
1835  hname = TString::Format("%s/fHistNEvents", param->GetName());
1836  htitle = hname + ";Event status;counts";
1837  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
1838  h->GetXaxis()->SetBinLabel(1, "Accepted");
1839  h->GetXaxis()->SetBinLabel(2, "Rejected");
1840 
1841  hname = TString::Format("%s/fHistEventRejectionReasons", param->GetName());
1842  htitle = hname + ";Rejection reason;counts";
1843  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
1844 
1845  hname = TString::Format("%s/fHistRejectedDMesonPt", param->GetName());
1846  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
1847  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
1848 
1849  hname = TString::Format("%s/fHistRejectedDMesonEta", param->GetName());
1850  htitle = hname + ";#it{#eta}_{D};counts";
1851  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
1852 
1853  hname = TString::Format("%s/fHistRejectedDMesonPhi", param->GetName());
1854  htitle = hname + ";#it{#phi}_{D};counts";
1855  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
1856 
1857  if (param->fCandidateType == kD0toKpi) {
1858  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param->GetName());
1859  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
1860  fHistManager.CreateTH1(hname, htitle, param->fNMassBins, param->fMinMass, param->fMaxMass);
1861  }
1862  else if (param->fCandidateType == kDstartoKpipi) {
1863  Double_t min = 0;
1864  Double_t max = 0;
1865 
1866  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param->GetName());
1867  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
1868  CalculateMassLimits(param->fMaxMass - param->fMinMass, 421, param->fNMassBins, min, max);
1869  fHistManager.CreateTH1(hname, htitle, param->fNMassBins, min, max);
1870 
1871  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1872  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param->GetName());
1873  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
1874  CalculateMassLimits(0.20, 413, param->fNMassBins*6, min, max);
1875  fHistManager.CreateTH1(hname, htitle, param->fNMassBins*6, min-D0mass, max-D0mass);
1876  }
1877 
1878  for (std::vector<AliHFJetDefinition>::iterator itdef = param->fJetDefinitions.begin(); itdef != param->fJetDefinitions.end(); itdef++) {
1879  AliHFJetDefinition* jetDef = &(*itdef);
1880  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef->GetName());
1881 
1882  fHistManager.CreateHistoGroup(jetDef->GetName(), param->GetName());
1883 
1884  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param->GetName(), jetDef->GetName());
1885  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
1886  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
1887  SetRejectionReasonLabels(h->GetXaxis());
1888 
1889  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param->GetName(), jetDef->GetName());
1890  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
1891  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
1892  SetRejectionReasonLabels(h->GetXaxis());
1893 
1894  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param->GetName(), jetDef->GetName());
1895  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
1896  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
1897  SetRejectionReasonLabels(h->GetXaxis());
1898 
1899  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param->GetName(), jetDef->GetName());
1900  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
1901  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
1902 
1903  hname = TString::Format("%s/%s/fHistRejectedJetPt", param->GetName(), jetDef->GetName());
1904  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
1905  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
1906 
1907  hname = TString::Format("%s/%s/fHistRejectedJetEta", param->GetName(), jetDef->GetName());
1908  htitle = hname + ";#it{#eta}_{jet};counts";
1909  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
1910 
1911  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param->GetName(), jetDef->GetName());
1912  htitle = hname + ";#it{#phi}_{jet};counts";
1913  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
1914  }
1915  if (fTreeOutput) {
1916  TTree* tree = param->BuildTree();
1917  fOutput->Add(tree);
1918  }
1919  else {
1921  }
1922  }
1923 
1924  fOutput->Add(fHistManager.GetListOfHistograms());
1925 
1926  PostData(1, fOutput);
1927 }
1928 
1932 {
1934 
1935  // Load the event
1936  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
1937 
1938  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
1939 
1940  fFastJetWrapper->SetAreaType(fastjet::active_area);
1942 
1943  if (!fAodEvent) {
1944  AliError(Form("This task need an AOD event! Task '%s' will be disabled!", GetName()));
1945  return;
1946  }
1947 
1948  for (std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin(); it != fAnalysisEngines.end(); it++) {
1949  AnalysisEngine* params = &(*it);
1950 
1951  params->fAodEvent = fAodEvent;
1952  params->fFastJetWrapper = fFastJetWrapper;
1953  params->Init(fGeom, fAodEvent->GetRunNumber());
1954 
1955  if (params->fMCMode != kMCTruth) {
1956  params->fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params->fBranchName.Data()));
1957 
1958  if (params->fCandidateArray) {
1959  if (!params->fCandidateArray->GetClass()->InheritsFrom("AliAODRecoDecayHF2Prong")) {
1960  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
1961  "%s: Objects of type %s in %s are not inherited from AliAODRecoDecayHF2Prong! Task will be disabled!",
1962  GetName(), params->fCandidateArray->GetClass()->GetName(), params->fCandidateArray->GetName());
1963  params->fCandidateArray = 0;
1964  params->fInhibit = kTRUE;
1965  }
1966  }
1967  else {
1968  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
1969  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
1970  params->fBranchName.Data(), params->GetName());
1971  params->fInhibit = kTRUE;
1972  }
1973  }
1974 
1975  if (params->fMCMode != kNoMC) {
1976  params->fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(0));
1977 
1978  if (!params->fMCContainer) params->fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(1));
1979 
1980  if (!params->fMCContainer) {
1981  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
1982  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
1983  params->GetName());
1984  params->fInhibit = kTRUE;
1985  }
1986  }
1987 
1988  if (params->fMCMode != kMCTruth) {
1989  params->fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(0));
1990  if (!params->fTrackContainer) params->fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(1));
1991 
1993 
1994  if (!params->fTrackContainer && !params->fClusterContainer) {
1995  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
1996  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
1997  params->GetName());
1998  params->fInhibit = kTRUE;
1999  }
2000  }
2001  }
2002 }
2003 
2008 {
2009  if (!fAodEvent) return kFALSE;
2010 
2011  TString hname;
2012 
2013  // fix for temporary bug in ESDfilter
2014  // the AODs with null vertex pointer didn't pass the PhysSel
2015  if (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001) return kFALSE;
2016 
2017  for (std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin(); it != fAnalysisEngines.end(); it++) {
2018  AnalysisEngine* eng = &(*it);
2019 
2020  if (eng->fInhibit) continue;
2021 
2022  //Event selection
2023  hname = TString::Format("%s/fHistNEvents", eng->GetName());
2024  Bool_t iseventselected = eng->fRDHFCuts->IsEventSelected(fAodEvent);
2025  if (!iseventselected) {
2026  fHistManager.FillTH1(hname, "Rejected");
2027  hname = TString::Format("%s/fHistEventRejectionReasons", eng->GetName());
2028  UInt_t bitmap = eng->fRDHFCuts->GetEventRejectionBitMap();
2029  TString label;
2030  do {
2031  label = GetHFEventRejectionReasonLabel(bitmap);
2032  if (label.IsNull()) break;
2033  fHistManager.FillTH1(hname, label);
2034  } while (true);
2035 
2036  continue;
2037  }
2038 
2039  fHistManager.FillTH1(hname, "Accepted");
2040 
2041  AliDebug(2, "Event selected");
2042 
2043  eng->RunAnalysis();
2044  }
2045  return kTRUE;
2046 }
2047 
2052 {
2053  TString hname;
2054  for (std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin(); it != fAnalysisEngines.end(); it++) {
2055  AnalysisEngine* param = &(*it);
2056 
2057  if (param->fInhibit) continue;
2058 
2059  if (fTreeOutput) {
2060  param->FillTree(fApplyKinematicCuts);
2061  }
2062  else {
2064  }
2065  }
2066  return kTRUE;
2067 }
2068 
2076 void AliAnalysisTaskDmesonJets::CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t& minMass, Double_t& maxMass)
2077 {
2078  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
2079 
2080  Double_t mass = part->Mass();
2081 
2082  // To make sure that the PDG mass value is not at the edge of a bin
2083  if (nbins % 2 == 0) {
2084  minMass = mass - range / 2 - range / nbins / 2;
2085  maxMass = mass + range / 2 - range / nbins / 2;
2086  }
2087  else {
2088  minMass = mass - range / 2;
2089  maxMass = mass + range / 2;
2090  }
2091 }
2092 
2099 {
2100  static TString label;
2101  label = "";
2102 
2103  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
2104  label = "NotSelTrigger";
2105  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
2106  return label.Data();
2107  }
2108  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
2109  label = "NoVertex";
2110  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
2111  return label.Data();
2112  }
2113  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
2114  label = "TooFewVtxContrib";
2115  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
2116  return label.Data();
2117  }
2118  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
2119  label = "ZVtxOutFid";
2120  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
2121  return label.Data();
2122  }
2123  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
2124  label = "Pileup";
2125  bitmap &= ~BIT(AliRDHFCuts::kPileup);
2126  return label.Data();
2127  }
2128  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
2129  label = "OutsideCentrality";
2130  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
2131  return label.Data();
2132  }
2133  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
2134  label = "PhysicsSelection";
2135  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
2136  return label.Data();
2137  }
2138  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
2139  label = "BadSPDVertex";
2140  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
2141  return label.Data();
2142  }
2143  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
2144  label = "ZVtxSPDOutFid";
2145  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
2146  return label.Data();
2147  }
2148  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
2149  label = "CentralityFlattening";
2150  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
2151  return label.Data();
2152  }
2153  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
2154  label = "BadTrackV0Correl";
2155  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
2156  return label.Data();
2157  }
2158 
2159  return label.Data();
2160 }
Int_t pdg
void Print() const
Prints the content of this object in the standard output.
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
Analysis task for D meson jets.
virtual void Set(const AliDmesonJetInfo &source)
Invariant mass of the D0 meson candidate in GeV/c2.
virtual Int_t Run()
Bool_t IsJetInAcceptance(const AliJetInfo &jet) const
Double_t fSoftPionPt
! Transverse momentum of the soft pion of the D* candidate
const char * title
Definition: MakeQAPdf.C:26
AliRDHFCuts * fRDHFCuts
D meson candidates cuts.
Lightweight class that encapsulates D meson jets.
Bool_t fTreeOutput
If true, output will be posted in a TTree rather than a THnSparse.
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
ERecoScheme_t fRecoScheme
Jet recombination scheme (pt scheme, E scheme, ...)
Base task in the EMCAL framework.
Double_t fMaxNeutralPt
Transverse momentum of the leading neutral particle (or cluster)
void RunDetectorLevelAnalysis()
Run a detector level analysis.
Double_t mass
virtual void UserCreateOutputObjects()
Creates the output containers.
Lightweight class that encapsulates D meson jets.
Double_t fMinBinPt
min pt in histograms
Double_t fInvMass2Prong
! 2-prong mass of the D* candidate (w/o the soft pion)
static void CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t &minMass, Double_t &maxMass)
AliAODEvent * fAodEvent
! AOD event
Bool_t FindJet(AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, AliHFJetDefinition &jetDef)
UInt_t fEnabledAxis
Use bit defined in EAxis_t to enable axis in the THnSparse.
EMCMode_t fMCMode
MC mode: No MC (data and MC detector level), background-only (MC), signal-only (MC), MC truth (particle level)
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:82
AnalysisEngine()
This is the default constructor, used for ROOT I/O purposes.
TList * fOutput
!output list
void Init(const AliEMCALGeometry *const geom, Int_t runNumber)
Initialize the analysis engine.
static FJRecoScheme ConvertToFJRecoScheme(ERecoScheme_t reco)
static TString GenerateJetName(EJetType_t jetType, EJetAlgo_t jetAlgo, ERecoScheme_t recoScheme, Double_t radius, AliParticleContainer *partCont, AliClusterContainer *clusCont, TString tag)
void SetR(Double_t r)
Definition: AliFJWrapper.h:87
Select tracks based on specific prescriptions of HF analysis.
ECandidateType_t fCandidateType
Candidate type.
static EMesonOrigin_t CheckOrigin(AliAODMCParticle *part, TClonesArray *mcArray)
Double_t GetDistance(std::string n, Double_t &deta, Double_t &dphi) const
AliParticleContainer * GetParticleContainer(Int_t i=0) const
virtual Bool_t AcceptObject(Int_t i, UInt_t &rejectionReason) const =0
virtual void Set(const AliDmesonJetInfo &source)
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:32
EJetType_t fJetType
Jet type (charged, full, neutral)
virtual void Set(const AliDmesonJetInfo &source)
AliTLorentzVector fD
! 4-momentum of the D meson candidate
UShort_t GetRejectionReasonBitPosition(UInt_t rejectionReason) const
AliAnalysisTaskDmesonJets()
This is the default constructor, used for ROOT I/O purposes.
AliClusterContainer * fClusterContainer
! Cluster container
virtual void Reset()
Reset the current object.
AliEMCALGeometry * fGeom
!emcal geometry
AnalysisEngine & operator=(const AnalysisEngine &source)
Struct that encapsulates analysis parameters.
std::vector< AliHFJetDefinition >::iterator FindJetDefinition(const AliHFJetDefinition &eng)
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:81
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Base class for container structures within the EMCAL framework.
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)
virtual Bool_t GetMomentum(TLorentzVector &mom, Int_t i) const =0
static EMesonDecayChannel_t CheckDecayChannel(AliAODMCParticle *part, TClonesArray *mcArray)
EJetAcceptanceType_t fAcceptance
Jet acceptance.
Double_t fMaxMass
Max mass in histogram axis.
EJetAlgo_t fJetAlgo
Jet algorithm (kt, anti-kt,...)
static const char * GetHFEventRejectionReasonLabel(UInt_t &bitmap)
virtual void Clear(const Option_t *="")
UInt_t GetEventRejectionBitMap() const
Definition: AliRDHFCuts.h:295
Bool_t ExtractDstarAttributes(const AliAODRecoCascadeHF *DstarCand, AliDmesonJetInfo &DmesonJet, UInt_t i)
const AliJetInfo * GetJet(std::string n) const
void BuildHnSparse(UInt_t enabledAxis, Int_t nBins, Double_t minBinPt, Double_t maxBinPt)
TClonesArray * fCandidateArray
! D meson candidate array
vector< AliHFJetDefinition > fJetDefinitions
Jet definitions.
virtual ~AliAnalysisTaskDmesonJets()
This is the standard destructor.
Bool_t ExtractD0Attributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
Bool_t ExtractParticleLevelHFAttributes(const AliAODMCParticle *part, AliDmesonJetInfo &DmesonJet)
Select MC particles based on specific prescriptions of HF analysis.
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:85
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)
Bool_t IsEventSelected(AliVEvent *event)
std::map< std::string, AliJetInfo > fJets
! list of jets
Double_t fMaxBinPt
max pt in histograms
Double_t minMass
void RunParticleLevelAnalysis()
Run a particle level analysis.
void SetMakeGeneralHistograms(Bool_t g)
Bool_t fApplyKinematicCuts
Apply jet kinematic cuts.
const char * GetName() const
Generate a name for this jet definition.
void SetDetectorJetEtaPhiRange(const AliEMCALGeometry *const geom, Int_t run)
Bool_t ExtractRecoDecayAttributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
Double_t InvMassD0() const
list< AnalysisEngine > fAnalysisEngines
Array of analysis parameters.
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="")
void SetRejectionReasonLabels(TAxis *axis)
Int_t GetNEntries() const
const Int_t nbins
Double_t maxMass
TString fBranchName
AOD branch where the D meson candidate are found.
AliHFJetDefinition()
This is the default constructor, used for ROOT I/O purposes.
AliHFJetDefinition & operator=(const AliHFJetDefinition &source)
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:83
AliHFAODMCParticleContainer * fMCContainer
! MC particle container
AliFJWrapper * fFastJetWrapper
! Fastjet wrapper
AliHFJetDefinition * AddJetDefinition(EJetType_t type, Double_t r, EJetAlgo_t algo, ERecoScheme_t reco)
void AddInputVectors(AliEmcalContainer *cont, Int_t offset, TH2 *rejectHist)
Lightweight class that encapsulates D0.
AliHFTrackContainer * fTrackContainer
! Track container
Bool_t IsAnyJetInAcceptance(const AliDmesonJetInfo &dMesonJet) const
Int_t fNbins
no. of pt bins
AliFJWrapper * fFastJetWrapper
! Fastjet wrapper
THistManager fHistManager
Histogram manager.
void Reset()
Reset all fields to their default values.
virtual void Set(const AliDmesonJetInfo &source, std::string n)
Double_t fMinMass
Min mass in histogram axis.