AliPhysics  64f4410 (64f4410)
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 //C++
17 #include <sstream>
18 #include <array>
19 
20 // Root
21 #include <TClonesArray.h>
22 #include <TDatabasePDG.h>
23 #include <TParticlePDG.h>
24 #include <TVector3.h>
25 #include <THnSparse.h>
26 #include <TParticle.h>
27 #include <TMath.h>
28 #include <THashList.h>
29 #include <TFile.h>
30 #include <TRandom3.h>
31 
32 // Aliroot general
33 #include "AliLog.h"
34 #include "AliEMCALGeometry.h"
35 #include "AliAnalysisManager.h"
36 #include "AliVEventHandler.h"
37 
38 // Aliroot HF
39 #include "AliAODRecoCascadeHF.h"
41 #include "AliRDHFCutsD0toKpi.h"
44 #include "AliHFTrackContainer.h"
45 #include "AliAnalysisVertexingHF.h"
46 
47 // Aliroot EMCal jet framework
48 #include "AliEmcalJetTask.h"
49 #include "AliEmcalJet.h"
50 #include "AliJetContainer.h"
51 #include "AliParticleContainer.h"
52 #include "AliEmcalParticle.h"
53 #include "AliFJWrapper.h"
54 #include "AliRhoParameter.h"
55 
57 
58 AliAnalysisTaskDmesonJets::AliEventNotFound::AliEventNotFound(const std::string& class_name, const std::string& method_name) :
59  std::exception(),
60  fClassName(class_name),
61  fAccessMethodName(method_name)
62 {
63  std::stringstream what_str;
64  what_str << "ALICE event not found in class '" << fClassName << "' using method '" << method_name << "'.";
65  fWhat = what_str.str();
66 }
67 
68 #if !(defined(__CINT__) || defined(__MAKECINT__))
70 {
71  return fWhat.c_str();
72 }
73 #endif
74 
75 
76 // Definitions of class AliAnalysisTaskDmesonJets::AliEventInfoSummary
77 
81 
87  fWeight(1),
88  fPtHard(0)
89 {
90  Set(event);
91 }
92 
95 {
96  fWeight = 1;
97  fPtHard = 0;
98 }
99 
105 {
106  fWeight = event.fWeight;
107  fPtHard = event.fPtHard;
108 }
109 
110 
111 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfo
112 
116 
124 {
125  dphi = TVector2::Phi_mpi_pi(fMomentum.Phi() - jet.Phi());;
126  deta = fMomentum.Eta() - jet.Eta();
127  return TMath::Sqrt(dphi*dphi + deta*deta);
128 }
129 
135 {
136  Double_t deta = 0;
137  Double_t dphi = 0;
138  return GetDistance(jet, deta, dphi);
139 }
140 
141 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonJetInfo
142 
146 
149  fDmesonParticle(0),
150  fD(),
151  fSoftPionPt(0),
152  fInvMass2Prong(0),
153  fJets(),
154  fMCLabel(-1),
155  fReconstructed(kFALSE),
156  fParton(0),
157  fPartonType(0),
158  fAncestor(0),
159  fD0D0bar(kFALSE),
160  fSelectionType(0),
161  fEvent(nullptr)
162 {
163 }
164 
170  fD(source.fD),
171  fSoftPionPt(source.fSoftPionPt),
173  fJets(source.fJets),
174  fMCLabel(source.fMCLabel),
176  fParton(source.fParton),
177  fPartonType(source.fPartonType),
178  fAncestor(source.fAncestor),
179  fD0D0bar(source.fD0D0bar),
181  fEvent(source.fEvent)
182 {
183 }
184 
189 {
190  new (this) AliDmesonJetInfo(source);
191  return *this;
192 }
193 
196 {
197  fD.SetPtEtaPhiE(0,0,0,0);
198  fSoftPionPt = 0;
199  fInvMass2Prong = 0;
200  fDmesonParticle = 0;
201  fMCLabel = -1;
202  fReconstructed = kFALSE;
203  fParton = 0;
204  fPartonType = 0;
205  fAncestor = 0;
206  fD0D0bar = kFALSE;
207  for (auto &jet : fJets) {
208  jet.second.fMomentum.SetPtEtaPhiE(0,0,0,0);
209  jet.second.fNConstituents = 0;
210  jet.second.fNEF = 0;
211  jet.second.fMaxChargedPt = 0;
212  jet.second.fMaxNeutralPt = 0;
213  }
214 }
215 
218 {
219  Printf("Printing D Meson Jet object.");
220  Printf("D Meson: pT = %.3f, eta = %.3f, phi = %.3f, inv. mass = %.3f", fD.Pt(), fD.Eta(), fD.Phi_0_2pi(), fD.M());
221  Printf("Soft pion pT: %.3f. 2-Prong Invariant mass = %.3f", fSoftPionPt, fInvMass2Prong);
222  for (auto &jet : fJets) {
223  Printf("Jet %s: pT = %.3f, eta = %.3f, phi = %.3f", jet.first.c_str(), jet.second.Pt(), jet.second.Eta(), jet.second.Phi_0_2pi());
224  Printf("Jet N Consituents = %d", jet.second.fNConstituents);
225  }
226 }
227 
232 {
233  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
234  if (it == fJets.end()) return 0;
235 
236  Double_t z = 0;
237 
238  if ((*it).second.Pt() > 0) {
239  TVector3 dvect = fD.Vect();
240  TVector3 jvect = (*it).second.fMomentum.Vect();
241 
242  Double_t jetMom = jvect * jvect;
243 
244  if (jetMom < 1e-6) {
245  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetZ", "Zero jet momentum!");
246  z = 0.999;
247  }
248  else {
249  z = (dvect * jvect) / jetMom;
250  }
251 
252  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
253  }
254 
255  return z;
256 }
257 
262 {
263  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
264  if (it == fJets.end()) return 0;
265 
266  Double_t z = 0;
267 
268  if ((*it).second.Pt() > 0) {
269  TVector3 dvect = fD.Vect();
270  TVector3 jvect = (*it).second.fMomentum.Vect();
271  // If the corr pt is < 0, assign 0.
272  Double_t corrpt = (*it).second.fCorrPt > 0 ? (*it).second.fCorrPt : 0.;
273  jvect.SetPerp(corrpt);
274 
275  Double_t jetMom = jvect * jvect;
276 
277  if (jetMom < 1e-6) {
278  z = 1.0;
279  }
280  else {
281  z = (dvect * jvect) / jetMom;
282  }
283  }
284 
285  return z;
286 }
287 
295 {
296  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
297  if (it == fJets.end()) return 0;
298 
299  return GetDistance((*it).second, deta, dphi);
300 }
301 
307 {
308  Double_t deta = 0;
309  Double_t dphi = 0;
310  return GetDistance(n, deta, dphi);
311 }
312 
320 {
321  dphi = TVector2::Phi_mpi_pi(fD.Phi() - jet.Phi());;
322  deta = fD.Eta() - jet.Eta();
323  return TMath::Sqrt(dphi*dphi + deta*deta);
324 }
325 
331 {
332  Double_t deta = 0;
333  Double_t dphi = 0;
334  return GetDistance(jet, deta, dphi);
335 }
336 
342 {
343  std::map<std::string, AliJetInfo>::const_iterator it = fJets.find(n);
344  if (it == fJets.end()) {
345  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
346  n.c_str());
347  return 0;
348  }
349  return &((*it).second);
350 }
351 
357 {
358  std::map<std::string, AliJetInfo>::iterator it = fJets.find(n);
359  if (it == fJets.end()) {
360  ::Error("AliAnalysisTaskDmesonJets::AliDmesonJetInfo::GetJet", "Could not find jet info for the jet definition '%s'!",
361  n.c_str());
362  return 0;
363  }
364  return &((*it).second);
365 }
366 
367 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoSummary
368 
372 
378  fPt(0),
379  fEta(0),
380  fPhi(0),
381  fR(0),
382  fZ(0),
383  fN(0)
384 {
385  Set(source, n);
386 }
387 
390 {
391  fPt = 0;
392  fEta = 0;
393  fPhi = 0;
394  fR = 0;
395  fZ = 0;
396  fN = 0;
397 }
398 
404 {
405  std::map<std::string, AliJetInfo>::const_iterator it = source.fJets.find(n);
406  if (it == source.fJets.end()) return;
407 
408  Set((*it).second);
409 
410  fR = source.GetDistance(n);
411  fZ = source.GetZ(n);
412 }
413 
418 {
419  fPt = source.Pt();
420  fEta = source.Eta();
421  fPhi = source.Phi_0_2pi();
422  fN = source.GetNConstituents();
423  fR = 0;
424  fZ = 0;
425 }
426 
427 // Definitions of class AliAnalysisTaskDmesonJets::AliJetInfoPbPbSummary
428 
432 
439  fCorrPt(0),
440  fCorrZ(0),
441  fArea(0)
442 {
443  Set(source, n);
444 }
445 
448 {
450  fCorrPt = 0;
451  fCorrZ = 0;
452  fArea = 0;
453 }
454 
459 {
460  AliJetInfoSummary::Set(source);
461  fArea = source.fArea;
462  fCorrPt = source.fCorrPt;
463 }
464 
470 {
471  AliJetInfoSummary::Set(source, n);
472  fCorrZ = source.GetCorrZ(n);
473 }
474 
475 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonInfoSummary
476 
480 
485  fPt(0),
486  fEta(0),
487  fPhi(0)
488 {
489  Set(source);
490 }
491 
496 {
497  fPt = source.fD.Pt();
498  fEta = source.fD.Eta();
499  fPhi = source.fD.Phi_0_2pi();
500 }
501 
504 {
505  fPt = 0;
506  fEta = 0;
507  fPhi = 0;
508 }
509 
510 // Definitions of class AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary
511 
515 
520  AliDmesonInfoSummary(source),
521  fPartonType(0),
522  fPartonPt(0),
523  fAncestorPDG(0)
524 {
525  Set(source);
526 }
527 
532 {
534 
535  fPartonType = source.fPartonType;
536 
537  if (source.fParton) {
538  fPartonPt = source.fParton->Pt();
539  }
540  else {
541  fPartonPt = 0.;
542  }
543 
544  if (source.fAncestor) {
545  fAncestorPDG = (UShort_t)((UInt_t)(TMath::Abs(source.fAncestor->GetPdgCode())));
546  }
547 }
548 
551 {
553  fPartonType = 0,
554  fPartonPt = 0.;
555  fAncestorPDG = 0;
556 }
557 
558 // Definitions of class AliAnalysisTaskDmesonJets::AliD0InfoSummary
559 
563 
568  AliDmesonInfoSummary(source),
569  fInvMass(source.fD.M()),
570  fSelectionType(0)
571 {
572 }
573 
578 {
579  fInvMass = source.fD.M();
582 }
583 
586 {
588  fSelectionType = 0;
589  fInvMass = 0;
590 }
591 
592 // Definitions of class AliAnalysisTaskDmesonJets::AliD0ExtendedInfoSummary
593 
597 
602  AliD0InfoSummary(source),
603  fDCA(0),
604  fCosThetaStar(0),
605  fd0K(0),
606  fd0Pi(0),
607  fd0d0(0),
608  fCosPointing(0),
609  fMaxNormd0(0)
610 {
611  Set(source);
612 }
613 
618 {
619  AliD0InfoSummary::Set(source);
620 
621  AliAODRecoDecayHF2Prong* recoDecay = dynamic_cast<AliAODRecoDecayHF2Prong*>(source.fDmesonParticle);
622  if (recoDecay) {
623  fDCA = recoDecay->GetDCA();
624  if (source.fSelectionType == 1) { // D0
625  fCosThetaStar = recoDecay->CosThetaStarD0();
626  fPtK = recoDecay->PtProng(0);
627  fPtPi = recoDecay->PtProng(1);
628  fd0K = recoDecay->Getd0Prong(0);
629  fd0Pi = recoDecay->Getd0Prong(1);
630  }
631  else { //D0bar
632  fCosThetaStar = recoDecay->CosThetaStarD0bar();
633  fPtK = recoDecay->PtProng(1);
634  fPtPi = recoDecay->PtProng(0);
635  fd0K = recoDecay->Getd0Prong(1);
636  fd0Pi = recoDecay->Getd0Prong(0);
637  }
638 
639  fMaxNormd0 = 0.;
640  // Based on Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent* aod)
641  // Line 480 and following
642  if (source.fEvent) {
643  for (Int_t ipr=0; ipr < 2; ipr++) {
644  Double_t diffIP = 0., errdiffIP = 0.;
645  recoDecay->Getd0MeasMinusExpProng(ipr, source.fEvent->GetMagneticField(), diffIP, errdiffIP);
646  Double_t normdd0 = 0.;
647  if (errdiffIP > 0.) {
648  normdd0 = diffIP / errdiffIP;
649  }
650  else {
651  if (diffIP == 0) {
652  normdd0 = 0;
653  }
654  else {
655  normdd0 = diffIP > 0 ? 9999. : -9999.;
656  }
657  }
658  if (TMath::Abs(normdd0) > TMath::Abs(fMaxNormd0)) {
659  fMaxNormd0 = normdd0;
660  }
661  }
662  }
663  else {
664  throw AliAnalysisTaskDmesonJets::AliEventNotFound("AliAnalysisTaskDmesonJets::AliDmesonJetInfo", "fEvent");
665  }
666 
667  fd0d0 = recoDecay->Prodd0d0();
668  fCosPointing = recoDecay->CosPointingAngle();
669  }
670 }
671 
674 {
676  fDCA = 0;
677  fCosThetaStar = 0;
678  fd0K = 0;
679  fd0Pi = 0;
680  fd0d0 = 0;
681  fCosPointing = 0;
682  fMaxNormd0 = 0;
683 }
684 
685 // Definitions of class AliAnalysisTaskDmesonJets::AliDStarInfoSummary
686 
690 
695  AliDmesonInfoSummary(source),
696  f2ProngInvMass(source.fInvMass2Prong),
697  fDeltaInvMass(source.fD.M() - source.fInvMass2Prong)
698 {
699 }
700 
705 {
707  fDeltaInvMass = source.fD.M() - source.fInvMass2Prong;
709 }
710 
713 {
715 
716  f2ProngInvMass = 0;
717  fDeltaInvMass = 0;
718 }
719 
720 // Definitions of class AliAnalysisTaskDmesonJets::OutputHandler
721 
724  fCandidateType(kD0toKpi),
725  fMCMode(kNoMC),
726  fNMassBins(0),
727  fMinMass(0),
728  fMaxMass(0),
729  fJetDefinitions(nullptr),
730  fPtBinWidth(0.5),
731  fMaxPt(100),
732  fD0Extended(kFALSE),
733  fEventInfo(nullptr),
734  fDmesonJets(nullptr),
736  fName()
737 {
738 }
739 
743  fMCMode(eng->fMCMode),
744  fNMassBins(eng->fNMassBins),
745  fMinMass(eng->fMinMass),
746  fMaxMass(eng->fMaxMass),
748  fPtBinWidth(eng->fPtBinWidth),
749  fMaxPt(eng->fMaxPt),
750  fD0Extended(eng->fD0Extended),
751  fEventInfo(&eng->fEventInfo),
752  fDmesonJets(&eng->fDmesonJets),
754  fName(eng->GetName())
755 {
756 }
757 
763 {
764  TString hname;
765 
766  TH1* histAncestor = nullptr;
767  TH1* histPrompt = nullptr;
768 
769  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
770  hname = TString::Format("%s/fHistPrompt", fName.Data());
771  histPrompt = static_cast<TH1*>(fHistManager->FindObject(hname));
772 
773  hname = TString::Format("%s/fHistAncestor", fName.Data());
774  histAncestor = static_cast<TH1*>(fHistManager->FindObject(hname));
775  }
776 
777  std::map<AliAODMCParticle*, Short_t> partons ; // set of the partons in the shower that produced each D meson
778  for (auto& dmeson_pair : *fDmesonJets) {
779  Int_t accJets = 0;
780  for (UInt_t ij = 0; ij < fJetDefinitions->size(); ij++) {
781  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions->at(ij).GetName());
782  if (!jet) continue;
783  if (applyKinCuts && !fJetDefinitions->at(ij).IsJetInAcceptance(*jet)) {
784  hname = TString::Format("%s/%s/fHistRejectedJetPt", fName.Data(), fJetDefinitions->at(ij).GetName());
785  fHistManager->FillTH1(hname, jet->Pt());
786  hname = TString::Format("%s/%s/fHistRejectedJetPhi", fName.Data(), fJetDefinitions->at(ij).GetName());
787  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
788  hname = TString::Format("%s/%s/fHistRejectedJetEta", fName.Data(), fJetDefinitions->at(ij).GetName());
789  fHistManager->FillTH1(hname, jet->Eta());
790  continue;
791  }
792  accJets++;
793  }
794  if (accJets > 0) {
795  if (histPrompt) {
796  if (dmeson_pair.second.fParton) {
797  partons[dmeson_pair.second.fParton] = dmeson_pair.second.fPartonType;
798  UInt_t absPdgParton = TMath::Abs(dmeson_pair.second.fParton->GetPdgCode());
799  if (absPdgParton == 4) {
800  histPrompt->Fill("Prompt", 1);
801  }
802  else if (absPdgParton == 5) {
803  histPrompt->Fill("Non-Prompt", 1);
804  }
805  else {
806  histPrompt->Fill("Unknown", 1);
807  }
808  }
809  else {
810  histPrompt->Fill("Unknown", 1);
811  }
812  }
813 
814  if (histAncestor) {
815  if (dmeson_pair.second.fAncestor) {
816  UInt_t absPdgAncestor = TMath::Abs(dmeson_pair.second.fAncestor->GetPdgCode());
817  if (absPdgAncestor == 4) {
818  histAncestor->Fill("Charm", 1);
819  }
820  else if (absPdgAncestor == 5) {
821  histAncestor->Fill("Bottom", 1);
822  }
823  else if (absPdgAncestor == 2212) {
824  histAncestor->Fill("Proton", 1);
825  }
826  else {
827  histAncestor->Fill("Unknown", 1);
828  }
829  }
830  else {
831  histAncestor->Fill("Unknown", 1);
832  }
833  }
834  }
835  else {
836  hname = TString::Format("%s/fHistRejectedDMesonPt", fName.Data());
837  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
838  hname = TString::Format("%s/fHistRejectedDMesonPhi", fName.Data());
839  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
840  hname = TString::Format("%s/fHistRejectedDMesonEta", fName.Data());
841  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
842  if (fMCMode != kMCTruth) {
844  hname = TString::Format("%s/fHistRejectedDMesonInvMass", fName.Data());
845  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
846  }
847  else if (fCandidateType == kDstartoKpipi) {
848  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", fName.Data());
849  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
850 
851  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", fName.Data());
852  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
853  }
854  }
855  }
856  }
857 
858  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
859  hname = TString::Format("%s/fHistPartonPt", fName.Data());
860  TH1* histPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
861  hname = TString::Format("%s/fHistPartonEta", fName.Data());
862  TH1* histPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
863  hname = TString::Format("%s/fHistPartonPhi", fName.Data());
864  TH1* histPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
865  hname = TString::Format("%s/fHistPartonType", fName.Data());
866  TH1* histPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
867 
868  for (auto parton : partons) {
869  if (!parton.first) continue;
870  histPartonPt->Fill(parton.first->Pt());
871  histPartonEta->Fill(parton.first->Eta());
872  histPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
873  histPartonType->Fill(parton.second);
874  }
875  }
876 
877  return kTRUE;
878 }
879 
880 // Definitions of class AliAnalysisTaskDmesonJets::OutputHandlerTHnSparse
881 
884  OutputHandler(),
885  fEnabledAxis(0)
886 {
887 }
888 
891  OutputHandler(eng),
892  fEnabledAxis(0)
893 {
894 }
895 
896 
901 {
902  TString hname;
903 
904  Int_t nPtBins = TMath::CeilNint(fMaxPt / fPtBinWidth);
905 
906  for (auto &jetDef : *fJetDefinitions) {
907 
908  AliDebugGeneralStream("AliAnalysisTaskDmesonJets::OutputHandlerTHnSparse::BuildOutputObject", 2) << "Now working on '" << jetDef.GetName() << "'" << std::endl;
909 
910  Double_t radius = jetDef.fRadius;
911 
912  TString title[30] = {""};
913  Int_t nbins[30] = {0 };
914  Double_t min [30] = {0.};
915  Double_t max [30] = {0.};
916  Int_t dim = 0 ;
917 
918  title[dim] = "#it{p}_{T,D} (GeV/#it{c})";
919  nbins[dim] = nPtBins;
920  min[dim] = 0;
921  max[dim] = fMaxPt;
922  dim++;
923 
924  if ((fEnabledAxis & kPositionD) != 0) {
925  title[dim] = "#eta_{D}";
926  nbins[dim] = 50;
927  min[dim] = -1;
928  max[dim] = 1;
929  dim++;
930 
931  title[dim] = "#phi_{D} (rad)";
932  nbins[dim] = 150;
933  min[dim] = 0;
934  max[dim] = TMath::TwoPi();
935  dim++;
936  }
937 
938  if ((fEnabledAxis & kInvMass) != 0 && fCandidateType == kDstartoKpipi) {
939  title[dim] = "#it{M}_{K#pi#pi} (GeV/#it{c}^{2})";
940  nbins[dim] = fNMassBins;
941  min[dim] = fMinMass;
942  max[dim] = fMaxMass;
943  dim++;
944  }
945 
947  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
948  nbins[dim] = fNMassBins;
949  min[dim] = fMinMass;
950  max[dim] = fMaxMass;
951  dim++;
952  }
953 
955  title[dim] = "#it{M}_{K#pi} (GeV/#it{c}^{2})";
956  nbins[dim] = fNMassBins;
957  CalculateMassLimits(fMaxMass - fMinMass, 421, fNMassBins, min[dim], max[dim]);
958  dim++;
959  }
960 
961  if (fCandidateType == kDstartoKpipi) {
962  title[dim] = "#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})";
963  nbins[dim] = fNMassBins*6;
964  CalculateMassLimits(0.20, 413, nbins[dim], min[dim], max[dim]);
965 
966  // subtract mass of D0
967  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
968  min[dim] -= D0mass;
969  max[dim] -= D0mass;
970 
971  dim++;
972  }
973 
975  title[dim] = "#it{p}_{T,#pi} (GeV/#it{c})";
976  nbins[dim] = 100;
977  min[dim] = 0;
978  max[dim] = 25;
979  dim++;
980  }
981 
982  title[dim] = "#it{z}_{D}";
983  nbins[dim] = 110;
984  min[dim] = 0;
985  max[dim] = 1.10;
986  dim++;
987 
988  if ((fEnabledAxis & kDeltaR) != 0) {
989  title[dim] = "#Delta R_{D-jet}";
990  nbins[dim] = 100;
991  min[dim] = 0;
992  max[dim] = radius * 1.5;
993  dim++;
994  }
995 
996  if ((fEnabledAxis & kDeltaEta) != 0) {
997  title[dim] = "#eta_{D} - #eta_{jet}";
998  nbins[dim] = 100;
999  min[dim] = -radius * 1.2;
1000  max[dim] = radius * 1.2;
1001  dim++;
1002  }
1003 
1004  if ((fEnabledAxis & kDeltaPhi) != 0) {
1005  title[dim] = "#phi_{D} - #phi_{jet} (rad)";
1006  nbins[dim] = 100;
1007  min[dim] = -radius * 1.2;
1008  max[dim] = radius * 1.2;
1009  dim++;
1010  }
1011 
1012  title[dim] = "#it{p}_{T,jet} (GeV/#it{c})";
1013  nbins[dim] = nPtBins;
1014  min[dim] = 0;
1015  max[dim] = fMaxPt;
1016  dim++;
1017 
1018  if ((fEnabledAxis & kPositionJet) != 0) {
1019  title[dim] = "#eta_{jet}";
1020  nbins[dim] = 50;
1021  min[dim] = -1;
1022  max[dim] = 1;
1023  dim++;
1024 
1025  title[dim] = "#phi_{jet} (rad)";
1026  nbins[dim] = 150;
1027  min[dim] = 0;
1028  max[dim] = TMath::TwoPi();
1029  dim++;
1030  }
1031 
1032  if ((fEnabledAxis & kJetConstituents) != 0) {
1033  title[dim] = "No. of constituents";
1034  nbins[dim] = 50;
1035  min[dim] = -0.5;
1036  max[dim] = 49.5;
1037  dim++;
1038  }
1039 
1040  hname = TString::Format("%s/%s/fDmesonJets", fName.Data(), jetDef.GetName());
1041  THnSparse* h = fHistManager->CreateTHnSparse(hname,hname,dim,nbins,min,max);
1042  for (Int_t j = 0; j < dim; j++) {
1043  h->GetAxis(j)->SetTitle(title[j]);
1044  }
1045  }
1046 }
1047 
1052 {
1053  for (UInt_t i = 0; i < fJetDefinitions->size(); i++) {
1054  if (fJetDefinitions->at(i).IsJetInAcceptance(dMesonJet, fJetDefinitions->at(i).GetName())) return kTRUE;
1055  }
1056 
1057  return kFALSE;
1058 }
1059 
1064 {
1065  TString hname;
1066 
1067  for (auto& dmeson_pair : *fDmesonJets) {
1068  if (!IsAnyJetInAcceptance(dmeson_pair.second)) {
1069  hname = TString::Format("%s/fHistRejectedDMesonPt", fName.Data());
1070  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1071  hname = TString::Format("%s/fHistRejectedDMesonPhi", fName.Data());
1072  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1073  hname = TString::Format("%s/fHistRejectedDMesonEta", fName.Data());
1074  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1075  }
1076  }
1077 
1078  for (auto &jetDef : *fJetDefinitions) {
1079 
1080  hname = TString::Format("%s/%s/fDmesonJets", fName.Data(), jetDef.GetName());
1081  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
1082 
1083  for (auto& dmeson_pair : *fDmesonJets) {
1084  const AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
1085  if (!jet) continue;
1086  if (!jetDef.IsJetInAcceptance(*jet)) {
1087  hname = TString::Format("%s/%s/fHistRejectedJetPt", fName.Data(), jetDef.GetName());
1088  fHistManager->FillTH1(hname, jet->Pt());
1089  hname = TString::Format("%s/%s/fHistRejectedJetPhi", fName.Data(), jetDef.GetName());
1090  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1091  hname = TString::Format("%s/%s/fHistRejectedJetEta", fName.Data(), jetDef.GetName());
1092  fHistManager->FillTH1(hname, jet->Eta());
1093  continue;
1094  }
1095  FillHnSparse(h, dmeson_pair.second, jetDef.GetName());
1096  }
1097  }
1098 
1099  return kTRUE;
1100 }
1101 
1108 {
1109  // Fill the THnSparse histogram.
1110 
1111  Double_t contents[30] = {0.};
1112 
1113  Double_t z = DmesonJet.GetZ(n);
1114  Double_t deltaPhi = 0;
1115  Double_t deltaEta = 0;
1116  Double_t deltaR = DmesonJet.GetDistance(n, deltaEta, deltaPhi);
1117 
1118  std::map<std::string, AliJetInfo>::const_iterator it = DmesonJet.fJets.find(n);
1119  if (it == DmesonJet.fJets.end()) return kFALSE;
1120 
1121  for (Int_t i = 0; i < h->GetNdimensions(); i++) {
1122  TString title(h->GetAxis(i)->GetTitle());
1123  if (title=="#it{p}_{T,D} (GeV/#it{c})") contents[i] = DmesonJet.fD.Pt();
1124  else if (title=="#eta_{D}") contents[i] = DmesonJet.fD.Eta();
1125  else if (title=="#phi_{D} (rad)") contents[i] = DmesonJet.fD.Phi_0_2pi();
1126  else if (title=="#it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fInvMass2Prong > 0 ? DmesonJet.fInvMass2Prong : DmesonJet.fD.M();
1127  else if (title=="#it{M}_{K#pi#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M();
1128  else if (title=="#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2})") contents[i] = DmesonJet.fD.M() - DmesonJet.fInvMass2Prong;
1129  else if (title=="#it{p}_{T,#pi} (GeV/#it{c})") contents[i] = DmesonJet.fSoftPionPt;
1130  else if (title=="#it{z}_{D}") contents[i] = z;
1131  else if (title=="#Delta R_{D-jet}") contents[i] = deltaR;
1132  else if (title=="#eta_{D} - #eta_{jet}") contents[i] = deltaEta;
1133  else if (title=="#phi_{D} - #phi_{jet} (rad)") contents[i] = deltaPhi;
1134  else if (title=="#it{p}_{T,jet} (GeV/#it{c})") contents[i] = (*it).second.Pt();
1135  else if (title=="#eta_{jet}") contents[i] = (*it).second.Eta();
1136  else if (title=="#phi_{jet} (rad)") contents[i] = (*it).second.Phi_0_2pi();
1137  else if (title=="No. of constituents") contents[i] = (*it).second.fNConstituents;
1138  else {
1139  AliWarningGeneralStream("AliAnalysisTaskDmesonJets::OutputHandlerTHnSparse::FillHnSparse") << "Unable to fill dimension '" << title.Data() << "'!" << std::endl;
1140  }
1141  }
1142 
1143  h->Fill(contents);
1144 
1145  return kTRUE;
1146 }
1147 
1148 // Definitions of class AliAnalysisTaskDmesonJets::OutputHandlerTTree
1149 
1152  OutputHandler(),
1153  fDataSlotNumber(-1),
1154  fTree(0),
1155  fCurrentDmesonJetInfo(0),
1156  fCurrentJetInfo(0)
1157 {
1158 }
1159 
1162  OutputHandler(eng),
1163  fDataSlotNumber(-1),
1164  fTree(0),
1166  fCurrentJetInfo(0)
1167 {
1168 }
1169 
1174 {
1175  TString classname;
1176  if (fMCMode == kMCTruth) {
1177  classname = "AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary";
1179  }
1180  else {
1181  switch (fCandidateType) {
1182  case kD0toKpi:
1183  case kD0toKpiLikeSign:
1184  if (fD0Extended) {
1185  classname = "AliAnalysisTaskDmesonJets::AliD0ExtendedInfoSummary";
1187  }
1188  else {
1189  classname = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1191  }
1192  break;
1193  case kDstartoKpipi:
1194  classname = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1196  break;
1197  }
1198  }
1199  TString treeName = TString::Format("%s_%s", taskName, fName.Data());
1200  fTree = new TTree(treeName, treeName);
1201  fTree->Branch("DmesonJet", classname, &fCurrentDmesonJetInfo);
1203  for (Int_t i = 0; i < fJetDefinitions->size(); i++) {
1204  if (fJetDefinitions->at(i).fRhoName.IsNull()) {
1206  fTree->Branch(fJetDefinitions->at(i).GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoSummary", &fCurrentJetInfo[i]);
1207  }
1208  else {
1210  fTree->Branch(fJetDefinitions->at(i).GetName(), "AliAnalysisTaskDmesonJets::AliJetInfoPbPbSummary", &fCurrentJetInfo[i]);
1211  }
1212  }
1213 }
1214 
1219 {
1220  TString hname;
1221 
1222  std::map<AliAODMCParticle*, Short_t> partons ;
1223 
1224  TH1* histAncestor = nullptr;
1225  TH1* histPrompt = nullptr;
1226 
1227  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
1228  hname = TString::Format("%s/fHistPrompt", fName.Data());
1229  histPrompt = static_cast<TH1*>(fHistManager->FindObject(hname));
1230 
1231  hname = TString::Format("%s/fHistAncestor", fName.Data());
1232  histAncestor = static_cast<TH1*>(fHistManager->FindObject(hname));
1233  }
1234 
1235  for (auto& dmeson_pair : *fDmesonJets) {
1236  fCurrentDmesonJetInfo->Set(dmeson_pair.second);
1237  Int_t accJets = 0;
1238  for (UInt_t ij = 0; ij < fJetDefinitions->size(); ij++) {
1239  fCurrentJetInfo[ij]->Reset();
1240  AliJetInfo* jet = dmeson_pair.second.GetJet(fJetDefinitions->at(ij).GetName());
1241  if (!jet) continue;
1242  if (applyKinCuts && !fJetDefinitions->at(ij).IsJetInAcceptance(*jet)) {
1243  hname = TString::Format("%s/%s/fHistRejectedJetPt", fName.Data(), fJetDefinitions->at(ij).GetName());
1244  fHistManager->FillTH1(hname, jet->Pt());
1245  hname = TString::Format("%s/%s/fHistRejectedJetPhi", fName.Data(), fJetDefinitions->at(ij).GetName());
1246  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1247  hname = TString::Format("%s/%s/fHistRejectedJetEta", fName.Data(), fJetDefinitions->at(ij).GetName());
1248  fHistManager->FillTH1(hname, jet->Eta());
1249  continue;
1250  }
1251  fCurrentJetInfo[ij]->Set(dmeson_pair.second, fJetDefinitions->at(ij).GetName());
1252  accJets++;
1253  }
1254  if (accJets > 0) {
1255  if (histPrompt) {
1256  if (dmeson_pair.second.fParton) {
1257  partons[dmeson_pair.second.fParton] = dmeson_pair.second.fPartonType;
1258  UInt_t absPdgParton = TMath::Abs(dmeson_pair.second.fParton->GetPdgCode());
1259  if (absPdgParton == 4) {
1260  histPrompt->Fill("Prompt", 1);
1261  }
1262  else if (absPdgParton == 5) {
1263  histPrompt->Fill("Non-Prompt", 1);
1264  }
1265  else {
1266  histPrompt->Fill("Unknown", 1);
1267  }
1268  }
1269  else {
1270  histPrompt->Fill("Unknown", 1);
1271  }
1272  }
1273 
1274  if (histAncestor) {
1275  if (dmeson_pair.second.fAncestor) {
1276  UInt_t absPdgAncestor = TMath::Abs(dmeson_pair.second.fAncestor->GetPdgCode());
1277  if (absPdgAncestor == 4) {
1278  histAncestor->Fill("Charm", 1);
1279  }
1280  else if (absPdgAncestor == 5) {
1281  histAncestor->Fill("Bottom", 1);
1282  }
1283  else if (absPdgAncestor == 2212) {
1284  histAncestor->Fill("Proton", 1);
1285  }
1286  else {
1287  histAncestor->Fill("Unknown", 1);
1288  }
1289  }
1290  else {
1291  histAncestor->Fill("Unknown", 1);
1292  }
1293  }
1294 
1295  fTree->Fill();
1296  }
1297  else {
1298  hname = TString::Format("%s/fHistRejectedDMesonPt", fName.Data());
1299  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1300  hname = TString::Format("%s/fHistRejectedDMesonPhi", fName.Data());
1301  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1302  hname = TString::Format("%s/fHistRejectedDMesonEta", fName.Data());
1303  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1304  if (fMCMode != kMCTruth) {
1306  hname = TString::Format("%s/fHistRejectedDMesonInvMass", fName.Data());
1307  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1308  }
1309  else if (fCandidateType == kDstartoKpipi) {
1310  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", fName.Data());
1311  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1312 
1313  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", fName.Data());
1314  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1315  }
1316  }
1317  }
1318  }
1319 
1320  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
1321  hname = TString::Format("%s/fHistPartonPt", fName.Data());
1322  TH1* histPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
1323  hname = TString::Format("%s/fHistPartonEta", fName.Data());
1324  TH1* histPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
1325  hname = TString::Format("%s/fHistPartonPhi", fName.Data());
1326  TH1* histPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
1327  hname = TString::Format("%s/fHistPartonType", fName.Data());
1328  TH1* histPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
1329 
1330  for (auto parton : partons) {
1331  if (!parton.first) continue;
1332  histPartonPt->Fill(parton.first->Pt());
1333  histPartonEta->Fill(parton.first->Eta());
1334  histPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
1335  histPartonType->Fill(parton.second);
1336  }
1337  }
1338 
1339  return kTRUE;
1340 }
1341 
1342 // Definitions of class AliAnalysisTaskDmesonJets::OutputHandlerTTreeExtended
1343 
1346  OutputHandler(),
1347  fDataSlotNumber(-1),
1348  fTree(0),
1349  fEventClassName(),
1350  fDMesonClassName(),
1351  fJetClassName()
1352 {
1353 }
1354 
1357  OutputHandler(eng),
1358  fDataSlotNumber(-1),
1359  fTree(0),
1360  fEventClassName(),
1361  fDMesonClassName(),
1362  fJetClassName()
1363 {
1364 }
1365 
1370 {
1371  TString event_class_name = "AliAnalysisTaskDmesonJets::AliEventInfoSummary";
1372  TString d_meson_class_name;
1373 
1374  TString jet_class_name = "std::vector<AliAnalysisTaskDmesonJets::AliJetInfoSummary>";
1375  Bool_t RhoJet = kFALSE;
1376  for (auto jetDef : eng->GetJetDefinitions()) {
1377  if (!jetDef.fRhoName.IsNull()) {
1378  RhoJet = kTRUE;
1379  jet_class_name = "std::vector<AliAnalysisTaskDmesonJets::AliJetInfoPbPbSummary>";
1380  }
1381  }
1382 
1383  OutputHandlerTTreeExtendedBase* result = nullptr;
1384  if (eng->GetMCMode() == kMCTruth) {
1385  d_meson_class_name = "std::vector<AliAnalysisTaskDmesonJets::AliDmesonMCInfoSummary>";
1386  if (RhoJet) {
1388  }
1389  else {
1391  }
1392  }
1393  else {
1394  switch (eng->GetCandidateType()) {
1395  case kD0toKpi:
1396  case kD0toKpiLikeSign:
1397  if (eng->IsD0Extended()) {
1398  d_meson_class_name = "std::vector<AliAnalysisTaskDmesonJets::AliD0ExtendedInfoSummary>";
1399  if (RhoJet) {
1401  }
1402  else {
1404  }
1405  }
1406  else {
1407  d_meson_class_name = "AliAnalysisTaskDmesonJets::AliD0InfoSummary";
1408  if (RhoJet) {
1410  }
1411  else {
1413  }
1414  }
1415  break;
1416  case kDstartoKpipi:
1417  d_meson_class_name = "AliAnalysisTaskDmesonJets::AliDStarInfoSummary";
1418  if (RhoJet) {
1420  }
1421  else {
1423  }
1424  break;
1425  }
1426  }
1427 
1428  result->fEventClassName = event_class_name;
1429  result->fDMesonClassName = d_meson_class_name;
1430  result->fJetClassName = jet_class_name;
1431 
1432  return result;
1433 }
1434 
1435 // Definitions of class AliAnalysisTaskDmesonJets::OutputHandlerTTree
1436 
1438 template<class EVENTTYPE, class DMESONTYPE, class JETTYPE>
1441  fCurrentEventInfo(),
1442  fCurrentDmesonInfo(),
1443  fCurrentJetInfo()
1444 {
1445 }
1446 
1448 template<class EVENTTYPE, class DMESONTYPE, class JETTYPE>
1453  fCurrentJetInfo()
1454 {
1455 }
1456 
1460 template<class EVENTTYPE, class DMESONTYPE, class JETTYPE>
1462 {
1463  TString treeName = TString::Format("%s_%s", taskName, fName.Data());
1464  fTree = new TTree(treeName, treeName);
1465  fTree->Branch("Event", fEventClassName, &fCurrentEventInfo);
1466  fTree->Branch("Dmesons", fDMesonClassName, &fCurrentDmesonInfo);
1467 
1468  for (auto jetDef : *fJetDefinitions) {
1469  fCurrentJetInfo[jetDef.GetName()] = std::vector<JETTYPE>();
1470  fTree->Branch(jetDef.GetName(), fJetClassName, &fCurrentJetInfo[jetDef.GetName()]);
1471  }
1472 }
1473 
1477 template<class EVENTTYPE, class DMESONTYPE, class JETTYPE>
1479 {
1480  if (fDmesonJets->empty()) return kFALSE;
1481 
1482  TString hname;
1483 
1484  std::map<AliAODMCParticle*, Short_t> partons ;
1485 
1486  TH1* histAncestor = nullptr;
1487  TH1* histPrompt = nullptr;
1488 
1489  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
1490  hname = TString::Format("%s/fHistPrompt", fName.Data());
1491  histPrompt = static_cast<TH1*>(fHistManager->FindObject(hname));
1492 
1493  hname = TString::Format("%s/fHistAncestor", fName.Data());
1494  histAncestor = static_cast<TH1*>(fHistManager->FindObject(hname));
1495  }
1496 
1498 
1499  fCurrentDmesonInfo.clear();
1500  for (auto& jetInfo : fCurrentJetInfo) {
1501  jetInfo.second.clear();
1502  }
1503 
1504  for (auto& dmeson_pair : *fDmesonJets) {
1505  DMESONTYPE dmeson_tree;
1506  dmeson_tree.Set(dmeson_pair.second);
1507  Int_t accJets = 0;
1508  for (auto jetDef : *fJetDefinitions) {
1509  JETTYPE jet_tree;
1510  AliJetInfo* jet = dmeson_pair.second.GetJet(jetDef.GetName());
1511  if (jet) {
1512  if (applyKinCuts && !jetDef.IsJetInAcceptance(*jet)) {
1513  hname = TString::Format("%s/%s/fHistRejectedJetPt", fName.Data(), jetDef.GetName());
1514  fHistManager->FillTH1(hname, jet->Pt());
1515  hname = TString::Format("%s/%s/fHistRejectedJetPhi", fName.Data(), jetDef.GetName());
1516  fHistManager->FillTH1(hname, jet->Phi_0_2pi());
1517  hname = TString::Format("%s/%s/fHistRejectedJetEta", fName.Data(), jetDef.GetName());
1518  fHistManager->FillTH1(hname, jet->Eta());
1519  }
1520  else {
1521  jet_tree.Set(dmeson_pair.second, jetDef.GetName());
1522  accJets++;
1523  }
1524  }
1525  fCurrentJetInfo[jetDef.GetName()].push_back(jet_tree);
1526  }
1527  if (accJets > 0) {
1528  if (histPrompt) {
1529  if (dmeson_pair.second.fParton) {
1530  partons[dmeson_pair.second.fParton] = dmeson_pair.second.fPartonType;
1531  UInt_t absPdgParton = TMath::Abs(dmeson_pair.second.fParton->GetPdgCode());
1532  if (absPdgParton == 4) {
1533  histPrompt->Fill("Prompt", 1);
1534  }
1535  else if (absPdgParton == 5) {
1536  histPrompt->Fill("Non-Prompt", 1);
1537  }
1538  else {
1539  histPrompt->Fill("Unknown", 1);
1540  }
1541  }
1542  else {
1543  histPrompt->Fill("Unknown", 1);
1544  }
1545  }
1546 
1547  if (histAncestor) {
1548  if (dmeson_pair.second.fAncestor) {
1549  UInt_t absPdgAncestor = TMath::Abs(dmeson_pair.second.fAncestor->GetPdgCode());
1550  if (absPdgAncestor == 4) {
1551  histAncestor->Fill("Charm", 1);
1552  }
1553  else if (absPdgAncestor == 5) {
1554  histAncestor->Fill("Bottom", 1);
1555  }
1556  else if (absPdgAncestor == 2212) {
1557  histAncestor->Fill("Proton", 1);
1558  }
1559  else {
1560  histAncestor->Fill("Unknown", 1);
1561  }
1562  }
1563  else {
1564  histAncestor->Fill("Unknown", 1);
1565  }
1566  }
1567 
1568  fCurrentDmesonInfo.push_back(dmeson_tree);
1569  }
1570  else {
1571  hname = TString::Format("%s/fHistRejectedDMesonPt", fName.Data());
1572  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Pt());
1573  hname = TString::Format("%s/fHistRejectedDMesonPhi", fName.Data());
1574  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Phi_0_2pi());
1575  hname = TString::Format("%s/fHistRejectedDMesonEta", fName.Data());
1576  fHistManager->FillTH1(hname, dmeson_pair.second.fD.Eta());
1577  if (fMCMode != kMCTruth) {
1579  hname = TString::Format("%s/fHistRejectedDMesonInvMass", fName.Data());
1580  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M());
1581  }
1582  else if (fCandidateType == kDstartoKpipi) {
1583  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", fName.Data());
1584  fHistManager->FillTH1(hname, dmeson_pair.second.fInvMass2Prong);
1585 
1586  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", fName.Data());
1587  fHistManager->FillTH1(hname, dmeson_pair.second.fD.M() - dmeson_pair.second.fInvMass2Prong);
1588  }
1589  }
1590  for (auto& jetInfo : fCurrentJetInfo) {
1591  jetInfo.second.pop_back();
1592  }
1593  }
1594  }
1595 
1596  if (!fCurrentDmesonInfo.empty()) fTree->Fill();
1597 
1598  if (fMCMode == kSignalOnly || fMCMode == kMCTruth) {
1599  hname = TString::Format("%s/fHistPartonPt", fName.Data());
1600  TH1* histPartonPt = static_cast<TH1*>(fHistManager->FindObject(hname));
1601  hname = TString::Format("%s/fHistPartonEta", fName.Data());
1602  TH1* histPartonEta = static_cast<TH1*>(fHistManager->FindObject(hname));
1603  hname = TString::Format("%s/fHistPartonPhi", fName.Data());
1604  TH1* histPartonPhi = static_cast<TH1*>(fHistManager->FindObject(hname));
1605  hname = TString::Format("%s/fHistPartonType", fName.Data());
1606  TH1* histPartonType = static_cast<TH1*>(fHistManager->FindObject(hname));
1607 
1608  for (auto parton : partons) {
1609  if (!parton.first) continue;
1610  histPartonPt->Fill(parton.first->Pt());
1611  histPartonEta->Fill(parton.first->Eta());
1612  histPartonPhi->Fill(TVector2::Phi_0_2pi(parton.first->Phi()));
1613  histPartonType->Fill(parton.second);
1614  }
1615  }
1616 
1617  return kTRUE;
1618 }
1619 
1620 // Definitions of class AliAnalysisTaskDmesonJets::AliJetDefinition
1621 
1625 
1628  TObject(),
1629  fJetType(AliJetContainer::kChargedJet),
1630  fRadius(0),
1631  fJetAlgo(AliJetContainer::antikt_algorithm),
1632  fRecoScheme(AliJetContainer::pt_scheme),
1633  fMinJetPt(0.),
1634  fMaxJetPt(500.),
1635  fMinJetPhi(0.),
1636  fMaxJetPhi(0.),
1637  fMinJetEta(-1.),
1638  fMaxJetEta(1.),
1639  fMinChargedPt(0.),
1640  fMaxChargedPt(0.),
1641  fMinNeutralPt(0.),
1642  fMaxNeutralPt(0.),
1643  fRhoName(),
1644  fRho(0)
1645 {
1646 }
1647 
1655  TObject(),
1656  fJetType(type),
1657  fRadius(r),
1658  fJetAlgo(algo),
1659  fRecoScheme(reco),
1660  fMinJetPt(0.),
1661  fMaxJetPt(500.),
1662  fMinJetPhi(0.),
1663  fMaxJetPhi(0.),
1664  fMinJetEta(-1.),
1665  fMaxJetEta(1.),
1666  fMinChargedPt(0.),
1667  fMaxChargedPt(0.),
1668  fMinNeutralPt(0.),
1669  fMaxNeutralPt(0.),
1670  fRhoName(),
1671  fRho(0)
1672 {
1673 }
1674 
1682  TObject(),
1683  fJetType(type),
1684  fRadius(r),
1685  fJetAlgo(algo),
1686  fRecoScheme(reco),
1687  fMinJetPt(0.),
1688  fMaxJetPt(0.),
1689  fMinJetPhi(0.),
1690  fMaxJetPhi(0.),
1691  fMinJetEta(0.),
1692  fMaxJetEta(0.),
1693  fMinChargedPt(0.),
1694  fMaxChargedPt(0.),
1695  fMinNeutralPt(0.),
1696  fMaxNeutralPt(0.),
1697  fRhoName(rhoName),
1698  fRho(0)
1699 {
1700 }
1701 
1706  TObject(),
1707  fJetType(source.fJetType),
1708  fRadius(source.fRadius),
1709  fJetAlgo(source.fJetAlgo),
1710  fRecoScheme(source.fRecoScheme),
1711  fMinJetPt(source.fMinJetPt),
1712  fMaxJetPt(source.fMaxJetPt),
1713  fMinJetPhi(source.fMinJetPhi),
1714  fMaxJetPhi(source.fMaxJetPhi),
1715  fMinJetEta(source.fMinJetEta),
1716  fMaxJetEta(source.fMaxJetEta),
1717  fMinChargedPt(source.fMinChargedPt),
1718  fMaxChargedPt(source.fMaxChargedPt),
1719  fMinNeutralPt(source.fMinNeutralPt),
1720  fMaxNeutralPt(source.fMaxNeutralPt),
1721  fRhoName(source.fRhoName),
1722  fRho(0)
1723 {
1724 }
1725 
1730 {
1731  new (this) AliHFJetDefinition(source);
1732  return *this;
1733 }
1734 
1737 {
1738  static TString name;
1739 
1741 
1742  return name.Data();
1743 }
1744 
1750 {
1751  if (fMinJetEta < fMaxJetEta && (jet.Eta() < fMinJetEta || jet.Eta() > fMaxJetEta)) return kFALSE;
1752  if (fMinJetPhi < fMaxJetPhi && (jet.Phi() < fMinJetPhi || jet.Phi() > fMaxJetPhi)) return kFALSE;
1753  if (fMinJetPt < fMaxJetPt && (jet.Pt() > fMaxJetPt || jet.Pt() < fMinJetPt)) return kFALSE;
1754  if (fMinChargedPt < fMaxChargedPt && (jet.fMaxChargedPt < fMinChargedPt || jet.fMaxChargedPt > fMaxChargedPt)) return kFALSE;
1755  if (fMinNeutralPt < fMaxNeutralPt && (jet.fMaxNeutralPt < fMinNeutralPt || jet.fMaxNeutralPt > fMaxNeutralPt)) return kFALSE;
1756 
1757  return kTRUE;
1758 }
1759 
1765 {
1766  const AliJetInfo* jet = dMesonJet.GetJet(n);
1767  if (!jet) return kFALSE;
1768  return IsJetInAcceptance((*jet));
1769 }
1770 
1777 {
1778  if (lhs.fJetType > rhs.fJetType) return false;
1779  else if (lhs.fJetType < rhs.fJetType) return true;
1780  else {
1781  if (lhs.fRadius > rhs.fRadius) return false;
1782  else if (lhs.fRadius < rhs.fRadius) return true;
1783  else {
1784  if (lhs.fJetAlgo > rhs.fJetAlgo) return false;
1785  else if (lhs.fJetAlgo < rhs.fJetAlgo) return true;
1786  else {
1787  if (lhs.fRecoScheme < rhs.fRecoScheme) return true;
1788  else return false;
1789  }
1790  }
1791  }
1792 }
1793 
1800 {
1801  if (lhs.fJetType != rhs.fJetType) return false;
1802  if (lhs.fRadius != rhs.fRadius) return false;
1803  if (lhs.fJetAlgo != rhs.fJetAlgo) return false;
1804  if (lhs.fRecoScheme != rhs.fRecoScheme) return false;
1805  return true;
1806 }
1807 
1808 // Definitions of class AliAnalysisTaskDmesonJets::AnalysisEngine
1809 
1813 
1816  TObject(),
1817  fCandidateType(kD0toKpi),
1818  fCandidateName(),
1819  fCandidatePDG(0),
1820  fNDaughters(0),
1821  fPDGdaughters(),
1822  fBranchName(),
1823  fMCMode(kNoMC),
1824  fNMassBins(0),
1825  fMinMass(0),
1826  fMaxMass(0),
1827  fRDHFCuts(0),
1828  fRejectedOrigin(0),
1829  fAcceptedDecay(0),
1830  fInhibit(kFALSE),
1831  fJetDefinitions(),
1832  fPtBinWidth(0.5),
1833  fMaxPt(100),
1834  fD0Extended(kFALSE),
1835  fOutputHandler(nullptr),
1836  fRandomGen(0),
1837  fTrackEfficiency(0),
1838  fRejectISR(kFALSE),
1839  fDmesonJets(),
1840  fCandidateArray(0),
1841  fMCContainer(),
1842  fTrackContainers(),
1843  fClusterContainers(),
1844  fAodEvent(0),
1845  fFastJetWrapper(0),
1846  fHistManager(0),
1847  fEventInfo(),
1848  fName()
1849 {
1850 }
1851 
1860  TObject(),
1861  fCandidateType(type),
1862  fCandidateName(),
1863  fCandidatePDG(0),
1864  fNDaughters(0),
1865  fPDGdaughters(),
1866  fBranchName(),
1867  fMCMode(MCmode),
1868  fNMassBins(nMassBins),
1869  fMinMass(0),
1870  fMaxMass(0),
1871  fRDHFCuts(cuts),
1872  fRejectedOrigin(0),
1874  fInhibit(kFALSE),
1875  fJetDefinitions(),
1876  fPtBinWidth(0.5),
1877  fMaxPt(100),
1878  fD0Extended(kFALSE),
1880  fRandomGen(0),
1881  fTrackEfficiency(0),
1882  fDmesonJets(),
1883  fCandidateArray(0),
1884  fMCContainer(),
1885  fTrackContainers(),
1887  fAodEvent(0),
1888  fFastJetWrapper(0),
1889  fHistManager(0),
1890  fEventInfo(),
1891  fName()
1892 {
1893  SetCandidateProperties(range);
1894 }
1895 
1900  TObject(source),
1903  fCandidatePDG(source.fCandidatePDG),
1904  fNDaughters(source.fNDaughters),
1905  fPDGdaughters(source.fPDGdaughters),
1906  fBranchName(source.fBranchName),
1907  fMCMode(source.fMCMode),
1908  fNMassBins(source.fNMassBins),
1909  fMinMass(source.fMinMass),
1910  fMaxMass(source.fMaxMass),
1911  fRDHFCuts(),
1914  fInhibit(source.fInhibit),
1916  fPtBinWidth(source.fPtBinWidth),
1917  fMaxPt(source.fMaxPt),
1918  fD0Extended(source.fD0Extended),
1919  fRandomGen(source.fRandomGen),
1921  fDmesonJets(),
1923  fMCContainer(source.fMCContainer),
1926  fAodEvent(source.fAodEvent),
1928  fHistManager(source.fHistManager),
1929  fEventInfo(),
1930  fName()
1931 {
1932  SetRDHFCuts(source.fRDHFCuts);
1933 }
1934 
1935 // Destructor
1937 {
1938  delete fRDHFCuts;
1939 }
1940 
1945 {
1946  new (this) AnalysisEngine(source);
1947  return *this;
1948 }
1949 
1951 void AliAnalysisTaskDmesonJets::AnalysisEngine::Init(const AliEMCALGeometry* const /*geom*/, Int_t /*runNumber*/)
1952 {
1953 }
1954 
1959 {
1960  switch (fCandidateType) {
1961  case kD0toKpi:
1962  fCandidatePDG = 421;
1963  fCandidateName = "D0";
1964  fNDaughters = 2;
1966  fPDGdaughters.Reset();
1967  fPDGdaughters[0] = 211; // pi
1968  fPDGdaughters[1] = 321; // K
1969  fBranchName = "D0toKpi";
1971  break;
1972  case kD0toKpiLikeSign:
1973  fCandidatePDG = 421;
1974  fCandidateName = "2ProngLikeSign";
1975  fNDaughters = 2;
1977  fPDGdaughters.Reset();
1978  fPDGdaughters[0] = 211; // pi
1979  fPDGdaughters[1] = 321; // K
1980  fBranchName = "LikeSign2Prong";
1981  break;
1982  case kDstartoKpipi:
1983  fCandidatePDG = 413;
1984  fCandidateName = "DStar";
1985  fNDaughters = 3;
1987  fPDGdaughters.Reset();
1988  fPDGdaughters[0] = 211; // pi soft
1989  fPDGdaughters[1] = 211; // pi fromD0
1990  fPDGdaughters[2] = 321; // K from D0
1991  fBranchName = "Dstar";
1993  break;
1994  default:
1995  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::SetCandidateProperties","Candidate %d unknown!", fCandidateType);
1996  }
1997 
1999 }
2000 
2005 {
2006  if (fRDHFCuts) delete fRDHFCuts;
2007  fRDHFCuts = cuts;
2008 }
2009 
2014 {
2015  if (!cuts) return;
2016  if (fRDHFCuts) delete fRDHFCuts;
2017  fRDHFCuts = static_cast<AliRDHFCuts*>(cuts->Clone());
2018 }
2019 
2024 {
2025  static TString name;
2026 
2027  name = TString::Format("%s_%s", GetName(), jetDef.GetName());
2028 
2029  return name.Data();
2030 }
2031 
2036 {
2038  switch (fMCMode) {
2039  case kBackgroundOnly:
2040  fName += "_BackgroundOnly";
2041  break;
2042  case kSignalOnly:
2043  fName += "_SignalOnly";
2044  break;
2045  case kMCTruth:
2046  fName += "_MCTruth";
2047  break;
2048  case kD0Reflection:
2049  fName += "_D0Reflection";
2050  break;
2051  case kOnlyWrongPIDAccepted:
2052  fName += "_OnlyWrongPIDAccepted";
2053  break;
2054  default:
2055  break;
2056  }
2057 
2058  if (fRDHFCuts) fName += TString::Format("_%s", fRDHFCuts->GetName());
2059 
2060  return fName.Data();
2061 }
2062 
2070 {
2071  std::vector<AliHFJetDefinition>::iterator it = FindJetDefinition(def);
2072 
2073  if (it == fJetDefinitions.end() || *it != def) { // No jet definition was found, adding a new one
2074  fJetDefinitions.push_back(def);
2075  ::Info("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "Jet definition '%s' has been added to analysis engine '%s'."
2076  "Total number of jet definitions is now %lu.",
2077  def.GetName(), GetName(), fJetDefinitions.size());
2078  // For detector level set maximum track pt to 100 GeV/c
2079  if (fMCMode != kMCTruth) fJetDefinitions[fJetDefinitions.size()-1].SetChargedPtRange(0., 100.);
2080  }
2081  else {
2082  ::Warning("AliAnalysisTaskDmesonJets::AnalysisEngine::AddJetDefinition", "The same jet definition '%s' was already added in analysis engine '%s'.", def.GetName(), GetName());
2083  }
2084 
2085  return &(*it);
2086 }
2087 
2099 {
2100  AliHFJetDefinition def(type, r, algo, reco);
2101 
2102  return AddJetDefinition(def);
2103 }
2104 
2110 std::vector<AliAnalysisTaskDmesonJets::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJets::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJets::AliHFJetDefinition& def)
2111 {
2112  std::vector<AliHFJetDefinition>::iterator it = fJetDefinitions.begin();
2113  while (it != fJetDefinitions.end() && (*it) != def) it++;
2114  return it;
2115 }
2116 
2121 {
2122  for (auto &jetdef : fJetDefinitions) jetdef.SetJetPhiRange(min, max);
2123 }
2124 
2129 {
2130  for (auto &jetdef : fJetDefinitions) jetdef.SetJetEtaRange(min, max);
2131 }
2132 
2137 {
2138  for (auto &jetdef : fJetDefinitions) jetdef.SetJetPtRange(min, max);
2139 }
2140 
2145 {
2146  for (auto &jetdef : fJetDefinitions) jetdef.SetChargedPtRange(min, max);
2147 }
2148 
2153 {
2154  for (auto &jetdef : fJetDefinitions) jetdef.SetNeutralPtRange(min, max);
2155 }
2156 
2163 {
2164  if (lhs.fCandidateType < rhs.fCandidateType) {
2165  return true;
2166  }
2167  else if (lhs.fCandidateType > rhs.fCandidateType) {
2168  return false;
2169  }
2170  else if (lhs.fMCMode < rhs.fMCMode) {
2171  return true;
2172  }
2173  else if (lhs.fMCMode > rhs.fMCMode) {
2174  return false;
2175  }
2176  else if (lhs.fRDHFCuts && !rhs.fRDHFCuts) {
2177  return true;
2178  }
2179  else if (lhs.fRDHFCuts && rhs.fRDHFCuts && strcmp(lhs.fRDHFCuts->GetName(), rhs.fRDHFCuts->GetName()) < 0) {
2180  return true;
2181  }
2182  else {
2183  return false;
2184  }
2185 }
2186 
2193 {
2194  if (lhs.fCandidateType != rhs.fCandidateType) return false;
2195  if (lhs.fMCMode != rhs.fMCMode) return false;
2196  if (lhs.fRDHFCuts == nullptr && rhs.fRDHFCuts != nullptr) return false;
2197  if (lhs.fRDHFCuts != nullptr && rhs.fRDHFCuts == nullptr) return false;
2198  if (lhs.fRDHFCuts && rhs.fRDHFCuts && strcmp(lhs.fRDHFCuts->GetName(), rhs.fRDHFCuts->GetName()) != 0) return false;
2199  return true;
2200 }
2201 
2210 {
2211  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) { // D0 candidate
2212  return ExtractD0Attributes(Dcand, DmesonJet, i);
2213  }
2214  else if (fCandidateType == kDstartoKpipi) { // Dstar candidate
2215  return ExtractDstarAttributes(static_cast<const AliAODRecoCascadeHF*>(Dcand), DmesonJet, i);
2216  }
2217  else {
2218  return kFALSE;
2219  }
2220 }
2221 
2230 {
2231  AliDebug(10,"Checking if D0 meson is selected");
2232  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
2233  if (isSelected == 0) return kFALSE;
2234 
2235  Int_t MCtruthPdgCode = 0;
2236 
2237  Double_t invMassD = 0;
2238 
2239  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
2240  // Checks also the origin, and if it matches the rejected origin mask, return false
2241  if (fMCMode != kNoMC) {
2242  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
2243  DmesonJet.fMCLabel = mcLab;
2244 
2245  // Retrieve the generated particle (if exists) and its PDG code
2246  if (mcLab >= 0) {
2247  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
2248 
2249  if (aodMcPart) {
2250  // Check origin and return false if it matches the rejected origin mask
2251  if (fRejectedOrigin) {
2252  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
2253  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
2254  }
2255  MCtruthPdgCode = aodMcPart->PdgCode();
2256  }
2257  }
2258  }
2259 
2260  if (isSelected == 1) { // selected as a D0
2261  if (i != 0) return kFALSE; // only one mass hypothesis thanks to PID
2262 
2263  if (fMCMode == kNoMC ||
2264  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
2265  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly) ||
2266  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kD0Reflection || fMCMode == kOnlyWrongPIDAccepted))) {
2267  // 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)
2268  AliDebug(10,"Selected as D0");
2269  invMassD = Dcand->InvMassD0();
2270  }
2271  else { // conditions above not passed, so return FALSE
2272  return kFALSE;
2273  }
2274  }
2275  else if (isSelected == 2) { // selected as a D0bar
2276  if (i != 1) return kFALSE; // only one mass hypothesis thanks to PID
2277 
2278  if (fMCMode == kNoMC ||
2279  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
2280  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly) ||
2281  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kD0Reflection || fMCMode == kOnlyWrongPIDAccepted))) {
2282  // 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)
2283  AliDebug(10,"Selected as D0bar");
2284  invMassD = Dcand->InvMassD0bar();
2285  }
2286  else { // conditions above not passed, so return FALSE
2287  return kFALSE;
2288  }
2289  }
2290  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
2291  AliDebug(10,"Selected as either D0 or D0bar");
2292 
2293  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
2294  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
2295  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kD0Reflection))) {
2296  if (i != 0) return kFALSE;
2297  AliDebug(10, "MC truth is D0");
2298  invMassD = Dcand->InvMassD0();
2299  }
2300  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
2301  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kD0Reflection))) {
2302  if (i != 1) return kFALSE;
2303  AliDebug(10, "MC truth is D0bar");
2304  invMassD = Dcand->InvMassD0bar();
2305  }
2306  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
2307 
2308  // Only accept it if background-only OR background-and-signal was requested
2309  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
2310  // Select D0 or D0bar depending on the i-parameter
2311  if (i == 0) {
2312  AliDebug(10, "Returning invariant mass with D0 hypothesis");
2313  invMassD = Dcand->InvMassD0();
2314  }
2315  else if (i == 1) {
2316  AliDebug(10, "Returning invariant mass with D0bar hypothesis");
2317  invMassD = Dcand->InvMassD0bar();
2318  }
2319  else { // i > 1
2320  return kFALSE;
2321  }
2322  }
2323  else { // signal-only was requested but this is not a true D0
2324  return kFALSE;
2325  }
2326  }
2327  }
2328  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
2329  return kTRUE;
2330 }
2331 
2340 {
2341  AliDebug(10,"Checking if D* meson is selected");
2342  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoCascadeHF*>(DstarCand), AliRDHFCuts::kAll, fAodEvent);
2343  if (isSelected == 0) return kFALSE;
2344 
2345  if ((i == 1 && DstarCand->Charge()>0) || (i == 0 && DstarCand->Charge()<0) || i > 1) return kFALSE; // only one mass hypothesis for the D*
2346 
2347  Int_t MCtruthPdgCode = 0;
2348 
2349  Double_t invMassD = 0;
2350 
2351  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
2352  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
2353  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
2354 
2355  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
2356  AliDebug(10, Form("MC label is %d", mcLab));
2357  DmesonJet.fMCLabel = mcLab;
2358  if (mcLab >= 0) {
2359  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
2360 
2361  if (aodMcPart) {
2362  if (fRejectedOrigin) {
2363  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
2364  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
2365  }
2366 
2367  MCtruthPdgCode = aodMcPart->PdgCode();
2368  AliDebug(10, Form("MC truth pdg code is %d",MCtruthPdgCode));
2369  }
2370  }
2371  }
2372 
2373  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
2374  if (fMCMode == kNoMC ||
2375  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
2376  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
2377  // 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)
2378  invMassD = DstarCand->InvMassDstarKpipi();
2379  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
2380  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
2381  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
2382  return kTRUE;
2383  }
2384  else { // conditions above not passed, so return FALSE
2385  return kFALSE;
2386  }
2387 }
2388 
2396 {
2397  if (!part) return kUnknownDecay;
2398  if (!mcArray) return kUnknownDecay;
2399 
2401 
2402  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
2403 
2404  if (part->GetNDaughters() == 2) {
2405 
2406  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
2407  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
2408 
2409  if (!d1 || !d2) {
2410  return decay;
2411  }
2412 
2413  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
2414  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
2415 
2416  if (absPdgPart == 421) { // D0 -> K pi
2417  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
2418  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
2419  decay = kDecayD0toKpi;
2420  }
2421  }
2422 
2423  if (absPdgPart == 413) { // D* -> D0 pi
2424  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
2425  Int_t D0decay = CheckDecayChannel(d1, mcArray);
2426  if (D0decay == kDecayD0toKpi) {
2427  decay = kDecayDStartoKpipi;
2428  }
2429  }
2430  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
2431  Int_t D0decay = CheckDecayChannel(d2, mcArray);
2432  if (D0decay == kDecayD0toKpi) {
2433  decay = kDecayDStartoKpipi;
2434  }
2435  }
2436  }
2437  }
2438 
2439  return decay;
2440 }
2441 
2448 std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> AliAnalysisTaskDmesonJets::AnalysisEngine::IsPromptCharm(const AliAODMCParticle* part, TClonesArray* mcArray)
2449 {
2450  std::pair<AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle*> result(kUnknownQuark, 0);
2451 
2452  if (!part) return result;
2453  if (!mcArray) return result;
2454 
2455  static std::set<UInt_t> partons = { 4, 5 };
2456 
2457  AliAODMCParticle* parton = FindParticleOrigin(part, mcArray, kFindLast, partons);
2458  if (parton) {
2459  result.second = parton;
2460  UInt_t absPdgParton = TMath::Abs(parton->GetPdgCode());
2461  if (absPdgParton == 4) result.first = kFromCharm;
2462  else if (absPdgParton == 5) result.first = kFromBottom;
2463  }
2464 
2465  return result;
2466 }
2467 
2475 
2476 AliAODMCParticle* AliAnalysisTaskDmesonJets::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode)
2477 {
2478  static std::set<UInt_t> pdgSet;
2479 
2480  return FindParticleOrigin(part, mcArray, mode, pdgSet);
2481 }
2482 
2491 
2492 AliAODMCParticle* AliAnalysisTaskDmesonJets::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode, const std::set<UInt_t>& pdgSet)
2493 {
2494  AliAODMCParticle* result = nullptr;
2495 
2496  Int_t mother = part->GetMother();
2497  while (mother >= 0) {
2498  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
2499  if (mcGranma) {
2500  UInt_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
2501 
2502  // If the current particle is one of the particle types that is being searched assign it to the result pointer
2503  if (pdgSet.empty() || pdgSet.count(abspdgGranma) > 0) {
2504  result = mcGranma;
2505  // If the last particle in the fragmentation tree (first when going reverse) was requested then stop the loop
2506  if (mode == kFindLast) break;
2507  }
2508  if (mother == mcGranma->GetMother()) { // avoid infinite loop!
2509  AliWarningClassStream() << "Particle " << mother << " (PDG=" << mcGranma->PdgCode() << ") is the mother of itself!?" << std::endl;
2510  break;
2511  }
2512  mother = mcGranma->GetMother();
2513  }
2514  else {
2515  AliErrorClassStream() << "Could not retrieve mother particle " << mother << "!" << std::endl;
2516  break;
2517  }
2518  }
2519 
2520  return result;
2521 }
2522 
2525 {
2526  for (auto& jetDef : fJetDefinitions) {
2527  jetDef.fJets.clear();
2528  }
2529 
2530  if (fMCMode == kMCTruth) {
2532  }
2533  else {
2535  }
2536 }
2537 
2540 {
2541  // Fill the vertex info of the candidates
2542  // Needed for reduced delta AOD, where the vertex info has been deleted
2543  // to reduce the delta AOD file size
2545 
2546  const Int_t nD = fCandidateArray->GetEntriesFast();
2547 
2548  AliDmesonJetInfo DmesonJet;
2549  DmesonJet.fEvent = this->fAodEvent;
2550 
2551  std::map<AliHFJetDefinition*,Double_t> maxJetPt;
2552  for (auto& def : fJetDefinitions) maxJetPt[&def] = 0;
2553  Double_t maxDPt = 0;
2554 
2555  std::array<int, 3> nAccCharm = {0};
2556  std::array<std::array<int, 3>, 5> nAccCharmPt = {{{0}}};
2557 
2558  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
2559  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
2560  if (!charmCand) continue;
2561  if(!(vHF.FillRecoCand(fAodEvent,charmCand))) continue;
2562 
2563  // region of interest + cuts
2564  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
2565  Int_t nMassHypo = 0; // number of mass hypothesis accepted for this D meson
2566  if (charmCand->Pt() > maxDPt) maxDPt = charmCand->Pt();
2567  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
2568  DmesonJet.Reset();
2569  DmesonJet.fDmesonParticle = charmCand;
2570  DmesonJet.fSelectionType = im + 1;
2571  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
2572  for (auto& def : fJetDefinitions) {
2573  if (FindJet(charmCand, DmesonJet, def)) {
2574  Double_t jetPt = DmesonJet.fJets[def.GetName()].fMomentum.Pt();
2575  if (jetPt > maxJetPt[&def]) maxJetPt[&def] = jetPt;
2576  }
2577  else {
2578  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
2579  def.GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
2580  }
2581  }
2582  fDmesonJets[(icharm+1)*(1-(im*2))] = DmesonJet;
2583  nMassHypo++;
2584  nAccCharm[im]++;
2585 
2586  for (int i = 0; i < nAccCharmPt.size(); i++) {
2587  if (charmCand->Pt() < i) break;
2588  nAccCharmPt[i][im]++;
2589  }
2590  }
2591  }
2592  if (nMassHypo == 2) { // both mass hypothesis accepted
2593  nAccCharm[0]--;
2594  nAccCharm[1]--;
2595  nAccCharm[2]++;
2596 
2597  for (int i = 0; i < nAccCharmPt.size(); i++) {
2598  if (charmCand->Pt() < i) break;
2599  nAccCharmPt[i][0]--;
2600  nAccCharmPt[i][1]--;
2601  nAccCharmPt[i][2]++;
2602  }
2603 
2604  fDmesonJets[(icharm+1)].fD0D0bar = kTRUE;
2605  fDmesonJets[-(icharm+1)].fD0D0bar = kTRUE;
2606  }
2607  } // end of D cand loop
2608 
2609  TString hname;
2610 
2611  Int_t ntracks = 0;
2612 
2613  for (auto track_cont : fTrackContainers) {
2614  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
2615  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(nullptr);
2616  ntracks += track_cont->GetNAcceptEntries();
2617  }
2618 
2619  for (auto& def : fJetDefinitions) {
2620  if (!def.fRho) continue;
2621  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
2622  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
2623 
2624  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
2625  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
2626 
2627  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
2628  fHistManager->FillTH2(hname, fEventInfo.fCent, def.fRho->GetVal());
2629 
2630  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
2631  fHistManager->FillTH2(hname, fEventInfo.fCent, maxJetPt[&def]);
2632 
2633  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
2634  fHistManager->FillTH2(hname, fEventInfo.fCent, maxDPt);
2635 
2636  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
2637  fHistManager->FillTH2(hname, ntracks, def.fRho->GetVal());
2638 
2639  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
2640  fHistManager->FillTH2(hname, ntracks, maxJetPt[&def]);
2641 
2642  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
2643  fHistManager->FillTH2(hname, ntracks, maxDPt);
2644  }
2645 
2646  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
2647  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
2648  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
2649  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
2650 
2651  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
2652  fHistManager->FillTH2(hname, ntracks, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
2653 
2654  for (int i = 0; i < nAccCharmPt.size(); i++) {
2655  hname = TString::Format("%s/fHistNTotAcceptedDmesonsPt%d", GetName(), i);
2656  fHistManager->FillTH1(hname, "D", nAccCharmPt[i][0]);
2657  fHistManager->FillTH1(hname, "Anti-D", nAccCharmPt[i][1]);
2658  fHistManager->FillTH1(hname, "Both", nAccCharmPt[i][2]);
2659 
2660  hname = TString::Format("%s/fHistNAcceptedDmesonsPt%d", GetName(), i);
2661  fHistManager->FillTH1(hname, nAccCharmPt[i][0]+nAccCharmPt[i][1]+nAccCharmPt[i][2]);
2662  }
2663 
2664  hname = TString::Format("%s/fHistNDmesons", GetName());
2665  fHistManager->FillTH1(hname, nD);
2666 }
2667 
2678 {
2679  TString hname;
2680 
2681  Double_t rho = 0;
2682  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
2683 
2685  fFastJetWrapper->SetR(jetDef.fRadius);
2688 
2689  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
2690 
2691  if (jetDef.fJetType != AliJetContainer::kNeutralJet) {
2692  for (auto track_cont : fTrackContainers) {
2693  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
2694  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(Dcand);
2695  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
2696  AddInputVectors(track_cont, 100, static_cast<TH2*>(fHistManager->FindObject(hname)), fTrackEfficiency);
2697 
2698  if (hftrack_cont) {
2699  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
2700  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
2701  const TObjArray& daughters = hftrack_cont->GetDaughterList();
2702  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
2703  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
2704  if (!hftrack_cont->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
2705  }
2706  }
2707  }
2708  }
2709 
2710  if (jetDef.fJetType != AliJetContainer::kChargedJet) {
2711  for (auto clus_cont : fClusterContainers) {
2712  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
2713  AddInputVectors(clus_cont, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
2714  }
2715  }
2716 
2717  // run jet finder
2718  fFastJetWrapper->Run();
2719 
2720  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
2721 
2722  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
2723  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
2724 
2725  Bool_t isDmesonJet = kFALSE;
2726 
2727  Double_t maxChPt = 0;
2728  Double_t maxNePt = 0;
2729  Double_t totalNeutralPt = 0;
2730  Int_t nConst = 1;
2731 
2732  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
2733  if (constituents[ic].user_index() == 0) {
2734  isDmesonJet = kTRUE;
2735  }
2736  else if (constituents[ic].user_index() >= 100) {
2737  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
2738  nConst++;
2739  }
2740  else if (constituents[ic].user_index() <= -100) {
2741  totalNeutralPt += constituents[ic].pt();
2742  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
2743  nConst++;
2744  }
2745  }
2746 
2747  if (isDmesonJet) {
2748  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
2749  DmesonJet.fJets[jetDef.GetName()].fNConstituents = nConst;
2750  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
2751  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
2752  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
2753  DmesonJet.fJets[jetDef.GetName()].fArea = jets_incl[ijet].area();
2754  DmesonJet.fJets[jetDef.GetName()].fCorrPt = DmesonJet.fJets[jetDef.GetName()].fMomentum.Pt() - jets_incl[ijet].area() * rho;
2755 
2756  return kTRUE;
2757  }
2758  }
2759 
2760  return kFALSE;
2761 }
2762 
2766 void AliAnalysisTaskDmesonJets::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist, Double_t eff)
2767 {
2768  auto itcont = cont->all_momentum();
2769  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
2770  UInt_t rejectionReason = 0;
2771  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
2772  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
2773  continue;
2774  }
2775  if (fRandomGen && eff > 0 && eff < 1) {
2776  Double_t rnd = fRandomGen->Rndm();
2777  if (eff < rnd) {
2778  if (rejectHist) rejectHist->Fill(6, it->first.Pt());
2779  continue;
2780  }
2781  }
2782  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
2783  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
2784  }
2785 }
2786 
2789 {
2790  TString hname;
2791 
2796 
2797  if (!fMCContainer->IsSpecialPDGFound()) return;
2798 
2799  std::array<int,2> nAccCharm = {0};
2800  std::array<std::array<int, 2>, 5> nAccCharmPt = {{{0}}};
2801 
2802  std::map<AliHFJetDefinition*, Double_t> maxJetPt;
2803  Double_t maxDPt = 0;
2804 
2805  for (auto &jetDef : fJetDefinitions) {
2806  maxJetPt[&jetDef] = 0;
2807  Double_t rho = 0;
2808  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
2809  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", GetName(), jetDef.GetName());
2810  TH1* histNDmesonsVsNconstituents = static_cast<TH1*>(fHistManager->FindObject(hname));
2811 
2812  switch (jetDef.fJetType) {
2814  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
2815  break;
2817  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
2818  break;
2820  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
2821  break;
2822  }
2823 
2825  fFastJetWrapper->SetR(jetDef.fRadius);
2828 
2829  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
2830  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
2831 
2832  fFastJetWrapper->Run();
2833 
2834  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
2835 
2836  for (auto jet : jets_incl) {
2837  Int_t nDmesonsInJet = 0;
2838 
2839  for (auto constituent : jet.constituents()) {
2840  Int_t iPart = constituent.user_index() - 100;
2841  if (constituent.perp() < 1e-6) continue; // reject ghost particles
2842  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
2843  if (!part) {
2844  ::Error("AliAnalysisTaskDmesonJets::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
2845  continue;
2846  }
2847  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
2848  nDmesonsInJet++;
2849  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
2850  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
2851  if (part->Pt() > maxDPt) maxDPt = part->Pt();
2852  std::pair<int, AliDmesonJetInfo> element;
2853  element.first = iPart;
2854  dMesonJetIt = fDmesonJets.insert(element).first;
2855  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
2856  (*dMesonJetIt).second.fDmesonParticle = part;
2857  (*dMesonJetIt).second.fSelectionType = part->PdgCode() > 0 ? 1 : 2;
2858 
2859  UShort_t p = 0;
2860  UInt_t rs = 0;
2861 
2862  auto origin = IsPromptCharm(part, fMCContainer->GetArray());
2863  p = 0;
2864  rs = origin.first;
2865  while (rs >>= 1) { p++; }
2866  (*dMesonJetIt).second.fPartonType = p;
2867  (*dMesonJetIt).second.fParton = origin.second;
2868 
2869  (*dMesonJetIt).second.fAncestor = FindParticleOrigin(part, fMCContainer->GetArray(), kFindFirst);
2870 
2871  Int_t im = -1;
2872  if (part->PdgCode() > 0) { // D0
2873  im = 0;
2874  }
2875  else { // D0bar
2876  im = 1;
2877  }
2878 
2879  nAccCharm[im]++;
2880  for (int i = 0; i < nAccCharmPt.size(); i++) {
2881  if (part->Pt() < i) break;
2882  nAccCharmPt[i][im]++;
2883  }
2884  }
2885 
2886  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
2887  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = jet.constituents().size();
2888  (*dMesonJetIt).second.fJets[jetDef.GetName()].fArea = jet.area();
2889  (*dMesonJetIt).second.fJets[jetDef.GetName()].fCorrPt = (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.Pt() - jet.area() * rho;
2890  if (jet.perp() > maxJetPt[&jetDef]) maxJetPt[&jetDef] = jet.perp();
2891  } // if constituent is a D meson
2892  } // for each constituent
2893  if (nDmesonsInJet > 0) histNDmesonsVsNconstituents->Fill(jet.constituents().size(), nDmesonsInJet);
2894  } // for each jet
2895  } // for each jet definition
2896 
2898 
2899  for (auto& def : fJetDefinitions) {
2900  if (!def.fRho) continue;
2901  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
2902  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
2903 
2904  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
2905  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
2906 
2907  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
2908  fHistManager->FillTH2(hname, fEventInfo.fCent, def.fRho->GetVal());
2909 
2910  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
2911  fHistManager->FillTH2(hname, fEventInfo.fCent, maxJetPt[&def]);
2912 
2913  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
2914  fHistManager->FillTH2(hname, fEventInfo.fCent, maxDPt);
2915 
2916  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
2917  fHistManager->FillTH2(hname, npart, def.fRho->GetVal());
2918 
2919  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
2920  fHistManager->FillTH2(hname, npart, maxJetPt[&def]);
2921 
2922  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
2923  fHistManager->FillTH2(hname, npart, maxDPt);
2924  }
2925 
2926  if (fDmesonJets.size() != nAccCharm[0]+nAccCharm[1]) AliError(Form("I found %lu mesons (%d)?", fDmesonJets.size(), nAccCharm[0]+nAccCharm[1]));
2927  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
2928  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
2929  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
2930 
2931  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
2932  fHistManager->FillTH2(hname, npart, nAccCharm[0]+nAccCharm[1]);
2933 
2934  for (int i = 0; i < nAccCharmPt.size(); i++) {
2935  hname = TString::Format("%s/fHistNTotAcceptedDmesonsPt%d", GetName(), i);
2936  fHistManager->FillTH1(hname, "D", nAccCharmPt[i][0]);
2937  fHistManager->FillTH1(hname, "Anti-D", nAccCharmPt[i][1]);
2938 
2939  hname = TString::Format("%s/fHistNAcceptedDmesonsPt%d", GetName(), i);
2940  fHistManager->FillTH1(hname, nAccCharmPt[i][0]+nAccCharmPt[i][1]);
2941  }
2942 
2943  hname = TString::Format("%s/fHistNDmesons", GetName());
2944  fHistManager->FillTH1(hname, nAccCharm[0]+nAccCharm[1]); // same as the number of accepted D mesons, since no selection is performed
2945 }
2946 
2947 
2948 
2949 // Definitions of class AliAnalysisTaskDmesonJets
2950 
2952 ClassImp(AliAnalysisTaskDmesonJets);
2954 
2958  fAnalysisEngines(),
2959  fEnabledAxis(0),
2961  fHistManager(),
2962  fApplyKinematicCuts(kTRUE),
2963  fNOutputTrees(0),
2964  fTrackEfficiency(0),
2965  fRejectISR(kFALSE),
2966  fJetAreaType(fastjet::active_area),
2967  fJetGhostArea(0.005),
2968  fMCContainer(0),
2969  fAodEvent(0),
2970  fFastJetWrapper(0)
2971 {
2972  SetMakeGeneralHistograms(kTRUE);
2973 }
2974 
2979  AliAnalysisTaskEmcalLight(name, kTRUE),
2980  fAnalysisEngines(),
2983  fHistManager(name),
2984  fApplyKinematicCuts(kTRUE),
2985  fNOutputTrees(nOutputTrees),
2986  fTrackEfficiency(0),
2987  fRejectISR(kFALSE),
2988  fJetAreaType(fastjet::active_area),
2989  fJetGhostArea(0.005),
2990  fMCContainer(0),
2991  fAodEvent(0),
2992  fFastJetWrapper(0)
2993 {
2994  SetMakeGeneralHistograms(kTRUE);
2995  for (Int_t i = 0; i < nOutputTrees; i++){
2996  DefineOutput(2+i, TTree::Class());
2997  }
2998 }
2999 
3002 {
3003  if (fFastJetWrapper) delete fFastJetWrapper;
3004 }
3005 
3013 {
3014  AliRDHFCuts* analysiscuts = 0;
3015  TFile* filecuts = TFile::Open(cutfname);
3016  if (!filecuts || filecuts->IsZombie()) {
3017  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Input file not found: will use std cuts.");
3018  filecuts = 0;
3019  }
3020 
3021  if (filecuts) analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
3022 
3023  if (!analysiscuts) {
3024  ::Error("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Could not find analysis cuts '%s' in '%s'.", cutsname.Data(), cutfname.Data());
3025  if (filecuts) {
3026  filecuts->ls();
3027  }
3028  }
3029  else {
3030  ::Info("AliAnalysisTaskDmesonJets::LoadDMesonCutsFromFile", "Cuts '%s' loaded from file '%s'", cutsname.Data(), cutfname.Data());
3031  }
3032 
3033  return analysiscuts;
3034 }
3035 
3048 {
3050  return AddAnalysisEngine(type, cutfname, cuttype, MCmode, jetDef, rhoName);
3051 }
3052 
3064 {
3065  AliRDHFCuts* cuts = 0;
3066 
3067  if (!cutfname.IsNull()) {
3068  TString cutsname;
3069 
3070  switch (type) {
3071  case kD0toKpi:
3072  case kD0toKpiLikeSign:
3073  cutsname = "D0toKpiCuts";
3074  break;
3075  case kDstartoKpipi:
3076  cutsname = "DStartoKpipiCuts";
3077  break;
3078  default:
3079  return 0;
3080  }
3081 
3082  if (!cuttype.IsNull()) {
3083  cutsname += TString::Format("_%s", cuttype.Data());
3084  }
3085 
3086  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
3087  if (cuts) cuts->PrintAll();
3088  }
3089 
3090  AnalysisEngine eng(type, MCmode, cuts);
3091 
3092  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
3093 
3094  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
3095  eng.AddJetDefinition(jetDef);
3096  it = fAnalysisEngines.insert(it, eng);
3097  ::Info("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
3098  }
3099  else {
3100  AnalysisEngine* found_eng = &(*it);
3101  ::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());
3102  found_eng->AddJetDefinition(jetDef);
3103 
3104  if (cuts) {
3105  if (found_eng->fRDHFCuts != 0) ::Warning("AliAnalysisTaskDmesonJets::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
3106  found_eng->SetRDHFCuts(cuts);
3107  }
3108  }
3109 
3110  return &(*it);
3111 }
3112 
3113 std::list<AliAnalysisTaskDmesonJets::AnalysisEngine>::iterator AliAnalysisTaskDmesonJets::FindAnalysisEngine(const AliAnalysisTaskDmesonJets::AnalysisEngine& eng)
3114 {
3115  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
3116  while (it != fAnalysisEngines.end() && (*it) != eng) it++;
3117  return it;
3118 }
3119 
3122 {
3123  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
3124 
3126 
3127  // Define histograms
3128  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
3129 
3130  TString hname;
3131  TString htitle;
3132  TH1* h = 0;
3133  Int_t treeSlot = 0;
3134 
3135  Int_t maxTracks = 6000;
3136  Double_t maxRho = 500;
3137  if (fForceBeamType == kpp) {
3138  maxRho = 50;
3139  maxTracks = 200;
3140  }
3141  else if (fForceBeamType == kpA) {
3142  maxRho = 200;
3143  maxTracks = 500;
3144  }
3145 
3146  hname = "fHistCharmPt";
3147  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
3148  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3149 
3150  hname = "fHistCharmEta";
3151  htitle = hname + ";#eta_{charm};counts";
3152  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
3153 
3154  hname = "fHistCharmPhi";
3155  htitle = hname + ";#phi_{charm};counts";
3156  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
3157 
3158  hname = "fHistBottomPt";
3159  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
3160  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3161 
3162  hname = "fHistBottomEta";
3163  htitle = hname + ";#eta_{bottom};counts";
3164  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
3165 
3166  hname = "fHistBottomPhi";
3167  htitle = hname + ";#phi_{bottom};counts";
3168  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
3169 
3170  hname = "fHistHighestPartonPt";
3171  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
3172  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3173 
3174  hname = "fHistHighestPartonType";
3175  htitle = hname + ";type;counts";
3176  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
3177 
3178  hname = "fHistNHeavyQuarks";
3179  htitle = hname + ";number of heavy-quarks;counts";
3180  fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
3181 
3182  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
3183  for (auto &param : fAnalysisEngines) {
3184  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
3185 
3186  param.fHistManager = &fHistManager;
3187 
3188  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", param.GetName());
3189  htitle = hname + ";#it{N}_{tracks};#it{N}_{D};events";
3190  h = fHistManager.CreateTH2(hname, htitle, 251, -0.5, 250.5, 21, -0.5, 20.5);
3191 
3192  hname = TString::Format("%s/fHistNTotAcceptedDmesons", param.GetName());
3193  htitle = hname + ";;#it{N}_{D}";
3194  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
3195 
3196  for (int i = 0 ; i < 5; i++) {
3197  hname = TString::Format("%s/fHistNAcceptedDmesonsPt%d", param.GetName(), i);
3198  htitle = hname + ";#it{N}_{D};events";
3199  h = fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
3200 
3201  hname = TString::Format("%s/fHistNTotAcceptedDmesonsPt%d", param.GetName(), i);
3202  htitle = hname + ";;#it{N}_{D}";
3203  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
3204  }
3205 
3206  hname = TString::Format("%s/fHistNDmesons", param.GetName());
3207  htitle = hname + ";#it{N}_{D};events";
3208  h = fHistManager.CreateTH1(hname, htitle, 501, -0.5, 500.5);
3209 
3210  hname = TString::Format("%s/fHistNEvents", param.GetName());
3211  htitle = hname + ";Event status;counts";
3212  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
3213  h->GetXaxis()->SetBinLabel(1, "Accepted");
3214  h->GetXaxis()->SetBinLabel(2, "Rejected");
3215 
3216  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
3217  htitle = hname + ";Rejection reason;counts";
3218  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
3219 
3220  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
3221  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
3222  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
3223 
3224  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
3225  htitle = hname + ";#it{#eta}_{D};counts";
3226  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
3227 
3228  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
3229  htitle = hname + ";#it{#phi}_{D};counts";
3230  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
3231 
3232  if (param.fMCMode != kMCTruth) {
3233  if (param.fCandidateType == kD0toKpi || param.fCandidateType == kD0toKpiLikeSign) {
3234  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
3235  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
3236  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
3237  }
3238  else if (param.fCandidateType == kDstartoKpipi) {
3239  Double_t min = 0;
3240  Double_t max = 0;
3241 
3242  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
3243  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
3244  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
3245  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
3246 
3247  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
3248  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
3249  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
3250  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
3251  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
3252  }
3253  }
3254 
3255  if (param.fMCMode == kMCTruth) {
3256  hname = TString::Format("%s/fHistPartonPt", param.GetName());
3257  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
3258  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3259 
3260  hname = TString::Format("%s/fHistPartonEta", param.GetName());
3261  htitle = hname + ";#eta_{parton};counts";
3262  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
3263 
3264  hname = TString::Format("%s/fHistPartonPhi", param.GetName());
3265  htitle = hname + ";#phi_{parton};counts";
3266  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
3267 
3268  hname = TString::Format("%s/fHistPartonType", param.GetName());
3269  htitle = hname + ";type;counts";
3270  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
3271 
3272  hname = TString::Format("%s/fHistPrompt", param.GetName());
3273  htitle = hname + ";Type;counts";
3274  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
3275  h->GetXaxis()->SetBinLabel(1, "Unknown");
3276  h->GetXaxis()->SetBinLabel(2, "Prompt");
3277  h->GetXaxis()->SetBinLabel(3, "Non-Prompt");
3278 
3279  hname = TString::Format("%s/fHistAncestor", param.GetName());
3280  htitle = hname + ";Ancestor;counts";
3281  h = fHistManager.CreateTH1(hname, htitle, 4, 0, 4);
3282  h->GetXaxis()->SetBinLabel(1, "Unknown");
3283  h->GetXaxis()->SetBinLabel(2, "Charm");
3284  h->GetXaxis()->SetBinLabel(3, "Bottom");
3285  h->GetXaxis()->SetBinLabel(4, "Proton");
3286  }
3287 
3288  for (auto& jetDef : param.fJetDefinitions) {
3289  ::Info("AliAnalysisTaskDmesonJets::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef.GetName());
3290 
3291  if (param.fMCMode == kMCTruth) {
3292  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", param.GetName(), jetDef.GetName());
3293  htitle = hname + ";#it{N}_{constituents};#it{N}_{D};counts";
3294  h = fHistManager.CreateTH2(hname, htitle, 51, -0.5, 50.5, 10, 0.5, 10.5);
3295  }
3296 
3297  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef.GetName());
3298  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
3299  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
3300  SetRejectionReasonLabels(h->GetXaxis());
3301 
3302  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef.GetName());
3303  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
3304  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
3305  SetRejectionReasonLabels(h->GetXaxis());
3306 
3307  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef.GetName());
3308  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
3309  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
3310  SetRejectionReasonLabels(h->GetXaxis());
3311 
3312  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef.GetName());
3313  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
3314  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
3315 
3316  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef.GetName());
3317  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
3318  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
3319 
3320  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef.GetName());
3321  htitle = hname + ";#it{#eta}_{jet};counts";
3322  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
3323 
3324  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef.GetName());
3325  htitle = hname + ";#it{#phi}_{jet};counts";
3326  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
3327 
3328  if (!jetDef.fRhoName.IsNull()) {
3329  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", param.GetName(), jetDef.GetName());
3330  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});#rho (GeV/#it{c} #times rad^{-1});counts";
3331  fHistManager.CreateTH2(hname, htitle, 300, 0, 150, 1000, 0, maxRho);
3332 
3333  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", param.GetName(), jetDef.GetName());
3334  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});#rho (GeV/#it{c} #times rad^{-1});counts";
3335  fHistManager.CreateTH2(hname, htitle, 300, 0, 150, 1000, 0, maxRho);
3336 
3337  hname = TString::Format("%s/%s/fHistRhoVsCent", param.GetName(), jetDef.GetName());
3338  htitle = hname + ";Centrality (%);#rho (GeV/#it{c} #times rad^{-1});counts";
3339  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 1000, 0, maxRho);
3340 
3341  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", param.GetName(), jetDef.GetName());
3342  htitle = hname + ";Centrality (%);#it{p}_{T,jet} (GeV/#it{c});counts";
3343  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
3344 
3345  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", param.GetName(), jetDef.GetName());
3346  htitle = hname + ";Centrality (%);#it{p}_{T,D} (GeV/#it{c});counts";
3347  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
3348 
3349  hname = TString::Format("%s/%s/fHistRhoVsNTracks", param.GetName(), jetDef.GetName());
3350  htitle = hname + ";no. of tracks;#rho (GeV/#it{c} #times rad^{-1});counts";
3351  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 1000, 0, maxRho);
3352 
3353  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", param.GetName(), jetDef.GetName());
3354  htitle = hname + ";no. of tracks;#it{p}_{T,jet} (GeV/#it{c});counts";
3355  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
3356 
3357  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", param.GetName(), jetDef.GetName());
3358  htitle = hname + ";no. of tracks;#it{p}_{T,D} (GeV/#it{c});counts";
3359  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
3360  }
3361  }
3362  switch (fOutputType) {
3363  case kTreeOutput:
3364  {
3365  OutputHandlerTTree* tree_handler = new OutputHandlerTTree(&param);
3366  param.fOutputHandler = tree_handler;
3367  tree_handler->BuildOutputObject(GetName());
3368  if (treeSlot < fNOutputTrees) {
3369  tree_handler->AssignDataSlot(treeSlot+2);
3370  treeSlot++;
3371  PostDataFromAnalysisEngine(tree_handler);
3372  }
3373  else {
3374  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
3375  }
3376  }
3377  break;
3378  case kTHnOutput:
3379  {
3380  OutputHandlerTHnSparse* thnsparse_handler = new OutputHandlerTHnSparse(&param);
3381  param.fOutputHandler = thnsparse_handler;
3382  thnsparse_handler->SetEnabledAxis(fEnabledAxis);
3383  thnsparse_handler->BuildOutputObject(GetName());
3384  }
3385  break;
3386  case kOnlyQAOutput:
3387  {
3388  OutputHandler* qa_handler = new OutputHandler(&param);
3389  param.fOutputHandler = qa_handler;
3390  }
3391  break;
3392  case kNoOutput:
3393  break;
3394  case kTreeExtendedOutput:
3395  {
3397  param.fOutputHandler = tree_handler;
3398  tree_handler->BuildOutputObject(GetName());
3399  if (treeSlot < fNOutputTrees) {
3400  tree_handler->AssignDataSlot(treeSlot+2);
3401  treeSlot++;
3402  PostDataFromAnalysisEngine(tree_handler);
3403  }
3404  else {
3405  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
3406  }
3407  }
3408  break;
3409  }
3410  }
3411 
3413 
3414  PostData(1, fOutput);
3415 }
3416 
3420 {
3422 
3423  // Load the event
3424  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
3425 
3426  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
3427 
3428  fFastJetWrapper->SetAreaType((fastjet::AreaType)fJetAreaType);
3430 
3431  if (!fAodEvent) {
3432  AliError(Form("This task need an AOD event (Task '%s'). Expect troubles...", GetName()));
3433  //return;
3434  }
3435 
3436  TRandom* rnd = 0;
3437  if (fTrackEfficiency > 0 && fTrackEfficiency < 1) rnd = new TRandom3(0);
3438 
3439  for (auto cont_it : fParticleCollArray) {
3440  AliHFAODMCParticleContainer* part_cont = dynamic_cast<AliHFAODMCParticleContainer*>(cont_it.second);
3441  if (part_cont) fMCContainer = part_cont;
3442  }
3443 
3444  for (auto &params : fAnalysisEngines) {
3445 
3446  params.fAodEvent = fAodEvent;
3447  params.fFastJetWrapper = fFastJetWrapper;
3448  params.fTrackEfficiency = fTrackEfficiency;
3449  params.fRejectISR = fRejectISR;
3450  params.fRandomGen = rnd;
3451 
3452  for (auto &jetdef: params.fJetDefinitions) {
3453  if (!jetdef.fRhoName.IsNull()) {
3454  jetdef.fRho = dynamic_cast<AliRhoParameter*>(fInputEvent->FindListObject(jetdef.fRhoName));
3455  if (!jetdef.fRho) {
3456  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3457  "%s: Could not find rho object '%s' for engine '%s'",
3458  GetName(), jetdef.fRhoName.Data(), params.GetName());
3459  }
3460  }
3461  }
3462 
3463  if (!params.fRDHFCuts) {
3464  if (params.fMCMode == kMCTruth) {
3465  ::Warning("AliAnalysisTaskDmesonJets::ExecOnce",
3466  "%s: RDHF cuts not provided for engine '%s'.",
3467  GetName(), params.GetName());
3468  }
3469  else {
3470  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3471  "%s: RDHF cuts not provided. Engine '%s' disabled.",
3472  GetName(), params.GetName());
3473  params.fInhibit = kTRUE;
3474  }
3475  }
3476 
3477  params.fMCContainer = fMCContainer;
3478 
3479  for (auto cont_it : fParticleCollArray) {
3480  AliTrackContainer* track_cont = dynamic_cast<AliTrackContainer*>(cont_it.second);
3481  if (track_cont) params.fTrackContainers.push_back(track_cont);
3482  }
3483 
3484  for (auto cont_it : fClusterCollArray) params.fClusterContainers.push_back(cont_it.second);
3485 
3486  if (fAodEvent) params.Init(fGeom, fAodEvent->GetRunNumber());
3487 
3488  if (params.fMCMode != kMCTruth && fAodEvent) {
3489  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
3490 
3491  if (params.fCandidateArray) {
3492  TString className;
3493  if (params.fCandidateType == kD0toKpi || params.fCandidateType == kD0toKpiLikeSign) {
3494  className = "AliAODRecoDecayHF2Prong";
3495  }
3496  else if (params.fCandidateType == kDstartoKpipi) {
3497  className = "AliAODRecoCascadeHF";
3498  }
3499  if (!params.fCandidateArray->GetClass()->InheritsFrom(className)) {
3500  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3501  "%s: Objects of type %s in %s are not inherited from %s! Task will be disabled!",
3502  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName(), className.Data()); // @suppress("Ambiguous problem")
3503  params.fCandidateArray = 0;
3504  params.fInhibit = kTRUE;
3505  }
3506  }
3507  else {
3508  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3509  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
3510  params.fBranchName.Data(), params.GetName());
3511  params.fInhibit = kTRUE;
3512  }
3513  }
3514 
3515  if (params.fMCMode != kNoMC) {
3516  if (!params.fMCContainer) {
3517  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3518  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
3519  params.GetName());
3520  params.fInhibit = kTRUE;
3521  }
3522  }
3523 
3524  if (params.fMCMode != kMCTruth) {
3525  if (params.fTrackContainers.size() == 0 && params.fClusterContainers.size() == 0) {
3526  ::Error("AliAnalysisTaskDmesonJets::ExecOnce",
3527  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
3528  params.GetName());
3529  params.fInhibit = kTRUE;
3530  }
3531  }
3532  }
3533 }
3534 
3539 {
3540  TString hname;
3541 
3542  // Fix for temporary bug in ESDfilter
3543  // The AODs with null vertex pointer didn't pass the PhysSel
3544  // Now adding an entry in the histogram so as to check that this is actually cutting anything out
3545  if (fAodEvent && (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001)) {
3546  for (auto &eng : fAnalysisEngines) {
3547  if (eng.fInhibit) continue;
3548  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
3549  fHistManager.FillTH1(hname, "ESDfilterBug");
3550  }
3551  return kFALSE;
3552  }
3553 
3554  if (fAodEvent) {
3555  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
3556  if (matchingAODdeltaAODlevel <= 0) {
3557  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
3558  for (auto &eng : fAnalysisEngines) {
3559  if (eng.fInhibit) continue;
3560  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
3561  fHistManager.FillTH1(hname, "MismatchDeltaAOD");
3562  }
3563  return kFALSE;
3564  }
3565  }
3566 
3567  for (auto &eng : fAnalysisEngines) {
3568  eng.fDmesonJets.clear();
3569  if (eng.fInhibit) continue;
3570 
3571  eng.fEventInfo = EventInfo(fCent, fEPV0, fEventWeight, fPtHard);
3572 
3573  //Event selection
3574  hname = TString::Format("%s/fHistNEvents", eng.GetName());
3575  if (fAodEvent) {
3576  Bool_t iseventselected = kTRUE;
3577  if (eng.fRDHFCuts) iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
3578  if (!iseventselected) {
3579  fHistManager.FillTH1(hname, "Rejected");
3580  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
3581  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
3582  TString label;
3583  do {
3584  label = GetHFEventRejectionReasonLabel(bitmap);
3585  if (label.IsNull()) break;
3586  fHistManager.FillTH1(hname, label);
3587  } while (true);
3588 
3589  continue;
3590  }
3591  }
3592 
3593  fHistManager.FillTH1(hname, "Accepted");
3594 
3595  AliDebug(2, "Event selected");
3596 
3597  eng.RunAnalysis();
3598  }
3599  return kTRUE;
3600 }
3601 
3606 {
3607  for (auto &param : fAnalysisEngines) {
3608  if (param.fInhibit) continue;
3609  param.fOutputHandler->FillOutput(fApplyKinematicCuts);
3610  PostDataFromAnalysisEngine(param.fOutputHandler);
3611  }
3613  return kTRUE;
3614 }
3615 
3618 {
3619  auto itcont = fMCContainer->all_momentum();
3620  Int_t nHQ = 0;
3621  Double_t highestPartonPt = 0;
3622  Int_t PdgHighParton = 0;
3623  for (auto it = itcont.begin(); it != itcont.end(); it++) {
3624  auto part = *it;
3625  if (part.first.Pt() == 0) continue;
3626 
3627  Int_t PdgCode = part.second->GetPdgCode();
3628 
3629  // Skip all particles that are not either quarks or gluons
3630  if ((PdgCode < -9 || PdgCode > 9) && PdgCode != 21 && PdgCode != -21) continue;
3631 
3632  AliDebugStream(5) << "Parton " << it.current_index() <<
3633  " with pdg=" << PdgCode <<
3634  ", px=" << part.first.Px() <<
3635  ", py=" << part.first.Py() <<
3636  ", pz=" << part.first.Pz() <<
3637  ", n daughters = " << part.second->GetNDaughters() <<
3638  std::endl;
3639 
3640  // Skip partons that do not have any children
3641  // Unclear how this can happen, it would seem that this parton were not fragmented by the generator
3642  if (part.second->GetNDaughters() == 0) continue;
3643 
3644  // Look for highest momentum parton
3645  if (highestPartonPt < part.first.Pt()) {
3646  highestPartonPt = part.first.Pt();
3647  PdgHighParton = PdgCode;
3648  }
3649 
3650  // Skip partons that are not HF
3651  if (PdgCode != 4 && PdgCode != 5 && PdgCode != -4 && PdgCode != -5) continue;
3652 
3653  Bool_t lastInPartonShower = kTRUE;
3654  Bool_t hadronDaughter = kFALSE;
3655  for (Int_t daughterIndex = part.second->GetFirstDaughter(); daughterIndex <= part.second->GetLastDaughter(); daughterIndex++){
3656  if (daughterIndex < 0) {
3657  AliDebugStream(5) << "Could not find daughter index!" << std::endl;
3658  continue;
3659  }
3660  AliAODMCParticle *daughter = fMCContainer->GetMCParticle(daughterIndex);
3661  if (!daughter) {
3662  AliDebugStream(5) << "Could not find particle with index " << daughterIndex << "!" << std::endl;
3663  continue;
3664  }
3665  Int_t daughterPdgCode = daughter->GetPdgCode();
3666  if (daughter->GetMother() != it.current_index()) {
3667  AliDebugStream(5) << "Particle " << daughterIndex << " with pdg=" << daughterPdgCode <<
3668  ", px=" << daughter->Px() <<
3669  ", py=" << daughter->Py() <<
3670  ", pz=" << daughter->Pz() <<
3671  ", is not a daughter of " << it.current_index() <<
3672  "!" << std::endl;
3673  continue;
3674  }
3675 
3676  AliDebugStream(5) << "Found daughter " << daughterIndex <<
3677  " with pdg=" << daughterPdgCode <<
3678  ", px=" << daughter->Px() <<
3679  ", py=" << daughter->Py() <<
3680  ", pz=" << daughter->Pz() <<
3681  std::endl;
3682  // Codes between 81 and 100 are for internal MC code use, they may be intermediate states used e.g. in hadronizaion models
3683  if (daughterPdgCode == PdgCode) lastInPartonShower = kFALSE; // this parton is not the last parton in the shower
3684  if (TMath::Abs(daughterPdgCode) >= 111 || (daughterPdgCode >= 81 && daughterPdgCode <= 100)) hadronDaughter = kTRUE;
3685  }
3686  if (hadronDaughter) {
3687  AliDebugStream(5) << "This particle has at least a hadron among its daughters!" << std::endl;
3688  if (!lastInPartonShower) AliDebugStream(2) << "Odly, quark " << it.current_index() << " with PDG " << PdgCode << " (pt = " << part.first.Pt() << ", eta = " << part.first.Eta() << ") is not the last in the parton shower but at least a hadron found among its daughters?!" << std::endl;
3689  }
3690  else {
3691  AliDebugStream(5) << "This particle does not have hadrons among its daughters!" << std::endl;
3692  if (lastInPartonShower) AliDebugStream(2) << "Odly, quark " << it.current_index() << " with PDG " << PdgCode << " (pt = " << part.first.Pt() << ", eta = " << part.first.Eta() << ") is the last in the parton shower but no hadron found among its daughters?!" << std::endl;
3693  continue;
3694  }
3695 
3696  if (PdgCode == 4 || PdgCode == -4) {
3697  fHistManager.FillTH1("fHistCharmPt", part.first.Pt());
3698  fHistManager.FillTH1("fHistCharmEta", part.first.Eta());
3699  fHistManager.FillTH1("fHistCharmPhi", part.first.Phi_0_2pi());
3700  }
3701  else if (PdgCode == 5 || PdgCode == -5) {
3702  fHistManager.FillTH1("fHistBottomPt", part.first.Pt());
3703  fHistManager.FillTH1("fHistBottomEta", part.first.Eta());
3704  fHistManager.FillTH1("fHistBottomPhi", part.first.Phi_0_2pi());
3705  }
3706  nHQ++;
3707  }
3708  fHistManager.FillTH1("fHistNHeavyQuarks", nHQ);
3709  fHistManager.FillTH1("fHistHighestPartonPt",highestPartonPt);
3710  Int_t partonType = 0;
3711  Int_t absPdgHighParton = TMath::Abs(PdgHighParton);
3712  if (absPdgHighParton == 9 || absPdgHighParton == 21) {
3713  partonType = 7;
3714  }
3715  else {
3716  partonType = absPdgHighParton;
3717  }
3718  fHistManager.FillTH1("fHistHighestPartonType",partonType);
3719 }
3720 
3729 {
3730  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
3731 
3732  Double_t mass = part->Mass();
3733 
3734  // To make sure that the PDG mass value is not at the edge of a bin
3735  if (nbins % 2 == 0) {
3736  minMass = mass - range / 2 - range / nbins / 2;
3737  maxMass = mass + range / 2 - range / nbins / 2;
3738  }
3739  else {
3740  minMass = mass - range / 2;
3741  maxMass = mass + range / 2;
3742  }
3743 }
3744 
3751 {
3752  static TString label;
3753  label = "";
3754 
3755  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
3756  label = "NotSelTrigger";
3757  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
3758  return label.Data();
3759  }
3760  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
3761  label = "NoVertex";
3762  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
3763  return label.Data();
3764  }
3765  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
3766  label = "TooFewVtxContrib";
3767  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
3768  return label.Data();
3769  }
3770  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
3771  label = "ZVtxOutFid";
3772  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
3773  return label.Data();
3774  }
3775  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
3776  label = "Pileup";
3777  bitmap &= ~BIT(AliRDHFCuts::kPileup);
3778  return label.Data();
3779  }
3780  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
3781  label = "OutsideCentrality";
3782  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
3783  return label.Data();
3784  }
3785  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
3786  label = "PhysicsSelection";
3787  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
3788  return label.Data();
3789  }
3790  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
3791  label = "BadSPDVertex";
3792  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
3793  return label.Data();
3794  }
3795  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
3796  label = "ZVtxSPDOutFid";
3797  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
3798  return label.Data();
3799  }
3800  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
3801  label = "CentralityFlattening";
3802  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
3803  return label.Data();
3804  }
3805  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
3806  label = "BadTrackV0Correl";
3807  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
3808  return label.Data();
3809  }
3810 
3811  return label.Data();
3812 }
3813 
3820 {
3821  if (handler->GetDataSlotNumber() >= 0 && handler->GetOutputObject()) {
3822  PostData(handler->GetDataSlotNumber(), handler->GetOutputObject());
3823  return handler->GetDataSlotNumber();
3824  }
3825  else {
3826  return -1;
3827  }
3828 }
3829 
3839 {
3840  // Get the pointer to the existing analysis manager via the static access method
3841  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
3842  if (!mgr) {
3843  ::Error("AddTaskDmesonJets", "No analysis manager to connect to.");
3844  return NULL;
3845  }
3846 
3847  // Check the analysis type using the event handlers connected to the analysis manager
3848  AliVEventHandler* handler = mgr->GetInputEventHandler();
3849  if (!handler) {
3850  ::Error("AddTaskDmesonJets", "This task requires an input event handler");
3851  return NULL;
3852  }
3853 
3855 
3856  if (handler->InheritsFrom("AliESDInputHandler")) {
3857  dataType = kESD;
3858  }
3859  else if (handler->InheritsFrom("AliAODInputHandler")) {
3860  dataType = kAOD;
3861  }
3862 
3863  // Init the task and do settings
3864  if (ntracks == "usedefault") {
3865  if (dataType == kESD) {
3866  ntracks = "Tracks";
3867  }
3868  else if (dataType == kAOD) {
3869  ntracks = "tracks";
3870  }
3871  else {
3872  ntracks = "";
3873  }
3874  }
3875 
3876  if (nclusters == "usedefault") {
3877  if (dataType == kESD) {
3878  nclusters = "CaloClusters";
3879  }
3880  else if (dataType == kAOD) {
3881  nclusters = "caloClusters";
3882  }
3883  else {
3884  nclusters = "";
3885  }
3886  }
3887 
3888  if (nMCpart == "usedefault") {
3889  nMCpart = "mcparticles"; // Always needs AliAODMCParticle objects
3890  }
3891 
3892  TString name("AliAnalysisTaskDmesonJets");
3893  if (strcmp(suffix, "") != 0) {
3894  name += TString::Format("_%s", suffix.Data());
3895  }
3896 
3897  AliAnalysisTaskDmesonJets* jetTask = new AliAnalysisTaskDmesonJets(name, nMaxTrees);
3898 
3899  if (!ntracks.IsNull()) {
3900  AliHFTrackContainer* trackCont = new AliHFTrackContainer(ntracks);
3901  jetTask->AdoptParticleContainer(trackCont);
3902  }
3903 
3904  if (!nMCpart.IsNull()) {
3905  AliMCParticleContainer* partCont = new AliHFAODMCParticleContainer(nMCpart);
3906  partCont->SetEtaLimits(-1.5, 1.5);
3907  partCont->SetPtLimits(0, 1000);
3908  jetTask->AdoptParticleContainer(partCont);
3909  }
3910 
3911  jetTask->AddClusterContainer(nclusters.Data());
3912 
3913  // Final settings, pass to manager and set the containers
3914  mgr->AddTask(jetTask);
3915 
3916  // Create containers for input/output
3917  AliAnalysisDataContainer* cinput1 = mgr->GetCommonInputContainer();
3918  TString contname1(name);
3919  contname1 += "_histos";
3920  AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(contname1.Data(),
3921  TList::Class(), AliAnalysisManager::kOutputContainer,
3922  Form("%s", AliAnalysisManager::GetCommonFileName()));
3923 
3924  mgr->ConnectInput(jetTask, 0, cinput1);
3925  mgr->ConnectOutput(jetTask, 1, coutput1);
3926 
3927  for (Int_t i = 0; i < nMaxTrees; i++) {
3928  TString contname = TString::Format("%s_tree_%d", name.Data(), i);
3929  AliAnalysisDataContainer *coutput = mgr->CreateContainer(contname.Data(),
3930  TTree::Class(),AliAnalysisManager::kOutputContainer,
3931  Form("%s", AliAnalysisManager::GetCommonFileName()));
3932  mgr->ConnectOutput(jetTask, 2+i, coutput);
3933  }
3934  return jetTask;
3935 }
Bool_t IsAnyJetInAcceptance(const AliDmesonJetInfo &dMesonJet) const
Double32_t fCorrZ
Z of the D meson after subtracting average background.
Lightweight class that encapsulates D meson jets.
Int_t pdg
Double_t fMaxMass
! Max mass in histogram axis
void Print() const
Prints the content of this object in the standard output.
AliAODMCParticle * fParton
! pointer to the parton in the shower tree of the D meson (only for particle level D mesons) ...
std::list< AnalysisEngine > fAnalysisEngines
Array of analysis parameters.
std::string fClassName
Class name where the event was not found.
AliEMCALGeometry * fGeom
!emcal geometry
void SetDMesonCandidate(AliAODRecoDecay *c)
Analysis task for D meson jets.
UInt_t fRejectedOrigin
Bit mask with D meson origins that are rejected (used for MC analysis, i.e. signal-only, background-only and particle-level)
virtual void Set(const AliDmesonJetInfo &source)
Invariant mass of the D0 meson candidate in GeV/c2.
Bool_t fInhibit
!inhibit execution of the task
virtual Int_t Run()
double Double_t
Definition: External.C:58
Bool_t IsJetInAcceptance(const AliJetInfo &jet) const
Double_t fSoftPionPt
! Transverse momentum of the soft pion of the D* candidate
EDataType_t
Switch for the data type.
virtual AliAODMCParticle * GetMCParticle(Int_t i=-1) const
Int_t fDataSlotNumber
! Data slot where the tree output is posted
void Getd0MeasMinusExpProng(Int_t ip, Double_t magf, Double_t &d0diff, Double_t &errd0diff) const
const char * title
Definition: MakeQAPdf.C:27
AliRDHFCuts * fRDHFCuts
D meson candidates cuts.
virtual void Set(const AliDmesonJetInfo &source)
TString fRhoName
Name of the object that holds the average background value.
Lightweight class that encapsulates D meson jets.
std::vector< AliTrackContainer * > fTrackContainers
! Track containers
AnalysisEngine * AddAnalysisEngine(ECandidateType_t type, TString cutfname, TString cuttype, EMCMode_t bkgMode, EJetType_t jettype, Double_t jetradius, TString rhoName="")
AliJetInfoSummary ** fCurrentJetInfo
! Current jet info
Double_t fTrackEfficiency
Artificial tracking inefficiency (0...1)
std::vector< AliClusterContainer * > fClusterContainers
! Cluster containers
EMCMode_t fMCMode
! MC mode: No MC (data and MC detector level), background-only (MC), signal-only (MC), MC truth (particle level)
Double_t fEPV0
!event plane V0
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
AliHFAODMCParticleContainer * fMCContainer
! MC particle container
ERecoScheme_t fRecoScheme
Jet recombination scheme (pt scheme, E scheme, ...)
const TObjArray & GetDaughterList() const
Double32_t fWeight
Centrality of the collision.
Double_t fMaxNeutralPt
Transverse momentum of the leading neutral particle (or cluster)
Class that encapsulates event properties in a very compact structure (useful for pp simulation analys...
void RunDetectorLevelAnalysis()
Run a detector level analysis.
Double_t mass
Container with name, TClonesArray and cuts for particles.
Double32_t fPartonPt
Transverse momentum of the parton.
friend bool operator<(const AnalysisEngine &lhs, const AnalysisEngine &rhs)
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
virtual void UserCreateOutputObjects()
Creates the output containers.
void FillPartonLevelHistograms()
Fill histograms with parton-level information.
static Int_t CheckMatchingAODdeltaAODevents()
Double32_t fInvMass
Invariant mass of the D0 meson candidate in GeV/c2.
Lightweight class that encapsulates D meson jets.
Double_t fMinMass
! Min mass in histogram axis
Double_t fInvMass2Prong
! 2-prong mass of the D* candidate (w/o the soft pion)
friend bool operator==(const AnalysisEngine &lhs, const AnalysisEngine &rhs)
Double_t fMinChargedPt
Minimum pt of the leading charged particle (or track)
virtual Bool_t FillOutput(Bool_t applyKinCuts)
UInt_t fEnabledAxis
! Use bit defined in EAxis_t to enable axis in the THnSparse
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.
Double32_t fPtK
Transverse momentum of the kaon.
TRandom * fRandomGen
! Random number generator
std::vector< AliAnalysisTaskDmesonJets::AliHFJetDefinition > & GetJetDefinitions()
Double32_t fCosPointing
Cosine of the pointing angle.
AliRhoParameter * fRho
Object that holds the average background value.
EMCMode_t fMCMode
MC mode: No MC (data and MC detector level), background-only (MC), signal-only (MC), MC truth (particle level)
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
Double_t fJetGhostArea
Area of the ghost particles.
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:122
AnalysisEngine()
This is the default constructor, used for ROOT I/O purposes.
Output handler for D meson jet analysis.
Int_t PostDataFromAnalysisEngine(OutputHandler const *handler)
Double32_t fPtPi
Transverse momentum of the pion.
void Init(const AliEMCALGeometry *const geom, Int_t runNumber)
Initialize the analysis engine.
Int_t fDataSlotNumber
! Data slot where the tree output is posted
std::vector< AliHFJetDefinition > * fJetDefinitions
! Jet definitions
static FJRecoScheme ConvertToFJRecoScheme(ERecoScheme_t reco)
AliVParticle * fDmesonParticle
! pointer to the particle object
static TString GenerateJetName(EJetType_t jetType, EJetAlgo_t jetAlgo, ERecoScheme_t recoScheme, Double_t radius, AliParticleContainer *partCont, AliClusterContainer *clusCont, TString tag)
static std::pair< AliAnalysisTaskDmesonJets::EMesonOrigin_t, AliAODMCParticle * > IsPromptCharm(const AliAODMCParticle *part, TClonesArray *mcArray)
std::vector< DMESONTYPE > fCurrentDmesonInfo
! Current D meson jet info
virtual void Set(const AliDmesonJetInfo &source)
Double32_t fN
Number of jet constituents.
void SetR(Double_t r)
Definition: AliFJWrapper.h:127
Select tracks based on specific prescriptions of HF analysis.
Double32_t fPt
Transverse momentum of the jet in GeV/c.
ECandidateType_t fCandidateType
Candidate type.
Double_t GetDistance(std::string n, Double_t &deta, Double_t &dphi) const
Int_t fMCLabel
! MC label, i.e. index of the generator level D meson (only for detector level D meson candidates) ...
EBeamType_t fForceBeamType
forced beam type
const Int_t nPtBins
Double32_t fR
Distance between D meson and jet axis.
virtual void Set(const AliDmesonJetInfo &source)
void AddInputVectors(AliEmcalContainer *cont, Int_t offset, TH2 *rejectHist=0, Double_t eff=0.)
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:34
EventInfo * fEventInfo
! Object conatining the event information (centrality, pt hard, weight, etc.)
EJetType_t fJetType
Jet type (charged, full, neutral)
Double32_t fPt
Transverse momentum of the D meson in GeV/c.
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
Create a new TH2 within the container.
virtual void Set(const AliDmesonJetInfo &source)
AliTLorentzVector fD
! 4-momentum of the D meson candidate
Double_t GetDistance(const AliJetInfo &jet, Double_t &deta, Double_t &dphi) const
TObject * FindObject(const char *name) const
Find an object inside the container.
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
Double_t CosThetaStarD0bar() const
angle of K
AliAnalysisTaskDmesonJets()
This is the default constructor, used for ROOT I/O purposes.
virtual void Reset()
Reset the current object.
Output handler for D meson jet analysis.
AnalysisEngine & operator=(const AnalysisEngine &source)
Lightweight class that encapsulates D meson jets for PbPb analysis.
Struct that encapsulates analysis parameters.
std::vector< AliHFJetDefinition >::iterator FindJetDefinition(const AliHFJetDefinition &eng)
Base task in the EMCAL framework (lighter version of AliAnalysisTaskEmcal)
Double_t Phi_0_2pi() const
static AliAODMCParticle * FindParticleOrigin(const AliAODMCParticle *part, TClonesArray *mcArray, EFindParticleOriginMode_t mode, const std::set< UInt_t > &pdgSet)
Double_t fMinNeutralPt
Minimum pt of the leading neutral particle (or cluster)
Bool_t fD0D0bar
! kTRUE if selected both as D0 and D0bar
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:121
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
Bool_t fD0Extended
! Store extended information in the tree (only for D0 mesons)
AliAODTrack * GetBachelor() const
Bool_t fD0Extended
Store extended information in the tree (only for D0 mesons)
AliRDHFCuts * LoadDMesonCutsFromFile(TString cutfname, TString cutsname)
std::list< AnalysisEngine >::iterator FindAnalysisEngine(const AnalysisEngine &eng)
std::vector< AliHFJetDefinition > fJetDefinitions
Jet definitions.
Double_t fMaxNeutralPt
Maximum pt of the leading neutral particle (or cluster)
Double_t fMaxMass
Max mass in histogram axis.
EJetAlgo_t fJetAlgo
Jet algorithm (kt, anti-kt,...)
static const char * GetHFEventRejectionReasonLabel(UInt_t &bitmap)
Lightweight class that encapsulates D*.
virtual void Clear(const Option_t *="")
friend bool operator<(const AliHFJetDefinition &lhs, const AliHFJetDefinition &rhs)
Int_t mode
Definition: anaM.C:41
Bool_t ExtractDstarAttributes(const AliAODRecoCascadeHF *DstarCand, AliDmesonJetInfo &DmesonJet, UInt_t i)
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
const AliJetInfo * GetJet(std::string n) const
Short_t fPartonType
! type of the parton in the shower tree (only for particle level D mesons)
TClonesArray * fCandidateArray
! D meson candidate array
const AliMCParticleIterableMomentumContainer all_momentum() const
virtual ~AliAnalysisTaskDmesonJets()
This is the standard destructor.
Byte_t fSelectionType
! for D0: 0=not selected, 1=D0, 2=D0bar
Bool_t ExtractD0Attributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
Select MC particles based on specific prescriptions of HF analysis.
std::map< int, AliDmesonJetInfo > * fDmesonJets
! Array containing the D meson jets
static EMesonDecayChannel_t CheckDecayChannel(const AliAODMCParticle *part, TClonesArray *mcArray)
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:125
AliDmesonJetInfo & operator=(const AliDmesonJetInfo &source)
Double_t fMaxChargedPt
Transverse momentum of the leading charged particle (or track)
Double32_t fCorrPt
Transverse momentum of the jet in GeV/c after subtracting average background.
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
EventInfo fEventInfo
! Event info (centrality, weight, pt hard etc.)
friend bool operator==(const AliHFJetDefinition &lhs, const AliHFJetDefinition &rhs)
static FJJetAlgo ConvertToFJAlgo(EJetAlgo_t algo)
Double_t fMaxChargedPt
Maximum pt of the leading charged particle (or track)
ECandidateType_t fCandidateType
! Candidate type
decay
Definition: HFPtSpectrum.C:41
std::map< std::string, AliJetInfo > fJets
! list of jets
Look for the very first particle in the fragmentation tree.
AliAODMCParticle * fAncestor
! pointer to the ancestor particle in the shower tree of the D meson (only for particle level D meson...
Double32_t fDeltaInvMass
< Difference between the Kpipi and the Kpi invariant masses in GeV/c2
OutputHandler * fOutputHandler
! Output handler
Definition: External.C:220
Double_t minMass
Bool_t fRejectISR
! Reject initial state radiation
void RunParticleLevelAnalysis()
Run a particle level analysis.
Bool_t FillHnSparse(THnSparse *h, const AliDmesonJetInfo &DmesonJet, std::string n)
std::map< std::string, AliParticleContainer * > fParticleCollArray
particle/track collection array
std::map< int, AliDmesonJetInfo > fDmesonJets
! Array containing the D meson jets
Double_t fTrackEfficiency
! Artificial tracking inefficiency (0...1) -> set automatically at ExecOnce by AliAnalysisTaskDmesonJ...
Bool_t fApplyKinematicCuts
Apply jet kinematic cuts.
Double32_t fSelectionType
Selection type: D0, D0bar, both.
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:292
const char * GetName() const
Generate a name for this jet definition.
AliEventNotFound(const std::string &class_name, const std::string &method_name)
UInt_t fAcceptedDecay
Bit mask with D meson decays that are accepted (only used for particle-level analysis) ...
virtual void Set(const AliDmesonJetInfo &source, std::string n)
unsigned short UShort_t
Definition: External.C:28
Bool_t fReconstructed
! Whether this D meson was reconstructed (only for particle level D mesons)
Bool_t ExtractRecoDecayAttributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
Double_t InvMassD0() const
Bool_t fRejectISR
Reject initial state radiation.
AliDmesonInfoSummary * fCurrentDmesonJetInfo
! Current D meson jet info
EOutputType_t fOutputType
Output type: none, TTree or THnSparse.
virtual void AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
std::map< std::string, std::vector< JETTYPE > > fCurrentJetInfo
! Current jet info
Double_t InvMassDstarKpipi() const
const AliVEvent * fEvent
! pointer to the ESD/AOD event
TArrayI fPDGdaughters
List of the PDG code of the daughters.
const Int_t nbins
Double_t maxMass
TString fBranchName
AOD branch where the D meson candidate are found.
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
Int_t fNOutputTrees
Maximum number of output trees.
static OutputHandlerTTreeExtendedBase * GenerateOutputHandler(AnalysisEngine *eng)
AliHFJetDefinition()
This is the default constructor, used for ROOT I/O purposes.
AliHFJetDefinition & operator=(const AliHFJetDefinition &source)
Int_t GetNAcceptedParticles() const
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:123
std::map< std::string, AliClusterContainer * > fClusterCollArray
cluster collection array
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
Create a new THnSparse within the container.
Double_t fCorrPt
Transverse momentum of the jet after subtracting the average background.
AliHFAODMCParticleContainer * fMCContainer
! MC particle container
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:312
AliFJWrapper * fFastJetWrapper
! Fastjet wrapper
Float_t fPtBinWidth
! Histogram pt bin width
AliHFJetDefinition * AddJetDefinition(EJetType_t type, Double_t r, EJetAlgo_t algo, ERecoScheme_t reco)
Container for MC-true particles within the EMCAL framework.
Lightweight class that encapsulates D0.
Double_t fCent
!event centrality
void SetCharge(EChargeCut_t c)
Container for jet within the EMCAL jet framework.
Definition: External.C:196
virtual void BuildOutputObject(const char *taskName)
AliFJWrapper * fFastJetWrapper
! Fastjet wrapper
THistManager fHistManager
Histogram manager.
void Reset()
Reset all fields to their default values.
static AliAnalysisTaskDmesonJets * AddTaskDmesonJets(TString ntracks="usedefault", TString nclusters="usedefault", TString nMCpart="", Int_t nMaxTrees=2, TString suffix="")
virtual void Set(const AliDmesonJetInfo &source, std::string n)
Double_t fMinMass
Min mass in histogram axis.