AliPhysics  master (3d17d9d)
AliAnalysisTaskDmesonJetsSub.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 //copy paste the task of Salvatore with some additions for substructure
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 AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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("AliAnalysisTaskDmesonJetsSub::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("AliAnalysisTaskDmesonJetsSub::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("AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::AliEventNotFound("AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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("AliAnalysisTaskDmesonJetsSub::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("AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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 = "AliAnalysisTaskDmesonJetsSub::AliDmesonMCInfoSummary";
1179  }
1180  else {
1181  switch (fCandidateType) {
1182  case kD0toKpi:
1183  case kD0toKpiLikeSign:
1184  if (fD0Extended) {
1185  classname = "AliAnalysisTaskDmesonJetsSub::AliD0ExtendedInfoSummary";
1187  }
1188  else {
1189  classname = "AliAnalysisTaskDmesonJetsSub::AliD0InfoSummary";
1191  }
1192  break;
1193  case kDstartoKpipi:
1194  classname = "AliAnalysisTaskDmesonJetsSub::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(), "AliAnalysisTaskDmesonJetsSub::AliJetInfoSummary", &fCurrentJetInfo[i]);
1207  }
1208  else {
1210  fTree->Branch(fJetDefinitions->at(i).GetName(), "AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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 = "AliAnalysisTaskDmesonJetsSub::AliEventInfoSummary";
1372  TString d_meson_class_name;
1373 
1374  TString jet_class_name = "std::vector<AliAnalysisTaskDmesonJetsSub::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<AliAnalysisTaskDmesonJetsSub::AliJetInfoPbPbSummary>";
1380  }
1381  }
1382 
1383  OutputHandlerTTreeExtendedBase* result = nullptr;
1384  if (eng->GetMCMode() == kMCTruth) {
1385  d_meson_class_name = "std::vector<AliAnalysisTaskDmesonJetsSub::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<AliAnalysisTaskDmesonJetsSub::AliD0ExtendedInfoSummary>";
1399  if (RhoJet) {
1401  }
1402  else {
1404  }
1405  }
1406  else {
1407  d_meson_class_name = "AliAnalysisTaskDmesonJetsSub::AliD0InfoSummary";
1408  if (RhoJet) {
1410  }
1411  else {
1413  }
1414  }
1415  break;
1416  case kDstartoKpipi:
1417  d_meson_class_name = "AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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 AliAnalysisTaskDmesonJetsSub::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("AliAnalysisTaskDmesonJetsSub::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("AliAnalysisTaskDmesonJetsSub::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("AliAnalysisTaskDmesonJetsSub::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<AliAnalysisTaskDmesonJetsSub::AliHFJetDefinition>::iterator AliAnalysisTaskDmesonJetsSub::AnalysisEngine::FindJetDefinition(const AliAnalysisTaskDmesonJetsSub::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 
2223 {
2224  if (fCandidateType == kD0toKpi || fCandidateType == kD0toKpiLikeSign) { // D0 candidate
2225 
2226  return ExtractD0Efficiencies(Dcand, DmesonJet, jetDef, i);
2227  }
2228  else {
2229  return kFALSE;
2230  }
2231 }
2232 
2234 {
2235  AliDebug(10,"Checking if D0 meson is selected");
2236  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
2237  if (isSelected == 0) return kFALSE;
2238  TString hname1;
2239  TString hname2;
2240 
2241  Int_t MCtruthPdgCode = 0;
2242 
2243  Int_t myflag=0;
2244  Double_t jeteta=0;
2245  Double_t jetpt=0;
2246  hname1 = TString::Format("%s/EfficiencyMatchesPrompt", fName.Data());
2247  TH2* EfficiencyMatchesPrompt = static_cast<TH2*>(fHistManager->FindObject(hname1));
2248  AliAODMCParticle* aodMcPart;
2249 
2250 
2251  hname2 = TString::Format("%s/EfficiencyMatchesNonPrompt", fName.Data());
2252  TH2* EfficiencyMatchesNonPrompt = static_cast<TH2*>(fHistManager->FindObject(hname2));
2253 
2254  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
2255  // Checks also the origin, and if it matches the rejected origin mask, return false
2256  Double_t jetPtdet = DmesonJet.fJets[jetDef.GetName()].fMomentum.Pt();
2257  Double_t jetEtadet= DmesonJet.fJets[jetDef.GetName()].fMomentum.Eta();
2258  if(TMath::Abs(jetEtadet)>0.5) return kFALSE;
2259  if(jetPtdet>100) return kFALSE;
2260 
2261 
2262  if (fMCMode != kNoMC) {
2263  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
2264  DmesonJet.fMCLabel = mcLab;
2265 
2266  // Retrieve the generated particle (if exists) and its PDG code
2267  if (mcLab >= 0) {
2268  aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
2269 
2270  if (aodMcPart) {
2271  // Check origin and return false if it matches the rejected origin mask
2272  if (fRejectedOrigin) {
2273  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
2274  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
2275  }
2276  MCtruthPdgCode = aodMcPart->PdgCode();
2277  }
2278  }
2279  }
2280 
2281 
2282  if(fMCMode==kSignalOnly){
2283  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
2285  fFastJetWrapper->SetR(jetDef.fRadius);
2289  fFastJetWrapper->Run();
2290  std::vector<fastjet::PseudoJet> jets_incl = sorted_by_pt(fFastJetWrapper->GetInclusiveJets());
2291 
2292  for (auto jet : jets_incl) {
2293  std::vector<fastjet::PseudoJet> constituents = jet.constituents();
2294  // if(constituents.size()<2) continue;
2295  for (auto constituent : jet.constituents()) {
2296  Int_t iPart = constituent.user_index() - 100;
2297  if (constituent.perp() < 1e-6) continue; // reject ghost particles
2298  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
2299  if (!part) {
2300  ::Error("AliAnalysisTaskDmesonJetsSub::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
2301  continue;
2302  }
2303  if(part==aodMcPart){myflag=1;
2304  jetpt=jet.perp();
2305  jeteta=jet.eta();
2306  break;}
2307 
2308  }}
2309  if(myflag==0) return kFALSE;
2310  if(TMath::Abs(jeteta)>0.5) return kFALSE;
2311 
2312  Double_t maxFiducialY=-0.2/15*aodMcPart->Pt()*aodMcPart->Pt()+1.9/15*aodMcPart->Pt()+0.5;
2313  Double_t minFiducialY = 0.2/15*aodMcPart->Pt()*aodMcPart->Pt()-1.9/15*aodMcPart->Pt()-0.5;
2314 
2315  if (isSelected == 1) { // selected as a D0
2316  if (i != 0) return kFALSE; // only one mass hypothesis thanks to PID
2317 
2318  if(MCtruthPdgCode == fCandidatePDG){
2319  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
2320 
2321  if(origin.first == kFromCharm) {
2322  if(aodMcPart->Pt()<=5) if(aodMcPart->Y()>=minFiducialY && aodMcPart->Y()<=maxFiducialY) EfficiencyMatchesPrompt->Fill(aodMcPart->Pt(),jetpt);
2323  if(aodMcPart->Pt()>5) if(TMath::Abs(aodMcPart->Y())<0.8)EfficiencyMatchesPrompt->Fill(aodMcPart->Pt(),jetpt);}
2324  if(origin.first == kFromBottom) {
2325  if(aodMcPart->Pt()<=5) if(aodMcPart->Y()>=minFiducialY && aodMcPart->Y()<=maxFiducialY)EfficiencyMatchesNonPrompt->Fill(aodMcPart->Pt(),jetpt);
2326  if(aodMcPart->Pt()>5) if(TMath::Abs(aodMcPart->Y())<0.8)EfficiencyMatchesNonPrompt->Fill(aodMcPart->Pt(),jetpt);}
2327  }
2328 
2329  }
2330  else if (isSelected == 2) { // selected as a D0bar
2331  if (i != 1) return kFALSE; // only one mass hypothesis thanks to PID
2332 
2333 
2334  if(MCtruthPdgCode == -fCandidatePDG){
2335  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
2336  if(origin.first == kFromCharm){
2337 
2338  if(aodMcPart->Pt()<=5) if(aodMcPart->Y()>=minFiducialY && aodMcPart->Y()<=maxFiducialY)EfficiencyMatchesPrompt->Fill(aodMcPart->Pt(),jetpt);
2339  if(aodMcPart->Pt()>5) if(TMath::Abs(aodMcPart->Y())<0.8)EfficiencyMatchesPrompt->Fill(aodMcPart->Pt(),jetpt);}
2340  if(origin.first == kFromBottom){
2341  if(aodMcPart->Pt()<=5) if(aodMcPart->Y()>=minFiducialY && aodMcPart->Y()<=maxFiducialY)EfficiencyMatchesNonPrompt->Fill(aodMcPart->Pt(),jetpt);
2342  if(aodMcPart->Pt()>5) if(TMath::Abs(aodMcPart->Y())<0.8)EfficiencyMatchesNonPrompt->Fill(aodMcPart->Pt(),jetpt);
2343  }
2344  }
2345 
2346  }
2347 
2348  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
2349 
2350 
2351  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
2352  if (MCtruthPdgCode == fCandidatePDG){
2353 
2354  if (i == 0){
2355  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
2356  if(origin.first == kFromCharm){
2357  if(aodMcPart->Pt()<=5) if(aodMcPart->Y()>=minFiducialY && aodMcPart->Y()<=maxFiducialY)EfficiencyMatchesPrompt->Fill(aodMcPart->Pt(),jetpt);
2358  if(aodMcPart->Pt()>5) if(TMath::Abs(aodMcPart->Y())<0.8)EfficiencyMatchesPrompt->Fill(aodMcPart->Pt(),jetpt);
2359  }
2360  if(origin.first == kFromBottom){
2361  if(aodMcPart->Pt()<=5) if(aodMcPart->Y()>=minFiducialY && aodMcPart->Y()<=maxFiducialY)EfficiencyMatchesNonPrompt->Fill(aodMcPart->Pt(),jetpt);
2362  if(aodMcPart->Pt()>5) if(TMath::Abs(aodMcPart->Y())<0.8)EfficiencyMatchesNonPrompt->Fill(aodMcPart->Pt(),jetpt);
2363 
2364 
2365  }}
2366  }
2367  else if (MCtruthPdgCode == -fCandidatePDG){
2368  if (i == 1){
2369 
2370  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
2371  if(origin.first == kFromCharm){
2372  if(aodMcPart->Pt()<=5) if(aodMcPart->Y()>=minFiducialY && aodMcPart->Y()<=maxFiducialY)EfficiencyMatchesPrompt->Fill(aodMcPart->Pt(),jetpt);
2373  if(aodMcPart->Pt()>5) if(TMath::Abs(aodMcPart->Y())<0.8)EfficiencyMatchesPrompt->Fill(aodMcPart->Pt(),jetpt);}
2374  if(origin.first == kFromBottom){
2375  if(aodMcPart->Pt()<=5) if(aodMcPart->Y()>=minFiducialY && aodMcPart->Y()<=maxFiducialY)EfficiencyMatchesNonPrompt->Fill(aodMcPart->Pt(),jetpt);
2376  if(aodMcPart->Pt()>5) if(TMath::Abs(aodMcPart->Y())<0.8)EfficiencyMatchesNonPrompt->Fill(aodMcPart->Pt(),jetpt);}}
2377 
2378  }
2379  }
2380 
2381 
2382 
2383 
2384 
2385 
2386 
2387 
2388 
2389 
2390  }
2391  return kTRUE;
2392 }
2393 
2394 
2396 {
2397 
2398  TString hname;
2399  TString hname1;
2400  TString hname2;
2401 
2402  Double_t jeteta=0;
2403  Double_t jetpt=0;
2404  Int_t TheTrueCode = 0;
2405 
2411  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
2412 
2413  if (!fMCContainer->IsSpecialPDGFound()) return kFALSE;
2414 
2415  hname1 = TString::Format("%s/EfficiencyGeneratorPrompt", fName.Data());
2416  TH2* EfficiencyGeneratorPrompt = static_cast<TH2*>(fHistManager->FindObject(hname1));
2417 
2418 
2419 
2420  hname2 = TString::Format("%s/EfficiencyGeneratorNonPrompt", fName.Data());
2421  TH2* EfficiencyGeneratorNonPrompt = static_cast<TH2*>(fHistManager->FindObject(hname2));
2422  if(fMCMode==kSignalOnly){
2423 
2425  fFastJetWrapper->SetR(jetDef.fRadius);
2428  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
2429  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
2430 
2431  fFastJetWrapper->Run();
2432  std::vector<fastjet::PseudoJet> jets_incl = sorted_by_pt(fFastJetWrapper->GetInclusiveJets());
2433 
2434 
2435 
2436  for (auto jet : jets_incl) {
2437  jetpt=0;
2438  jeteta=0;
2439  std::vector<fastjet::PseudoJet> constituents = jet.constituents();
2440 
2441  for (auto constituent : jet.constituents()) {
2442  Int_t iPart = constituent.user_index() - 100;
2443  if (constituent.perp() < 1e-6) continue; // reject ghost particles
2444  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
2445  if (!part) {
2446  ::Error("AliAnalysisTaskDmesonJetsSub::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
2447  continue;
2448  }
2449 
2450  TheTrueCode=part->PdgCode();
2451  // Int_t mother=part->GetMother();
2452  // AliAODMCParticle* mypart = fMCContainer->GetMCParticle(mother);
2453  if(TMath::Abs(TheTrueCode)==fCandidatePDG){
2454 
2455  jetpt=jet.perp();
2456  jeteta=jet.eta();
2457  if(TMath::Abs(jeteta)>0.5) continue;
2458 
2459  auto origin = IsPromptCharm(part, fMCContainer->GetArray());
2460  Double_t maxFiducialY=-0.2/15*part->Pt()*part->Pt()+1.9/15*part->Pt()+0.5;
2461  Double_t minFiducialY = 0.2/15*part->Pt()*part->Pt()-1.9/15*part->Pt()-0.5;
2462  if(origin.first == kFromCharm){
2463  if(part->Pt()>5) if(TMath::Abs(part->Y())<0.8) EfficiencyGeneratorPrompt->Fill(part->Pt(),jetpt);
2464  if(part->Pt()<=5) if(part->Y()>=minFiducialY && part->Y()<=maxFiducialY) EfficiencyGeneratorPrompt->Fill(part->Pt(),jetpt);
2465  }
2466  if(origin.first == kFromBottom){
2467  if(part->Pt()>5)if(TMath::Abs(part->Y())<0.8) EfficiencyGeneratorNonPrompt->Fill(part->Pt(),jetpt);
2468  if(part->Pt()<=5) if(part->Y()>=minFiducialY && part->Y()<=maxFiducialY) EfficiencyGeneratorNonPrompt->Fill(part->Pt(),jetpt);
2469  }
2470 
2471  }
2472 
2473 
2474 
2475  }}
2476 
2477  }
2478 
2479  return kTRUE;
2480 }
2481 
2483 {
2484 
2485  TString hname;
2486  TString hname1;
2487  TString hname2;
2488 
2489  Double_t jeteta=0;
2490  Double_t jetpt=0;
2491  Int_t TheTrueCode = 0;
2492 
2493  vector<int> dlabel;
2500  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
2501  if (!fMCContainer->IsSpecialPDGFound()) return kFALSE;
2502 
2503 
2504  hname1 = TString::Format("%s/EfficiencyGeneratorPrompt", fName.Data());
2505  TH2* EfficiencyGeneratorPrompt = static_cast<TH2*>(fHistManager->FindObject(hname1));
2506 
2507 
2508 
2509  hname2 = TString::Format("%s/EfficiencyGeneratorNonPrompt", fName.Data());
2510  TH2* EfficiencyGeneratorNonPrompt = static_cast<TH2*>(fHistManager->FindObject(hname2));
2511  if(fMCMode==kSignalOnly){
2512 
2513  //here I loop over the container and count the D mesons and store their indexes
2514  auto cont = fMCContainer->all();
2515  for (auto it = cont.begin(); it != cont.end(); ++it) {
2516  UInt_t rejectionReason = 0;
2517  if((*it)->PdgCode()==fCandidatePDG){
2518 
2519  dlabel.push_back(it.current_index());}
2520 
2521 
2522  if (!fMCContainer->AcceptObject(it.current_index(), rejectionReason)) {
2523 
2524  continue;
2525  }
2526 
2527 
2528  }
2529 
2530  // then, for each D meson I replace only its decays (not other D decays) and I only keep for the jet finding the given D meson
2531  for(Int_t j=0;j<dlabel.size();j++){
2532 
2534  fMCContainer->SetSpecialIndex(dlabel[j]);
2535 
2537  fFastJetWrapper->SetR(jetDef.fRadius);
2540  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
2541  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
2542  fFastJetWrapper->Run();
2543  std::vector<fastjet::PseudoJet> jets_incl = sorted_by_pt(fFastJetWrapper->GetInclusiveJets());
2544 
2545 
2546 
2547  for (auto jet : jets_incl) {
2548  jetpt=0;
2549  jeteta=0;
2550  std::vector<fastjet::PseudoJet> constituents = jet.constituents();
2551 
2552  for (auto constituent : jet.constituents()) {
2553  Int_t iPart = constituent.user_index() - 100;
2554  if (constituent.perp() < 1e-6) continue; // reject ghost particles
2555  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
2556  if (!part) {
2557  ::Error("AliAnalysisTaskDmesonJetsSub::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
2558  continue;
2559  }
2560 
2561  TheTrueCode=part->PdgCode();
2562 
2563  if(TMath::Abs(TheTrueCode)==fCandidatePDG){
2564 
2565  jetpt=jet.perp();
2566  jeteta=jet.eta();
2567  if(TMath::Abs(jeteta)>0.5) continue;
2568 
2569  auto origin = IsPromptCharm(part, fMCContainer->GetArray());
2570  Double_t maxFiducialY=-0.2/15*part->Pt()*part->Pt()+1.9/15*part->Pt()+0.5;
2571  Double_t minFiducialY = 0.2/15*part->Pt()*part->Pt()-1.9/15*part->Pt()-0.5;
2572  if(origin.first == kFromCharm){
2573  if(part->Pt()>5) if(TMath::Abs(part->Y())<0.8) EfficiencyGeneratorPrompt->Fill(part->Pt(),jetpt);
2574  if(part->Pt()<=5) if(part->Y()>=minFiducialY && part->Y()<=maxFiducialY) EfficiencyGeneratorPrompt->Fill(part->Pt(),jetpt);}
2575  if(origin.first == kFromBottom){
2576  if(part->Pt()>5) if(TMath::Abs(part->Y())<0.8) EfficiencyGeneratorNonPrompt->Fill(part->Pt(),jetpt);
2577  if(part->Pt()<=5) if(part->Y()>=minFiducialY && part->Y()<=maxFiducialY) EfficiencyGeneratorNonPrompt->Fill(part->Pt(),jetpt);
2578 
2579  }
2580 
2581  }
2582 
2583 
2584 
2585  }}
2586 
2587  }
2588  dlabel.clear();
2589  }
2590 
2591  return kTRUE;
2592 }
2593 
2594 
2595 
2596 
2605 {
2606  AliDebug(10,"Checking if D0 meson is selected");
2607  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
2608  if (isSelected == 0) return kFALSE;
2609  TString hname;
2610  TString hname2;
2611  Int_t MCtruthPdgCode = 0;
2612 
2613  Double_t invMassD = 0;
2614 
2615  AliAODMCParticle* aodMcPart;
2616 
2617 
2618  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
2619  // Checks also the origin, and if it matches the rejected origin mask, return false
2620  if (fMCMode != kNoMC) {
2621  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
2622  DmesonJet.fMCLabel = mcLab;
2623 
2624  // Retrieve the generated particle (if exists) and its PDG code
2625  if (mcLab >= 0) {
2626  aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
2627 
2628  if (aodMcPart) {
2629  // Check origin and return false if it matches the rejected origin mask
2630  if (fRejectedOrigin) {
2631  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
2632  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
2633  }
2634  MCtruthPdgCode = aodMcPart->PdgCode();
2635  }
2636  }
2637  }
2638 
2639 
2640 
2641  if (isSelected == 1) { // selected as a D0
2642  if (i != 0) return kFALSE; // only one mass hypothesis thanks to PID
2643 
2644  if (fMCMode == kNoMC ||
2645  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
2646  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly) ||
2647  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kD0Reflection || fMCMode == kOnlyWrongPIDAccepted))) {
2648  // 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)
2649  AliDebug(10,"Selected as D0");
2650  invMassD = Dcand->InvMassD0();
2651  }
2652  else { // conditions above not passed, so return FALSE
2653  return kFALSE;
2654  }
2655 
2656 
2657 
2658 
2659  }
2660  else if (isSelected == 2) { // selected as a D0bar
2661  if (i != 1) return kFALSE; // only one mass hypothesis thanks to PID
2662 
2663  if (fMCMode == kNoMC ||
2664  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
2665  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly) ||
2666  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kD0Reflection || fMCMode == kOnlyWrongPIDAccepted))) {
2667  // 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)
2668  AliDebug(10,"Selected as D0bar");
2669  invMassD = Dcand->InvMassD0bar();
2670  }
2671  else { // conditions above not passed, so return FALSE
2672  return kFALSE;
2673  }
2674 
2675  }
2676  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
2677  AliDebug(10,"Selected as either D0 or D0bar");
2678 
2679  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
2680  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
2681  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kD0Reflection))) {
2682  if (i != 0) return kFALSE;
2683  AliDebug(10, "MC truth is D0");
2684  invMassD = Dcand->InvMassD0();
2685  }
2686  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
2687  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kD0Reflection))) {
2688  if (i != 1) return kFALSE;
2689  AliDebug(10, "MC truth is D0bar");
2690  invMassD = Dcand->InvMassD0bar();
2691  }
2692  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
2693 
2694  // Only accept it if background-only OR background-and-signal was requested
2695  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
2696  // Select D0 or D0bar depending on the i-parameter
2697  if (i == 0) {
2698  AliDebug(10, "Returning invariant mass with D0 hypothesis");
2699  invMassD = Dcand->InvMassD0();
2700  }
2701  else if (i == 1) {
2702  AliDebug(10, "Returning invariant mass with D0bar hypothesis");
2703  invMassD = Dcand->InvMassD0bar();
2704  }
2705  else { // i > 1
2706  return kFALSE;
2707  }
2708  }
2709  else { // signal-only was requested but this is not a true D0
2710  return kFALSE;
2711  }
2712  }
2713  }
2714  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
2715  return kTRUE;
2716 }
2717 
2726 {
2727  AliDebug(10,"Checking if D* meson is selected");
2728  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoCascadeHF*>(DstarCand), AliRDHFCuts::kAll, fAodEvent);
2729  if (isSelected == 0) return kFALSE;
2730 
2731  if ((i == 1 && DstarCand->Charge()>0) || (i == 0 && DstarCand->Charge()<0) || i > 1) return kFALSE; // only one mass hypothesis for the D*
2732 
2733  Int_t MCtruthPdgCode = 0;
2734 
2735  Double_t invMassD = 0;
2736 
2737  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
2738  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
2739  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
2740 
2741  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
2742  AliDebug(10, Form("MC label is %d", mcLab));
2743  DmesonJet.fMCLabel = mcLab;
2744  if (mcLab >= 0) {
2745  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
2746 
2747  if (aodMcPart) {
2748  if (fRejectedOrigin) {
2749  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
2750  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
2751  }
2752 
2753  MCtruthPdgCode = aodMcPart->PdgCode();
2754  AliDebug(10, Form("MC truth pdg code is %d",MCtruthPdgCode));
2755  }
2756  }
2757  }
2758 
2759  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
2760  if (fMCMode == kNoMC ||
2761  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
2762  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
2763  // 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)
2764  invMassD = DstarCand->InvMassDstarKpipi();
2765  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
2766  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
2767  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
2768  return kTRUE;
2769  }
2770  else { // conditions above not passed, so return FALSE
2771  return kFALSE;
2772  }
2773 }
2774 
2782 {
2783  if (!part) return kUnknownDecay;
2784  if (!mcArray) return kUnknownDecay;
2785 
2787 
2788  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
2789 
2790  if (part->GetNDaughters() == 2) {
2791 
2792  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughterLabel(0)));
2793  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughterLabel(1)));
2794 
2795  if (!d1 || !d2) {
2796  return decay;
2797  }
2798 
2799  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
2800  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
2801 
2802  if (absPdgPart == 421) { // D0 -> K pi
2803  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
2804  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
2805  decay = kDecayD0toKpi;
2806  }
2807  }
2808 
2809  if (absPdgPart == 413) { // D* -> D0 pi
2810  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
2811  Int_t D0decay = CheckDecayChannel(d1, mcArray);
2812  if (D0decay == kDecayD0toKpi) {
2813  decay = kDecayDStartoKpipi;
2814  }
2815  }
2816  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
2817  Int_t D0decay = CheckDecayChannel(d2, mcArray);
2818  if (D0decay == kDecayD0toKpi) {
2819  decay = kDecayDStartoKpipi;
2820  }
2821  }
2822  }
2823  }
2824 
2825  return decay;
2826 }
2827 
2834 std::pair<AliAnalysisTaskDmesonJetsSub::EMesonOrigin_t, AliAODMCParticle*> AliAnalysisTaskDmesonJetsSub::AnalysisEngine::IsPromptCharm(const AliAODMCParticle* part, TClonesArray* mcArray)
2835 {
2836  std::pair<AliAnalysisTaskDmesonJetsSub::EMesonOrigin_t, AliAODMCParticle*> result(kUnknownQuark, 0);
2837 
2838  if (!part) return result;
2839  if (!mcArray) return result;
2840 
2841  static std::set<UInt_t> partons = { 4, 5 };
2842 
2843  AliAODMCParticle* parton = FindParticleOrigin(part, mcArray, kFindLast, partons);
2844  if (parton) {
2845  result.second = parton;
2846  UInt_t absPdgParton = TMath::Abs(parton->GetPdgCode());
2847  if (absPdgParton == 4) result.first = kFromCharm;
2848  else if (absPdgParton == 5) result.first = kFromBottom;
2849  }
2850 
2851  return result;
2852 }
2853 
2861 
2862 AliAODMCParticle* AliAnalysisTaskDmesonJetsSub::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode)
2863 {
2864  static std::set<UInt_t> pdgSet;
2865 
2866  return FindParticleOrigin(part, mcArray, mode, pdgSet);
2867 }
2868 
2877 
2878 AliAODMCParticle* AliAnalysisTaskDmesonJetsSub::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode, const std::set<UInt_t>& pdgSet)
2879 {
2880  AliAODMCParticle* result = nullptr;
2881 
2882  Int_t mother = part->GetMother();
2883  while (mother >= 0) {
2884  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
2885  if (mcGranma) {
2886  UInt_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
2887 
2888  // If the current particle is one of the particle types that is being searched assign it to the result pointer
2889  if (pdgSet.empty() || pdgSet.count(abspdgGranma) > 0) {
2890  result = mcGranma;
2891  // If the last particle in the fragmentation tree (first when going reverse) was requested then stop the loop
2892  if (mode == kFindLast) break;
2893  }
2894  if (mother == mcGranma->GetMother()) { // avoid infinite loop!
2895  AliWarningClassStream() << "Particle " << mother << " (PDG=" << mcGranma->PdgCode() << ") is the mother of itself!?" << std::endl;
2896  break;
2897  }
2898  mother = mcGranma->GetMother();
2899  }
2900  else {
2901  AliErrorClassStream() << "Could not retrieve mother particle " << mother << "!" << std::endl;
2902  break;
2903  }
2904  }
2905 
2906  return result;
2907 }
2908 
2911 {
2912  for (auto& jetDef : fJetDefinitions) {
2913  jetDef.fJets.clear();
2914  }
2915 
2916  if (fMCMode == kMCTruth) {
2918  }
2919  else {
2921  }
2922 }
2923 
2926 {
2927  // Fill the vertex info of the candidates
2928  // Needed for reduced delta AOD, where the vertex info has been deleted
2929  // to reduce the delta AOD file size
2931 
2932  const Int_t nD = fCandidateArray->GetEntriesFast();
2933 
2934  AliDmesonJetInfo DmesonJet;
2935  DmesonJet.fEvent = this->fAodEvent;
2936 
2937  std::map<AliHFJetDefinition*,Double_t> maxJetPt;
2938  for (auto& def : fJetDefinitions) maxJetPt[&def] = 0;
2939  Double_t maxDPt = 0;
2940 
2941  std::array<int, 3> nAccCharm = {0};
2942  std::array<std::array<int, 3>, 5> nAccCharmPt = {{{0}}};
2943 
2944 
2945  //fill the mc efficiency//
2946  for (auto& def : fJetDefinitions)GetEfficiencyDenominator(def);
2947 
2948 
2949  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
2950  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
2951  if (!charmCand) continue;
2952  if(!(vHF.FillRecoCand(fAodEvent,charmCand))) continue;
2953 
2954 
2955 
2956 
2957  //region of interest + cuts
2958  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
2959  Int_t nMassHypo = 0; // number of mass hypothesis accepted for this D meson
2960  if (charmCand->Pt() > maxDPt) maxDPt = charmCand->Pt();
2961  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
2962  DmesonJet.Reset();
2963  DmesonJet.fDmesonParticle = charmCand;
2964  DmesonJet.fSelectionType = im + 1;
2965  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
2966  for (auto& def : fJetDefinitions) {
2967  if (FindJet(charmCand, DmesonJet, def,im)) {
2968  Double_t jetPt = DmesonJet.fJets[def.GetName()].fMomentum.Pt();
2969  ExtractEfficiencies(charmCand,DmesonJet,def,im);
2970  if (jetPt > maxJetPt[&def]) maxJetPt[&def] = jetPt;
2971  }
2972  else {
2973  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
2974  def.GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
2975  }
2976  }
2977  fDmesonJets[(icharm+1)*(1-(im*2))] = DmesonJet;
2978  nMassHypo++;
2979  nAccCharm[im]++;
2980 
2981  for (int i = 0; i < nAccCharmPt.size(); i++) {
2982  if (charmCand->Pt() < i) break;
2983  nAccCharmPt[i][im]++;
2984  }
2985  }
2986  }
2987  if (nMassHypo == 2) { // both mass hypothesis accepted
2988  nAccCharm[0]--;
2989  nAccCharm[1]--;
2990  nAccCharm[2]++;
2991 
2992  for (int i = 0; i < nAccCharmPt.size(); i++) {
2993  if (charmCand->Pt() < i) break;
2994  nAccCharmPt[i][0]--;
2995  nAccCharmPt[i][1]--;
2996  nAccCharmPt[i][2]++;
2997  }
2998 
2999  fDmesonJets[(icharm+1)].fD0D0bar = kTRUE;
3000  fDmesonJets[-(icharm+1)].fD0D0bar = kTRUE;
3001  }
3002  } // end of D cand loop
3003 
3004  TString hname;
3005 
3006  Int_t ntracks = 0;
3007 
3008  for (auto track_cont : fTrackContainers) {
3009  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
3010  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(nullptr);
3011  ntracks += track_cont->GetNAcceptEntries();
3012  }
3013 
3014  for (auto& def : fJetDefinitions) {
3015  if (!def.fRho) continue;
3016  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
3017  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
3018 
3019  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
3020  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
3021 
3022  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
3023  fHistManager->FillTH2(hname, fEventInfo.fCent, def.fRho->GetVal());
3024 
3025  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
3026  fHistManager->FillTH2(hname, fEventInfo.fCent, maxJetPt[&def]);
3027 
3028  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
3029  fHistManager->FillTH2(hname, fEventInfo.fCent, maxDPt);
3030 
3031  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
3032  fHistManager->FillTH2(hname, ntracks, def.fRho->GetVal());
3033 
3034  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
3035  fHistManager->FillTH2(hname, ntracks, maxJetPt[&def]);
3036 
3037  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
3038  fHistManager->FillTH2(hname, ntracks, maxDPt);
3039  }
3040 
3041  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
3042  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
3043  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
3044  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
3045 
3046  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
3047  fHistManager->FillTH2(hname, ntracks, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
3048 
3049  for (int i = 0; i < nAccCharmPt.size(); i++) {
3050  hname = TString::Format("%s/fHistNTotAcceptedDmesonsPt%d", GetName(), i);
3051  fHistManager->FillTH1(hname, "D", nAccCharmPt[i][0]);
3052  fHistManager->FillTH1(hname, "Anti-D", nAccCharmPt[i][1]);
3053  fHistManager->FillTH1(hname, "Both", nAccCharmPt[i][2]);
3054 
3055  hname = TString::Format("%s/fHistNAcceptedDmesonsPt%d", GetName(), i);
3056  fHistManager->FillTH1(hname, nAccCharmPt[i][0]+nAccCharmPt[i][1]+nAccCharmPt[i][2]);
3057  }
3058 
3059  hname = TString::Format("%s/fHistNDmesons", GetName());
3060  fHistManager->FillTH1(hname, nD);
3061 }
3062 
3073 {
3074  TString hname;
3075 
3076  Double_t rho = 0;
3077  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
3078 
3080  fFastJetWrapper->SetR(jetDef.fRadius);
3083 
3084  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
3085 
3086  if (jetDef.fJetType != AliJetContainer::kNeutralJet) {
3087  for (auto track_cont : fTrackContainers) {
3088  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
3089  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(Dcand);
3090  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
3091  AddInputVectors(track_cont, 100, static_cast<TH2*>(fHistManager->FindObject(hname)), fTrackEfficiency);
3092 
3093  if (hftrack_cont) {
3094  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
3095  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
3096  const TObjArray& daughters = hftrack_cont->GetDaughterList();
3097  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
3098  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
3099  if (!hftrack_cont->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
3100  }
3101  }
3102  }
3103  }
3104 
3105  if (jetDef.fJetType != AliJetContainer::kChargedJet) {
3106  for (auto clus_cont : fClusterContainers) {
3107  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
3108  AddInputVectors(clus_cont, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
3109  }
3110  }
3111 
3112  // run jet finder
3113  fFastJetWrapper->Run();
3114 
3115  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
3116 
3117  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
3118  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
3119 
3120  Bool_t isDmesonJet = kFALSE;
3121 
3122  Double_t maxChPt = 0;
3123  Double_t maxNePt = 0;
3124  Double_t totalNeutralPt = 0;
3125  Int_t nConst = 1;
3126 
3127  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
3128  if (constituents[ic].user_index() == 0) {
3129  isDmesonJet = kTRUE;
3130  }
3131  else if (constituents[ic].user_index() >= 100) {
3132  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
3133  nConst++;
3134  }
3135  else if (constituents[ic].user_index() <= -100) {
3136  totalNeutralPt += constituents[ic].pt();
3137  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
3138  nConst++;
3139  }
3140  }
3141 
3142  if (isDmesonJet) {
3143  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
3144  DmesonJet.fJets[jetDef.GetName()].fNConstituents = nConst;
3145  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
3146  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
3147  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
3148  DmesonJet.fJets[jetDef.GetName()].fArea = jets_incl[ijet].area();
3149  DmesonJet.fJets[jetDef.GetName()].fCorrPt = DmesonJet.fJets[jetDef.GetName()].fMomentum.Pt() - jets_incl[ijet].area() * rho;
3150  IterativeDeclustering(ijet,1,jetDef, DmesonJet.fD.M());
3151  return kTRUE;
3152  }
3153  if(!isDmesonJet && numcand!=1) IterativeDeclustering(ijet,0,jetDef,0.);
3154  }
3155 
3156  return kFALSE;
3157 }
3158 
3163 {
3164  auto itcont = cont->all_momentum();
3165  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
3166  UInt_t rejectionReason = 0;
3167  //only physical primaries are accepted for the jet finding
3168  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
3169  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
3170  continue;
3171  }
3172  if (fRandomGen && eff > 0 && eff < 1) {
3173  Double_t rnd = fRandomGen->Rndm();
3174  if (eff < rnd) {
3175  if (rejectHist) rejectHist->Fill(6, it->first.Pt());
3176  continue;
3177  }
3178  }
3179  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
3180 
3181  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
3182  }
3183 }
3184 
3185 
3187 {
3188 
3189  double nall = 0;
3190  double zg = 0.;
3191  double flagSubjet=0;
3192  double xconstperp=0;
3193  TString hname;
3194  TString hname2;
3195  fastjet::JetAlgorithm jet_algo(fastjet::cambridge_algorithm);
3196  double jet_radius_ca = 1.0;
3197  fastjet::JetDefinition jet_def(jet_algo, jet_radius_ca,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
3198  hname = TString::Format("%s/%s/LundIterative", GetName(), jetDef.GetName());
3199  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
3200  hname2 = TString::Format("%s/%s/AngleDifference", GetName(), jetDef.GetName());
3201  TH2* hdiffangle = static_cast<TH2*>(fHistManager->FindObject(hname2));
3202  try{
3203  std::vector<fastjet::PseudoJet> particles(fFastJetWrapper->GetJetConstituents(ijet));
3204  fastjet::ClusterSequence cs_ca(particles, jet_def);
3205  std::vector<fastjet::PseudoJet> output_jets = cs_ca.inclusive_jets(0);
3206  output_jets = sorted_by_pt(output_jets);
3207 
3208 
3209  fastjet::PseudoJet jj = output_jets[0];
3210  fastjet::PseudoJet j1;
3211  fastjet::PseudoJet j2;
3212 
3213  while(jj.has_parents(j1,j2)){
3214  nall = nall + 1;
3215  if(j1.perp() < j2.perp()) std::swap(j1,j2);
3216  flagSubjet=0;
3217  vector < fastjet::PseudoJet > constitj1 = sorted_by_pt(j1.constituents());
3218  if(type==1){
3219  for(Int_t j=0;j<constitj1.size();j++){
3220  if(constitj1[j].user_index()==0){
3221  xconstperp=constitj1[j].perp();
3222  flagSubjet=1; }}}
3223 
3224 
3225  double delta_R = j1.delta_R(j2);
3226  double delta_Raxis=j2.delta_R(output_jets[0]);
3227  zg = j2.perp()/(j1.perp()+j2.perp());
3228  double yh=j1.e()+j2.e();
3229  double y = log(1.0/delta_R);
3230  double lnpt_rel = log(j2.perp()*delta_R);
3231 
3232  double lundEntries[10] = {y, lnpt_rel, output_jets[0].perp(), nall, type, flagSubjet, xconstperp, invmass,yh,TMath::Abs(output_jets[0].eta())};
3233  h->Fill(lundEntries);
3234  hdiffangle->Fill(delta_R, delta_Raxis);
3235  jj=j1;
3236  }
3237 
3238  if(nall==0){ double lundEntrieszero[10]={0,0,output_jets[0].perp(),0,type,0,0,invmass,0,TMath::Abs(output_jets[0].eta())};
3239  h->Fill(lundEntrieszero);}
3240 
3241  } catch (fastjet::Error) { /*return -1;*/ }
3242 
3243 }
3244 
3247 {
3248  TString hname;
3249 
3254 
3255  if (!fMCContainer->IsSpecialPDGFound()) return;
3256 
3257  std::array<int,2> nAccCharm = {0};
3258  std::array<std::array<int, 2>, 5> nAccCharmPt = {{{0}}};
3259 
3260  std::map<AliHFJetDefinition*, Double_t> maxJetPt;
3261  Double_t maxDPt = 0;
3262 
3263  for (auto &jetDef : fJetDefinitions) {
3264  maxJetPt[&jetDef] = 0;
3265  Double_t rho = 0;
3266  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
3267  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", GetName(), jetDef.GetName());
3268  TH1* histNDmesonsVsNconstituents = static_cast<TH1*>(fHistManager->FindObject(hname));
3269 
3270  switch (jetDef.fJetType) {
3272  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
3273  break;
3275  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
3276  break;
3278  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
3279  break;
3280  }
3281 
3283  fFastJetWrapper->SetR(jetDef.fRadius);
3286 
3287  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
3288  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
3289 
3290  fFastJetWrapper->Run();
3291 
3292  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
3293 
3294  for (auto jet : jets_incl) {
3295  Int_t nDmesonsInJet = 0;
3296 
3297  for (auto constituent : jet.constituents()) {
3298  Int_t iPart = constituent.user_index() - 100;
3299  if (constituent.perp() < 1e-6) continue; // reject ghost particles
3300  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
3301  if (!part) {
3302  ::Error("AliAnalysisTaskDmesonJetsSub::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
3303  continue;
3304  }
3305  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
3306  nDmesonsInJet++;
3307  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
3308  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
3309  if (part->Pt() > maxDPt) maxDPt = part->Pt();
3310  std::pair<int, AliDmesonJetInfo> element;
3311  element.first = iPart;
3312  dMesonJetIt = fDmesonJets.insert(element).first;
3313  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
3314  (*dMesonJetIt).second.fDmesonParticle = part;
3315  (*dMesonJetIt).second.fSelectionType = part->PdgCode() > 0 ? 1 : 2;
3316 
3317  UShort_t p = 0;
3318  UInt_t rs = 0;
3319 
3320  auto origin = IsPromptCharm(part, fMCContainer->GetArray());
3321  p = 0;
3322  rs = origin.first;
3323  while (rs >>= 1) { p++; }
3324  (*dMesonJetIt).second.fPartonType = p;
3325  (*dMesonJetIt).second.fParton = origin.second;
3326 
3327  (*dMesonJetIt).second.fAncestor = FindParticleOrigin(part, fMCContainer->GetArray(), kFindFirst);
3328 
3329  Int_t im = -1;
3330  if (part->PdgCode() > 0) { // D0
3331  im = 0;
3332  }
3333  else { // D0bar
3334  im = 1;
3335  }
3336 
3337  nAccCharm[im]++;
3338  for (int i = 0; i < nAccCharmPt.size(); i++) {
3339  if (part->Pt() < i) break;
3340  nAccCharmPt[i][im]++;
3341  }
3342  }
3343 
3344  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
3345  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = jet.constituents().size();
3346  (*dMesonJetIt).second.fJets[jetDef.GetName()].fArea = jet.area();
3347  (*dMesonJetIt).second.fJets[jetDef.GetName()].fCorrPt = (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.Pt() - jet.area() * rho;
3348  if (jet.perp() > maxJetPt[&jetDef]) maxJetPt[&jetDef] = jet.perp();
3349  } // if constituent is a D meson
3350  } // for each constituent
3351  if (nDmesonsInJet > 0) histNDmesonsVsNconstituents->Fill(jet.constituents().size(), nDmesonsInJet);
3352  } // for each jet
3353  } // for each jet definition
3354 
3356 
3357  for (auto& def : fJetDefinitions) {
3358  if (!def.fRho) continue;
3359  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
3360  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
3361 
3362  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
3363  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
3364 
3365  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
3366  fHistManager->FillTH2(hname, fEventInfo.fCent, def.fRho->GetVal());
3367 
3368  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
3369  fHistManager->FillTH2(hname, fEventInfo.fCent, maxJetPt[&def]);
3370 
3371  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
3372  fHistManager->FillTH2(hname, fEventInfo.fCent, maxDPt);
3373 
3374  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
3375  fHistManager->FillTH2(hname, npart, def.fRho->GetVal());
3376 
3377  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
3378  fHistManager->FillTH2(hname, npart, maxJetPt[&def]);
3379 
3380  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
3381  fHistManager->FillTH2(hname, npart, maxDPt);
3382  }
3383 
3384  if (fDmesonJets.size() != nAccCharm[0]+nAccCharm[1]) AliError(Form("I found %lu mesons (%d)?", fDmesonJets.size(), nAccCharm[0]+nAccCharm[1]));
3385  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
3386  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
3387  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
3388 
3389  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
3390  fHistManager->FillTH2(hname, npart, nAccCharm[0]+nAccCharm[1]);
3391 
3392  for (int i = 0; i < nAccCharmPt.size(); i++) {
3393  hname = TString::Format("%s/fHistNTotAcceptedDmesonsPt%d", GetName(), i);
3394  fHistManager->FillTH1(hname, "D", nAccCharmPt[i][0]);
3395  fHistManager->FillTH1(hname, "Anti-D", nAccCharmPt[i][1]);
3396 
3397  hname = TString::Format("%s/fHistNAcceptedDmesonsPt%d", GetName(), i);
3398  fHistManager->FillTH1(hname, nAccCharmPt[i][0]+nAccCharmPt[i][1]);
3399  }
3400 
3401  hname = TString::Format("%s/fHistNDmesons", GetName());
3402  fHistManager->FillTH1(hname, nAccCharm[0]+nAccCharm[1]); // same as the number of accepted D mesons, since no selection is performed
3403 }
3404 
3405 
3406 
3407 // Definitions of class AliAnalysisTaskDmesonJetsSub
3408 
3412 
3416  fAnalysisEngines(),
3417  fEnabledAxis(0),
3419  fHistManager(),
3420  fApplyKinematicCuts(kTRUE),
3421  fNOutputTrees(0),
3422  fTrackEfficiency(0),
3423  fRejectISR(kFALSE),
3424  fJetAreaType(fastjet::active_area),
3425  fJetGhostArea(0.005),
3426  fMCContainer(0),
3427  fAodEvent(0),
3428  fFastJetWrapper(0)
3429 {
3430  SetMakeGeneralHistograms(kTRUE);
3431 }
3432 
3437  AliAnalysisTaskEmcalLight(name, kTRUE),
3438  fAnalysisEngines(),
3441  fHistManager(name),
3442  fApplyKinematicCuts(kTRUE),
3443  fNOutputTrees(nOutputTrees),
3444  fTrackEfficiency(0),
3445  fRejectISR(kFALSE),
3446  fJetAreaType(fastjet::active_area),
3447  fJetGhostArea(0.005),
3448  fMCContainer(0),
3449  fAodEvent(0),
3450  fFastJetWrapper(0)
3451 {
3452  SetMakeGeneralHistograms(kTRUE);
3453  for (Int_t i = 0; i < nOutputTrees; i++){
3454  DefineOutput(2+i, TTree::Class());
3455  }
3456 }
3457 
3460 {
3461  if (fFastJetWrapper) delete fFastJetWrapper;
3462 }
3463 
3471 {
3472  AliRDHFCuts* analysiscuts = 0;
3473  TFile* filecuts = TFile::Open(cutfname);
3474  if (!filecuts || filecuts->IsZombie()) {
3475  ::Error("AliAnalysisTaskDmesonJetsSub::LoadDMesonCutsFromFile", "Input file not found: will use std cuts.");
3476  filecuts = 0;
3477  }
3478 
3479  if (filecuts) analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
3480 
3481  if (!analysiscuts) {
3482  ::Error("AliAnalysisTaskDmesonJetsSub::LoadDMesonCutsFromFile", "Could not find analysis cuts '%s' in '%s'.", cutsname.Data(), cutfname.Data());
3483  if (filecuts) {
3484  filecuts->ls();
3485  }
3486  }
3487  else {
3488  ::Info("AliAnalysisTaskDmesonJetsSub::LoadDMesonCutsFromFile", "Cuts '%s' loaded from file '%s'", cutsname.Data(), cutfname.Data());
3489  }
3490 
3491  return analysiscuts;
3492 }
3493 
3506 {
3508  return AddAnalysisEngine(type, cutfname, cuttype, MCmode, jetDef, rhoName);
3509 }
3510 
3522 {
3523  AliRDHFCuts* cuts = 0;
3524 
3525  if (!cutfname.IsNull()) {
3526  TString cutsname;
3527 
3528  switch (type) {
3529  case kD0toKpi:
3530  case kD0toKpiLikeSign:
3531  cutsname = "D0toKpiCuts";
3532  break;
3533  case kDstartoKpipi:
3534  cutsname = "DStartoKpipiCuts";
3535  break;
3536  default:
3537  return 0;
3538  }
3539 
3540  if (!cuttype.IsNull()) {
3541  cutsname += TString::Format("_%s", cuttype.Data());
3542  }
3543 
3544  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
3545  if (cuts) cuts->PrintAll();
3546  }
3547 
3548  AnalysisEngine eng(type, MCmode, cuts);
3549 
3550  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
3551 
3552  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
3553  eng.AddJetDefinition(jetDef);
3554  it = fAnalysisEngines.insert(it, eng);
3555  ::Info("AliAnalysisTaskDmesonJetsSub::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
3556  }
3557  else {
3558  AnalysisEngine* found_eng = &(*it);
3559  ::Info("AliAnalysisTaskDmesonJetsSub::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());
3560  found_eng->AddJetDefinition(jetDef);
3561 
3562  if (cuts) {
3563  if (found_eng->fRDHFCuts != 0) ::Warning("AliAnalysisTaskDmesonJetsSub::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
3564  found_eng->SetRDHFCuts(cuts);
3565  }
3566  }
3567 
3568  return &(*it);
3569 }
3570 
3571 std::list<AliAnalysisTaskDmesonJetsSub::AnalysisEngine>::iterator AliAnalysisTaskDmesonJetsSub::FindAnalysisEngine(const AliAnalysisTaskDmesonJetsSub::AnalysisEngine& eng)
3572 {
3573  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
3574  while (it != fAnalysisEngines.end() && (*it) != eng) it++;
3575  return it;
3576 }
3577 
3580 {
3581  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
3582 
3584 
3585  // Define histograms
3586  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
3587 
3588  TString hname;
3589  TString htitle;
3590  TH1* h = 0;
3591  Int_t treeSlot = 0;
3592 
3593  Int_t maxTracks = 6000;
3594  Double_t maxRho = 500;
3595  if (fForceBeamType == kpp) {
3596  maxRho = 50;
3597  maxTracks = 200;
3598  }
3599  else if (fForceBeamType == kpA) {
3600  maxRho = 200;
3601  maxTracks = 500;
3602  }
3603 
3604  Int_t dimx = 10;
3605  Int_t nbinsx[10] = {50,100,10,20,2,2,200,150,100,9};
3606  Double_t minx[10] = {0,-10,0,0,0,0,0,1.6,0,0};
3607  Double_t maxx[10] = {5,10,100,20,2,2,100,2.3,100,0.9};
3608  TString titlex[10]={"log(1/deltaR)","log(zteta)","jet pt","n","type","flagSubjet","ptD","invmass","frac","abs(eta)"};
3609 
3610 
3611 
3612 
3613  hname = "fHistCharmPt";
3614  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
3615  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3616 
3617  hname = "fHistCharmEta";
3618  htitle = hname + ";#eta_{charm};counts";
3619  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
3620 
3621  hname = "fHistCharmPhi";
3622  htitle = hname + ";#phi_{charm};counts";
3623  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
3624 
3625  hname = "fHistBottomPt";
3626  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
3627  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3628 
3629  hname = "fHistBottomEta";
3630  htitle = hname + ";#eta_{bottom};counts";
3631  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
3632 
3633  hname = "fHistBottomPhi";
3634  htitle = hname + ";#phi_{bottom};counts";
3635  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
3636 
3637  hname = "fHistHighestPartonPt";
3638  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
3639  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3640 
3641  hname = "fHistHighestPartonType";
3642  htitle = hname + ";type;counts";
3643  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
3644 
3645  hname = "fHistNHeavyQuarks";
3646  htitle = hname + ";number of heavy-quarks;counts";
3647  fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
3648 
3649  ::Info("AliAnalysisTaskDmesonJetsSub::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
3650  for (auto &param : fAnalysisEngines) {
3651  ::Info("AliAnalysisTaskDmesonJetsSub::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
3652 
3653  param.fHistManager = &fHistManager;
3654 
3655  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", param.GetName());
3656  htitle = hname + ";#it{N}_{tracks};#it{N}_{D};events";
3657  h = fHistManager.CreateTH2(hname, htitle, 251, -0.5, 250.5, 21, -0.5, 20.5);
3658 
3659  hname = TString::Format("%s/fHistNTotAcceptedDmesons", param.GetName());
3660  htitle = hname + ";;#it{N}_{D}";
3661  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
3662 
3663  for (int i = 0 ; i < 5; i++) {
3664  hname = TString::Format("%s/fHistNAcceptedDmesonsPt%d", param.GetName(), i);
3665  htitle = hname + ";#it{N}_{D};events";
3666  h = fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
3667 
3668  hname = TString::Format("%s/fHistNTotAcceptedDmesonsPt%d", param.GetName(), i);
3669  htitle = hname + ";;#it{N}_{D}";
3670  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
3671  }
3672 
3673  hname = TString::Format("%s/fHistNDmesons", param.GetName());
3674  htitle = hname + ";#it{N}_{D};events";
3675  h = fHistManager.CreateTH1(hname, htitle, 501, -0.5, 500.5);
3676 
3677  hname = TString::Format("%s/fHistNEvents", param.GetName());
3678  htitle = hname + ";Event status;counts";
3679  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
3680  h->GetXaxis()->SetBinLabel(1, "Accepted");
3681  h->GetXaxis()->SetBinLabel(2, "Rejected");
3682 
3683  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
3684  htitle = hname + ";Rejection reason;counts";
3685  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
3686 
3687  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
3688  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
3689  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
3690 
3691  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
3692  htitle = hname + ";#it{#eta}_{D};counts";
3693  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
3694 
3695  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
3696  htitle = hname + ";#it{#phi}_{D};counts";
3697  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
3698 
3699  hname = TString::Format("%s/EfficiencyMatchesPrompt",param.GetName());
3700  htitle = hname + ";D meson matches prompt";
3701  fHistManager.CreateTH2(hname,htitle,100,0,50,20,0,100);
3702 
3703  hname = TString::Format("%s/EfficiencyGeneratorPrompt",param.GetName());
3704  htitle = hname + ";D meson part level prompt";
3705  fHistManager.CreateTH2(hname,htitle,100,0,50,20,0,100);
3706 
3707  hname = TString::Format("%s/EfficiencyMatchesNonPrompt",param.GetName());
3708  htitle = hname + ";D meson matches non prompt";
3709  fHistManager.CreateTH2(hname,htitle,100,0,50,20,0,100);
3710 
3711  hname = TString::Format("%s/EfficiencyGeneratorNonPrompt",param.GetName());
3712  htitle = hname + ";D meson part level non prompt";
3713  fHistManager.CreateTH2(hname,htitle,100,0,50,20,0,100);
3714 
3715  if (param.fMCMode != kMCTruth) {
3716  if (param.fCandidateType == kD0toKpi || param.fCandidateType == kD0toKpiLikeSign) {
3717  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
3718  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
3719  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
3720  }
3721  else if (param.fCandidateType == kDstartoKpipi) {
3722  Double_t min = 0;
3723  Double_t max = 0;
3724 
3725  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
3726  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
3727  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
3728  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
3729 
3730  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
3731  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
3732  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
3733  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
3734  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
3735  }
3736  }
3737 
3738  if (param.fMCMode == kMCTruth) {
3739  hname = TString::Format("%s/fHistPartonPt", param.GetName());
3740  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
3741  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3742 
3743  hname = TString::Format("%s/fHistPartonEta", param.GetName());
3744  htitle = hname + ";#eta_{parton};counts";
3745  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
3746 
3747  hname = TString::Format("%s/fHistPartonPhi", param.GetName());
3748  htitle = hname + ";#phi_{parton};counts";
3749  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
3750 
3751  hname = TString::Format("%s/fHistPartonType", param.GetName());
3752  htitle = hname + ";type;counts";
3753  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
3754 
3755  hname = TString::Format("%s/fHistPrompt", param.GetName());
3756  htitle = hname + ";Type;counts";
3757  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
3758  h->GetXaxis()->SetBinLabel(1, "Unknown");
3759  h->GetXaxis()->SetBinLabel(2, "Prompt");
3760  h->GetXaxis()->SetBinLabel(3, "Non-Prompt");
3761 
3762  hname = TString::Format("%s/fHistAncestor", param.GetName());
3763  htitle = hname + ";Ancestor;counts";
3764  h = fHistManager.CreateTH1(hname, htitle, 4, 0, 4);
3765  h->GetXaxis()->SetBinLabel(1, "Unknown");
3766  h->GetXaxis()->SetBinLabel(2, "Charm");
3767  h->GetXaxis()->SetBinLabel(3, "Bottom");
3768  h->GetXaxis()->SetBinLabel(4, "Proton");
3769  }
3770 
3771  for (auto& jetDef : param.fJetDefinitions) {
3772  ::Info("AliAnalysisTaskDmesonJetsSub::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef.GetName());
3773 
3774  if (param.fMCMode == kMCTruth) {
3775  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", param.GetName(), jetDef.GetName());
3776  htitle = hname + ";#it{N}_{constituents};#it{N}_{D};counts";
3777  h = fHistManager.CreateTH2(hname, htitle, 51, -0.5, 50.5, 10, 0.5, 10.5);
3778  }
3779 
3780  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef.GetName());
3781  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
3782  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
3783  SetRejectionReasonLabels(h->GetXaxis());
3784 
3785  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef.GetName());
3786  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
3787  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
3788  SetRejectionReasonLabels(h->GetXaxis());
3789 
3790  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef.GetName());
3791  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
3792  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
3793  SetRejectionReasonLabels(h->GetXaxis());
3794 
3795  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef.GetName());
3796  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
3797  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
3798 
3799  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef.GetName());
3800  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
3801  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
3802 
3803  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef.GetName());
3804  htitle = hname + ";#it{#eta}_{jet};counts";
3805  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
3806 
3807  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef.GetName());
3808  htitle = hname + ";#it{#phi}_{jet};counts";
3809  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
3810 
3811  if (!jetDef.fRhoName.IsNull()) {
3812  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", param.GetName(), jetDef.GetName());
3813  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});#rho (GeV/#it{c} #times rad^{-1});counts";
3814  fHistManager.CreateTH2(hname, htitle, 300, 0, 150, 1000, 0, maxRho);
3815 
3816  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", param.GetName(), jetDef.GetName());
3817  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});#rho (GeV/#it{c} #times rad^{-1});counts";
3818  fHistManager.CreateTH2(hname, htitle, 300, 0, 150, 1000, 0, maxRho);
3819 
3820  hname = TString::Format("%s/%s/fHistRhoVsCent", param.GetName(), jetDef.GetName());
3821  htitle = hname + ";Centrality (%);#rho (GeV/#it{c} #times rad^{-1});counts";
3822  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 1000, 0, maxRho);
3823 
3824  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", param.GetName(), jetDef.GetName());
3825  htitle = hname + ";Centrality (%);#it{p}_{T,jet} (GeV/#it{c});counts";
3826  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
3827 
3828  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", param.GetName(), jetDef.GetName());
3829  htitle = hname + ";Centrality (%);#it{p}_{T,D} (GeV/#it{c});counts";
3830  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
3831 
3832  hname = TString::Format("%s/%s/fHistRhoVsNTracks", param.GetName(), jetDef.GetName());
3833  htitle = hname + ";no. of tracks;#rho (GeV/#it{c} #times rad^{-1});counts";
3834  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 1000, 0, maxRho);
3835 
3836  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", param.GetName(), jetDef.GetName());
3837  htitle = hname + ";no. of tracks;#it{p}_{T,jet} (GeV/#it{c});counts";
3838  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
3839 
3840  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", param.GetName(), jetDef.GetName());
3841  htitle = hname + ";no. of tracks;#it{p}_{T,D} (GeV/#it{c});counts";
3842  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
3843  }
3844  hname = TString::Format("%s/%s/LundIterative",param.GetName(),jetDef.GetName());
3845  cout<<"at the begining"<<hname<<endl;
3846  THnSparse* h = fHistManager.CreateTHnSparse(hname,hname,dimx,nbinsx,minx,maxx);
3847  for (Int_t j = 0; j < dimx; j++) {
3848  h->GetAxis(j)->SetTitle(titlex[j]);}
3849 
3850 
3851 
3852  hname = TString::Format("%s/%s/AngleDifference",param.GetName(),jetDef.GetName());
3853  htitle = hname + ";angle iterative declustering;angle to axis";
3854  fHistManager.CreateTH2(hname,htitle,100,0.,2*3.1416,100,0,2*3.1416);
3855 
3856  }
3857  switch (fOutputType) {
3858  case kTreeOutput:
3859  {
3860  OutputHandlerTTree* tree_handler = new OutputHandlerTTree(&param);
3861  param.fOutputHandler = tree_handler;
3862  tree_handler->BuildOutputObject(GetName());
3863  if (treeSlot < fNOutputTrees) {
3864  tree_handler->AssignDataSlot(treeSlot+2);
3865  treeSlot++;
3866  PostDataFromAnalysisEngine(tree_handler);
3867  }
3868  else {
3869  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
3870  }
3871  }
3872  break;
3873  case kTHnOutput:
3874  {
3875  OutputHandlerTHnSparse* thnsparse_handler = new OutputHandlerTHnSparse(&param);
3876  param.fOutputHandler = thnsparse_handler;
3877  thnsparse_handler->SetEnabledAxis(fEnabledAxis);
3878  thnsparse_handler->BuildOutputObject(GetName());
3879  }
3880  break;
3881  case kOnlyQAOutput:
3882  {
3883  OutputHandler* qa_handler = new OutputHandler(&param);
3884  param.fOutputHandler = qa_handler;
3885  }
3886  break;
3887  case kNoOutput:
3888  break;
3889  case kTreeExtendedOutput:
3890  {
3892  param.fOutputHandler = tree_handler;
3893  tree_handler->BuildOutputObject(GetName());
3894  if (treeSlot < fNOutputTrees) {
3895  tree_handler->AssignDataSlot(treeSlot+2);
3896  treeSlot++;
3897  PostDataFromAnalysisEngine(tree_handler);
3898  }
3899  else {
3900  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
3901  }
3902  }
3903  break;
3904  }
3905  }
3906 
3908 
3909  PostData(1, fOutput);
3910 }
3911 
3915 {
3917 
3918  // Load the event
3919  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
3920 
3921  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
3922 
3923  fFastJetWrapper->SetAreaType((fastjet::AreaType)fJetAreaType);
3925 
3926  if (!fAodEvent) {
3927  AliError(Form("This task need an AOD event (Task '%s'). Expect troubles...", GetName()));
3928  //return;
3929  }
3930 
3931  TRandom* rnd = 0;
3932  if (fTrackEfficiency > 0 && fTrackEfficiency < 1) rnd = new TRandom3(0);
3933 
3934  for (auto cont_it : fParticleCollArray) {
3935  AliHFAODMCParticleContainer* part_cont = dynamic_cast<AliHFAODMCParticleContainer*>(cont_it.second);
3936  if (part_cont) fMCContainer = part_cont;
3937  }
3938 
3939  for (auto &params : fAnalysisEngines) {
3940 
3941  params.fAodEvent = fAodEvent;
3942  params.fFastJetWrapper = fFastJetWrapper;
3943  params.fTrackEfficiency = fTrackEfficiency;
3944  params.fRejectISR = fRejectISR;
3945  params.fRandomGen = rnd;
3946 
3947  for (auto &jetdef: params.fJetDefinitions) {
3948  if (!jetdef.fRhoName.IsNull()) {
3949  jetdef.fRho = dynamic_cast<AliRhoParameter*>(fInputEvent->FindListObject(jetdef.fRhoName));
3950  if (!jetdef.fRho) {
3951  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3952  "%s: Could not find rho object '%s' for engine '%s'",
3953  GetName(), jetdef.fRhoName.Data(), params.GetName());
3954  }
3955  }
3956  }
3957 
3958  if (!params.fRDHFCuts) {
3959  if (params.fMCMode == kMCTruth) {
3960  ::Warning("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3961  "%s: RDHF cuts not provided for engine '%s'.",
3962  GetName(), params.GetName());
3963  }
3964  else {
3965  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3966  "%s: RDHF cuts not provided. Engine '%s' disabled.",
3967  GetName(), params.GetName());
3968  params.fInhibit = kTRUE;
3969  }
3970  }
3971 
3972  params.fMCContainer = fMCContainer;
3973 
3974  for (auto cont_it : fParticleCollArray) {
3975  AliTrackContainer* track_cont = dynamic_cast<AliTrackContainer*>(cont_it.second);
3976  if (track_cont) params.fTrackContainers.push_back(track_cont);
3977  }
3978 
3979  for (auto cont_it : fClusterCollArray) params.fClusterContainers.push_back(cont_it.second);
3980 
3981  if (fAodEvent) params.Init(fGeom, fAodEvent->GetRunNumber());
3982 
3983  if (params.fMCMode != kMCTruth && fAodEvent) {
3984  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
3985 
3986  if (params.fCandidateArray) {
3987  TString className;
3988  if (params.fCandidateType == kD0toKpi || params.fCandidateType == kD0toKpiLikeSign) {
3989  className = "AliAODRecoDecayHF2Prong";
3990  }
3991  else if (params.fCandidateType == kDstartoKpipi) {
3992  className = "AliAODRecoCascadeHF";
3993  }
3994  if (!params.fCandidateArray->GetClass()->InheritsFrom(className)) {
3995  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3996  "%s: Objects of type %s in %s are not inherited from %s! Task will be disabled!",
3997  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName(), className.Data()); // @suppress("Ambiguous problem")
3998  params.fCandidateArray = 0;
3999  params.fInhibit = kTRUE;
4000  }
4001  }
4002  else {
4003  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
4004  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
4005  params.fBranchName.Data(), params.GetName());
4006  params.fInhibit = kTRUE;
4007  }
4008  }
4009 
4010  if (params.fMCMode != kNoMC) {
4011  if (!params.fMCContainer) {
4012  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
4013  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
4014  params.GetName());
4015  params.fInhibit = kTRUE;
4016  }
4017  }
4018 
4019  if (params.fMCMode != kMCTruth) {
4020  if (params.fTrackContainers.size() == 0 && params.fClusterContainers.size() == 0) {
4021  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
4022  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
4023  params.GetName());
4024  params.fInhibit = kTRUE;
4025  }
4026  }
4027  }
4028 }
4029 
4034 {
4035  TString hname;
4036 
4037  // Fix for temporary bug in ESDfilter
4038  // The AODs with null vertex pointer didn't pass the PhysSel
4039  // Now adding an entry in the histogram so as to check that this is actually cutting anything out
4040  if (fAodEvent && (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001)) {
4041  for (auto &eng : fAnalysisEngines) {
4042  if (eng.fInhibit) continue;
4043  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
4044  fHistManager.FillTH1(hname, "ESDfilterBug");
4045  }
4046  return kFALSE;
4047  }
4048 
4049  if (fAodEvent) {
4050  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
4051  if (matchingAODdeltaAODlevel <= 0) {
4052  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
4053  for (auto &eng : fAnalysisEngines) {
4054  if (eng.fInhibit) continue;
4055  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
4056  fHistManager.FillTH1(hname, "MismatchDeltaAOD");
4057  }
4058  return kFALSE;
4059  }
4060  }
4061 
4062  for (auto &eng : fAnalysisEngines) {
4063  eng.fDmesonJets.clear();
4064  if (eng.fInhibit) continue;
4065 
4066  eng.fEventInfo = EventInfo(fCent, fEPV0, fEventWeight, fPtHard);
4067 
4068  //Event selection
4069  hname = TString::Format("%s/fHistNEvents", eng.GetName());
4070  if (fAodEvent) {
4071  Bool_t iseventselected = kTRUE;
4072  if (eng.fRDHFCuts) iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
4073  if (!iseventselected) {
4074  fHistManager.FillTH1(hname, "Rejected");
4075  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
4076  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
4077  TString label;
4078  do {
4079  label = GetHFEventRejectionReasonLabel(bitmap);
4080  if (label.IsNull()) break;
4081  fHistManager.FillTH1(hname, label);
4082  } while (true);
4083 
4084  continue;
4085  }
4086  }
4087 
4088  fHistManager.FillTH1(hname, "Accepted");
4089 
4090  AliDebug(2, "Event selected");
4091 
4092  eng.RunAnalysis();
4093  }
4094  return kTRUE;
4095 }
4096 
4101 {
4102  for (auto &param : fAnalysisEngines) {
4103  if (param.fInhibit) continue;
4104  param.fOutputHandler->FillOutput(fApplyKinematicCuts);
4105  PostDataFromAnalysisEngine(param.fOutputHandler);
4106  }
4108  return kTRUE;
4109 }
4110 
4113 {
4114  auto itcont = fMCContainer->all_momentum();
4115  Int_t nHQ = 0;
4116  Double_t highestPartonPt = 0;
4117  Int_t PdgHighParton = 0;
4118  for (auto it = itcont.begin(); it != itcont.end(); it++) {
4119  auto part = *it;
4120  if (part.first.Pt() == 0) continue;
4121 
4122  Int_t PdgCode = part.second->GetPdgCode();
4123 
4124  // Skip all particles that are not either quarks or gluons
4125  if ((PdgCode < -9 || PdgCode > 9) && PdgCode != 21 && PdgCode != -21) continue;
4126 
4127  AliDebugStream(5) << "Parton " << it.current_index() <<
4128  " with pdg=" << PdgCode <<
4129  ", px=" << part.first.Px() <<
4130  ", py=" << part.first.Py() <<
4131  ", pz=" << part.first.Pz() <<
4132  ", n daughters = " << part.second->GetNDaughters() <<
4133  std::endl;
4134 
4135  // Skip partons that do not have any children
4136  // Unclear how this can happen, it would seem that this parton were not fragmented by the generator
4137  if (part.second->GetNDaughters() == 0) continue;
4138 
4139  // Look for highest momentum parton
4140  if (highestPartonPt < part.first.Pt()) {
4141  highestPartonPt = part.first.Pt();
4142  PdgHighParton = PdgCode;
4143  }
4144 
4145  // Skip partons that are not HF
4146  if (PdgCode != 4 && PdgCode != 5 && PdgCode != -4 && PdgCode != -5) continue;
4147 
4148  Bool_t lastInPartonShower = kTRUE;
4149  Bool_t hadronDaughter = kFALSE;
4150  for (Int_t daughterIndex = part.second->GetDaughterFirst(); daughterIndex <= part.second->GetDaughterLast(); daughterIndex++){
4151  if (daughterIndex < 0) {
4152  AliDebugStream(5) << "Could not find daughter index!" << std::endl;
4153  continue;
4154  }
4155  AliAODMCParticle *daughter = fMCContainer->GetMCParticle(daughterIndex);
4156  if (!daughter) {
4157  AliDebugStream(5) << "Could not find particle with index " << daughterIndex << "!" << std::endl;
4158  continue;
4159  }
4160  Int_t daughterPdgCode = daughter->GetPdgCode();
4161  if (daughter->GetMother() != it.current_index()) {
4162  AliDebugStream(5) << "Particle " << daughterIndex << " with pdg=" << daughterPdgCode <<
4163  ", px=" << daughter->Px() <<
4164  ", py=" << daughter->Py() <<
4165  ", pz=" << daughter->Pz() <<
4166  ", is not a daughter of " << it.current_index() <<
4167  "!" << std::endl;
4168  continue;
4169  }
4170 
4171  AliDebugStream(5) << "Found daughter " << daughterIndex <<
4172  " with pdg=" << daughterPdgCode <<
4173  ", px=" << daughter->Px() <<
4174  ", py=" << daughter->Py() <<
4175  ", pz=" << daughter->Pz() <<
4176  std::endl;
4177  // Codes between 81 and 100 are for internal MC code use, they may be intermediate states used e.g. in hadronizaion models
4178  if (daughterPdgCode == PdgCode) lastInPartonShower = kFALSE; // this parton is not the last parton in the shower
4179  if (TMath::Abs(daughterPdgCode) >= 111 || (daughterPdgCode >= 81 && daughterPdgCode <= 100)) hadronDaughter = kTRUE;
4180  }
4181  if (hadronDaughter) {
4182  AliDebugStream(5) << "This particle has at least a hadron among its daughters!" << std::endl;
4183  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;
4184  }
4185  else {
4186  AliDebugStream(5) << "This particle does not have hadrons among its daughters!" << std::endl;
4187  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;
4188  continue;
4189  }
4190 
4191  if (PdgCode == 4 || PdgCode == -4) {
4192  fHistManager.FillTH1("fHistCharmPt", part.first.Pt());
4193  fHistManager.FillTH1("fHistCharmEta", part.first.Eta());
4194  fHistManager.FillTH1("fHistCharmPhi", part.first.Phi_0_2pi());
4195  }
4196  else if (PdgCode == 5 || PdgCode == -5) {
4197  fHistManager.FillTH1("fHistBottomPt", part.first.Pt());
4198  fHistManager.FillTH1("fHistBottomEta", part.first.Eta());
4199  fHistManager.FillTH1("fHistBottomPhi", part.first.Phi_0_2pi());
4200  }
4201  nHQ++;
4202  }
4203  fHistManager.FillTH1("fHistNHeavyQuarks", nHQ);
4204  fHistManager.FillTH1("fHistHighestPartonPt",highestPartonPt);
4205  Int_t partonType = 0;
4206  Int_t absPdgHighParton = TMath::Abs(PdgHighParton);
4207  if (absPdgHighParton == 9 || absPdgHighParton == 21) {
4208  partonType = 7;
4209  }
4210  else {
4211  partonType = absPdgHighParton;
4212  }
4213  fHistManager.FillTH1("fHistHighestPartonType",partonType);
4214 }
4215 
4224 {
4225  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
4226 
4227  Double_t mass = part->Mass();
4228 
4229  // To make sure that the PDG mass value is not at the edge of a bin
4230  if (nbins % 2 == 0) {
4231  minMass = mass - range / 2 - range / nbins / 2;
4232  maxMass = mass + range / 2 - range / nbins / 2;
4233  }
4234  else {
4235  minMass = mass - range / 2;
4236  maxMass = mass + range / 2;
4237  }
4238 }
4239 
4246 {
4247  static TString label;
4248  label = "";
4249 
4250  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
4251  label = "NotSelTrigger";
4252  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
4253  return label.Data();
4254  }
4255  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
4256  label = "NoVertex";
4257  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
4258  return label.Data();
4259  }
4260  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
4261  label = "TooFewVtxContrib";
4262  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
4263  return label.Data();
4264  }
4265  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
4266  label = "ZVtxOutFid";
4267  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
4268  return label.Data();
4269  }
4270  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
4271  label = "Pileup";
4272  bitmap &= ~BIT(AliRDHFCuts::kPileup);
4273  return label.Data();
4274  }
4275  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
4276  label = "OutsideCentrality";
4277  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
4278  return label.Data();
4279  }
4280  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
4281  label = "PhysicsSelection";
4282  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
4283  return label.Data();
4284  }
4285  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
4286  label = "BadSPDVertex";
4287  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
4288  return label.Data();
4289  }
4290  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
4291  label = "ZVtxSPDOutFid";
4292  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
4293  return label.Data();
4294  }
4295  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
4296  label = "CentralityFlattening";
4297  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
4298  return label.Data();
4299  }
4300  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
4301  label = "BadTrackV0Correl";
4302  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
4303  return label.Data();
4304  }
4305 
4306  return label.Data();
4307 }
4308 
4315 {
4316  if (handler->GetDataSlotNumber() >= 0 && handler->GetOutputObject()) {
4317  PostData(handler->GetDataSlotNumber(), handler->GetOutputObject());
4318  return handler->GetDataSlotNumber();
4319  }
4320  else {
4321  return -1;
4322  }
4323 }
4324 
4334 {
4335  // Get the pointer to the existing analysis manager via the static access method
4336  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
4337  if (!mgr) {
4338  ::Error("AddTaskDmesonJetsSub", "No analysis manager to connect to.");
4339  return NULL;
4340  }
4341 
4342  // Check the analysis type using the event handlers connected to the analysis manager
4343  AliVEventHandler* handler = mgr->GetInputEventHandler();
4344  if (!handler) {
4345  ::Error("AddTaskDmesonJetsSub", "This task requires an input event handler");
4346  return NULL;
4347  }
4348 
4350 
4351  if (handler->InheritsFrom("AliESDInputHandler")) {
4352  dataType = kESD;
4353  }
4354  else if (handler->InheritsFrom("AliAODInputHandler")) {
4355  dataType = kAOD;
4356  }
4357 
4358  // Init the task and do settings
4359  if (ntracks == "usedefault") {
4360  if (dataType == kESD) {
4361  ntracks = "Tracks";
4362  }
4363  else if (dataType == kAOD) {
4364  ntracks = "tracks";
4365  }
4366  else {
4367  ntracks = "";
4368  }
4369  }
4370 
4371  if (nclusters == "usedefault") {
4372  if (dataType == kESD) {
4373  nclusters = "CaloClusters";
4374  }
4375  else if (dataType == kAOD) {
4376  nclusters = "caloClusters";
4377  }
4378  else {
4379  nclusters = "";
4380  }
4381  }
4382 
4383  if (nMCpart == "usedefault") {
4384  nMCpart = "mcparticles"; // Always needs AliAODMCParticle objects
4385  }
4386 
4387  TString name("AliAnalysisTaskDmesonJetsSub");
4388  if (strcmp(suffix, "") != 0) {
4389  name += TString::Format("_%s", suffix.Data());
4390  }
4391 
4392  AliAnalysisTaskDmesonJetsSub* jetTask = new AliAnalysisTaskDmesonJetsSub(name, nMaxTrees);
4393 
4394  if (!ntracks.IsNull()) {
4395  AliHFTrackContainer* trackCont = new AliHFTrackContainer(ntracks);
4396  jetTask->AdoptParticleContainer(trackCont);
4397  }
4398 
4399  if (!nMCpart.IsNull()) {
4400  AliMCParticleContainer* partCont = new AliHFAODMCParticleContainer(nMCpart);
4401  partCont->SetEtaLimits(-1.5, 1.5);
4402  partCont->SetPtLimits(0, 1000);
4403  jetTask->AdoptParticleContainer(partCont);
4404  }
4405 
4406  jetTask->AddClusterContainer(nclusters.Data());
4407 
4408  // Final settings, pass to manager and set the containers
4409  mgr->AddTask(jetTask);
4410 
4411  // Create containers for input/output
4412  AliAnalysisDataContainer* cinput1 = mgr->GetCommonInputContainer();
4413  TString contname1(name);
4414  contname1 += "_histos";
4415  AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(contname1.Data(),
4416  TList::Class(), AliAnalysisManager::kOutputContainer,
4417  Form("%s", AliAnalysisManager::GetCommonFileName()));
4418 
4419  mgr->ConnectInput(jetTask, 0, cinput1);
4420  mgr->ConnectOutput(jetTask, 1, coutput1);
4421 
4422  for (Int_t i = 0; i < nMaxTrees; i++) {
4423  TString contname = TString::Format("%s_tree_%d", name.Data(), i);
4424  AliAnalysisDataContainer *coutput = mgr->CreateContainer(contname.Data(),
4425  TTree::Class(),AliAnalysisManager::kOutputContainer,
4426  Form("%s", AliAnalysisManager::GetCommonFileName()));
4427  mgr->ConnectOutput(jetTask, 2+i, coutput);
4428  }
4429  return jetTask;
4430 }
Double_t fSoftPionPt
! Transverse momentum of the soft pion of the D* candidate
TString fBranchName
AOD branch where the D meson candidate are found.
Class that encapsulates event properties in a very compact structure (useful for pp simulation analys...
Int_t fDataSlotNumber
! Data slot where the tree output is posted
Int_t pdg
static std::pair< AliAnalysisTaskDmesonJetsSub::EMesonOrigin_t, AliAODMCParticle * > IsPromptCharm(const AliAODMCParticle *part, TClonesArray *mcArray)
Double32_t fDeltaInvMass
< Difference between the Kpipi and the Kpi invariant masses in GeV/c2
void Reset()
Reset all fields to their default values.
friend bool operator==(const AnalysisEngine &lhs, const AnalysisEngine &rhs)
AliEMCALGeometry * fGeom
!emcal geometry
void SetDMesonCandidate(AliAODRecoDecay *c)
virtual void Set(const AliDmesonJetInfo &source, std::string n)
virtual void Set(const AliDmesonJetInfo &source)
Bool_t fInhibit
!inhibit execution of the task
Double32_t fCorrZ
Z of the D meson after subtracting average background.
virtual Int_t Run()
double Double_t
Definition: External.C:58
EventInfo * fEventInfo
! Object conatining the event information (centrality, pt hard, weight, etc.)
Bool_t ExtractRecoDecayAttributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
EDataType_t
Switch for the data type.
UInt_t fEnabledAxis
Use bit defined in EAxis_t to enable axis in the THnSparse.
virtual ~AliAnalysisTaskDmesonJetsSub()
This is the standard destructor.
virtual AliAODMCParticle * GetMCParticle(Int_t i=-1) const
void Getd0MeasMinusExpProng(Int_t ip, Double_t magf, Double_t &d0diff, Double_t &errd0diff) const
const char * title
Definition: MakeQAPdf.C:27
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...
std::map< int, AliDmesonJetInfo > fDmesonJets
! Array containing the D meson jets
AliHFJetDefinition * AddJetDefinition(EJetType_t type, Double_t r, EJetAlgo_t algo, ERecoScheme_t reco)
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
const TObjArray & GetDaughterList() const
AliTLorentzVector fD
! 4-momentum of the D meson candidate
bidirectional stl iterator over the EMCAL iterable container
std::list< AnalysisEngine > fAnalysisEngines
Array of analysis parameters.
static AliAODMCParticle * FindParticleOrigin(const AliAODMCParticle *part, TClonesArray *mcArray, EFindParticleOriginMode_t mode, const std::set< UInt_t > &pdgSet)
Double_t fMinMass
! Min mass in histogram axis
Bool_t ExtractDstarAttributes(const AliAODRecoCascadeHF *DstarCand, AliDmesonJetInfo &DmesonJet, UInt_t i)
virtual void Set(const AliDmesonJetInfo &source)
Invariant mass of the D0 meson candidate in GeV/c2.
Double_t fJetGhostArea
Area of the ghost particles.
Double_t mass
Container with name, TClonesArray and cuts for particles.
Bool_t ExtractD0Attributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
TArrayI fPDGdaughters
List of the PDG code of the daughters.
EventInfo fEventInfo
! Event info (centrality, weight, pt hard etc.)
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
AliAODMCParticle * fParton
! pointer to the parton in the shower tree of the D meson (only for particle level D mesons) ...
static Int_t CheckMatchingAODdeltaAODevents()
AnalysisEngine()
This is the default constructor, used for ROOT I/O purposes.
virtual void UserCreateOutputObjects()
Creates the output containers.
Lightweight class that encapsulates D meson jets.
Double_t fTrackEfficiency
Artificial tracking inefficiency (0...1)
Bool_t GetEfficiencyDenominatorOneByOne(AliHFJetDefinition &jetDef)
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
Int_t fDataSlotNumber
! Data slot where the tree output is posted
static UShort_t GetRejectionReasonBitPosition(UInt_t rejectionReason)
Returns the highest bit in the rejection map as reason why the object was rejected.
Double_t fMinNeutralPt
Minimum pt of the leading neutral particle (or cluster)
const AliMCParticleIterableContainer all() const
AliAnalysisTaskDmesonJetsSub()
This is the default constructor, used for ROOT I/O purposes.
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:122
Int_t fNOutputTrees
Maximum number of output trees.
Bool_t FindJet(AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, AliHFJetDefinition &jetDef, Int_t numcand)
static AliAnalysisTaskDmesonJetsSub * AddTaskDmesonJetsSub(TString ntracks="usedefault", TString nclusters="usedefault", TString nMCpart="", Int_t nMaxTrees=2, TString suffix="")
Bool_t fRejectISR
Reject initial state radiation.
void SetEtaLimits(Double_t min, Double_t max)
const AliVEvent * fEvent
! pointer to the ESD/AOD event
friend bool operator<(const AnalysisEngine &lhs, const AnalysisEngine &rhs)
Double_t fInvMass2Prong
! 2-prong mass of the D* candidate (w/o the soft pion)
static FJRecoScheme ConvertToFJRecoScheme(ERecoScheme_t reco)
Bool_t fD0D0bar
! kTRUE if selected both as D0 and D0bar
static TString GenerateJetName(EJetType_t jetType, EJetAlgo_t jetAlgo, ERecoScheme_t recoScheme, Double_t radius, AliParticleContainer *partCont, AliClusterContainer *clusCont, TString tag)
static void CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t &minMass, Double_t &maxMass)
Double32_t fSelectionType
Selection type: D0, D0bar, both.
AliEventNotFound(const std::string &class_name, const std::string &method_name)
AliRDHFCuts * fRDHFCuts
D meson candidates cuts.
void SetR(Double_t r)
Definition: AliFJWrapper.h:127
Select tracks based on specific prescriptions of HF analysis.
ERecoScheme_t fRecoScheme
Jet recombination scheme (pt scheme, E scheme, ...)
AliHFJetDefinition & operator=(const AliHFJetDefinition &source)
Output handler for D meson jet analysis.
Lightweight class that encapsulates D meson jets.
EBeamType_t fForceBeamType
forced beam type
UInt_t fAcceptedDecay
Bit mask with D meson decays that are accepted (only used for particle-level analysis) ...
AliHFAODMCParticleContainer * fMCContainer
! MC particle container
virtual Bool_t AcceptObject(Int_t i, UInt_t &rejectionReason) const =0
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:34
virtual Bool_t AcceptObject(Int_t i, UInt_t &rejectionReason) const
std::vector< AliTrackContainer * > fTrackContainers
! Track containers
AliDmesonInfoSummary * fCurrentDmesonJetInfo
! Current D meson jet info
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.
TObject * FindObject(const char *name) const
Find an object inside the container.
int Int_t
Definition: External.C:63
Double_t GetDistance(const AliJetInfo &jet, Double_t &deta, Double_t &dphi) const
void RunParticleLevelAnalysis()
Run a particle level analysis.
void Init(const AliEMCALGeometry *const geom, Int_t runNumber)
Initialize the analysis engine.
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
AliAODMCParticle * fAncestor
! pointer to the ancestor particle in the shower tree of the D meson (only for particle level D meson...
AliDmesonJetInfo & operator=(const AliDmesonJetInfo &source)
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)
std::vector< DMESONTYPE > fCurrentDmesonInfo
! Current D meson jet info
Int_t fMCLabel
! MC label, i.e. index of the generator level D meson (only for detector level D meson candidates) ...
Double32_t fPt
Transverse momentum of the jet in GeV/c.
AnalysisEngine * AddAnalysisEngine(ECandidateType_t type, TString cutfname, TString cuttype, EMCMode_t bkgMode, EJetType_t jettype, Double_t jetradius, TString rhoName="")
void IterativeDeclustering(Int_t jetnum, Double_t type, AliHFJetDefinition &jetDef, Double_t invm)
Base task in the EMCAL framework (lighter version of AliAnalysisTaskEmcal)
friend bool operator<(const AliHFJetDefinition &lhs, const AliHFJetDefinition &rhs)
Double_t Phi_0_2pi() const
UInt_t fEnabledAxis
! Use bit defined in EAxis_t to enable axis in the THnSparse
Bool_t fApplyKinematicCuts
Apply jet kinematic cuts.
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 IsAnyJetInAcceptance(const AliDmesonJetInfo &dMesonJet) const
AliHFJetDefinition()
This is the default constructor, used for ROOT I/O purposes.
Bool_t fRejectISR
! Reject initial state radiation
AliRDHFCuts * LoadDMesonCutsFromFile(TString cutfname, TString cutsname)
Base class for container structures within the EMCAL framework.
AliAODTrack * GetBachelor() const
TClonesArray * fCandidateArray
! D meson candidate array
Bool_t fD0Extended
Store extended information in the tree (only for D0 mesons)
std::vector< AliClusterContainer * > fClusterContainers
! Cluster containers
static OutputHandlerTTreeExtendedBase * GenerateOutputHandler(AnalysisEngine *eng)
EJetType_t fJetType
Jet type (charged, full, neutral)
Bool_t fD0Extended
! Store extended information in the tree (only for D0 mesons)
TString fRhoName
Name of the object that holds the average background value.
AliFJWrapper * fFastJetWrapper
! Fastjet wrapper
AnalysisEngine & operator=(const AnalysisEngine &source)
TClonesArray * GetArray() const
void FillPartonLevelHistograms()
Fill histograms with parton-level information.
std::string fClassName
Class name where the event was not found.
virtual void Clear(const Option_t *="")
Double_t fMaxChargedPt
Transverse momentum of the leading charged particle (or track)
void RunDetectorLevelAnalysis()
Run a detector level analysis.
Int_t mode
Definition: anaM.C:41
std::vector< AliHFJetDefinition > fJetDefinitions
Jet definitions.
void SetPtLimits(Double_t min, Double_t max)
AliJetInfoSummary ** fCurrentJetInfo
! Current jet info
const char * GetName() const
Generate a name for this jet definition.
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
std::map< std::string, std::vector< JETTYPE > > fCurrentJetInfo
! Current jet info
void AddInputVectors(AliEmcalContainer *cont, Int_t offset, TH2 *rejectHist=0, Double_t eff=0.)
Double32_t fR
Distance between D meson and jet axis.
const AliMCParticleIterableMomentumContainer all_momentum() const
Select MC particles based on specific prescriptions of HF analysis.
virtual void Set(const AliDmesonJetInfo &source, std::string n)
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:125
virtual void BuildOutputObject(const char *taskName)
Int_t nbinsx
void Print() const
Prints the content of this object in the standard output.
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
static const char * GetHFEventRejectionReasonLabel(UInt_t &bitmap)
Double_t fMaxChargedPt
Maximum pt of the leading charged particle (or track)
static FJJetAlgo ConvertToFJAlgo(EJetAlgo_t algo)
std::map< std::string, AliJetInfo > fJets
! list of jets
decay
Definition: HFPtSpectrum.C:42
AliVParticle * fDmesonParticle
! pointer to the particle object
Bool_t ExtractEfficiencies(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, AliHFJetDefinition &jetDef, UInt_t i)
Double_t fMaxMass
! Max mass in histogram axis
Double_t fCorrPt
Transverse momentum of the jet after subtracting the average background.
Definition: External.C:220
Look for the very first particle in the fragmentation tree.
std::vector< AliAnalysisTaskDmesonJetsSub::AliHFJetDefinition >::iterator FindJetDefinition(const AliHFJetDefinition &eng)
Double_t minMass
friend bool operator==(const AliHFJetDefinition &lhs, const AliHFJetDefinition &rhs)
Lightweight class that encapsulates D meson jets for PbPb analysis.