AliPhysics  b69f57c (b69f57c)
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 
2230 {
2231  AliDebug(10,"Checking if D0 meson is selected");
2232  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoDecayHF2Prong*>(Dcand), AliRDHFCuts::kAll, fAodEvent);
2233  if (isSelected == 0) return kFALSE;
2234 
2235  Int_t MCtruthPdgCode = 0;
2236 
2237  Double_t invMassD = 0;
2238 
2239  // If the analysis require knowledge of the MC truth, look for generated D meson matched to reconstructed candidate
2240  // Checks also the origin, and if it matches the rejected origin mask, return false
2241  if (fMCMode != kNoMC) {
2242  Int_t mcLab = Dcand->MatchToMC(fCandidatePDG, fMCContainer->GetArray(), fNDaughters, fPDGdaughters.GetArray());
2243  DmesonJet.fMCLabel = mcLab;
2244 
2245  // Retrieve the generated particle (if exists) and its PDG code
2246  if (mcLab >= 0) {
2247  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
2248 
2249  if (aodMcPart) {
2250  // Check origin and return false if it matches the rejected origin mask
2251  if (fRejectedOrigin) {
2252  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
2253  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
2254  }
2255  MCtruthPdgCode = aodMcPart->PdgCode();
2256  }
2257  }
2258  }
2259 
2260  if (isSelected == 1) { // selected as a D0
2261  if (i != 0) return kFALSE; // only one mass hypothesis thanks to PID
2262 
2263  if (fMCMode == kNoMC ||
2264  (MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
2265  (MCtruthPdgCode != fCandidatePDG && fMCMode == kBackgroundOnly) ||
2266  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kD0Reflection || fMCMode == kOnlyWrongPIDAccepted))) {
2267  // both background and signal are requested OR (it is a true D0 AND signal is requested) OR (it is NOT a D0 and background is requested)
2268  AliDebug(10,"Selected as D0");
2269  invMassD = Dcand->InvMassD0();
2270  }
2271  else { // conditions above not passed, so return FALSE
2272  return kFALSE;
2273  }
2274  }
2275  else if (isSelected == 2) { // selected as a D0bar
2276  if (i != 1) return kFALSE; // only one mass hypothesis thanks to PID
2277 
2278  if (fMCMode == kNoMC ||
2279  (MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
2280  (MCtruthPdgCode != -fCandidatePDG && fMCMode == kBackgroundOnly) ||
2281  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kD0Reflection || fMCMode == kOnlyWrongPIDAccepted))) {
2282  // both background and signal are requested OR (it is a true D0bar AND signal is requested) OR (it is NOT a D0bar and background is requested)
2283  AliDebug(10,"Selected as D0bar");
2284  invMassD = Dcand->InvMassD0bar();
2285  }
2286  else { // conditions above not passed, so return FALSE
2287  return kFALSE;
2288  }
2289  }
2290  else if (isSelected == 3) { // selected as either a D0bar or a D0 (PID on K and pi undecisive)
2291  AliDebug(10,"Selected as either D0 or D0bar");
2292 
2293  // Accept the correct mass hypothesis for signal-only and the wrong one for background-only
2294  if ((MCtruthPdgCode == fCandidatePDG && fMCMode == kSignalOnly) ||
2295  (MCtruthPdgCode == -fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kD0Reflection))) {
2296  if (i != 0) return kFALSE;
2297  AliDebug(10, "MC truth is D0");
2298  invMassD = Dcand->InvMassD0();
2299  }
2300  else if ((MCtruthPdgCode == -fCandidatePDG && fMCMode == kSignalOnly) ||
2301  (MCtruthPdgCode == fCandidatePDG && (fMCMode == kBackgroundOnly || fMCMode == kD0Reflection))) {
2302  if (i != 1) return kFALSE;
2303  AliDebug(10, "MC truth is D0bar");
2304  invMassD = Dcand->InvMassD0bar();
2305  }
2306  else { // (This candidate is neither a D0 nor a D0bar) OR (background-and-signal was requested)
2307 
2308  // Only accept it if background-only OR background-and-signal was requested
2309  if (fMCMode == kBackgroundOnly || fMCMode == kNoMC) {
2310  // Select D0 or D0bar depending on the i-parameter
2311  if (i == 0) {
2312  AliDebug(10, "Returning invariant mass with D0 hypothesis");
2313  invMassD = Dcand->InvMassD0();
2314  }
2315  else if (i == 1) {
2316  AliDebug(10, "Returning invariant mass with D0bar hypothesis");
2317  invMassD = Dcand->InvMassD0bar();
2318  }
2319  else { // i > 1
2320  return kFALSE;
2321  }
2322  }
2323  else { // signal-only was requested but this is not a true D0
2324  return kFALSE;
2325  }
2326  }
2327  }
2328  DmesonJet.fD.SetPtEtaPhiM(Dcand->Pt(), Dcand->Eta(), Dcand->Phi(), invMassD);
2329  return kTRUE;
2330 }
2331 
2340 {
2341  AliDebug(10,"Checking if D* meson is selected");
2342  Int_t isSelected = fRDHFCuts->IsSelected(const_cast<AliAODRecoCascadeHF*>(DstarCand), AliRDHFCuts::kAll, fAodEvent);
2343  if (isSelected == 0) return kFALSE;
2344 
2345  if ((i == 1 && DstarCand->Charge()>0) || (i == 0 && DstarCand->Charge()<0) || i > 1) return kFALSE; // only one mass hypothesis for the D*
2346 
2347  Int_t MCtruthPdgCode = 0;
2348 
2349  Double_t invMassD = 0;
2350 
2351  if (fMCMode == kBackgroundOnly || fMCMode == kSignalOnly) {
2352  Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
2353  Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
2354 
2355  Int_t mcLab = DstarCand->MatchToMC(fCandidatePDG, 421, pdgDgDStartoD0pi, pdgDgD0toKpi, fMCContainer->GetArray());
2356  AliDebug(10, Form("MC label is %d", mcLab));
2357  DmesonJet.fMCLabel = mcLab;
2358  if (mcLab >= 0) {
2359  AliAODMCParticle* aodMcPart = static_cast<AliAODMCParticle*>(fMCContainer->GetArray()->At(mcLab));
2360 
2361  if (aodMcPart) {
2362  if (fRejectedOrigin) {
2363  auto origin = IsPromptCharm(aodMcPart, fMCContainer->GetArray());
2364  if ((origin.first & fRejectedOrigin) == origin.first) return kFALSE;
2365  }
2366 
2367  MCtruthPdgCode = aodMcPart->PdgCode();
2368  AliDebug(10, Form("MC truth pdg code is %d",MCtruthPdgCode));
2369  }
2370  }
2371  }
2372 
2373  Int_t absMCtruthPdgCode = TMath::Abs(MCtruthPdgCode);
2374  if (fMCMode == kNoMC ||
2375  (absMCtruthPdgCode == 413 && fMCMode == kSignalOnly) ||
2376  (absMCtruthPdgCode != 413 && fMCMode == kBackgroundOnly)) {
2377  // both background and signal are requested OR (it is a true D*/D*bar AND signal is requested) OR (it is NOT a D*/D*bar and background is requested)
2378  invMassD = DstarCand->InvMassDstarKpipi();
2379  DmesonJet.fSoftPionPt = DstarCand->GetBachelor()->Pt();
2380  DmesonJet.fInvMass2Prong = DstarCand->InvMassD0();
2381  DmesonJet.fD.SetPtEtaPhiM(DstarCand->Pt(), DstarCand->Eta(), DstarCand->Phi(), invMassD);
2382  return kTRUE;
2383  }
2384  else { // conditions above not passed, so return FALSE
2385  return kFALSE;
2386  }
2387 }
2388 
2396 {
2397  if (!part) return kUnknownDecay;
2398  if (!mcArray) return kUnknownDecay;
2399 
2401 
2402  Int_t absPdgPart = TMath::Abs(part->GetPdgCode());
2403 
2404  if (part->GetNDaughters() == 2) {
2405 
2406  AliAODMCParticle* d1 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(0)));
2407  AliAODMCParticle* d2 = static_cast<AliAODMCParticle*>(mcArray->At(part->GetDaughter(1)));
2408 
2409  if (!d1 || !d2) {
2410  return decay;
2411  }
2412 
2413  Int_t absPdg1 = TMath::Abs(d1->GetPdgCode());
2414  Int_t absPdg2 = TMath::Abs(d2->GetPdgCode());
2415 
2416  if (absPdgPart == 421) { // D0 -> K pi
2417  if ((absPdg1 == 211 && absPdg2 == 321) || // pi K
2418  (absPdg1 == 321 && absPdg2 == 211)) { // K pi
2419  decay = kDecayD0toKpi;
2420  }
2421  }
2422 
2423  if (absPdgPart == 413) { // D* -> D0 pi
2424  if (absPdg1 == 421 && absPdg2 == 211) { // D0 pi
2425  Int_t D0decay = CheckDecayChannel(d1, mcArray);
2426  if (D0decay == kDecayD0toKpi) {
2427  decay = kDecayDStartoKpipi;
2428  }
2429  }
2430  else if (absPdg1 == 211 && absPdg2 == 421) { // pi D0
2431  Int_t D0decay = CheckDecayChannel(d2, mcArray);
2432  if (D0decay == kDecayD0toKpi) {
2433  decay = kDecayDStartoKpipi;
2434  }
2435  }
2436  }
2437  }
2438 
2439  return decay;
2440 }
2441 
2448 std::pair<AliAnalysisTaskDmesonJetsSub::EMesonOrigin_t, AliAODMCParticle*> AliAnalysisTaskDmesonJetsSub::AnalysisEngine::IsPromptCharm(const AliAODMCParticle* part, TClonesArray* mcArray)
2449 {
2450  std::pair<AliAnalysisTaskDmesonJetsSub::EMesonOrigin_t, AliAODMCParticle*> result(kUnknownQuark, 0);
2451 
2452  if (!part) return result;
2453  if (!mcArray) return result;
2454 
2455  static std::set<UInt_t> partons = { 4, 5 };
2456 
2457  AliAODMCParticle* parton = FindParticleOrigin(part, mcArray, kFindLast, partons);
2458  if (parton) {
2459  result.second = parton;
2460  UInt_t absPdgParton = TMath::Abs(parton->GetPdgCode());
2461  if (absPdgParton == 4) result.first = kFromCharm;
2462  else if (absPdgParton == 5) result.first = kFromBottom;
2463  }
2464 
2465  return result;
2466 }
2467 
2475 
2476 AliAODMCParticle* AliAnalysisTaskDmesonJetsSub::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode)
2477 {
2478  static std::set<UInt_t> pdgSet;
2479 
2480  return FindParticleOrigin(part, mcArray, mode, pdgSet);
2481 }
2482 
2491 
2492 AliAODMCParticle* AliAnalysisTaskDmesonJetsSub::AnalysisEngine::FindParticleOrigin(const AliAODMCParticle* part, TClonesArray* mcArray, EFindParticleOriginMode_t mode, const std::set<UInt_t>& pdgSet)
2493 {
2494  AliAODMCParticle* result = nullptr;
2495 
2496  Int_t mother = part->GetMother();
2497  while (mother >= 0) {
2498  AliAODMCParticle* mcGranma = static_cast<AliAODMCParticle*>(mcArray->At(mother));
2499  if (mcGranma) {
2500  UInt_t abspdgGranma = TMath::Abs(mcGranma->GetPdgCode());
2501 
2502  // If the current particle is one of the particle types that is being searched assign it to the result pointer
2503  if (pdgSet.empty() || pdgSet.count(abspdgGranma) > 0) {
2504  result = mcGranma;
2505  // If the last particle in the fragmentation tree (first when going reverse) was requested then stop the loop
2506  if (mode == kFindLast) break;
2507  }
2508  if (mother == mcGranma->GetMother()) { // avoid infinite loop!
2509  AliWarningClassStream() << "Particle " << mother << " (PDG=" << mcGranma->PdgCode() << ") is the mother of itself!?" << std::endl;
2510  break;
2511  }
2512  mother = mcGranma->GetMother();
2513  }
2514  else {
2515  AliErrorClassStream() << "Could not retrieve mother particle " << mother << "!" << std::endl;
2516  break;
2517  }
2518  }
2519 
2520  return result;
2521 }
2522 
2525 {
2526  for (auto& jetDef : fJetDefinitions) {
2527  jetDef.fJets.clear();
2528  }
2529 
2530  if (fMCMode == kMCTruth) {
2532  }
2533  else {
2535  }
2536 }
2537 
2540 {
2541  // Fill the vertex info of the candidates
2542  // Needed for reduced delta AOD, where the vertex info has been deleted
2543  // to reduce the delta AOD file size
2545 
2546  const Int_t nD = fCandidateArray->GetEntriesFast();
2547 
2548  AliDmesonJetInfo DmesonJet;
2549  DmesonJet.fEvent = this->fAodEvent;
2550 
2551  std::map<AliHFJetDefinition*,Double_t> maxJetPt;
2552  for (auto& def : fJetDefinitions) maxJetPt[&def] = 0;
2553  Double_t maxDPt = 0;
2554 
2555  std::array<int, 3> nAccCharm = {0};
2556  std::array<std::array<int, 3>, 5> nAccCharmPt = {{{0}}};
2557 
2558  for (Int_t icharm = 0; icharm < nD; icharm++) { //loop over D candidates
2559  AliAODRecoDecayHF2Prong* charmCand = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(icharm)); // D candidates
2560  if (!charmCand) continue;
2561  if(!(vHF.FillRecoCand(fAodEvent,charmCand))) continue;
2562 
2563  // region of interest + cuts
2564  if (!fRDHFCuts->IsInFiducialAcceptance(charmCand->Pt(), charmCand->Y(fCandidatePDG))) continue;
2565  Int_t nMassHypo = 0; // number of mass hypothesis accepted for this D meson
2566  if (charmCand->Pt() > maxDPt) maxDPt = charmCand->Pt();
2567  for (Int_t im = 0; im < 2; im++) { // 2 mass hypothesis (when available)
2568  DmesonJet.Reset();
2569  DmesonJet.fDmesonParticle = charmCand;
2570  DmesonJet.fSelectionType = im + 1;
2571  if (ExtractRecoDecayAttributes(charmCand, DmesonJet, im)) {
2572  for (auto& def : fJetDefinitions) {
2573  if (FindJet(charmCand, DmesonJet, def,im)) {
2574  Double_t jetPt = DmesonJet.fJets[def.GetName()].fMomentum.Pt();
2575  if (jetPt > maxJetPt[&def]) maxJetPt[&def] = jetPt;
2576  }
2577  else {
2578  AliWarning(Form("Could not find jet '%s' for D meson '%s': pT = %.3f, eta = %.3f, phi = %.3f",
2579  def.GetName(), GetName(), DmesonJet.fD.Pt(), DmesonJet.fD.Eta(), DmesonJet.fD.Phi_0_2pi()));
2580  }
2581  }
2582  fDmesonJets[(icharm+1)*(1-(im*2))] = DmesonJet;
2583  nMassHypo++;
2584  nAccCharm[im]++;
2585 
2586  for (int i = 0; i < nAccCharmPt.size(); i++) {
2587  if (charmCand->Pt() < i) break;
2588  nAccCharmPt[i][im]++;
2589  }
2590  }
2591  }
2592  if (nMassHypo == 2) { // both mass hypothesis accepted
2593  nAccCharm[0]--;
2594  nAccCharm[1]--;
2595  nAccCharm[2]++;
2596 
2597  for (int i = 0; i < nAccCharmPt.size(); i++) {
2598  if (charmCand->Pt() < i) break;
2599  nAccCharmPt[i][0]--;
2600  nAccCharmPt[i][1]--;
2601  nAccCharmPt[i][2]++;
2602  }
2603 
2604  fDmesonJets[(icharm+1)].fD0D0bar = kTRUE;
2605  fDmesonJets[-(icharm+1)].fD0D0bar = kTRUE;
2606  }
2607  } // end of D cand loop
2608 
2609  TString hname;
2610 
2611  Int_t ntracks = 0;
2612 
2613  for (auto track_cont : fTrackContainers) {
2614  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
2615  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(nullptr);
2616  ntracks += track_cont->GetNAcceptEntries();
2617  }
2618 
2619  for (auto& def : fJetDefinitions) {
2620  if (!def.fRho) continue;
2621  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
2622  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
2623 
2624  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
2625  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
2626 
2627  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
2628  fHistManager->FillTH2(hname, fEventInfo.fCent, def.fRho->GetVal());
2629 
2630  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
2631  fHistManager->FillTH2(hname, fEventInfo.fCent, maxJetPt[&def]);
2632 
2633  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
2634  fHistManager->FillTH2(hname, fEventInfo.fCent, maxDPt);
2635 
2636  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
2637  fHistManager->FillTH2(hname, ntracks, def.fRho->GetVal());
2638 
2639  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
2640  fHistManager->FillTH2(hname, ntracks, maxJetPt[&def]);
2641 
2642  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
2643  fHistManager->FillTH2(hname, ntracks, maxDPt);
2644  }
2645 
2646  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
2647  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
2648  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
2649  fHistManager->FillTH1(hname, "Both", nAccCharm[2]);
2650 
2651  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
2652  fHistManager->FillTH2(hname, ntracks, nAccCharm[0]+nAccCharm[1]+nAccCharm[2]);
2653 
2654  for (int i = 0; i < nAccCharmPt.size(); i++) {
2655  hname = TString::Format("%s/fHistNTotAcceptedDmesonsPt%d", GetName(), i);
2656  fHistManager->FillTH1(hname, "D", nAccCharmPt[i][0]);
2657  fHistManager->FillTH1(hname, "Anti-D", nAccCharmPt[i][1]);
2658  fHistManager->FillTH1(hname, "Both", nAccCharmPt[i][2]);
2659 
2660  hname = TString::Format("%s/fHistNAcceptedDmesonsPt%d", GetName(), i);
2661  fHistManager->FillTH1(hname, nAccCharmPt[i][0]+nAccCharmPt[i][1]+nAccCharmPt[i][2]);
2662  }
2663 
2664  hname = TString::Format("%s/fHistNDmesons", GetName());
2665  fHistManager->FillTH1(hname, nD);
2666 }
2667 
2678 {
2679  TString hname;
2680 
2681  Double_t rho = 0;
2682  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
2683 
2685  fFastJetWrapper->SetR(jetDef.fRadius);
2688 
2689  fFastJetWrapper->AddInputVector(DmesonJet.fD.Px(), DmesonJet.fD.Py(), DmesonJet.fD.Pz(), DmesonJet.fD.E(), 0);
2690 
2691  if (jetDef.fJetType != AliJetContainer::kNeutralJet) {
2692  for (auto track_cont : fTrackContainers) {
2693  AliHFTrackContainer* hftrack_cont = dynamic_cast<AliHFTrackContainer*>(track_cont);
2694  if (hftrack_cont) hftrack_cont->SetDMesonCandidate(Dcand);
2695  hname = TString::Format("%s/%s/fHistTrackRejectionReason", GetName(), jetDef.GetName());
2696  AddInputVectors(track_cont, 100, static_cast<TH2*>(fHistManager->FindObject(hname)), fTrackEfficiency);
2697 
2698  if (hftrack_cont) {
2699  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", GetName(), jetDef.GetName());
2700  TH1* histDaughterNotInJet = static_cast<TH1*>(fHistManager->FindObject(hname));
2701  const TObjArray& daughters = hftrack_cont->GetDaughterList();
2702  for (Int_t i = 0; i < daughters.GetEntriesFast(); i++) {
2703  AliVParticle* daughter = static_cast<AliVParticle*>(daughters.At(i));
2704  if (!hftrack_cont->GetArray()->FindObject(daughter)) histDaughterNotInJet->Fill(daughter->Pt());
2705  }
2706  }
2707  }
2708  }
2709 
2710  if (jetDef.fJetType != AliJetContainer::kChargedJet) {
2711  for (auto clus_cont : fClusterContainers) {
2712  hname = TString::Format("%s/%s/fHistClusterRejectionReason", GetName(), jetDef.GetName());
2713  AddInputVectors(clus_cont, -100, static_cast<TH2*>(fHistManager->FindObject(hname)));
2714  }
2715  }
2716 
2717  // run jet finder
2718  fFastJetWrapper->Run();
2719 
2720  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
2721 
2722  for (UInt_t ijet = 0; ijet < jets_incl.size(); ++ijet) {
2723  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper->GetJetConstituents(ijet));
2724 
2725  Bool_t isDmesonJet = kFALSE;
2726 
2727  Double_t maxChPt = 0;
2728  Double_t maxNePt = 0;
2729  Double_t totalNeutralPt = 0;
2730  Int_t nConst = 1;
2731 
2732  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
2733  if (constituents[ic].user_index() == 0) {
2734  isDmesonJet = kTRUE;
2735  }
2736  else if (constituents[ic].user_index() >= 100) {
2737  if (constituents[ic].pt() > maxChPt) maxChPt = constituents[ic].pt();
2738  nConst++;
2739  }
2740  else if (constituents[ic].user_index() <= -100) {
2741  totalNeutralPt += constituents[ic].pt();
2742  if (constituents[ic].pt() > maxNePt) maxChPt = constituents[ic].pt();
2743  nConst++;
2744  }
2745  }
2746 
2747  if (isDmesonJet) {
2748  DmesonJet.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jets_incl[ijet].px(), jets_incl[ijet].py(), jets_incl[ijet].pz(), jets_incl[ijet].E());
2749  DmesonJet.fJets[jetDef.GetName()].fNConstituents = nConst;
2750  DmesonJet.fJets[jetDef.GetName()].fMaxChargedPt = maxChPt;
2751  DmesonJet.fJets[jetDef.GetName()].fMaxNeutralPt = maxNePt;
2752  DmesonJet.fJets[jetDef.GetName()].fNEF = totalNeutralPt / jets_incl[ijet].pt();
2753  DmesonJet.fJets[jetDef.GetName()].fArea = jets_incl[ijet].area();
2754  DmesonJet.fJets[jetDef.GetName()].fCorrPt = DmesonJet.fJets[jetDef.GetName()].fMomentum.Pt() - jets_incl[ijet].area() * rho;
2755  IterativeDeclustering(ijet,1,jetDef, DmesonJet.fD.M());
2756  return kTRUE;
2757  }
2758  if(!isDmesonJet && numcand!=1) IterativeDeclustering(ijet,0,jetDef,0.);
2759  }
2760 
2761  return kFALSE;
2762 }
2763 
2767 void AliAnalysisTaskDmesonJetsSub::AnalysisEngine::AddInputVectors(AliEmcalContainer* cont, Int_t offset, TH2* rejectHist, Double_t eff)
2768 {
2769  auto itcont = cont->all_momentum();
2770  for (AliEmcalIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
2771  UInt_t rejectionReason = 0;
2772  if (!cont->AcceptObject(it.current_index(), rejectionReason)) {
2773  if (rejectHist) rejectHist->Fill(AliEmcalContainer::GetRejectionReasonBitPosition(rejectionReason), it->first.Pt());
2774  continue;
2775  }
2776  if (fRandomGen && eff > 0 && eff < 1) {
2777  Double_t rnd = fRandomGen->Rndm();
2778  if (eff < rnd) {
2779  if (rejectHist) rejectHist->Fill(6, it->first.Pt());
2780  continue;
2781  }
2782  }
2783  Int_t uid = offset >= 0 ? it.current_index() + offset: -it.current_index() - offset;
2784  fFastJetWrapper->AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
2785  }
2786 }
2787 
2789 {
2790  double nsd = 0;
2791  double nall = 0;
2792  double zg = 0.;
2793  double flagSubjet=0;
2794  double xconstperp=0;
2795  TString hname;
2796  fastjet::JetAlgorithm jet_algo(fastjet::cambridge_algorithm);
2797  double jet_radius_ca = 1.0;
2798  fastjet::JetDefinition jet_def(jet_algo, jet_radius_ca,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
2799 
2800  try{
2801  std::vector<fastjet::PseudoJet> particles(fFastJetWrapper->GetJetConstituents(ijet));
2802  fastjet::ClusterSequence cs_ca(particles, jet_def);
2803  std::vector<fastjet::PseudoJet> output_jets = cs_ca.inclusive_jets(0);
2804  output_jets = sorted_by_pt(output_jets);
2805 
2806 
2807  fastjet::PseudoJet jj = output_jets[0];
2808  fastjet::PseudoJet j1;
2809  fastjet::PseudoJet j2;
2810 
2811  while(jj.has_parents(j1,j2)){
2812  nall = nall + 1;
2813  if(j1.perp() < j2.perp()) std::swap(j1,j2);
2814  flagSubjet=0;
2815  vector < fastjet::PseudoJet > constitj1 = sorted_by_pt(j1.constituents());
2816  if(type==1){
2817  for(Int_t j=0;j<constitj1.size();j++){
2818  if(constitj1[j].user_index()==0){
2819  xconstperp=constitj1[j].perp();
2820  flagSubjet=1; }}}
2821 
2822 
2823  double delta_R = j1.delta_R(j2);
2824  zg = j2.perp()/(j1.perp()+j2.perp());
2825  double yh=j1.e()+j2.e();
2826  double y = log(1.0/delta_R);
2827  double lnpt_rel = log(j2.perp()*delta_R);
2828  double frac=j1.perp()/output_jets[0].perp();
2829  double lundEntries[10] = {y, lnpt_rel, output_jets[0].perp(), nall, type, flagSubjet, xconstperp, invmass,yh,TMath::Abs(output_jets[0].eta())};
2830 
2831  hname = TString::Format("%s/%s/LundIterative", GetName(), jetDef.GetName());
2832  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
2833 
2834  h->Fill(lundEntries);
2835  jj=j1;
2836  }
2837 
2838 
2839  } catch (fastjet::Error) { /*return -1;*/ }
2840 
2841 }
2842 
2845 {
2846  TString hname;
2847 
2852 
2853  if (!fMCContainer->IsSpecialPDGFound()) return;
2854 
2855  std::array<int,2> nAccCharm = {0};
2856  std::array<std::array<int, 2>, 5> nAccCharmPt = {{{0}}};
2857 
2858  std::map<AliHFJetDefinition*, Double_t> maxJetPt;
2859  Double_t maxDPt = 0;
2860 
2861  for (auto &jetDef : fJetDefinitions) {
2862  maxJetPt[&jetDef] = 0;
2863  Double_t rho = 0;
2864  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
2865  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", GetName(), jetDef.GetName());
2866  TH1* histNDmesonsVsNconstituents = static_cast<TH1*>(fHistManager->FindObject(hname));
2867 
2868  switch (jetDef.fJetType) {
2870  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
2871  break;
2873  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
2874  break;
2876  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
2877  break;
2878  }
2879 
2881  fFastJetWrapper->SetR(jetDef.fRadius);
2884 
2885  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
2886  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
2887 
2888  fFastJetWrapper->Run();
2889 
2890  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
2891 
2892  for (auto jet : jets_incl) {
2893  Int_t nDmesonsInJet = 0;
2894 
2895  for (auto constituent : jet.constituents()) {
2896  Int_t iPart = constituent.user_index() - 100;
2897  if (constituent.perp() < 1e-6) continue; // reject ghost particles
2898  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
2899  if (!part) {
2900  ::Error("AliAnalysisTaskDmesonJetsSub::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
2901  continue;
2902  }
2903  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
2904  nDmesonsInJet++;
2905  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
2906  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
2907  if (part->Pt() > maxDPt) maxDPt = part->Pt();
2908  std::pair<int, AliDmesonJetInfo> element;
2909  element.first = iPart;
2910  dMesonJetIt = fDmesonJets.insert(element).first;
2911  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
2912  (*dMesonJetIt).second.fDmesonParticle = part;
2913  (*dMesonJetIt).second.fSelectionType = part->PdgCode() > 0 ? 1 : 2;
2914 
2915  UShort_t p = 0;
2916  UInt_t rs = 0;
2917 
2918  auto origin = IsPromptCharm(part, fMCContainer->GetArray());
2919  p = 0;
2920  rs = origin.first;
2921  while (rs >>= 1) { p++; }
2922  (*dMesonJetIt).second.fPartonType = p;
2923  (*dMesonJetIt).second.fParton = origin.second;
2924 
2925  (*dMesonJetIt).second.fAncestor = FindParticleOrigin(part, fMCContainer->GetArray(), kFindFirst);
2926 
2927  Int_t im = -1;
2928  if (part->PdgCode() > 0) { // D0
2929  im = 0;
2930  }
2931  else { // D0bar
2932  im = 1;
2933  }
2934 
2935  nAccCharm[im]++;
2936  for (int i = 0; i < nAccCharmPt.size(); i++) {
2937  if (part->Pt() < i) break;
2938  nAccCharmPt[i][im]++;
2939  }
2940  }
2941 
2942  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
2943  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = jet.constituents().size();
2944  (*dMesonJetIt).second.fJets[jetDef.GetName()].fArea = jet.area();
2945  (*dMesonJetIt).second.fJets[jetDef.GetName()].fCorrPt = (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.Pt() - jet.area() * rho;
2946  if (jet.perp() > maxJetPt[&jetDef]) maxJetPt[&jetDef] = jet.perp();
2947  } // if constituent is a D meson
2948  } // for each constituent
2949  if (nDmesonsInJet > 0) histNDmesonsVsNconstituents->Fill(jet.constituents().size(), nDmesonsInJet);
2950  } // for each jet
2951  } // for each jet definition
2952 
2954 
2955  for (auto& def : fJetDefinitions) {
2956  if (!def.fRho) continue;
2957  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
2958  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
2959 
2960  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
2961  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
2962 
2963  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
2964  fHistManager->FillTH2(hname, fEventInfo.fCent, def.fRho->GetVal());
2965 
2966  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
2967  fHistManager->FillTH2(hname, fEventInfo.fCent, maxJetPt[&def]);
2968 
2969  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
2970  fHistManager->FillTH2(hname, fEventInfo.fCent, maxDPt);
2971 
2972  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
2973  fHistManager->FillTH2(hname, npart, def.fRho->GetVal());
2974 
2975  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
2976  fHistManager->FillTH2(hname, npart, maxJetPt[&def]);
2977 
2978  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
2979  fHistManager->FillTH2(hname, npart, maxDPt);
2980  }
2981 
2982  if (fDmesonJets.size() != nAccCharm[0]+nAccCharm[1]) AliError(Form("I found %lu mesons (%d)?", fDmesonJets.size(), nAccCharm[0]+nAccCharm[1]));
2983  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
2984  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
2985  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
2986 
2987  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
2988  fHistManager->FillTH2(hname, npart, nAccCharm[0]+nAccCharm[1]);
2989 
2990  for (int i = 0; i < nAccCharmPt.size(); i++) {
2991  hname = TString::Format("%s/fHistNTotAcceptedDmesonsPt%d", GetName(), i);
2992  fHistManager->FillTH1(hname, "D", nAccCharmPt[i][0]);
2993  fHistManager->FillTH1(hname, "Anti-D", nAccCharmPt[i][1]);
2994 
2995  hname = TString::Format("%s/fHistNAcceptedDmesonsPt%d", GetName(), i);
2996  fHistManager->FillTH1(hname, nAccCharmPt[i][0]+nAccCharmPt[i][1]);
2997  }
2998 
2999  hname = TString::Format("%s/fHistNDmesons", GetName());
3000  fHistManager->FillTH1(hname, nAccCharm[0]+nAccCharm[1]); // same as the number of accepted D mesons, since no selection is performed
3001 }
3002 
3003 
3004 
3005 // Definitions of class AliAnalysisTaskDmesonJetsSub
3006 
3010 
3014  fAnalysisEngines(),
3015  fEnabledAxis(0),
3017  fHistManager(),
3018  fApplyKinematicCuts(kTRUE),
3019  fNOutputTrees(0),
3020  fTrackEfficiency(0),
3021  fRejectISR(kFALSE),
3022  fJetAreaType(fastjet::active_area),
3023  fJetGhostArea(0.005),
3024  fMCContainer(0),
3025  fAodEvent(0),
3026  fFastJetWrapper(0)
3027 {
3028  SetMakeGeneralHistograms(kTRUE);
3029 }
3030 
3035  AliAnalysisTaskEmcalLight(name, kTRUE),
3036  fAnalysisEngines(),
3039  fHistManager(name),
3040  fApplyKinematicCuts(kTRUE),
3041  fNOutputTrees(nOutputTrees),
3042  fTrackEfficiency(0),
3043  fRejectISR(kFALSE),
3044  fJetAreaType(fastjet::active_area),
3045  fJetGhostArea(0.005),
3046  fMCContainer(0),
3047  fAodEvent(0),
3048  fFastJetWrapper(0)
3049 {
3050  SetMakeGeneralHistograms(kTRUE);
3051  for (Int_t i = 0; i < nOutputTrees; i++){
3052  DefineOutput(2+i, TTree::Class());
3053  }
3054 }
3055 
3058 {
3059  if (fFastJetWrapper) delete fFastJetWrapper;
3060 }
3061 
3069 {
3070  AliRDHFCuts* analysiscuts = 0;
3071  TFile* filecuts = TFile::Open(cutfname);
3072  if (!filecuts || filecuts->IsZombie()) {
3073  ::Error("AliAnalysisTaskDmesonJetsSub::LoadDMesonCutsFromFile", "Input file not found: will use std cuts.");
3074  filecuts = 0;
3075  }
3076 
3077  if (filecuts) analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
3078 
3079  if (!analysiscuts) {
3080  ::Error("AliAnalysisTaskDmesonJetsSub::LoadDMesonCutsFromFile", "Could not find analysis cuts '%s' in '%s'.", cutsname.Data(), cutfname.Data());
3081  if (filecuts) {
3082  filecuts->ls();
3083  }
3084  }
3085  else {
3086  ::Info("AliAnalysisTaskDmesonJetsSub::LoadDMesonCutsFromFile", "Cuts '%s' loaded from file '%s'", cutsname.Data(), cutfname.Data());
3087  }
3088 
3089  return analysiscuts;
3090 }
3091 
3104 {
3106  return AddAnalysisEngine(type, cutfname, cuttype, MCmode, jetDef, rhoName);
3107 }
3108 
3120 {
3121  AliRDHFCuts* cuts = 0;
3122 
3123  if (!cutfname.IsNull()) {
3124  TString cutsname;
3125 
3126  switch (type) {
3127  case kD0toKpi:
3128  case kD0toKpiLikeSign:
3129  cutsname = "D0toKpiCuts";
3130  break;
3131  case kDstartoKpipi:
3132  cutsname = "DStartoKpipiCuts";
3133  break;
3134  default:
3135  return 0;
3136  }
3137 
3138  if (!cuttype.IsNull()) {
3139  cutsname += TString::Format("_%s", cuttype.Data());
3140  }
3141 
3142  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
3143  if (cuts) cuts->PrintAll();
3144  }
3145 
3146  AnalysisEngine eng(type, MCmode, cuts);
3147 
3148  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
3149 
3150  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
3151  eng.AddJetDefinition(jetDef);
3152  it = fAnalysisEngines.insert(it, eng);
3153  ::Info("AliAnalysisTaskDmesonJetsSub::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
3154  }
3155  else {
3156  AnalysisEngine* found_eng = &(*it);
3157  ::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());
3158  found_eng->AddJetDefinition(jetDef);
3159 
3160  if (cuts) {
3161  if (found_eng->fRDHFCuts != 0) ::Warning("AliAnalysisTaskDmesonJetsSub::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
3162  found_eng->SetRDHFCuts(cuts);
3163  }
3164  }
3165 
3166  return &(*it);
3167 }
3168 
3169 std::list<AliAnalysisTaskDmesonJetsSub::AnalysisEngine>::iterator AliAnalysisTaskDmesonJetsSub::FindAnalysisEngine(const AliAnalysisTaskDmesonJetsSub::AnalysisEngine& eng)
3170 {
3171  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
3172  while (it != fAnalysisEngines.end() && (*it) != eng) it++;
3173  return it;
3174 }
3175 
3178 {
3179  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
3180 
3182 
3183  // Define histograms
3184  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
3185 
3186  TString hname;
3187  TString htitle;
3188  TH1* h = 0;
3189  Int_t treeSlot = 0;
3190 
3191  Int_t maxTracks = 6000;
3192  Double_t maxRho = 500;
3193  if (fForceBeamType == kpp) {
3194  maxRho = 50;
3195  maxTracks = 200;
3196  }
3197  else if (fForceBeamType == kpA) {
3198  maxRho = 200;
3199  maxTracks = 500;
3200  }
3201 
3202  Int_t dimx = 10;
3203  Int_t nbinsx[10] = {50,100,10,20,2,2,20,100,100,9};
3204  Double_t minx[10] = {0,-10,0,0,0,0,0,1.6,0,0};
3205  Double_t maxx[10] = {5,10,100,20,2,2,20,2,100,0.9};
3206  TString titlex[10]={"log(1/deltaR)","log(zteta)","jet pt","n","type","flagSubjet","ptD","invmass","frac","abs(eta)"};
3207 
3208 
3209 
3210 
3211  hname = "fHistCharmPt";
3212  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
3213  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3214 
3215  hname = "fHistCharmEta";
3216  htitle = hname + ";#eta_{charm};counts";
3217  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
3218 
3219  hname = "fHistCharmPhi";
3220  htitle = hname + ";#phi_{charm};counts";
3221  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
3222 
3223  hname = "fHistBottomPt";
3224  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
3225  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3226 
3227  hname = "fHistBottomEta";
3228  htitle = hname + ";#eta_{bottom};counts";
3229  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
3230 
3231  hname = "fHistBottomPhi";
3232  htitle = hname + ";#phi_{bottom};counts";
3233  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
3234 
3235  hname = "fHistHighestPartonPt";
3236  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
3237  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3238 
3239  hname = "fHistHighestPartonType";
3240  htitle = hname + ";type;counts";
3241  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
3242 
3243  hname = "fHistNHeavyQuarks";
3244  htitle = hname + ";number of heavy-quarks;counts";
3245  fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
3246 
3247  ::Info("AliAnalysisTaskDmesonJetsSub::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
3248  for (auto &param : fAnalysisEngines) {
3249  ::Info("AliAnalysisTaskDmesonJetsSub::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
3250 
3251  param.fHistManager = &fHistManager;
3252 
3253  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", param.GetName());
3254  htitle = hname + ";#it{N}_{tracks};#it{N}_{D};events";
3255  h = fHistManager.CreateTH2(hname, htitle, 251, -0.5, 250.5, 21, -0.5, 20.5);
3256 
3257  hname = TString::Format("%s/fHistNTotAcceptedDmesons", param.GetName());
3258  htitle = hname + ";;#it{N}_{D}";
3259  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
3260 
3261  for (int i = 0 ; i < 5; i++) {
3262  hname = TString::Format("%s/fHistNAcceptedDmesonsPt%d", param.GetName(), i);
3263  htitle = hname + ";#it{N}_{D};events";
3264  h = fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
3265 
3266  hname = TString::Format("%s/fHistNTotAcceptedDmesonsPt%d", param.GetName(), i);
3267  htitle = hname + ";;#it{N}_{D}";
3268  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
3269  }
3270 
3271  hname = TString::Format("%s/fHistNDmesons", param.GetName());
3272  htitle = hname + ";#it{N}_{D};events";
3273  h = fHistManager.CreateTH1(hname, htitle, 501, -0.5, 500.5);
3274 
3275  hname = TString::Format("%s/fHistNEvents", param.GetName());
3276  htitle = hname + ";Event status;counts";
3277  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
3278  h->GetXaxis()->SetBinLabel(1, "Accepted");
3279  h->GetXaxis()->SetBinLabel(2, "Rejected");
3280 
3281  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
3282  htitle = hname + ";Rejection reason;counts";
3283  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
3284 
3285  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
3286  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
3287  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
3288 
3289  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
3290  htitle = hname + ";#it{#eta}_{D};counts";
3291  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
3292 
3293  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
3294  htitle = hname + ";#it{#phi}_{D};counts";
3295  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
3296 
3297  if (param.fMCMode != kMCTruth) {
3298  if (param.fCandidateType == kD0toKpi || param.fCandidateType == kD0toKpiLikeSign) {
3299  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
3300  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
3301  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
3302  }
3303  else if (param.fCandidateType == kDstartoKpipi) {
3304  Double_t min = 0;
3305  Double_t max = 0;
3306 
3307  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
3308  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
3309  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
3310  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
3311 
3312  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
3313  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
3314  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
3315  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
3316  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
3317  }
3318  }
3319 
3320  if (param.fMCMode == kMCTruth) {
3321  hname = TString::Format("%s/fHistPartonPt", param.GetName());
3322  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
3323  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3324 
3325  hname = TString::Format("%s/fHistPartonEta", param.GetName());
3326  htitle = hname + ";#eta_{parton};counts";
3327  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
3328 
3329  hname = TString::Format("%s/fHistPartonPhi", param.GetName());
3330  htitle = hname + ";#phi_{parton};counts";
3331  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
3332 
3333  hname = TString::Format("%s/fHistPartonType", param.GetName());
3334  htitle = hname + ";type;counts";
3335  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
3336 
3337  hname = TString::Format("%s/fHistPrompt", param.GetName());
3338  htitle = hname + ";Type;counts";
3339  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
3340  h->GetXaxis()->SetBinLabel(1, "Unknown");
3341  h->GetXaxis()->SetBinLabel(2, "Prompt");
3342  h->GetXaxis()->SetBinLabel(3, "Non-Prompt");
3343 
3344  hname = TString::Format("%s/fHistAncestor", param.GetName());
3345  htitle = hname + ";Ancestor;counts";
3346  h = fHistManager.CreateTH1(hname, htitle, 4, 0, 4);
3347  h->GetXaxis()->SetBinLabel(1, "Unknown");
3348  h->GetXaxis()->SetBinLabel(2, "Charm");
3349  h->GetXaxis()->SetBinLabel(3, "Bottom");
3350  h->GetXaxis()->SetBinLabel(4, "Proton");
3351  }
3352 
3353  for (auto& jetDef : param.fJetDefinitions) {
3354  ::Info("AliAnalysisTaskDmesonJetsSub::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef.GetName());
3355 
3356  if (param.fMCMode == kMCTruth) {
3357  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", param.GetName(), jetDef.GetName());
3358  htitle = hname + ";#it{N}_{constituents};#it{N}_{D};counts";
3359  h = fHistManager.CreateTH2(hname, htitle, 51, -0.5, 50.5, 10, 0.5, 10.5);
3360  }
3361 
3362  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef.GetName());
3363  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
3364  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
3365  SetRejectionReasonLabels(h->GetXaxis());
3366 
3367  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef.GetName());
3368  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
3369  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
3370  SetRejectionReasonLabels(h->GetXaxis());
3371 
3372  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef.GetName());
3373  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
3374  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
3375  SetRejectionReasonLabels(h->GetXaxis());
3376 
3377  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef.GetName());
3378  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
3379  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
3380 
3381  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef.GetName());
3382  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
3383  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
3384 
3385  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef.GetName());
3386  htitle = hname + ";#it{#eta}_{jet};counts";
3387  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
3388 
3389  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef.GetName());
3390  htitle = hname + ";#it{#phi}_{jet};counts";
3391  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
3392 
3393  if (!jetDef.fRhoName.IsNull()) {
3394  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", param.GetName(), jetDef.GetName());
3395  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});#rho (GeV/#it{c} #times rad^{-1});counts";
3396  fHistManager.CreateTH2(hname, htitle, 300, 0, 150, 1000, 0, maxRho);
3397 
3398  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", param.GetName(), jetDef.GetName());
3399  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});#rho (GeV/#it{c} #times rad^{-1});counts";
3400  fHistManager.CreateTH2(hname, htitle, 300, 0, 150, 1000, 0, maxRho);
3401 
3402  hname = TString::Format("%s/%s/fHistRhoVsCent", param.GetName(), jetDef.GetName());
3403  htitle = hname + ";Centrality (%);#rho (GeV/#it{c} #times rad^{-1});counts";
3404  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 1000, 0, maxRho);
3405 
3406  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", param.GetName(), jetDef.GetName());
3407  htitle = hname + ";Centrality (%);#it{p}_{T,jet} (GeV/#it{c});counts";
3408  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
3409 
3410  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", param.GetName(), jetDef.GetName());
3411  htitle = hname + ";Centrality (%);#it{p}_{T,D} (GeV/#it{c});counts";
3412  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
3413 
3414  hname = TString::Format("%s/%s/fHistRhoVsNTracks", param.GetName(), jetDef.GetName());
3415  htitle = hname + ";no. of tracks;#rho (GeV/#it{c} #times rad^{-1});counts";
3416  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 1000, 0, maxRho);
3417 
3418  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", param.GetName(), jetDef.GetName());
3419  htitle = hname + ";no. of tracks;#it{p}_{T,jet} (GeV/#it{c});counts";
3420  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
3421 
3422  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", param.GetName(), jetDef.GetName());
3423  htitle = hname + ";no. of tracks;#it{p}_{T,D} (GeV/#it{c});counts";
3424  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
3425  }
3426  hname = TString::Format("%s/%s/LundIterative",param.GetName(),jetDef.GetName());
3427  cout<<"at the begining"<<hname<<endl;
3428  THnSparse* h = fHistManager.CreateTHnSparse(hname,hname,dimx,nbinsx,minx,maxx);
3429  for (Int_t j = 0; j < dimx; j++) {
3430  h->GetAxis(j)->SetTitle(titlex[j]);}
3431 
3432  }
3433  switch (fOutputType) {
3434  case kTreeOutput:
3435  {
3436  OutputHandlerTTree* tree_handler = new OutputHandlerTTree(&param);
3437  param.fOutputHandler = tree_handler;
3438  tree_handler->BuildOutputObject(GetName());
3439  if (treeSlot < fNOutputTrees) {
3440  tree_handler->AssignDataSlot(treeSlot+2);
3441  treeSlot++;
3442  PostDataFromAnalysisEngine(tree_handler);
3443  }
3444  else {
3445  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
3446  }
3447  }
3448  break;
3449  case kTHnOutput:
3450  {
3451  OutputHandlerTHnSparse* thnsparse_handler = new OutputHandlerTHnSparse(&param);
3452  param.fOutputHandler = thnsparse_handler;
3453  thnsparse_handler->SetEnabledAxis(fEnabledAxis);
3454  thnsparse_handler->BuildOutputObject(GetName());
3455  }
3456  break;
3457  case kOnlyQAOutput:
3458  {
3459  OutputHandler* qa_handler = new OutputHandler(&param);
3460  param.fOutputHandler = qa_handler;
3461  }
3462  break;
3463  case kNoOutput:
3464  break;
3465  case kTreeExtendedOutput:
3466  {
3468  param.fOutputHandler = tree_handler;
3469  tree_handler->BuildOutputObject(GetName());
3470  if (treeSlot < fNOutputTrees) {
3471  tree_handler->AssignDataSlot(treeSlot+2);
3472  treeSlot++;
3473  PostDataFromAnalysisEngine(tree_handler);
3474  }
3475  else {
3476  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
3477  }
3478  }
3479  break;
3480  }
3481  }
3482 
3484 
3485  PostData(1, fOutput);
3486 }
3487 
3491 {
3493 
3494  // Load the event
3495  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
3496 
3497  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
3498 
3499  fFastJetWrapper->SetAreaType((fastjet::AreaType)fJetAreaType);
3501 
3502  if (!fAodEvent) {
3503  AliError(Form("This task need an AOD event (Task '%s'). Expect troubles...", GetName()));
3504  //return;
3505  }
3506 
3507  TRandom* rnd = 0;
3508  if (fTrackEfficiency > 0 && fTrackEfficiency < 1) rnd = new TRandom3(0);
3509 
3510  for (auto cont_it : fParticleCollArray) {
3511  AliHFAODMCParticleContainer* part_cont = dynamic_cast<AliHFAODMCParticleContainer*>(cont_it.second);
3512  if (part_cont) fMCContainer = part_cont;
3513  }
3514 
3515  for (auto &params : fAnalysisEngines) {
3516 
3517  params.fAodEvent = fAodEvent;
3518  params.fFastJetWrapper = fFastJetWrapper;
3519  params.fTrackEfficiency = fTrackEfficiency;
3520  params.fRejectISR = fRejectISR;
3521  params.fRandomGen = rnd;
3522 
3523  for (auto &jetdef: params.fJetDefinitions) {
3524  if (!jetdef.fRhoName.IsNull()) {
3525  jetdef.fRho = dynamic_cast<AliRhoParameter*>(fInputEvent->FindListObject(jetdef.fRhoName));
3526  if (!jetdef.fRho) {
3527  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3528  "%s: Could not find rho object '%s' for engine '%s'",
3529  GetName(), jetdef.fRhoName.Data(), params.GetName());
3530  }
3531  }
3532  }
3533 
3534  if (!params.fRDHFCuts) {
3535  if (params.fMCMode == kMCTruth) {
3536  ::Warning("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3537  "%s: RDHF cuts not provided for engine '%s'.",
3538  GetName(), params.GetName());
3539  }
3540  else {
3541  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3542  "%s: RDHF cuts not provided. Engine '%s' disabled.",
3543  GetName(), params.GetName());
3544  params.fInhibit = kTRUE;
3545  }
3546  }
3547 
3548  params.fMCContainer = fMCContainer;
3549 
3550  for (auto cont_it : fParticleCollArray) {
3551  AliTrackContainer* track_cont = dynamic_cast<AliTrackContainer*>(cont_it.second);
3552  if (track_cont) params.fTrackContainers.push_back(track_cont);
3553  }
3554 
3555  for (auto cont_it : fClusterCollArray) params.fClusterContainers.push_back(cont_it.second);
3556 
3557  if (fAodEvent) params.Init(fGeom, fAodEvent->GetRunNumber());
3558 
3559  if (params.fMCMode != kMCTruth && fAodEvent) {
3560  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
3561 
3562  if (params.fCandidateArray) {
3563  TString className;
3564  if (params.fCandidateType == kD0toKpi || params.fCandidateType == kD0toKpiLikeSign) {
3565  className = "AliAODRecoDecayHF2Prong";
3566  }
3567  else if (params.fCandidateType == kDstartoKpipi) {
3568  className = "AliAODRecoCascadeHF";
3569  }
3570  if (!params.fCandidateArray->GetClass()->InheritsFrom(className)) {
3571  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3572  "%s: Objects of type %s in %s are not inherited from %s! Task will be disabled!",
3573  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName(), className.Data()); // @suppress("Ambiguous problem")
3574  params.fCandidateArray = 0;
3575  params.fInhibit = kTRUE;
3576  }
3577  }
3578  else {
3579  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3580  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
3581  params.fBranchName.Data(), params.GetName());
3582  params.fInhibit = kTRUE;
3583  }
3584  }
3585 
3586  if (params.fMCMode != kNoMC) {
3587  if (!params.fMCContainer) {
3588  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3589  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
3590  params.GetName());
3591  params.fInhibit = kTRUE;
3592  }
3593  }
3594 
3595  if (params.fMCMode != kMCTruth) {
3596  if (params.fTrackContainers.size() == 0 && params.fClusterContainers.size() == 0) {
3597  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3598  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
3599  params.GetName());
3600  params.fInhibit = kTRUE;
3601  }
3602  }
3603  }
3604 }
3605 
3610 {
3611  TString hname;
3612 
3613  // Fix for temporary bug in ESDfilter
3614  // The AODs with null vertex pointer didn't pass the PhysSel
3615  // Now adding an entry in the histogram so as to check that this is actually cutting anything out
3616  if (fAodEvent && (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001)) {
3617  for (auto &eng : fAnalysisEngines) {
3618  if (eng.fInhibit) continue;
3619  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
3620  fHistManager.FillTH1(hname, "ESDfilterBug");
3621  }
3622  return kFALSE;
3623  }
3624 
3625  if (fAodEvent) {
3626  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
3627  if (matchingAODdeltaAODlevel <= 0) {
3628  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
3629  for (auto &eng : fAnalysisEngines) {
3630  if (eng.fInhibit) continue;
3631  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
3632  fHistManager.FillTH1(hname, "MismatchDeltaAOD");
3633  }
3634  return kFALSE;
3635  }
3636  }
3637 
3638  for (auto &eng : fAnalysisEngines) {
3639  eng.fDmesonJets.clear();
3640  if (eng.fInhibit) continue;
3641 
3642  eng.fEventInfo = EventInfo(fCent, fEPV0, fEventWeight, fPtHard);
3643 
3644  //Event selection
3645  hname = TString::Format("%s/fHistNEvents", eng.GetName());
3646  if (fAodEvent) {
3647  Bool_t iseventselected = kTRUE;
3648  if (eng.fRDHFCuts) iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
3649  if (!iseventselected) {
3650  fHistManager.FillTH1(hname, "Rejected");
3651  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
3652  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
3653  TString label;
3654  do {
3655  label = GetHFEventRejectionReasonLabel(bitmap);
3656  if (label.IsNull()) break;
3657  fHistManager.FillTH1(hname, label);
3658  } while (true);
3659 
3660  continue;
3661  }
3662  }
3663 
3664  fHistManager.FillTH1(hname, "Accepted");
3665 
3666  AliDebug(2, "Event selected");
3667 
3668  eng.RunAnalysis();
3669  }
3670  return kTRUE;
3671 }
3672 
3677 {
3678  for (auto &param : fAnalysisEngines) {
3679  if (param.fInhibit) continue;
3680  param.fOutputHandler->FillOutput(fApplyKinematicCuts);
3681  PostDataFromAnalysisEngine(param.fOutputHandler);
3682  }
3684  return kTRUE;
3685 }
3686 
3689 {
3690  auto itcont = fMCContainer->all_momentum();
3691  Int_t nHQ = 0;
3692  Double_t highestPartonPt = 0;
3693  Int_t PdgHighParton = 0;
3694  for (auto it = itcont.begin(); it != itcont.end(); it++) {
3695  auto part = *it;
3696  if (part.first.Pt() == 0) continue;
3697 
3698  Int_t PdgCode = part.second->GetPdgCode();
3699 
3700  // Skip all particles that are not either quarks or gluons
3701  if ((PdgCode < -9 || PdgCode > 9) && PdgCode != 21 && PdgCode != -21) continue;
3702 
3703  AliDebugStream(5) << "Parton " << it.current_index() <<
3704  " with pdg=" << PdgCode <<
3705  ", px=" << part.first.Px() <<
3706  ", py=" << part.first.Py() <<
3707  ", pz=" << part.first.Pz() <<
3708  ", n daughters = " << part.second->GetNDaughters() <<
3709  std::endl;
3710 
3711  // Skip partons that do not have any children
3712  // Unclear how this can happen, it would seem that this parton were not fragmented by the generator
3713  if (part.second->GetNDaughters() == 0) continue;
3714 
3715  // Look for highest momentum parton
3716  if (highestPartonPt < part.first.Pt()) {
3717  highestPartonPt = part.first.Pt();
3718  PdgHighParton = PdgCode;
3719  }
3720 
3721  // Skip partons that are not HF
3722  if (PdgCode != 4 && PdgCode != 5 && PdgCode != -4 && PdgCode != -5) continue;
3723 
3724  Bool_t lastInPartonShower = kTRUE;
3725  Bool_t hadronDaughter = kFALSE;
3726  for (Int_t daughterIndex = part.second->GetFirstDaughter(); daughterIndex <= part.second->GetLastDaughter(); daughterIndex++){
3727  if (daughterIndex < 0) {
3728  AliDebugStream(5) << "Could not find daughter index!" << std::endl;
3729  continue;
3730  }
3731  AliAODMCParticle *daughter = fMCContainer->GetMCParticle(daughterIndex);
3732  if (!daughter) {
3733  AliDebugStream(5) << "Could not find particle with index " << daughterIndex << "!" << std::endl;
3734  continue;
3735  }
3736  Int_t daughterPdgCode = daughter->GetPdgCode();
3737  if (daughter->GetMother() != it.current_index()) {
3738  AliDebugStream(5) << "Particle " << daughterIndex << " with pdg=" << daughterPdgCode <<
3739  ", px=" << daughter->Px() <<
3740  ", py=" << daughter->Py() <<
3741  ", pz=" << daughter->Pz() <<
3742  ", is not a daughter of " << it.current_index() <<
3743  "!" << std::endl;
3744  continue;
3745  }
3746 
3747  AliDebugStream(5) << "Found daughter " << daughterIndex <<
3748  " with pdg=" << daughterPdgCode <<
3749  ", px=" << daughter->Px() <<
3750  ", py=" << daughter->Py() <<
3751  ", pz=" << daughter->Pz() <<
3752  std::endl;
3753  // Codes between 81 and 100 are for internal MC code use, they may be intermediate states used e.g. in hadronizaion models
3754  if (daughterPdgCode == PdgCode) lastInPartonShower = kFALSE; // this parton is not the last parton in the shower
3755  if (TMath::Abs(daughterPdgCode) >= 111 || (daughterPdgCode >= 81 && daughterPdgCode <= 100)) hadronDaughter = kTRUE;
3756  }
3757  if (hadronDaughter) {
3758  AliDebugStream(5) << "This particle has at least a hadron among its daughters!" << std::endl;
3759  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;
3760  }
3761  else {
3762  AliDebugStream(5) << "This particle does not have hadrons among its daughters!" << std::endl;
3763  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;
3764  continue;
3765  }
3766 
3767  if (PdgCode == 4 || PdgCode == -4) {
3768  fHistManager.FillTH1("fHistCharmPt", part.first.Pt());
3769  fHistManager.FillTH1("fHistCharmEta", part.first.Eta());
3770  fHistManager.FillTH1("fHistCharmPhi", part.first.Phi_0_2pi());
3771  }
3772  else if (PdgCode == 5 || PdgCode == -5) {
3773  fHistManager.FillTH1("fHistBottomPt", part.first.Pt());
3774  fHistManager.FillTH1("fHistBottomEta", part.first.Eta());
3775  fHistManager.FillTH1("fHistBottomPhi", part.first.Phi_0_2pi());
3776  }
3777  nHQ++;
3778  }
3779  fHistManager.FillTH1("fHistNHeavyQuarks", nHQ);
3780  fHistManager.FillTH1("fHistHighestPartonPt",highestPartonPt);
3781  Int_t partonType = 0;
3782  Int_t absPdgHighParton = TMath::Abs(PdgHighParton);
3783  if (absPdgHighParton == 9 || absPdgHighParton == 21) {
3784  partonType = 7;
3785  }
3786  else {
3787  partonType = absPdgHighParton;
3788  }
3789  fHistManager.FillTH1("fHistHighestPartonType",partonType);
3790 }
3791 
3800 {
3801  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
3802 
3803  Double_t mass = part->Mass();
3804 
3805  // To make sure that the PDG mass value is not at the edge of a bin
3806  if (nbins % 2 == 0) {
3807  minMass = mass - range / 2 - range / nbins / 2;
3808  maxMass = mass + range / 2 - range / nbins / 2;
3809  }
3810  else {
3811  minMass = mass - range / 2;
3812  maxMass = mass + range / 2;
3813  }
3814 }
3815 
3822 {
3823  static TString label;
3824  label = "";
3825 
3826  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
3827  label = "NotSelTrigger";
3828  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
3829  return label.Data();
3830  }
3831  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
3832  label = "NoVertex";
3833  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
3834  return label.Data();
3835  }
3836  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
3837  label = "TooFewVtxContrib";
3838  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
3839  return label.Data();
3840  }
3841  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
3842  label = "ZVtxOutFid";
3843  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
3844  return label.Data();
3845  }
3846  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
3847  label = "Pileup";
3848  bitmap &= ~BIT(AliRDHFCuts::kPileup);
3849  return label.Data();
3850  }
3851  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
3852  label = "OutsideCentrality";
3853  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
3854  return label.Data();
3855  }
3856  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
3857  label = "PhysicsSelection";
3858  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
3859  return label.Data();
3860  }
3861  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
3862  label = "BadSPDVertex";
3863  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
3864  return label.Data();
3865  }
3866  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
3867  label = "ZVtxSPDOutFid";
3868  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
3869  return label.Data();
3870  }
3871  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
3872  label = "CentralityFlattening";
3873  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
3874  return label.Data();
3875  }
3876  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
3877  label = "BadTrackV0Correl";
3878  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
3879  return label.Data();
3880  }
3881 
3882  return label.Data();
3883 }
3884 
3891 {
3892  if (handler->GetDataSlotNumber() >= 0 && handler->GetOutputObject()) {
3893  PostData(handler->GetDataSlotNumber(), handler->GetOutputObject());
3894  return handler->GetDataSlotNumber();
3895  }
3896  else {
3897  return -1;
3898  }
3899 }
3900 
3910 {
3911  // Get the pointer to the existing analysis manager via the static access method
3912  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
3913  if (!mgr) {
3914  ::Error("AddTaskDmesonJetsSub", "No analysis manager to connect to.");
3915  return NULL;
3916  }
3917 
3918  // Check the analysis type using the event handlers connected to the analysis manager
3919  AliVEventHandler* handler = mgr->GetInputEventHandler();
3920  if (!handler) {
3921  ::Error("AddTaskDmesonJetsSub", "This task requires an input event handler");
3922  return NULL;
3923  }
3924 
3926 
3927  if (handler->InheritsFrom("AliESDInputHandler")) {
3928  dataType = kESD;
3929  }
3930  else if (handler->InheritsFrom("AliAODInputHandler")) {
3931  dataType = kAOD;
3932  }
3933 
3934  // Init the task and do settings
3935  if (ntracks == "usedefault") {
3936  if (dataType == kESD) {
3937  ntracks = "Tracks";
3938  }
3939  else if (dataType == kAOD) {
3940  ntracks = "tracks";
3941  }
3942  else {
3943  ntracks = "";
3944  }
3945  }
3946 
3947  if (nclusters == "usedefault") {
3948  if (dataType == kESD) {
3949  nclusters = "CaloClusters";
3950  }
3951  else if (dataType == kAOD) {
3952  nclusters = "caloClusters";
3953  }
3954  else {
3955  nclusters = "";
3956  }
3957  }
3958 
3959  if (nMCpart == "usedefault") {
3960  nMCpart = "mcparticles"; // Always needs AliAODMCParticle objects
3961  }
3962 
3963  TString name("AliAnalysisTaskDmesonJetsSub");
3964  if (strcmp(suffix, "") != 0) {
3965  name += TString::Format("_%s", suffix.Data());
3966  }
3967 
3968  AliAnalysisTaskDmesonJetsSub* jetTask = new AliAnalysisTaskDmesonJetsSub(name, nMaxTrees);
3969 
3970  if (!ntracks.IsNull()) {
3971  AliHFTrackContainer* trackCont = new AliHFTrackContainer(ntracks);
3972  jetTask->AdoptParticleContainer(trackCont);
3973  }
3974 
3975  if (!nMCpart.IsNull()) {
3976  AliMCParticleContainer* partCont = new AliHFAODMCParticleContainer(nMCpart);
3977  partCont->SetEtaLimits(-1.5, 1.5);
3978  partCont->SetPtLimits(0, 1000);
3979  jetTask->AdoptParticleContainer(partCont);
3980  }
3981 
3982  jetTask->AddClusterContainer(nclusters.Data());
3983 
3984  // Final settings, pass to manager and set the containers
3985  mgr->AddTask(jetTask);
3986 
3987  // Create containers for input/output
3988  AliAnalysisDataContainer* cinput1 = mgr->GetCommonInputContainer();
3989  TString contname1(name);
3990  contname1 += "_histos";
3991  AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(contname1.Data(),
3992  TList::Class(), AliAnalysisManager::kOutputContainer,
3993  Form("%s", AliAnalysisManager::GetCommonFileName()));
3994 
3995  mgr->ConnectInput(jetTask, 0, cinput1);
3996  mgr->ConnectOutput(jetTask, 1, coutput1);
3997 
3998  for (Int_t i = 0; i < nMaxTrees; i++) {
3999  TString contname = TString::Format("%s_tree_%d", name.Data(), i);
4000  AliAnalysisDataContainer *coutput = mgr->CreateContainer(contname.Data(),
4001  TTree::Class(),AliAnalysisManager::kOutputContainer,
4002  Form("%s", AliAnalysisManager::GetCommonFileName()));
4003  mgr->ConnectOutput(jetTask, 2+i, coutput);
4004  }
4005  return jetTask;
4006 }
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
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 FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
Int_t fDataSlotNumber
! Data slot where the tree output is posted
Double_t fMinNeutralPt
Minimum pt of the leading neutral particle (or cluster)
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.
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) ...
const Int_t nPtBins
AliHFAODMCParticleContainer * fMCContainer
! MC particle container
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:34
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)
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)
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.
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:41
AliVParticle * fDmesonParticle
! pointer to the particle object
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.
static EMesonDecayChannel_t CheckDecayChannel(const AliAODMCParticle *part, TClonesArray *mcArray)
std::map< std::string, AliParticleContainer * > fParticleCollArray
particle/track collection array
Double_t fMinChargedPt
Minimum pt of the leading charged particle (or track)
Short_t fPartonType
! type of the parton in the shower tree (only for particle level D mesons)
Bool_t FillHnSparse(THnSparse *h, const AliDmesonJetInfo &DmesonJet, std::string n)
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:292
AliHFAODMCParticleContainer * fMCContainer
! MC particle container
Double32_t fCorrPt
Transverse momentum of the jet in GeV/c after subtracting average background.
unsigned short UShort_t
Definition: External.C:28
Double_t InvMassD0() const
virtual void AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
Int_t PostDataFromAnalysisEngine(OutputHandler const *handler)
Double_t InvMassDstarKpipi() const
Struct that encapsulates analysis parameters.
Bool_t fReconstructed
! Whether this D meson was reconstructed (only for particle level D mesons)
Double_t GetDistance(std::string n, Double_t &deta, Double_t &dphi) const
const Int_t nbins
Double_t maxMass
EOutputType_t fOutputType
Output type: none, TTree or THnSparse.
bool Bool_t
Definition: External.C:53
Byte_t fSelectionType
! for D0: 0=not selected, 1=D0, 2=D0bar
Double_t CosPointingAngle() const
Double_t fMaxNeutralPt
Maximum pt of the leading neutral particle (or cluster)
std::vector< AliAnalysisTaskDmesonJetsSub::AliHFJetDefinition > & GetJetDefinitions()
Double32_t fInvMass
Invariant mass of the D0 meson candidate in GeV/c2.
Int_t GetNAcceptedParticles() const
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:123
std::map< std::string, AliClusterContainer * > fClusterCollArray
cluster collection array
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
Create a new THnSparse within the container.
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:312
std::vector< AliHFJetDefinition > * fJetDefinitions
! Jet definitions
EMCMode_t fMCMode
! MC mode: No MC (data and MC detector level), background-only (MC), signal-only (MC), MC truth (particle level)
void swap(PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &first, PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &second)
EJetAlgo_t fJetAlgo
Jet algorithm (kt, anti-kt,...)
Container for MC-true particles within the EMCAL framework.
Lightweight class that encapsulates D meson jets.
AliRhoParameter * fRho
Object that holds the average background value.
EMCMode_t fMCMode
MC mode: No MC (data and MC detector level), background-only (MC), signal-only (MC), MC truth (particle level)
std::list< AliAnalysisTaskDmesonJetsSub::AnalysisEngine >::iterator FindAnalysisEngine(const AnalysisEngine &eng)
Double_t fMaxNeutralPt
Transverse momentum of the leading neutral particle (or cluster)
Double32_t fPt
Transverse momentum of the D meson in GeV/c.
Double32_t fPartonPt
Transverse momentum of the parton.
Double_t fCent
!event centrality
void SetCharge(EChargeCut_t c)
Container for jet within the EMCAL jet framework.
THistManager fHistManager
Histogram manager.
Definition: External.C:196