AliPhysics  a9863a5 (a9863a5)
 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  fPtBinWidth(0.5),
520  fMaxPt(100),
521  fTree(0),
522  fCurrentDmesonJetInfo(0),
523  fCurrentJetInfo(0),
524  fCandidateArray(0),
525  fMCContainer(0),
526  fTrackContainer(0),
527  fClusterContainer(0),
528  fAodEvent(0),
529  fFastJetWrapper(0),
530  fHistManager(0)
531 {
532 }
533 
542  TObject(),
543  fCandidateType(type),
544  fCandidateName(),
545  fCandidatePDG(0),
546  fNDaughters(0),
547  fBranchName(),
548  fMCMode(MCmode),
549  fNMassBins(nMassBins),
550  fMinMass(0),
551  fMaxMass(0),
552  fRDHFCuts(cuts),
553  fRejectedOrigin(kUnknownQuark | kFromBottom),
554  fAcceptedDecay(kAnyDecay),
555  fInhibit(kFALSE),
556  fJetDefinitions(),
557  fPtBinWidth(0.5),
558  fMaxPt(100),
559  fTree(0),
560  fCurrentDmesonJetInfo(0),
561  fCurrentJetInfo(0),
562  fCandidateArray(0),
563  fMCContainer(0),
564  fTrackContainer(0),
565  fClusterContainer(0),
566  fAodEvent(0),
567  fFastJetWrapper(0),
568  fHistManager(0)
569 {
570  SetCandidateProperties(range);
571 }
572 
577  TObject(source),
578  fCandidateType(source.fCandidateType),
579  fCandidateName(source.fCandidateName),
580  fCandidatePDG(source.fCandidatePDG),
581  fNDaughters(source.fNDaughters),
582  fBranchName(source.fBranchName),
583  fMCMode(source.fMCMode),
584  fNMassBins(source.fNMassBins),
585  fMinMass(source.fMinMass),
586  fMaxMass(source.fMaxMass),
587  fRDHFCuts(),
588  fRejectedOrigin(source.fRejectedOrigin),
589  fAcceptedDecay(source.fAcceptedDecay),
590  fInhibit(source.fInhibit),
591  fJetDefinitions(source.fJetDefinitions),
592  fPtBinWidth(source.fPtBinWidth),
593  fMaxPt(source.fMaxPt),
594  fTree(0),
595  fCurrentDmesonJetInfo(0),
596  fCurrentJetInfo(0),
597  fCandidateArray(source.fCandidateArray),
598  fMCContainer(source.fMCContainer),
599  fTrackContainer(source.fTrackContainer),
600  fClusterContainer(source.fClusterContainer),
601  fAodEvent(source.fAodEvent),
603  fHistManager(source.fHistManager)
604 {
605  SetRDHFCuts(source.fRDHFCuts);
606 }
607 
608 // Destructor
610 {
611  if (fRDHFCuts) delete fRDHFCuts;
612 }
613 
618 {
619  new (this) AnalysisEngine(source);
620  return *this;
621 }
622 
627 {
628  for (UInt_t i = 0; i < fJetDefinitions.size(); i++) {
629  if (fJetDefinitions[i].IsJetInAcceptance(dMesonJet, fJetDefinitions[i].GetName())) return kTRUE;
630  }
631 
632  return kFALSE;
633 }
634 
636 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const geom, Int_t runNumber)
637 {
638  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
639  fJetDefinitions[i].SetDetectorJetEtaPhiRange(geom, runNumber);
640  }
641 }
642 
647 {
648  switch (fCandidateType) {
649  case kD0toKpi:
650  fCandidatePDG = 421;
651  fCandidateName = "D0";
652  fNDaughters = 2;
653  fPDGdaughters.Set(fNDaughters);
654  fPDGdaughters.Reset();
655  fPDGdaughters[0] = 211; // pi
656  fPDGdaughters[1] = 321; // K
657  fBranchName = "D0toKpi";
658  fAcceptedDecay = kD0toKpi;
659  if (!fRDHFCuts) {
660  fRDHFCuts = new AliRDHFCutsD0toKpi();
661  fRDHFCuts->SetStandardCutsPP2010();
662  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
663  fRDHFCuts->SetTriggerClass("","");
664  }
665  break;
666  case kDstartoKpipi:
667  fCandidatePDG = 413;
668  fCandidateName = "DStar";
669  fNDaughters = 3;
670  fPDGdaughters.Set(fNDaughters);
671  fPDGdaughters.Reset();
672  fPDGdaughters[0] = 211; // pi soft
673  fPDGdaughters[1] = 211; // pi fromD0
674  fPDGdaughters[2] = 321; // K from D0
675  fBranchName = "Dstar";
676  fAcceptedDecay = kDstartoKpipi;
677  if (!fRDHFCuts) {
678  fRDHFCuts = new AliRDHFCutsDStartoKpipi();
679  fRDHFCuts->SetStandardCutsPP2010();
680  fRDHFCuts->SetUsePhysicsSelection(kFALSE);
681  fRDHFCuts->SetTriggerClass("","");
682  }
683  break;
684  default:
685  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
686  }
687 
688  CalculateMassLimits(range, fCandidatePDG, fNMassBins, fMinMass, fMaxMass);
689 }
690 
695 {
696  if (fRDHFCuts) delete fRDHFCuts;
697  fRDHFCuts = cuts;
698 }
699 
704 {
705  if (!cuts) return;
706  if (fRDHFCuts) delete fRDHFCuts;
707  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
708 }
709 
714 {
715  static TString name;
716 
717  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
718 
719  return name.Data();
720 }
721 
726 {
727  static TString name;
728 
729  name = fCandidateName;
730  switch (fMCMode) {
731  case kBackgroundOnly:
732  name += "_kBackgroundOnly";
733  break;
734  case kSignalOnly:
735  name += "_kSignalOnly";
736  break;
737  case kMCTruth:
738  name += "_MCTruth";
739  break;
740  default:
741  break;
742  }
743 
744  return name.Data();
745 }
746 
754 {
755  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
756 
757  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
758  fJetDefinitions.push_back(def);
759  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
760  "Total number of jet definitions is now %lu.",
761  def.GetName(), GetName(), fJetDefinitions.size());
762  // For detector level set maximum track pt to 100 GeV/c
763  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
764  }
765  else {
766  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
767  }
768 
769  return &(*it);
770 }
771 
783 {
784  AliHFJetDefinition def(type, r, algo, reco);
785 
786  return AddJetDefinition(def);
787 }
788 
794 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
795 {
796  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
797  while (it != fJetDefinitions.end() && (*it) != def) it++;
798  return it;
799 }
800 
807 {
808  if (lhs.fCandidateType > rhs.fCandidateType) return false;
809  else if (lhs.fCandidateType < rhs.fCandidateType) return true;
810  else {
811  if (lhs.fMCMode < rhs.fMCMode) return true;
812  else return false;
813  }
814 }
815 
822 {
823  if (lhs.fCandidateType != rhs.fCandidateType) return false;
824  if (lhs.fMCMode != rhs.fMCMode) return false;
825  return true;
826 }
827 
836 {
837  DmesonJet.fD.SetPtEtaPhiM(part->Pt(), part->Eta(), part->Phi(), part->M());
838  return kTRUE;
839 }
840 
849 {
850  if (fCandidateType == kD0toKpi) { // D0 candidate
851  return ExtractD0Attributes(Dcand, DmesonJet, i);
852  }
853  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
854  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
855  }
856  else {
857  return kFALSE;
858  }
859 }
860 
869 {
870  Int_t MCtruthPdgCode = 0;
871 
872  Double_t invMassD = 0;
873 
874  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
875  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
876  if (mcLab >= 0) {
877  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
878 
879  if (aodMcPart) {
880  if (fRejectedOrigin && fMCMode == kSignalOnly) {
881  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
882 
883  if ((origin & fRejectedOrigin) == origin) return kFALSE;
884  }
885  MCtruthPdgCode = aodMcPart->PdgCode();
886  }
887  }
888  }
889 
890  //AliDebug(2,"Checking if D0 meson is selected");
891  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
892  if (isSelected == 1) { // selected as a D0
893  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
894 
895  if (fMCMode == kNoMC ||
896  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
897  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly)) {
898  // 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)
899  //AliDebug(2,"Selected as D0");
900  invMassD = Dcand->InvMassD0();
901  }
902  else { // conditions above not passed, so return FALSE
903  return kFALSE;
904  }
905  }
906  else if (isSelected == 2) { // selected as a D0bar
907  if (i > 0) return kFALSE; // only one mass hypothesis thanks to PID
908 
909  if (fMCMode == kNoMC ||
910  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
911  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly)) {
912  // 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)
913  //AliDebug(2,"Selected as D0bar");
914  invMassD = Dcand->InvMassD0bar();
915  }
916  else { // conditions above not passed, so return FALSE
917  return kFALSE;
918  }
919  }
920  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
921  //AliDebug(2,"Selected as either D0 or D0bar");
922 
923  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
924  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
925  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kBackgroundOnly)) {
926  if (i > 0) return kFALSE;
927  //AliDebug(2, "MC truth is D0");
928  invMassD = Dcand->InvMassD0();
929  }
930  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
931  (MCtruthPdgCode == fCandidatePDG && fMCMode == kBackgroundOnly)) {
932  if (i > 0) return kFALSE;
933  //AliDebug(2, "MC truth is D0bar");
934  invMassD = Dcand->InvMassD0bar();
935  }
936  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
937 
938  // Only accept it if background-only OR background-and-signal was requested
939  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
940  // Select D0 or D0bar depending on the i-parameter
941  if (i == 0) {
942  //AliDebug(2, "Returning invariant mass with D0 hypothesis");
943  invMassD = Dcand->InvMassD0();
944  }
945  else if (i == 1) {
946  //AliDebug(2, "Returning invariant mass with D0bar hypothesis");
947  invMassD = Dcand->InvMassD0bar();
948  }
949  else { // i > 1
950  return kFALSE;
951  }
952  }
953  else { // signal-only was requested but this is not a true D0
954  return kFALSE;
955  }
956  }
957  }
958  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
959  return kTRUE;
960 }
961 
970 {
971  if (i > 0) return kFALSE; // only one mass hypothesis for the D*
972 
973  Int_t MCtruthPdgCode = 0;
974 
975  Double_t invMassD = 0;
976 
977  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
978  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
979  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
980 
981  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
982  //AliDebug(2, Form("MC label is %d", mcLab));
983  if (mcLab >= 0) {
984  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
985 
986  if (aodMcPart) {
987  if (fRejectedOrigin && fMCMode == kSignalOnly) {
988  EMesonOrigin_t origin = CheckOrigin(aodMcPart, fMCContainer->GetArray());
989 
990  if ((origin & fRejectedOrigin) == origin) return kFALSE;
991  }
992 
993  MCtruthPdgCode = aodMcPart->PdgCode();
994  //AliDebug(2, Form("MC truth pdg code is %d",MCtruthPdgCode));
995  }
996  }
997  }
998 
999  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
1000  if (fMCMode == kNoMC ||
1001  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
1002  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
1003  // 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)
1004  invMassD = DstarCand->InvMassDstarKpipi();
1005  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
1006  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
1007  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
1008  return kTRUE;
1009  }
1010  else { // conditions above not passed, so return FALSE
1011  return kFALSE;
1012  }
1013 }
1014 
1022 {
1023  if (!part) return kDecayOther;
1024  if (!mcArray) return kDecayOther;
1025 
1027 
1028  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
1029 
1030  if (part->GetNDaughters() == 2) {
1031 
1032  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
1033  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
1034 
1035  if (!d1 || !d2) {
1036  return decay;
1037  }
1038 
1039  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
1040  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
1041 
1042  if (absPdgPart == 421) { // D0 -> K pi
1043 
1044  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
1045  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
1046  decay = kDecayD0toKpi;
1047  }
1048  }
1049 
1050  if (absPdgPart == 413) { // D* -> D0 pi
1051 
1052  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
1053  Int_t D0decay = CheckDecayChannel(d1, mcArray);
1054  if (D0decay == kDecayD0toKpi) {
1055  decay = kDecayDStartoKpipi;
1056  }
1057  }
1058 
1059  if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
1060  Int_t D0decay = CheckDecayChannel(d2, mcArray);
1061  if (D0decay == kDecayD0toKpi) {
1062  decay = kDecayDStartoKpipi;
1063  }
1064  }
1065  }
1066  }
1067 
1068  return decay;
1069 }
1070 
1078 {
1079  // Checks whether the mother of the particle comes from a charm or a bottom quark.
1080 
1081  if (!part) return kUnknownQuark;
1082  if (!mcArray) return kUnknownQuark;
1083 
1084  Int_t pdgGranma = 0;
1085  Int_t mother = part->GetMother();
1086  Int_t istep = 0;
1087  Int_t abspdgGranma = 0;
1088  Bool_t isFromB = kFALSE;
1089  Bool_t isQuarkFound = kFALSE;
1090 
1091  while (mother >= 0) {
1092  istep++;
1093  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
1094  if (mcGranma >= 0) {
1095  pdgGranma = mcGranma->GetPdgCode();
1096  abspdgGranma = TMath::Abs(pdgGranma);
1097  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1098  isFromB = kTRUE;
1099  }
1100 
1101  if (abspdgGranma == 4 || abspdgGranma == 5) isQuarkFound = kTRUE;
1102  mother = mcGranma->GetMother();
1103  }
1104  else {
1105  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::CheckOrigin", "Could not retrieve mother particle %d!", mother);
1106  break;
1107  }
1108  }
1109 
1110  if (isQuarkFound) {
1111  if (isFromB) {
1112  return kFromBottom;
1113  }
1114  else {
1115  return kFromCharm;
1116  }
1117  }
1118  else {
1119  return kUnknownQuark;
1120  }
1121 }
1122 
1125 {
1126  fDmesonJets.clear();
1127 
1128  if (fMCMode == kMCTruth) {
1129  RunParticleLevelAnalysis();
1130  }
1131  else {
1132  RunDetectorLevelAnalysis();
1133  }
1134 }
1135 
1138 {
1139  const Int_t nD = fCandidateArray->GetEntriesFast();
1140 
1141  AliDmesonJetInfo DmesonJet;
1142 
1143  Int_t nAccCharm = 0;
1144  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
1145  Int_t isSelected = 0;
1146 
1147  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
1148  if (!charmCand) continue;
1149 
1150  Int_t nprongs = charmCand->GetNProngs();
1151 
1152  if (fCandidateType == kDstartoKpipi) {
1153  if (!charmCand->InheritsFrom("AliAODRecoCascadeHF")) {
1154  ::Error("AliAnalysisTaskDmesonJets::AnalysisParams::RunDetectorLevelAnalysis","Candidate type is D* but object type is wrong (should be AliAODRecoCascadeHF)");
1155  continue;
1156  }
1157  }
1158 
1159  // region of interest + cuts
1160  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
1161 
1162  //candidate selected by cuts and PID
1163  isSelected = fRDHFCuts->IsSelected(charmCand, AliRDHFCuts::kAll, fAodEvent); //selected
1164 
1165  if (!isSelected) continue;
1166 
1167  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
1168  DmesonJet.Reset();
1169  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
1170  for (std::vector<AliHFJetDefinition>::iterator itdef = fJetDefinitions.begin(); itdef != fJetDefinitions.end(); itdef++) {
1171  if (!FindJet(charmCand, DmesonJet, *itdef)) {
1172  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
1173  (*itdef).GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
1174  }
1175  }
1176  fDmesonJets.push_back(DmesonJet);
1177  }
1178  }
1179  nAccCharm++;
1180  } // end of D cand loop
1181 
1182  TString hname;
1183 
1184  hname = TString::Format("%s/fHistNAcceptedDmesons", GetName());
1185  fHistManager->FillTH1(hname, nAccCharm);
1186 
1187  hname = TString::Format("%s/fHistNDmesons", GetName());
1188  fHistManager->FillTH1(hname, nD);
1189 }
1190 
1201 {
1202  TString hname;
1203 
1205  fFastJetWrapper->SetR(jetDef.fRadius);
1208 
1209  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
1210 
1211  if (fTrackContainer && jetDef.fJetType != AliJetContainer::kNeutralJet) {
1212  fTrackContainer->SetDMesonCandidate(Dcand);
1213  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
1214  AddInputVectors(fTrackContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1215 
1216  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
1217  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
1218  const TObjArray daughterNotInJet = fTrackContainer->GetDaughterList();
1219  for (Int_t i = 0; i < daughterNotInJet.GetEntriesFast(); i++) {
1220  AliVParticle* daughter = static_cast<AliVParticle*>(daughterNotInJet.At(i));
1221  if (!daughter) continue;
1222  histDaughterNotInJet->Fill(daughter->Pt());
1223  }
1224  }
1225 
1226  if (fClusterContainer && jetDef.fJetType != AliJetContainer::kChargedJet) {
1227  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
1228  AddInputVectors(fClusterContainer, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1229  }
1230 
1231  // run jet finder
1232  fFastJetWrapper->Run();
1233 
1234  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1235 
1236  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1237  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1238 
1239  Bool_t isDmesonJet = kFALSE;
1240 
1241  Double_t maxChPt = 0;
1242  Double_t maxNePt = 0;
1243  Double_t totalNeutralPt = 0;
1244 
1245  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1246  if (constituents[ic].user_index() == 0) {
1247  isDmesonJet = kTRUE;
1248  }
1249  else if (constituents[ic].user_index() >= 100) {
1250  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
1251  }
1252  else if (constituents[ic].user_index() <= -100) {
1253  totalNeutralPt += constituents[ic].pt();
1254  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
1255  }
1256  }
1257 
1258  if (isDmesonJet) {
1259  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1260  DmesonJet.fJets[jetDef.GetName()].fNConstituents = constituents.size();
1261  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
1262  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
1263  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
1264 
1265  return kTRUE;
1266  }
1267  }
1268 
1269  return kFALSE;
1270 }
1271 
1276 {
1277  AliTLorentzVector part;
1278  for (Int_t i = 0; i < cont->GetNEntries(); i++) {
1279  cont->GetMomentum(part, i);
1280  UInt_t rejectionReason = 0;
1281  if (!cont->AcceptObject(i, rejectionReason)) {
1282  rejectHist->Fill(cont->GetRejectionReasonBitPosition(rejectionReason), part.Pt());
1283  continue;
1284  }
1285  Int_t uid = offset >= 0 ? i : -i;
1286  uid += offset;
1287  fFastJetWrapper->AddInputVector(part.Px(), part.Py(), part.Pz(), part.E(), uid);
1288  }
1289 }
1290 
1293 {
1294  TString hname;
1295 
1296  fMCContainer->SetSpecialPDG(fCandidatePDG);
1297  fMCContainer->SetRejectedOriginMap(fRejectedOrigin);
1298  fMCContainer->SetAcceptedDecayMap(fAcceptedDecay);
1299 
1300  std::map<int, AliDmesonJetInfo> dMesonJets;
1301 
1302  for (std::vector<AliHFJetDefinition>::iterator itdef = fJetDefinitions.begin(); itdef != fJetDefinitions.end(); itdef++) {
1303  AliHFJetDefinition* jetDef = &(*itdef);
1304 
1306  fFastJetWrapper->SetR(jetDef->fRadius);
1309 
1310  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef->GetName());
1311  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
1312 
1313  fFastJetWrapper->Run();
1314 
1315  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
1316 
1317  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
1318  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
1319 
1320  Bool_t isDmesonJet = kFALSE;
1321 
1322  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
1323  Int_t iPart = constituents[ic].user_index() - 100;
1324  AliVParticle* part = fMCContainer->GetParticle(iPart);
1325  if (!part) {
1326  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
1327  continue;
1328  }
1329  if (part->PdgCode() == fCandidatePDG) {
1330  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = dMesonJets.find(iPart);
1331  if (dMesonJetIt == dMesonJets.end()) {
1332  std::pair<int, AliDmesonJetInfo> element;
1333  element.first = iPart;
1334  element.second = AliDmesonJetInfo();
1335  dMesonJetIt = dMesonJets.insert(element).first;
1336  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
1337  }
1338 
1339  (*dMesonJetIt).second.fJets[jetDef->GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
1340  break;
1341  }
1342  }
1343  }
1344  }
1345 
1346  for (std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = dMesonJets.begin();
1347  dMesonJetIt != dMesonJets.end();
1348  dMesonJetIt++) {
1349  fDmesonJets.push_back((*dMesonJetIt).second);
1350  }
1351 }
1352 
1357 {
1358  TString classname;
1359  switch (fCandidateType) {
1360  case kD0toKpi:
1361  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1362  fCurrentDmesonJetInfo = new AliD0InfoSummary();
1363  break;
1364  case kDstartoKpipi:
1365  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1366  fCurrentDmesonJetInfo = new AliDStarInfoSummary();
1367  break;
1368  }
1369  fTree = new TTree(GetName(), GetName());
1370  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo, 32000, 0);
1371  fCurrentJetInfo = new AliJetInfoSummary*[fJetDefinitions.size()];
1372  for (Int_t i = 0; i < fJetDefinitions.size(); i++) {
1373  fCurrentJetInfo[i] = new AliJetInfoSummary();
1374  fTree->Branch(fJetDefinitions[i].GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1375  }
1376 
1377  return fTree;
1378 }
1379 
1384 {
1385  TString hname;
1386 
1387  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
1388 
1389  for (std::vector<AliHFJetDefinition>::const_iterator it = fJetDefinitions.begin(); it != fJetDefinitions.end(); it++) {
1390  const AliHFJetDefinition* jetDef = &(*it);
1391 
1392  AliDebug(2,Form("Now working on '%s'", jetDef->GetName()));
1393 
1394  Double_t radius = jetDef->fRadius;
1395 
1396  TString title[30] = {""};
1397  Int_t nbins[30] = {0 };
1398  Double_t min [30] = {0.};
1399  Double_t max [30] = {0.};
1400  Int_t dim = 0 ;
1401 
1402  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
1403  nbins[dim] = nPtBins;
1404  min[dim] = 0;
1405  max[dim] = fMaxPt;
1406  dim++;
1407 
1408  if ((enabledAxis & kPositionD) != 0) {
1409  title[dim] = "#eta_{D}";
1410  nbins[dim] = 50;
1411  min[dim] = -1;
1412  max[dim] = 1;
1413  dim++;
1414 
1415  title[dim] = "#phi_{D} (rad)";
1416  nbins[dim] = 150;
1417  min[dim] = 0;
1418  max[dim] = TMath::TwoPi();
1419  dim++;
1420  }
1421 
1422  if ((enabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1423  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
1424  nbins[dim] = fNMassBins;
1425  min[dim] = fMinMass;
1426  max[dim] = fMaxMass;
1427  dim++;
1428  }
1429 
1430  if (fCandidateType == kD0toKpi) {
1431  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1432  nbins[dim] = fNMassBins;
1433  min[dim] = fMinMass;
1434  max[dim] = fMaxMass;
1435  dim++;
1436  }
1437 
1438  if ((enabledAxis & k2ProngInvMass) != 0 && fCandidateType == kDstartoKpipi) {
1439  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
1440  nbins[dim] = fNMassBins;
1441  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
1442  dim++;
1443  }
1444 
1445  if (fCandidateType == kDstartoKpipi) {
1446  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
1447  nbins[dim] = fNMassBins*6;
1448  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
1449 
1450  // subtract mass of D0
1451  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1452  min[dim] -= D0mass;
1453  max[dim] -= D0mass;
1454 
1455  dim++;
1456  }
1457 
1458  if ((enabledAxis & kSoftPionPt) != 0 && fCandidateType == kDstartoKpipi) {
1459  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
1460  nbins[dim] = 100;
1461  min[dim] = 0;
1462  max[dim] = 25;
1463  dim++;
1464  }
1465 
1466  title[dim] = "#it{z}_{D}";
1467  nbins[dim] = 110;
1468  min[dim] = 0;
1469  max[dim] = 1.10;
1470  dim++;
1471 
1472  if ((enabledAxis & kDeltaR) != 0) {
1473  title[dim] = "#Delta R_{D-jet}";
1474  nbins[dim] = 100;
1475  min[dim] = 0;
1476  max[dim] = radius * 1.5;
1477  dim++;
1478  }
1479 
1480  if ((enabledAxis & kDeltaEta) != 0) {
1481  title[dim] = "#eta_{D} - #eta_{jet}";
1482  nbins[dim] = 100;
1483  min[dim] = -radius * 1.2;
1484  max[dim] = radius * 1.2;
1485  dim++;
1486  }
1487 
1488  if ((enabledAxis & kDeltaPhi) != 0) {
1489  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1490  nbins[dim] = 100;
1491  min[dim] = -radius * 1.2;
1492  max[dim] = radius * 1.2;
1493  dim++;
1494  }
1495 
1496  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1497  nbins[dim] = nPtBins;
1498  min[dim] = 0;
1499  max[dim] = fMaxPt;
1500  dim++;
1501 
1502  if ((enabledAxis & kPositionJet) != 0) {
1503  title[dim] = "#eta_{jet}";
1504  nbins[dim] = 50;
1505  min[dim] = -1;
1506  max[dim] = 1;
1507  dim++;
1508 
1509  title[dim] = "#phi_{jet} (rad)";
1510  nbins[dim] = 150;
1511  min[dim] = 0;
1512  max[dim] = TMath::TwoPi();
1513  dim++;
1514  }
1515 
1516  if ((enabledAxis & kJetConstituents) != 0) {
1517  title[dim] = "No. of constituents";
1518  nbins[dim] = 50;
1519  min[dim] = -0.5;
1520  max[dim] = 49.5;
1521  dim++;
1522  }
1523 
1524  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef->GetName());
1525  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1526  for (Int_t j = 0; j < dim; j++) {
1527  h->GetAxis(j)->SetTitle(title[j]);
1528  }
1529  }
1530 }
1531 
1536 {
1537  TString hname;
1538 
1539  for (Int_t id = 0; id < fDmesonJets.size(); id++) {
1540  fCurrentDmesonJetInfo->Set(fDmesonJets[id]);
1541  Int_t accJets = 0;
1542  for (UInt_t ij = 0; ij < fJetDefinitions.size(); ij++) {
1543  fCurrentJetInfo[ij]->Reset();
1544  AliJetInfo* jet = fDmesonJets[id].GetJet(fJetDefinitions[ij].GetName());
1545  if (!jet) continue;
1546  if (applyKinCuts && !fJetDefinitions[ij].IsJetInAcceptance(*jet)) {
1547  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), fJetDefinitions[ij].GetName());
1548  fHistManager->FillTH1(hname, jet->Pt());
1549  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), fJetDefinitions[ij].GetName());
1550  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1551  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), fJetDefinitions[ij].GetName());
1552  fHistManager->FillTH1(hname, jet->Eta());
1553  continue;
1554  }
1555  fCurrentJetInfo[ij]->Set(fDmesonJets[id], fJetDefinitions[ij].GetName());
1556  accJets++;
1557  }
1558  if (accJets > 0) {
1559  fTree->Fill();
1560  }
1561  else {
1562  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1563  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Pt());
1564  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1565  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Phi_0_2pi());
1566  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1567  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Eta());
1568  if (fCandidateType == kD0toKpi) {
1569  hname = TString::Format("%s/fHistRejectedDMesonInvMass", GetName());
1570  fHistManager->FillTH1(hname, fDmesonJets[id].fD.M());
1571  }
1572  else if (fCandidateType == kDstartoKpipi) {
1573  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", GetName());
1574  fHistManager->FillTH1(hname, fDmesonJets[id].fInvMass2Prong);
1575 
1576  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", GetName());
1577  fHistManager->FillTH1(hname, fDmesonJets[id].fD.M() - fDmesonJets[id].fInvMass2Prong);
1578  }
1579  }
1580  }
1581  return kTRUE;
1582 }
1583 
1588 {
1589  TString hname;
1590 
1591  for (Int_t id = 0; id < fDmesonJets.size(); id++) {
1592  if (!IsAnyJetInAcceptance(fDmesonJets[id])) {
1593  hname = TString::Format("%s/fHistRejectedDMesonPt", GetName());
1594  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Pt());
1595  hname = TString::Format("%s/fHistRejectedDMesonPhi", GetName());
1596  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Phi_0_2pi());
1597  hname = TString::Format("%s/fHistRejectedDMesonEta", GetName());
1598  fHistManager->FillTH1(hname, fDmesonJets[id].fD.Eta());
1599  }
1600  }
1601 
1602  for (std::vector<AliHFJetDefinition>::iterator itdef = fJetDefinitions.begin(); itdef != fJetDefinitions.end(); itdef++) {
1603 
1604  AliHFJetDefinition* jetDef = &(*itdef);
1605 
1606  hname = TString::Format("%s/%s/fDmesonJets", GetName(), jetDef->GetName());
1607  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
1608 
1609  for (Int_t id = 0; id < fDmesonJets.size(); id++) {
1610  const AliJetInfo* jet = fDmesonJets[id].GetJet(jetDef->GetName());
1611  if (!jet) continue;
1612  if (!jetDef->IsJetInAcceptance(*jet)) {
1613  hname = TString::Format("%s/%s/fHistRejectedJetPt", GetName(), jetDef->GetName());
1614  fHistManager->FillTH1(hname, jet->Pt());
1615  hname = TString::Format("%s/%s/fHistRejectedJetPhi", GetName(), jetDef->GetName());
1616  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1617  hname = TString::Format("%s/%s/fHistRejectedJetEta", GetName(), jetDef->GetName());
1618  fHistManager->FillTH1(hname, jet->Eta());
1619  continue;
1620  }
1621  FillHnSparse(h, fDmesonJets[id], (*itdef).GetName());
1622  }
1623  }
1624 
1625  return kTRUE;
1626 }
1627 
1633 Bool_t AliAnalysisTaskDmesonJets::AnalysisEngine::FillHnSparse(THnSparse* h, const AliDmesonJetInfo& DmesonJet, std::string n)
1634 {
1635  // Fill the THnSparse histogram.
1636 
1637  Double_t contents[30] = {0.};
1638 
1639  Double_t z = DmesonJet.GetZ(n);
1640  Double_t deltaPhi = 0;
1641  Double_t deltaEta = 0;
1642  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
1643 
1644  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
1645  if (it == DmesonJet.fJets.end()) return kFALSE;
1646 
1647  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
1648  TString title(h->GetAxis(i)->GetTitle());
1649  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
1650  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
1651  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
1652  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
1653  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
1654  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
1655  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
1656  else if (title=="#it{z}_{D}") contents[i] = z;
1657  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
1658  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
1659  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
1660  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
1661  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
1662  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
1663  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
1664  else AliWarning(Form("Unable to fill dimension '%s'!",title.Data()));
1665  }
1666 
1667  h->Fill(contents);
1668 
1669  return kTRUE;
1670 }
1671 
1672 // Definitions of class AliAnalysisTaskDmesonJets
1673 
1677 
1681  fAnalysisEngines(),
1682  fEnabledAxis(0),
1683  fTreeOutput(kFALSE),
1684  fHistManager(),
1685  fApplyKinematicCuts(kTRUE),
1686  fAodEvent(0),
1687  fFastJetWrapper(0)
1688 {
1689  SetMakeGeneralHistograms(kTRUE);
1690 }
1691 
1696  AliAnalysisTaskEmcalLight(name, kTRUE),
1697  fAnalysisEngines(),
1698  fEnabledAxis(k2ProngInvMass),
1699  fTreeOutput(kFALSE),
1700  fHistManager(name),
1701  fApplyKinematicCuts(kTRUE),
1702  fAodEvent(0),
1703  fFastJetWrapper(0)
1704 {
1705  SetMakeGeneralHistograms(kTRUE);
1706 }
1707 
1710 {
1711  if (fFastJetWrapper) delete fFastJetWrapper;
1712 }
1713 
1721 {
1722  AliRDHFCuts* analysiscuts = 0;
1723  TFile* filecuts = TFile::Open(cutfname);
1724  if (!filecuts || filecuts->IsZombie()) {
1725  ::Warning("AddTaskDmesonJets", "Input file not found: will use std cuts.");
1726  filecuts = 0;
1727  }
1728 
1729  if (filecuts) {
1730  analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
1731  if (!analysiscuts) {
1732  ::Warning("AddTaskDmesonJetCorr", "Could not find analysis cuts '%s' in '%s'. Using std cuts.", cutsname.Data(), cutfname.Data());
1733  }
1734  }
1735 
1736  return analysiscuts;
1737 }
1738 
1748 {
1750  return AddAnalysisEngine(type, MCmode, jetDef, cutfname);
1751 }
1752 
1762 {
1763  AliRDHFCuts* cuts = 0;
1764 
1765  if (!cutfname.IsNull()) {
1766  TString cutsname;
1767 
1768  switch (type) {
1769  case kD0toKpi :
1770  cutsname = "D0toKpiCuts";
1771  break;
1772  case kDstartoKpipi :
1773  cutsname = "DStartoKpipiCuts";
1774  break;
1775  default:
1776  return 0;
1777  }
1778 
1779  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
1780  }
1781 
1782  AnalysisEngine eng(type, MCmode, cuts);
1783 
1784  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
1785 
1786  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
1787  eng.AddJetDefinition(jetDef);
1788  it = fAnalysisEngines.insert(it, eng);
1789  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(jetDef), fAnalysisEngines.size());
1790  }
1791  else {
1792  AnalysisEngine* found_eng = &(*it);
1793  ::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());
1794  found_eng->AddJetDefinition(jetDef);
1795 
1796  if (cuts && found_eng->fRDHFCuts != 0) {
1797  ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
1798  found_eng->SetRDHFCuts(cuts);
1799  }
1800  }
1801 
1802  return &(*it);
1803 }
1804 
1805 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
1806 {
1807  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
1808  while (it != fAnalysisEngines.end() && (*it) < eng) it++;
1809  return it;
1810 }
1811 
1814 {
1815  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
1816 
1818 
1819  // Define histograms
1820  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
1821 
1822  TString hname;
1823  TString htitle;
1824  TH1* h = 0;
1825 
1826  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
1827  for (std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin(); it != fAnalysisEngines.end(); it++) {
1828  AnalysisEngine* param = &(*it);
1829  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param->GetName(), param->fJetDefinitions.size());
1830 
1831  fHistManager.CreateHistoGroup(param->GetName());
1832 
1833  param->fHistManager = &fHistManager;
1834 
1835  hname = TString::Format("%s/fHistNAcceptedDmesons", param->GetName());
1836  htitle = hname + ";Number of D accepted meson candidates;counts";
1837  h = fHistManager.CreateTH1(hname, htitle, 51, -0.5, 50.5);
1838 
1839  hname = TString::Format("%s/fHistNDmesons", param->GetName());
1840  htitle = hname + ";Number of D meson candidates;counts";
1841  h = fHistManager.CreateTH1(hname, htitle, 101, -0.5, 100.5);
1842 
1843  hname = TString::Format("%s/fHistNEvents", param->GetName());
1844  htitle = hname + ";Event status;counts";
1845  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
1846  h->GetXaxis()->SetBinLabel(1, "Accepted");
1847  h->GetXaxis()->SetBinLabel(2, "Rejected");
1848 
1849  hname = TString::Format("%s/fHistEventRejectionReasons", param->GetName());
1850  htitle = hname + ";Rejection reason;counts";
1851  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
1852 
1853  hname = TString::Format("%s/fHistRejectedDMesonPt", param->GetName());
1854  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
1855  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
1856 
1857  hname = TString::Format("%s/fHistRejectedDMesonEta", param->GetName());
1858  htitle = hname + ";#it{#eta}_{D};counts";
1859  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
1860 
1861  hname = TString::Format("%s/fHistRejectedDMesonPhi", param->GetName());
1862  htitle = hname + ";#it{#phi}_{D};counts";
1863  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
1864 
1865  if (param->fCandidateType == kD0toKpi) {
1866  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param->GetName());
1867  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
1868  fHistManager.CreateTH1(hname, htitle, param->fNMassBins, param->fMinMass, param->fMaxMass);
1869  }
1870  else if (param->fCandidateType == kDstartoKpipi) {
1871  Double_t min = 0;
1872  Double_t max = 0;
1873 
1874  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param->GetName());
1875  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
1876  CalculateMassLimits(param->fMaxMass - param->fMinMass, 421, param->fNMassBins, min, max);
1877  fHistManager.CreateTH1(hname, htitle, param->fNMassBins, min, max);
1878 
1879  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
1880  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param->GetName());
1881  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
1882  CalculateMassLimits(0.20, 413, param->fNMassBins*6, min, max);
1883  fHistManager.CreateTH1(hname, htitle, param->fNMassBins*6, min-D0mass, max-D0mass);
1884  }
1885 
1886  for (std::vector<AliHFJetDefinition>::iterator itdef = param->fJetDefinitions.begin(); itdef != param->fJetDefinitions.end(); itdef++) {
1887  AliHFJetDefinition* jetDef = &(*itdef);
1888  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef->GetName());
1889 
1890  fHistManager.CreateHistoGroup(jetDef->GetName(), param->GetName());
1891 
1892  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param->GetName(), jetDef->GetName());
1893  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
1894  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
1895  SetRejectionReasonLabels(h->GetXaxis());
1896 
1897  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param->GetName(), jetDef->GetName());
1898  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
1899  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
1900  SetRejectionReasonLabels(h->GetXaxis());
1901 
1902  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param->GetName(), jetDef->GetName());
1903  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
1904  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
1905  SetRejectionReasonLabels(h->GetXaxis());
1906 
1907  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param->GetName(), jetDef->GetName());
1908  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
1909  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
1910 
1911  hname = TString::Format("%s/%s/fHistRejectedJetPt", param->GetName(), jetDef->GetName());
1912  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
1913  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
1914 
1915  hname = TString::Format("%s/%s/fHistRejectedJetEta", param->GetName(), jetDef->GetName());
1916  htitle = hname + ";#it{#eta}_{jet};counts";
1917  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
1918 
1919  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param->GetName(), jetDef->GetName());
1920  htitle = hname + ";#it{#phi}_{jet};counts";
1921  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
1922  }
1923  if (fTreeOutput) {
1924  TTree* tree = param->BuildTree();
1925  fOutput->Add(tree);
1926  }
1927  else {
1928  param->BuildHnSparse(fEnabledAxis);
1929  }
1930  }
1931 
1932  fOutput->Add(fHistManager.GetListOfHistograms());
1933 
1934  PostData(1, fOutput);
1935 }
1936 
1940 {
1942 
1943  // Load the event
1944  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
1945 
1946  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
1947 
1948  fFastJetWrapper->SetAreaType(fastjet::active_area);
1950 
1951  if (!fAodEvent) {
1952  AliError(Form("This task need an AOD event! Task '%s' will be disabled!", GetName()));
1953  return;
1954  }
1955 
1956  for (std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin(); it != fAnalysisEngines.end(); it++) {
1957  AnalysisEngine* params = &(*it);
1958 
1959  params->fAodEvent = fAodEvent;
1960  params->fFastJetWrapper = fFastJetWrapper;
1961  params->Init(fGeom, fAodEvent->GetRunNumber());
1962 
1963  if (params->fMCMode != kMCTruth) {
1964  params->fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params->fBranchName.Data()));
1965 
1966  if (params->fCandidateArray) {
1967  if (!params->fCandidateArray->GetClass()->InheritsFrom("AliAODRecoDecayHF2Prong")) {
1968  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
1969  "%s: Objects of type %s in %s are not inherited from AliAODRecoDecayHF2Prong! Task will be disabled!",
1970  GetName(), params->fCandidateArray->GetClass()->GetName(), params->fCandidateArray->GetName());
1971  params->fCandidateArray = 0;
1972  params->fInhibit = kTRUE;
1973  }
1974  }
1975  else {
1976  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
1977  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
1978  params->fBranchName.Data(), params->GetName());
1979  params->fInhibit = kTRUE;
1980  }
1981  }
1982 
1983  if (params->fMCMode != kNoMC) {
1984  params->fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(0));
1985 
1986  if (!params->fMCContainer) params->fMCContainer = dynamic_cast<AliHFAODMCParticleContainer*>(GetParticleContainer(1));
1987 
1988  if (!params->fMCContainer) {
1989  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
1990  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
1991  params->GetName());
1992  params->fInhibit = kTRUE;
1993  }
1994  }
1995 
1996  if (params->fMCMode != kMCTruth) {
1997  params->fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(0));
1998  if (!params->fTrackContainer) params->fTrackContainer = dynamic_cast<AliHFTrackContainer*>(GetParticleContainer(1));
1999 
2001 
2002  if (!params->fTrackContainer && !params->fClusterContainer) {
2003  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
2004  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
2005  params->GetName());
2006  params->fInhibit = kTRUE;
2007  }
2008  }
2009  }
2010 }
2011 
2016 {
2017  if (!fAodEvent) return kFALSE;
2018 
2019  TString hname;
2020 
2021  // fix for temporary bug in ESDfilter
2022  // the AODs with null vertex pointer didn't pass the PhysSel
2023  if (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001) return kFALSE;
2024 
2025  for (std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin(); it != fAnalysisEngines.end(); it++) {
2026  AnalysisEngine* eng = &(*it);
2027 
2028  if (eng->fInhibit) continue;
2029 
2030  //Event selection
2031  hname = TString::Format("%s/fHistNEvents", eng->GetName());
2032  Bool_t iseventselected = eng->fRDHFCuts->IsEventSelected(fAodEvent);
2033  if (!iseventselected) {
2034  fHistManager.FillTH1(hname, "Rejected");
2035  hname = TString::Format("%s/fHistEventRejectionReasons", eng->GetName());
2036  UInt_t bitmap = eng->fRDHFCuts->GetEventRejectionBitMap();
2037  TString label;
2038  do {
2039  label = GetHFEventRejectionReasonLabel(bitmap);
2040  if (label.IsNull()) break;
2041  fHistManager.FillTH1(hname, label);
2042  } while (true);
2043 
2044  continue;
2045  }
2046 
2047  fHistManager.FillTH1(hname, "Accepted");
2048 
2049  AliDebug(2, "Event selected");
2050 
2051  eng->RunAnalysis();
2052  }
2053  return kTRUE;
2054 }
2055 
2060 {
2061  TString hname;
2062  for (std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin(); it != fAnalysisEngines.end(); it++) {
2063  AnalysisEngine* param = &(*it);
2064 
2065  if (param->fInhibit) continue;
2066 
2067  if (fTreeOutput) {
2068  param->FillTree(fApplyKinematicCuts);
2069  }
2070  else {
2072  }
2073  }
2074  return kTRUE;
2075 }
2076 
2084 void AliAnalysisTaskDmesonJets::CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t& minMass, Double_t& maxMass)
2085 {
2086  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
2087 
2088  Double_t mass = part->Mass();
2089 
2090  // To make sure that the PDG mass value is not at the edge of a bin
2091  if (nbins % 2 == 0) {
2092  minMass = mass - range / 2 - range / nbins / 2;
2093  maxMass = mass + range / 2 - range / nbins / 2;
2094  }
2095  else {
2096  minMass = mass - range / 2;
2097  maxMass = mass + range / 2;
2098  }
2099 }
2100 
2107 {
2108  static TString label;
2109  label = "";
2110 
2111  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
2112  label = "NotSelTrigger";
2113  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
2114  return label.Data();
2115  }
2116  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
2117  label = "NoVertex";
2118  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
2119  return label.Data();
2120  }
2121  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
2122  label = "TooFewVtxContrib";
2123  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
2124  return label.Data();
2125  }
2126  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
2127  label = "ZVtxOutFid";
2128  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
2129  return label.Data();
2130  }
2131  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
2132  label = "Pileup";
2133  bitmap &= ~BIT(AliRDHFCuts::kPileup);
2134  return label.Data();
2135  }
2136  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
2137  label = "OutsideCentrality";
2138  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
2139  return label.Data();
2140  }
2141  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
2142  label = "PhysicsSelection";
2143  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
2144  return label.Data();
2145  }
2146  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
2147  label = "BadSPDVertex";
2148  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
2149  return label.Data();
2150  }
2151  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
2152  label = "ZVtxSPDOutFid";
2153  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
2154  return label.Data();
2155  }
2156  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
2157  label = "CentralityFlattening";
2158  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
2159  return label.Data();
2160  }
2161  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
2162  label = "BadTrackV0Correl";
2163  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
2164  return label.Data();
2165  }
2166 
2167  return label.Data();
2168 }
Int_t pdg
void Print() const
Prints the content of this object in the standard output.
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
AliEMCALGeometry * fGeom
!emcal geometry
Analysis task for D meson jets.
virtual void Set(const AliDmesonJetInfo &source)
Invariant mass of the D0 meson candidate in GeV/c2.
virtual Int_t Run()
Bool_t IsJetInAcceptance(const AliJetInfo &jet) const
Double_t fSoftPionPt
! Transverse momentum of the soft pion of the D* candidate
const char * title
Definition: MakeQAPdf.C:26
AliRDHFCuts * fRDHFCuts
D meson candidates cuts.
Lightweight class that encapsulates D meson jets.
Bool_t fTreeOutput
If true, output will be posted in a TTree rather than a THnSparse.
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
ERecoScheme_t fRecoScheme
Jet recombination scheme (pt scheme, E scheme, ...)
Double_t fMaxNeutralPt
Transverse momentum of the leading neutral particle (or cluster)
void RunDetectorLevelAnalysis()
Run a detector level analysis.
Double_t mass
virtual void UserCreateOutputObjects()
Creates the output containers.
Lightweight class that encapsulates D meson jets.
Double_t fInvMass2Prong
! 2-prong mass of the D* candidate (w/o the soft pion)
static void CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t &minMass, Double_t &maxMass)
AliAODEvent * fAodEvent
! AOD event
Bool_t FindJet(AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, AliHFJetDefinition &jetDef)
UInt_t fEnabledAxis
Use bit defined in EAxis_t to enable axis in the THnSparse.
EMCMode_t fMCMode
MC mode: No MC (data and MC detector level), background-only (MC), signal-only (MC), MC truth (particle level)
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:93
AnalysisEngine()
This is the default constructor, used for ROOT I/O purposes.
void Init(const AliEMCALGeometry *const geom, Int_t runNumber)
Initialize the analysis engine.
static FJRecoScheme ConvertToFJRecoScheme(ERecoScheme_t reco)
static TString GenerateJetName(EJetType_t jetType, EJetAlgo_t jetAlgo, ERecoScheme_t recoScheme, Double_t radius, AliParticleContainer *partCont, AliClusterContainer *clusCont, TString tag)
void SetR(Double_t r)
Definition: AliFJWrapper.h:98
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
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:33
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.
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)
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:92
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
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:96
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
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Double_t minMass
void RunParticleLevelAnalysis()
Run a particle level analysis.
Bool_t fApplyKinematicCuts
Apply jet kinematic cuts.
const char * GetName() const
Generate a name for this jet definition.
void SetDetectorJetEtaPhiRange(const AliEMCALGeometry *const geom, Int_t run)
Bool_t ExtractRecoDecayAttributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
Double_t InvMassD0() const
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="")
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:94
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.
const Int_t nPtBins
AliHFTrackContainer * fTrackContainer
! Track container
Bool_t IsAnyJetInAcceptance(const AliDmesonJetInfo &dMesonJet) const
Container for jet within the EMCAL jet framework.
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.