AliPhysics  9b6b435 (9b6b435)
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 
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 
2791  double nall = 0;
2792  double zg = 0.;
2793  double flagSubjet=0;
2794  double xconstperp=0;
2795  TString hname;
2796  TString hname2;
2797  fastjet::JetAlgorithm jet_algo(fastjet::cambridge_algorithm);
2798  double jet_radius_ca = 1.0;
2799  fastjet::JetDefinition jet_def(jet_algo, jet_radius_ca,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
2800  hname = TString::Format("%s/%s/LundIterative", GetName(), jetDef.GetName());
2801  THnSparse* h = static_cast<THnSparse*>(fHistManager->FindObject(hname));
2802  hname2 = TString::Format("%s/%s/AngleDifference", GetName(), jetDef.GetName());
2803  TH2* hdiffangle = static_cast<TH2*>(fHistManager->FindObject(hname2));
2804  try{
2805  std::vector<fastjet::PseudoJet> particles(fFastJetWrapper->GetJetConstituents(ijet));
2806  fastjet::ClusterSequence cs_ca(particles, jet_def);
2807  std::vector<fastjet::PseudoJet> output_jets = cs_ca.inclusive_jets(0);
2808  output_jets = sorted_by_pt(output_jets);
2809 
2810 
2811  fastjet::PseudoJet jj = output_jets[0];
2812  fastjet::PseudoJet j1;
2813  fastjet::PseudoJet j2;
2814 
2815  while(jj.has_parents(j1,j2)){
2816  nall = nall + 1;
2817  if(j1.perp() < j2.perp()) std::swap(j1,j2);
2818  flagSubjet=0;
2819  vector < fastjet::PseudoJet > constitj1 = sorted_by_pt(j1.constituents());
2820  if(type==1){
2821  for(Int_t j=0;j<constitj1.size();j++){
2822  if(constitj1[j].user_index()==0){
2823  xconstperp=constitj1[j].perp();
2824  flagSubjet=1; }}}
2825 
2826 
2827  double delta_R = j1.delta_R(j2);
2828  double delta_Raxis=j2.delta_R(output_jets[0]);
2829  zg = j2.perp()/(j1.perp()+j2.perp());
2830  double yh=j1.e()+j2.e();
2831  double y = log(1.0/delta_R);
2832  double lnpt_rel = log(j2.perp()*delta_R);
2833 
2834  double lundEntries[10] = {y, lnpt_rel, output_jets[0].perp(), nall, type, flagSubjet, xconstperp, invmass,yh,TMath::Abs(output_jets[0].eta())};
2835  h->Fill(lundEntries);
2836  hdiffangle->Fill(delta_R, delta_Raxis);
2837  jj=j1;
2838  }
2839 
2840  if(nall==0){ double lundEntrieszero[10]={0,0,output_jets[0].perp(),0,type,0,0,invmass,0,TMath::Abs(output_jets[0].eta())};
2841  h->Fill(lundEntrieszero);}
2842 
2843  } catch (fastjet::Error) { /*return -1;*/ }
2844 
2845 }
2846 
2849 {
2850  TString hname;
2851 
2856 
2857  if (!fMCContainer->IsSpecialPDGFound()) return;
2858 
2859  std::array<int,2> nAccCharm = {0};
2860  std::array<std::array<int, 2>, 5> nAccCharmPt = {{{0}}};
2861 
2862  std::map<AliHFJetDefinition*, Double_t> maxJetPt;
2863  Double_t maxDPt = 0;
2864 
2865  for (auto &jetDef : fJetDefinitions) {
2866  maxJetPt[&jetDef] = 0;
2867  Double_t rho = 0;
2868  if (jetDef.fRho) rho = jetDef.fRho->GetVal();
2869  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", GetName(), jetDef.GetName());
2870  TH1* histNDmesonsVsNconstituents = static_cast<TH1*>(fHistManager->FindObject(hname));
2871 
2872  switch (jetDef.fJetType) {
2874  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNoChargeCut);
2875  break;
2877  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
2878  break;
2880  fMCContainer->SetCharge(AliParticleContainer::EChargeCut_t::kNeutral);
2881  break;
2882  }
2883 
2885  fFastJetWrapper->SetR(jetDef.fRadius);
2888 
2889  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", GetName(), jetDef.GetName());
2890  AddInputVectors(fMCContainer, 100, static_cast<TH2*>(fHistManager->FindObject(hname)));
2891 
2892  fFastJetWrapper->Run();
2893 
2894  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper->GetInclusiveJets();
2895 
2896  for (auto jet : jets_incl) {
2897  Int_t nDmesonsInJet = 0;
2898 
2899  for (auto constituent : jet.constituents()) {
2900  Int_t iPart = constituent.user_index() - 100;
2901  if (constituent.perp() < 1e-6) continue; // reject ghost particles
2902  AliAODMCParticle* part = fMCContainer->GetMCParticle(iPart);
2903  if (!part) {
2904  ::Error("AliAnalysisTaskDmesonJetsSub::AnalysisEngine::RunParticleLevelAnalysis", "Could not find jet constituent %d!", iPart);
2905  continue;
2906  }
2907  if (TMath::Abs(part->PdgCode()) == fCandidatePDG) {
2908  nDmesonsInJet++;
2909  std::map<int, AliDmesonJetInfo>::iterator dMesonJetIt = fDmesonJets.find(iPart);
2910  if (dMesonJetIt == fDmesonJets.end()) { // This D meson does not exist yet
2911  if (part->Pt() > maxDPt) maxDPt = part->Pt();
2912  std::pair<int, AliDmesonJetInfo> element;
2913  element.first = iPart;
2914  dMesonJetIt = fDmesonJets.insert(element).first;
2915  (*dMesonJetIt).second.fD.SetPxPyPzE(part->Px(), part->Py(), part->Pz(), part->E());
2916  (*dMesonJetIt).second.fDmesonParticle = part;
2917  (*dMesonJetIt).second.fSelectionType = part->PdgCode() > 0 ? 1 : 2;
2918 
2919  UShort_t p = 0;
2920  UInt_t rs = 0;
2921 
2922  auto origin = IsPromptCharm(part, fMCContainer->GetArray());
2923  p = 0;
2924  rs = origin.first;
2925  while (rs >>= 1) { p++; }
2926  (*dMesonJetIt).second.fPartonType = p;
2927  (*dMesonJetIt).second.fParton = origin.second;
2928 
2929  (*dMesonJetIt).second.fAncestor = FindParticleOrigin(part, fMCContainer->GetArray(), kFindFirst);
2930 
2931  Int_t im = -1;
2932  if (part->PdgCode() > 0) { // D0
2933  im = 0;
2934  }
2935  else { // D0bar
2936  im = 1;
2937  }
2938 
2939  nAccCharm[im]++;
2940  for (int i = 0; i < nAccCharmPt.size(); i++) {
2941  if (part->Pt() < i) break;
2942  nAccCharmPt[i][im]++;
2943  }
2944  }
2945 
2946  (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.SetPxPyPzE(jet.px(), jet.py(), jet.pz(), jet.E());
2947  (*dMesonJetIt).second.fJets[jetDef.GetName()].fNConstituents = jet.constituents().size();
2948  (*dMesonJetIt).second.fJets[jetDef.GetName()].fArea = jet.area();
2949  (*dMesonJetIt).second.fJets[jetDef.GetName()].fCorrPt = (*dMesonJetIt).second.fJets[jetDef.GetName()].fMomentum.Pt() - jet.area() * rho;
2950  if (jet.perp() > maxJetPt[&jetDef]) maxJetPt[&jetDef] = jet.perp();
2951  } // if constituent is a D meson
2952  } // for each constituent
2953  if (nDmesonsInJet > 0) histNDmesonsVsNconstituents->Fill(jet.constituents().size(), nDmesonsInJet);
2954  } // for each jet
2955  } // for each jet definition
2956 
2958 
2959  for (auto& def : fJetDefinitions) {
2960  if (!def.fRho) continue;
2961  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", GetName(), def.GetName());
2962  fHistManager->FillTH2(hname, maxJetPt[&def], def.fRho->GetVal());
2963 
2964  hname = TString::Format("%s/%s/fHistRhoVsLeadDPt", GetName(), def.GetName());
2965  fHistManager->FillTH2(hname, maxDPt, def.fRho->GetVal());
2966 
2967  hname = TString::Format("%s/%s/fHistRhoVsCent", GetName(), def.GetName());
2968  fHistManager->FillTH2(hname, fEventInfo.fCent, def.fRho->GetVal());
2969 
2970  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", GetName(), def.GetName());
2971  fHistManager->FillTH2(hname, fEventInfo.fCent, maxJetPt[&def]);
2972 
2973  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", GetName(), def.GetName());
2974  fHistManager->FillTH2(hname, fEventInfo.fCent, maxDPt);
2975 
2976  hname = TString::Format("%s/%s/fHistRhoVsNTracks", GetName(), def.GetName());
2977  fHistManager->FillTH2(hname, npart, def.fRho->GetVal());
2978 
2979  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", GetName(), def.GetName());
2980  fHistManager->FillTH2(hname, npart, maxJetPt[&def]);
2981 
2982  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", GetName(), def.GetName());
2983  fHistManager->FillTH2(hname, npart, maxDPt);
2984  }
2985 
2986  if (fDmesonJets.size() != nAccCharm[0]+nAccCharm[1]) AliError(Form("I found %lu mesons (%d)?", fDmesonJets.size(), nAccCharm[0]+nAccCharm[1]));
2987  hname = TString::Format("%s/fHistNTotAcceptedDmesons", GetName());
2988  fHistManager->FillTH1(hname, "D", nAccCharm[0]);
2989  fHistManager->FillTH1(hname, "Anti-D", nAccCharm[1]);
2990 
2991  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", GetName());
2992  fHistManager->FillTH2(hname, npart, nAccCharm[0]+nAccCharm[1]);
2993 
2994  for (int i = 0; i < nAccCharmPt.size(); i++) {
2995  hname = TString::Format("%s/fHistNTotAcceptedDmesonsPt%d", GetName(), i);
2996  fHistManager->FillTH1(hname, "D", nAccCharmPt[i][0]);
2997  fHistManager->FillTH1(hname, "Anti-D", nAccCharmPt[i][1]);
2998 
2999  hname = TString::Format("%s/fHistNAcceptedDmesonsPt%d", GetName(), i);
3000  fHistManager->FillTH1(hname, nAccCharmPt[i][0]+nAccCharmPt[i][1]);
3001  }
3002 
3003  hname = TString::Format("%s/fHistNDmesons", GetName());
3004  fHistManager->FillTH1(hname, nAccCharm[0]+nAccCharm[1]); // same as the number of accepted D mesons, since no selection is performed
3005 }
3006 
3007 
3008 
3009 // Definitions of class AliAnalysisTaskDmesonJetsSub
3010 
3014 
3018  fAnalysisEngines(),
3019  fEnabledAxis(0),
3021  fHistManager(),
3022  fApplyKinematicCuts(kTRUE),
3023  fNOutputTrees(0),
3024  fTrackEfficiency(0),
3025  fRejectISR(kFALSE),
3026  fJetAreaType(fastjet::active_area),
3027  fJetGhostArea(0.005),
3028  fMCContainer(0),
3029  fAodEvent(0),
3030  fFastJetWrapper(0)
3031 {
3032  SetMakeGeneralHistograms(kTRUE);
3033 }
3034 
3039  AliAnalysisTaskEmcalLight(name, kTRUE),
3040  fAnalysisEngines(),
3043  fHistManager(name),
3044  fApplyKinematicCuts(kTRUE),
3045  fNOutputTrees(nOutputTrees),
3046  fTrackEfficiency(0),
3047  fRejectISR(kFALSE),
3048  fJetAreaType(fastjet::active_area),
3049  fJetGhostArea(0.005),
3050  fMCContainer(0),
3051  fAodEvent(0),
3052  fFastJetWrapper(0)
3053 {
3054  SetMakeGeneralHistograms(kTRUE);
3055  for (Int_t i = 0; i < nOutputTrees; i++){
3056  DefineOutput(2+i, TTree::Class());
3057  }
3058 }
3059 
3062 {
3063  if (fFastJetWrapper) delete fFastJetWrapper;
3064 }
3065 
3073 {
3074  AliRDHFCuts* analysiscuts = 0;
3075  TFile* filecuts = TFile::Open(cutfname);
3076  if (!filecuts || filecuts->IsZombie()) {
3077  ::Error("AliAnalysisTaskDmesonJetsSub::LoadDMesonCutsFromFile", "Input file not found: will use std cuts.");
3078  filecuts = 0;
3079  }
3080 
3081  if (filecuts) analysiscuts = dynamic_cast<AliRDHFCuts*>(filecuts->Get(cutsname));
3082 
3083  if (!analysiscuts) {
3084  ::Error("AliAnalysisTaskDmesonJetsSub::LoadDMesonCutsFromFile", "Could not find analysis cuts '%s' in '%s'.", cutsname.Data(), cutfname.Data());
3085  if (filecuts) {
3086  filecuts->ls();
3087  }
3088  }
3089  else {
3090  ::Info("AliAnalysisTaskDmesonJetsSub::LoadDMesonCutsFromFile", "Cuts '%s' loaded from file '%s'", cutsname.Data(), cutfname.Data());
3091  }
3092 
3093  return analysiscuts;
3094 }
3095 
3108 {
3110  return AddAnalysisEngine(type, cutfname, cuttype, MCmode, jetDef, rhoName);
3111 }
3112 
3124 {
3125  AliRDHFCuts* cuts = 0;
3126 
3127  if (!cutfname.IsNull()) {
3128  TString cutsname;
3129 
3130  switch (type) {
3131  case kD0toKpi:
3132  case kD0toKpiLikeSign:
3133  cutsname = "D0toKpiCuts";
3134  break;
3135  case kDstartoKpipi:
3136  cutsname = "DStartoKpipiCuts";
3137  break;
3138  default:
3139  return 0;
3140  }
3141 
3142  if (!cuttype.IsNull()) {
3143  cutsname += TString::Format("_%s", cuttype.Data());
3144  }
3145 
3146  cuts = LoadDMesonCutsFromFile(cutfname, cutsname);
3147  if (cuts) cuts->PrintAll();
3148  }
3149 
3150  AnalysisEngine eng(type, MCmode, cuts);
3151 
3152  std::list<AnalysisEngine>::iterator it = FindAnalysisEngine(eng);
3153 
3154  if (it == fAnalysisEngines.end() || *it != eng) { // No analysis engine was found, adding a new one
3155  eng.AddJetDefinition(jetDef);
3156  it = fAnalysisEngines.insert(it, eng);
3157  ::Info("AliAnalysisTaskDmesonJetsSub::AddAnalysisEngine", "A new analysis engine '%s' has been added. The total number of analysis engines is %lu.", eng.GetName(), fAnalysisEngines.size());
3158  }
3159  else {
3160  AnalysisEngine* found_eng = &(*it);
3161  ::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());
3162  found_eng->AddJetDefinition(jetDef);
3163 
3164  if (cuts) {
3165  if (found_eng->fRDHFCuts != 0) ::Warning("AliAnalysisTaskDmesonJetsSub::AddAnalysisEngine", "D meson cuts were already defined for this D meson type. They will be overwritten.");
3166  found_eng->SetRDHFCuts(cuts);
3167  }
3168  }
3169 
3170  return &(*it);
3171 }
3172 
3173 std::list<AliAnalysisTaskDmesonJetsSub::AnalysisEngine>::iterator AliAnalysisTaskDmesonJetsSub::FindAnalysisEngine(const AliAnalysisTaskDmesonJetsSub::AnalysisEngine& eng)
3174 {
3175  std::list<AnalysisEngine>::iterator it = fAnalysisEngines.begin();
3176  while (it != fAnalysisEngines.end() && (*it) != eng) it++;
3177  return it;
3178 }
3179 
3182 {
3183  ::Info("UserCreateOutputObjects", "CreateOutputObjects of task %s", GetName());
3184 
3186 
3187  // Define histograms
3188  // the TList fOutput is already defined in AliAnalysisTaskEmcalLight::UserCreateOutputObjects()
3189 
3190  TString hname;
3191  TString htitle;
3192  TH1* h = 0;
3193  Int_t treeSlot = 0;
3194 
3195  Int_t maxTracks = 6000;
3196  Double_t maxRho = 500;
3197  if (fForceBeamType == kpp) {
3198  maxRho = 50;
3199  maxTracks = 200;
3200  }
3201  else if (fForceBeamType == kpA) {
3202  maxRho = 200;
3203  maxTracks = 500;
3204  }
3205 
3206  Int_t dimx = 10;
3207  Int_t nbinsx[10] = {50,100,10,20,2,2,200,150,100,9};
3208  Double_t minx[10] = {0,-10,0,0,0,0,0,1.6,0,0};
3209  Double_t maxx[10] = {5,10,100,20,2,2,100,2.3,100,0.9};
3210  TString titlex[10]={"log(1/deltaR)","log(zteta)","jet pt","n","type","flagSubjet","ptD","invmass","frac","abs(eta)"};
3211 
3212 
3213 
3214 
3215  hname = "fHistCharmPt";
3216  htitle = hname + ";#it{p}_{T,charm} (GeV/#it{c});counts";
3217  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3218 
3219  hname = "fHistCharmEta";
3220  htitle = hname + ";#eta_{charm};counts";
3221  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
3222 
3223  hname = "fHistCharmPhi";
3224  htitle = hname + ";#phi_{charm};counts";
3225  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
3226 
3227  hname = "fHistBottomPt";
3228  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
3229  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3230 
3231  hname = "fHistBottomEta";
3232  htitle = hname + ";#eta_{bottom};counts";
3233  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
3234 
3235  hname = "fHistBottomPhi";
3236  htitle = hname + ";#phi_{bottom};counts";
3237  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
3238 
3239  hname = "fHistHighestPartonPt";
3240  htitle = hname + ";#it{p}_{T,bottom} (GeV/#it{c});counts";
3241  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3242 
3243  hname = "fHistHighestPartonType";
3244  htitle = hname + ";type;counts";
3245  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
3246 
3247  hname = "fHistNHeavyQuarks";
3248  htitle = hname + ";number of heavy-quarks;counts";
3249  fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
3250 
3251  ::Info("AliAnalysisTaskDmesonJetsSub::UserCreateOutputObjects", "Allocating histograms for task '%s' (%lu analysis engines)", GetName(), fAnalysisEngines.size());
3252  for (auto &param : fAnalysisEngines) {
3253  ::Info("AliAnalysisTaskDmesonJetsSub::UserCreateOutputObjects", "Allocating histograms for analysis engine '%s' (%lu jet definitions)", param.GetName(), param.fJetDefinitions.size());
3254 
3255  param.fHistManager = &fHistManager;
3256 
3257  hname = TString::Format("%s/fHistNAcceptedDmesonsVsNtracks", param.GetName());
3258  htitle = hname + ";#it{N}_{tracks};#it{N}_{D};events";
3259  h = fHistManager.CreateTH2(hname, htitle, 251, -0.5, 250.5, 21, -0.5, 20.5);
3260 
3261  hname = TString::Format("%s/fHistNTotAcceptedDmesons", param.GetName());
3262  htitle = hname + ";;#it{N}_{D}";
3263  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
3264 
3265  for (int i = 0 ; i < 5; i++) {
3266  hname = TString::Format("%s/fHistNAcceptedDmesonsPt%d", param.GetName(), i);
3267  htitle = hname + ";#it{N}_{D};events";
3268  h = fHistManager.CreateTH1(hname, htitle, 21, -0.5, 20.5);
3269 
3270  hname = TString::Format("%s/fHistNTotAcceptedDmesonsPt%d", param.GetName(), i);
3271  htitle = hname + ";;#it{N}_{D}";
3272  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
3273  }
3274 
3275  hname = TString::Format("%s/fHistNDmesons", param.GetName());
3276  htitle = hname + ";#it{N}_{D};events";
3277  h = fHistManager.CreateTH1(hname, htitle, 501, -0.5, 500.5);
3278 
3279  hname = TString::Format("%s/fHistNEvents", param.GetName());
3280  htitle = hname + ";Event status;counts";
3281  h = fHistManager.CreateTH1(hname, htitle, 2, 0, 2);
3282  h->GetXaxis()->SetBinLabel(1, "Accepted");
3283  h->GetXaxis()->SetBinLabel(2, "Rejected");
3284 
3285  hname = TString::Format("%s/fHistEventRejectionReasons", param.GetName());
3286  htitle = hname + ";Rejection reason;counts";
3287  h = fHistManager.CreateTH1(hname, htitle, 32, 0, 32);
3288 
3289  hname = TString::Format("%s/fHistRejectedDMesonPt", param.GetName());
3290  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});counts";
3291  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
3292 
3293  hname = TString::Format("%s/fHistRejectedDMesonEta", param.GetName());
3294  htitle = hname + ";#it{#eta}_{D};counts";
3295  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
3296 
3297  hname = TString::Format("%s/fHistRejectedDMesonPhi", param.GetName());
3298  htitle = hname + ";#it{#phi}_{D};counts";
3299  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
3300 
3301  if (param.fMCMode != kMCTruth) {
3302  if (param.fCandidateType == kD0toKpi || param.fCandidateType == kD0toKpiLikeSign) {
3303  hname = TString::Format("%s/fHistRejectedDMesonInvMass", param.GetName());
3304  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
3305  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, param.fMinMass, param.fMaxMass);
3306  }
3307  else if (param.fCandidateType == kDstartoKpipi) {
3308  Double_t min = 0;
3309  Double_t max = 0;
3310 
3311  hname = TString::Format("%s/fHistRejectedDMeson2ProngInvMass", param.GetName());
3312  htitle = hname + ";#it{M}_{K#pi} (GeV/#it{c}^{2});counts";
3313  CalculateMassLimits(param.fMaxMass - param.fMinMass, 421, param.fNMassBins, min, max);
3314  fHistManager.CreateTH1(hname, htitle, param.fNMassBins, min, max);
3315 
3316  Double_t D0mass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
3317  hname = TString::Format("%s/fHistRejectedDMesonDeltaInvMass", param.GetName());
3318  htitle = hname + ";#it{M}_{K#pi#pi} - #it{M}_{K#pi} (GeV/#it{c}^{2});counts";
3319  CalculateMassLimits(0.20, 413, param.fNMassBins*6, min, max);
3320  fHistManager.CreateTH1(hname, htitle, param.fNMassBins*6, min-D0mass, max-D0mass);
3321  }
3322  }
3323 
3324  if (param.fMCMode == kMCTruth) {
3325  hname = TString::Format("%s/fHistPartonPt", param.GetName());
3326  htitle = hname + ";#it{p}_{T,parton} (GeV/#it{c});counts";
3327  fHistManager.CreateTH1(hname, htitle, 500, 0, 1000);
3328 
3329  hname = TString::Format("%s/fHistPartonEta", param.GetName());
3330  htitle = hname + ";#eta_{parton};counts";
3331  fHistManager.CreateTH1(hname, htitle, 400, -10, 10);
3332 
3333  hname = TString::Format("%s/fHistPartonPhi", param.GetName());
3334  htitle = hname + ";#phi_{parton};counts";
3335  fHistManager.CreateTH1(hname, htitle, 125, 0, TMath::TwoPi());
3336 
3337  hname = TString::Format("%s/fHistPartonType", param.GetName());
3338  htitle = hname + ";type;counts";
3339  fHistManager.CreateTH1(hname, htitle, 10, 0, 10);
3340 
3341  hname = TString::Format("%s/fHistPrompt", param.GetName());
3342  htitle = hname + ";Type;counts";
3343  h = fHistManager.CreateTH1(hname, htitle, 3, 0, 3);
3344  h->GetXaxis()->SetBinLabel(1, "Unknown");
3345  h->GetXaxis()->SetBinLabel(2, "Prompt");
3346  h->GetXaxis()->SetBinLabel(3, "Non-Prompt");
3347 
3348  hname = TString::Format("%s/fHistAncestor", param.GetName());
3349  htitle = hname + ";Ancestor;counts";
3350  h = fHistManager.CreateTH1(hname, htitle, 4, 0, 4);
3351  h->GetXaxis()->SetBinLabel(1, "Unknown");
3352  h->GetXaxis()->SetBinLabel(2, "Charm");
3353  h->GetXaxis()->SetBinLabel(3, "Bottom");
3354  h->GetXaxis()->SetBinLabel(4, "Proton");
3355  }
3356 
3357  for (auto& jetDef : param.fJetDefinitions) {
3358  ::Info("AliAnalysisTaskDmesonJetsSub::UserCreateOutputObjects", "Allocating histograms for jet definition '%s'", jetDef.GetName());
3359 
3360  if (param.fMCMode == kMCTruth) {
3361  hname = TString::Format("%s/%s/fHistNDmesonsVsNconstituents", param.GetName(), jetDef.GetName());
3362  htitle = hname + ";#it{N}_{constituents};#it{N}_{D};counts";
3363  h = fHistManager.CreateTH2(hname, htitle, 51, -0.5, 50.5, 10, 0.5, 10.5);
3364  }
3365 
3366  hname = TString::Format("%s/%s/fHistMCParticleRejectionReason", param.GetName(), jetDef.GetName());
3367  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
3368  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
3369  SetRejectionReasonLabels(h->GetXaxis());
3370 
3371  hname = TString::Format("%s/%s/fHistTrackRejectionReason", param.GetName(), jetDef.GetName());
3372  htitle = hname + ";Track rejection reason;#it{p}_{T,track} (GeV/#it{c});counts";
3373  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
3374  SetRejectionReasonLabels(h->GetXaxis());
3375 
3376  hname = TString::Format("%s/%s/fHistClusterRejectionReason", param.GetName(), jetDef.GetName());
3377  htitle = hname + ";Cluster rejection reason;#it{p}_{T,cluster} (GeV/#it{c});counts";
3378  h = fHistManager.CreateTH2(hname, htitle, 32, 0, 32, 150, 0, 150);
3379  SetRejectionReasonLabels(h->GetXaxis());
3380 
3381  hname = TString::Format("%s/%s/fHistDMesonDaughterNotInJet", param.GetName(), jetDef.GetName());
3382  htitle = hname + ";#it{p}_{T,track} (GeV/#it{c});counts";
3383  fHistManager.CreateTH1(hname, htitle, 200, 0, 100);
3384 
3385  hname = TString::Format("%s/%s/fHistRejectedJetPt", param.GetName(), jetDef.GetName());
3386  htitle = hname + ";#it{p}_{T,jet} (GeV/#it{c});counts";
3387  fHistManager.CreateTH1(hname, htitle, 150, 0, 150);
3388 
3389  hname = TString::Format("%s/%s/fHistRejectedJetEta", param.GetName(), jetDef.GetName());
3390  htitle = hname + ";#it{#eta}_{jet};counts";
3391  fHistManager.CreateTH1(hname, htitle, 100, -2, 2);
3392 
3393  hname = TString::Format("%s/%s/fHistRejectedJetPhi", param.GetName(), jetDef.GetName());
3394  htitle = hname + ";#it{#phi}_{jet};counts";
3395  fHistManager.CreateTH1(hname, htitle, 200, 0, TMath::TwoPi());
3396 
3397  if (!jetDef.fRhoName.IsNull()) {
3398  hname = TString::Format("%s/%s/fHistRhoVsLeadJetPt", param.GetName(), jetDef.GetName());
3399  htitle = hname + ";#it{p}_{T,jet} (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/fHistRhoVsLeadDPt", param.GetName(), jetDef.GetName());
3403  htitle = hname + ";#it{p}_{T,D} (GeV/#it{c});#rho (GeV/#it{c} #times rad^{-1});counts";
3404  fHistManager.CreateTH2(hname, htitle, 300, 0, 150, 1000, 0, maxRho);
3405 
3406  hname = TString::Format("%s/%s/fHistRhoVsCent", param.GetName(), jetDef.GetName());
3407  htitle = hname + ";Centrality (%);#rho (GeV/#it{c} #times rad^{-1});counts";
3408  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 1000, 0, maxRho);
3409 
3410  hname = TString::Format("%s/%s/fHistLeadJetPtVsCent", param.GetName(), jetDef.GetName());
3411  htitle = hname + ";Centrality (%);#it{p}_{T,jet} (GeV/#it{c});counts";
3412  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
3413 
3414  hname = TString::Format("%s/%s/fHistLeadDPtVsCent", param.GetName(), jetDef.GetName());
3415  htitle = hname + ";Centrality (%);#it{p}_{T,D} (GeV/#it{c});counts";
3416  fHistManager.CreateTH2(hname, htitle, 100, 0, 100, 300, 0, 150);
3417 
3418  hname = TString::Format("%s/%s/fHistRhoVsNTracks", param.GetName(), jetDef.GetName());
3419  htitle = hname + ";no. of tracks;#rho (GeV/#it{c} #times rad^{-1});counts";
3420  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 1000, 0, maxRho);
3421 
3422  hname = TString::Format("%s/%s/fHistLeadJetPtVsNTracks", param.GetName(), jetDef.GetName());
3423  htitle = hname + ";no. of tracks;#it{p}_{T,jet} (GeV/#it{c});counts";
3424  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
3425 
3426  hname = TString::Format("%s/%s/fHistLeadDPtVsNTracks", param.GetName(), jetDef.GetName());
3427  htitle = hname + ";no. of tracks;#it{p}_{T,D} (GeV/#it{c});counts";
3428  fHistManager.CreateTH2(hname, htitle, 200, 0, maxTracks, 300, 0, 150);
3429  }
3430  hname = TString::Format("%s/%s/LundIterative",param.GetName(),jetDef.GetName());
3431  cout<<"at the begining"<<hname<<endl;
3432  THnSparse* h = fHistManager.CreateTHnSparse(hname,hname,dimx,nbinsx,minx,maxx);
3433  for (Int_t j = 0; j < dimx; j++) {
3434  h->GetAxis(j)->SetTitle(titlex[j]);}
3435 
3436 
3437  hname = TString::Format("%s/%s/AngleDifference",param.GetName(),jetDef.GetName());
3438  htitle = hname + ";angle iterative declustering;angle to axis";
3439  fHistManager.CreateTH2(hname,htitle,100,0.,2*3.1416,100,0,2*3.1416);
3440 
3441  }
3442  switch (fOutputType) {
3443  case kTreeOutput:
3444  {
3445  OutputHandlerTTree* tree_handler = new OutputHandlerTTree(&param);
3446  param.fOutputHandler = tree_handler;
3447  tree_handler->BuildOutputObject(GetName());
3448  if (treeSlot < fNOutputTrees) {
3449  tree_handler->AssignDataSlot(treeSlot+2);
3450  treeSlot++;
3451  PostDataFromAnalysisEngine(tree_handler);
3452  }
3453  else {
3454  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
3455  }
3456  }
3457  break;
3458  case kTHnOutput:
3459  {
3460  OutputHandlerTHnSparse* thnsparse_handler = new OutputHandlerTHnSparse(&param);
3461  param.fOutputHandler = thnsparse_handler;
3462  thnsparse_handler->SetEnabledAxis(fEnabledAxis);
3463  thnsparse_handler->BuildOutputObject(GetName());
3464  }
3465  break;
3466  case kOnlyQAOutput:
3467  {
3468  OutputHandler* qa_handler = new OutputHandler(&param);
3469  param.fOutputHandler = qa_handler;
3470  }
3471  break;
3472  case kNoOutput:
3473  break;
3474  case kTreeExtendedOutput:
3475  {
3477  param.fOutputHandler = tree_handler;
3478  tree_handler->BuildOutputObject(GetName());
3479  if (treeSlot < fNOutputTrees) {
3480  tree_handler->AssignDataSlot(treeSlot+2);
3481  treeSlot++;
3482  PostDataFromAnalysisEngine(tree_handler);
3483  }
3484  else {
3485  AliError(Form("Number of data output slots %d not sufficient. Tree of analysis engine %s will not be posted!", fNOutputTrees, param.GetName()));
3486  }
3487  }
3488  break;
3489  }
3490  }
3491 
3493 
3494  PostData(1, fOutput);
3495 }
3496 
3500 {
3502 
3503  // Load the event
3504  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
3505 
3506  fFastJetWrapper = new AliFJWrapper(fName, fTitle);
3507 
3508  fFastJetWrapper->SetAreaType((fastjet::AreaType)fJetAreaType);
3510 
3511  if (!fAodEvent) {
3512  AliError(Form("This task need an AOD event (Task '%s'). Expect troubles...", GetName()));
3513  //return;
3514  }
3515 
3516  TRandom* rnd = 0;
3517  if (fTrackEfficiency > 0 && fTrackEfficiency < 1) rnd = new TRandom3(0);
3518 
3519  for (auto cont_it : fParticleCollArray) {
3520  AliHFAODMCParticleContainer* part_cont = dynamic_cast<AliHFAODMCParticleContainer*>(cont_it.second);
3521  if (part_cont) fMCContainer = part_cont;
3522  }
3523 
3524  for (auto &params : fAnalysisEngines) {
3525 
3526  params.fAodEvent = fAodEvent;
3527  params.fFastJetWrapper = fFastJetWrapper;
3528  params.fTrackEfficiency = fTrackEfficiency;
3529  params.fRejectISR = fRejectISR;
3530  params.fRandomGen = rnd;
3531 
3532  for (auto &jetdef: params.fJetDefinitions) {
3533  if (!jetdef.fRhoName.IsNull()) {
3534  jetdef.fRho = dynamic_cast<AliRhoParameter*>(fInputEvent->FindListObject(jetdef.fRhoName));
3535  if (!jetdef.fRho) {
3536  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3537  "%s: Could not find rho object '%s' for engine '%s'",
3538  GetName(), jetdef.fRhoName.Data(), params.GetName());
3539  }
3540  }
3541  }
3542 
3543  if (!params.fRDHFCuts) {
3544  if (params.fMCMode == kMCTruth) {
3545  ::Warning("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3546  "%s: RDHF cuts not provided for engine '%s'.",
3547  GetName(), params.GetName());
3548  }
3549  else {
3550  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3551  "%s: RDHF cuts not provided. Engine '%s' disabled.",
3552  GetName(), params.GetName());
3553  params.fInhibit = kTRUE;
3554  }
3555  }
3556 
3557  params.fMCContainer = fMCContainer;
3558 
3559  for (auto cont_it : fParticleCollArray) {
3560  AliTrackContainer* track_cont = dynamic_cast<AliTrackContainer*>(cont_it.second);
3561  if (track_cont) params.fTrackContainers.push_back(track_cont);
3562  }
3563 
3564  for (auto cont_it : fClusterCollArray) params.fClusterContainers.push_back(cont_it.second);
3565 
3566  if (fAodEvent) params.Init(fGeom, fAodEvent->GetRunNumber());
3567 
3568  if (params.fMCMode != kMCTruth && fAodEvent) {
3569  params.fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject(params.fBranchName.Data()));
3570 
3571  if (params.fCandidateArray) {
3572  TString className;
3573  if (params.fCandidateType == kD0toKpi || params.fCandidateType == kD0toKpiLikeSign) {
3574  className = "AliAODRecoDecayHF2Prong";
3575  }
3576  else if (params.fCandidateType == kDstartoKpipi) {
3577  className = "AliAODRecoCascadeHF";
3578  }
3579  if (!params.fCandidateArray->GetClass()->InheritsFrom(className)) {
3580  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3581  "%s: Objects of type %s in %s are not inherited from %s! Task will be disabled!",
3582  GetName(), params.fCandidateArray->GetClass()->GetName(), params.fCandidateArray->GetName(), className.Data()); // @suppress("Ambiguous problem")
3583  params.fCandidateArray = 0;
3584  params.fInhibit = kTRUE;
3585  }
3586  }
3587  else {
3588  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3589  "Could not find candidate array '%s', skipping the event. Analysis engine '%s' will be disabled!",
3590  params.fBranchName.Data(), params.GetName());
3591  params.fInhibit = kTRUE;
3592  }
3593  }
3594 
3595  if (params.fMCMode != kNoMC) {
3596  if (!params.fMCContainer) {
3597  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3598  "No MC particle container was provided. Analysis engine '%s' will be disabled!",
3599  params.GetName());
3600  params.fInhibit = kTRUE;
3601  }
3602  }
3603 
3604  if (params.fMCMode != kMCTruth) {
3605  if (params.fTrackContainers.size() == 0 && params.fClusterContainers.size() == 0) {
3606  ::Error("AliAnalysisTaskDmesonJetsSub::ExecOnce",
3607  "No track container and no cluster container were provided. Analysis engine '%s' will be disabled!",
3608  params.GetName());
3609  params.fInhibit = kTRUE;
3610  }
3611  }
3612  }
3613 }
3614 
3619 {
3620  TString hname;
3621 
3622  // Fix for temporary bug in ESDfilter
3623  // The AODs with null vertex pointer didn't pass the PhysSel
3624  // Now adding an entry in the histogram so as to check that this is actually cutting anything out
3625  if (fAodEvent && (!fAodEvent->GetPrimaryVertex() || TMath::Abs(fAodEvent->GetMagneticField()) < 0.001)) {
3626  for (auto &eng : fAnalysisEngines) {
3627  if (eng.fInhibit) continue;
3628  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
3629  fHistManager.FillTH1(hname, "ESDfilterBug");
3630  }
3631  return kFALSE;
3632  }
3633 
3634  if (fAodEvent) {
3635  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
3636  if (matchingAODdeltaAODlevel <= 0) {
3637  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
3638  for (auto &eng : fAnalysisEngines) {
3639  if (eng.fInhibit) continue;
3640  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
3641  fHistManager.FillTH1(hname, "MismatchDeltaAOD");
3642  }
3643  return kFALSE;
3644  }
3645  }
3646 
3647  for (auto &eng : fAnalysisEngines) {
3648  eng.fDmesonJets.clear();
3649  if (eng.fInhibit) continue;
3650 
3651  eng.fEventInfo = EventInfo(fCent, fEPV0, fEventWeight, fPtHard);
3652 
3653  //Event selection
3654  hname = TString::Format("%s/fHistNEvents", eng.GetName());
3655  if (fAodEvent) {
3656  Bool_t iseventselected = kTRUE;
3657  if (eng.fRDHFCuts) iseventselected = eng.fRDHFCuts->IsEventSelected(fAodEvent);
3658  if (!iseventselected) {
3659  fHistManager.FillTH1(hname, "Rejected");
3660  hname = TString::Format("%s/fHistEventRejectionReasons", eng.GetName());
3661  UInt_t bitmap = eng.fRDHFCuts->GetEventRejectionBitMap();
3662  TString label;
3663  do {
3664  label = GetHFEventRejectionReasonLabel(bitmap);
3665  if (label.IsNull()) break;
3666  fHistManager.FillTH1(hname, label);
3667  } while (true);
3668 
3669  continue;
3670  }
3671  }
3672 
3673  fHistManager.FillTH1(hname, "Accepted");
3674 
3675  AliDebug(2, "Event selected");
3676 
3677  eng.RunAnalysis();
3678  }
3679  return kTRUE;
3680 }
3681 
3686 {
3687  for (auto &param : fAnalysisEngines) {
3688  if (param.fInhibit) continue;
3689  param.fOutputHandler->FillOutput(fApplyKinematicCuts);
3690  PostDataFromAnalysisEngine(param.fOutputHandler);
3691  }
3693  return kTRUE;
3694 }
3695 
3698 {
3699  auto itcont = fMCContainer->all_momentum();
3700  Int_t nHQ = 0;
3701  Double_t highestPartonPt = 0;
3702  Int_t PdgHighParton = 0;
3703  for (auto it = itcont.begin(); it != itcont.end(); it++) {
3704  auto part = *it;
3705  if (part.first.Pt() == 0) continue;
3706 
3707  Int_t PdgCode = part.second->GetPdgCode();
3708 
3709  // Skip all particles that are not either quarks or gluons
3710  if ((PdgCode < -9 || PdgCode > 9) && PdgCode != 21 && PdgCode != -21) continue;
3711 
3712  AliDebugStream(5) << "Parton " << it.current_index() <<
3713  " with pdg=" << PdgCode <<
3714  ", px=" << part.first.Px() <<
3715  ", py=" << part.first.Py() <<
3716  ", pz=" << part.first.Pz() <<
3717  ", n daughters = " << part.second->GetNDaughters() <<
3718  std::endl;
3719 
3720  // Skip partons that do not have any children
3721  // Unclear how this can happen, it would seem that this parton were not fragmented by the generator
3722  if (part.second->GetNDaughters() == 0) continue;
3723 
3724  // Look for highest momentum parton
3725  if (highestPartonPt < part.first.Pt()) {
3726  highestPartonPt = part.first.Pt();
3727  PdgHighParton = PdgCode;
3728  }
3729 
3730  // Skip partons that are not HF
3731  if (PdgCode != 4 && PdgCode != 5 && PdgCode != -4 && PdgCode != -5) continue;
3732 
3733  Bool_t lastInPartonShower = kTRUE;
3734  Bool_t hadronDaughter = kFALSE;
3735  for (Int_t daughterIndex = part.second->GetFirstDaughter(); daughterIndex <= part.second->GetLastDaughter(); daughterIndex++){
3736  if (daughterIndex < 0) {
3737  AliDebugStream(5) << "Could not find daughter index!" << std::endl;
3738  continue;
3739  }
3740  AliAODMCParticle *daughter = fMCContainer->GetMCParticle(daughterIndex);
3741  if (!daughter) {
3742  AliDebugStream(5) << "Could not find particle with index " << daughterIndex << "!" << std::endl;
3743  continue;
3744  }
3745  Int_t daughterPdgCode = daughter->GetPdgCode();
3746  if (daughter->GetMother() != it.current_index()) {
3747  AliDebugStream(5) << "Particle " << daughterIndex << " with pdg=" << daughterPdgCode <<
3748  ", px=" << daughter->Px() <<
3749  ", py=" << daughter->Py() <<
3750  ", pz=" << daughter->Pz() <<
3751  ", is not a daughter of " << it.current_index() <<
3752  "!" << std::endl;
3753  continue;
3754  }
3755 
3756  AliDebugStream(5) << "Found daughter " << daughterIndex <<
3757  " with pdg=" << daughterPdgCode <<
3758  ", px=" << daughter->Px() <<
3759  ", py=" << daughter->Py() <<
3760  ", pz=" << daughter->Pz() <<
3761  std::endl;
3762  // Codes between 81 and 100 are for internal MC code use, they may be intermediate states used e.g. in hadronizaion models
3763  if (daughterPdgCode == PdgCode) lastInPartonShower = kFALSE; // this parton is not the last parton in the shower
3764  if (TMath::Abs(daughterPdgCode) >= 111 || (daughterPdgCode >= 81 && daughterPdgCode <= 100)) hadronDaughter = kTRUE;
3765  }
3766  if (hadronDaughter) {
3767  AliDebugStream(5) << "This particle has at least a hadron among its daughters!" << std::endl;
3768  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;
3769  }
3770  else {
3771  AliDebugStream(5) << "This particle does not have hadrons among its daughters!" << std::endl;
3772  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;
3773  continue;
3774  }
3775 
3776  if (PdgCode == 4 || PdgCode == -4) {
3777  fHistManager.FillTH1("fHistCharmPt", part.first.Pt());
3778  fHistManager.FillTH1("fHistCharmEta", part.first.Eta());
3779  fHistManager.FillTH1("fHistCharmPhi", part.first.Phi_0_2pi());
3780  }
3781  else if (PdgCode == 5 || PdgCode == -5) {
3782  fHistManager.FillTH1("fHistBottomPt", part.first.Pt());
3783  fHistManager.FillTH1("fHistBottomEta", part.first.Eta());
3784  fHistManager.FillTH1("fHistBottomPhi", part.first.Phi_0_2pi());
3785  }
3786  nHQ++;
3787  }
3788  fHistManager.FillTH1("fHistNHeavyQuarks", nHQ);
3789  fHistManager.FillTH1("fHistHighestPartonPt",highestPartonPt);
3790  Int_t partonType = 0;
3791  Int_t absPdgHighParton = TMath::Abs(PdgHighParton);
3792  if (absPdgHighParton == 9 || absPdgHighParton == 21) {
3793  partonType = 7;
3794  }
3795  else {
3796  partonType = absPdgHighParton;
3797  }
3798  fHistManager.FillTH1("fHistHighestPartonType",partonType);
3799 }
3800 
3809 {
3810  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg));
3811 
3812  Double_t mass = part->Mass();
3813 
3814  // To make sure that the PDG mass value is not at the edge of a bin
3815  if (nbins % 2 == 0) {
3816  minMass = mass - range / 2 - range / nbins / 2;
3817  maxMass = mass + range / 2 - range / nbins / 2;
3818  }
3819  else {
3820  minMass = mass - range / 2;
3821  maxMass = mass + range / 2;
3822  }
3823 }
3824 
3831 {
3832  static TString label;
3833  label = "";
3834 
3835  if (bitmap & BIT(AliRDHFCuts::kNotSelTrigger)) {
3836  label = "NotSelTrigger";
3837  bitmap &= ~BIT(AliRDHFCuts::kNotSelTrigger);
3838  return label.Data();
3839  }
3840  if (bitmap & BIT(AliRDHFCuts::kNoVertex)) {
3841  label = "NoVertex";
3842  bitmap &= ~BIT(AliRDHFCuts::kNoVertex);
3843  return label.Data();
3844  }
3845  if (bitmap & BIT(AliRDHFCuts::kTooFewVtxContrib)) {
3846  label = "TooFewVtxContrib";
3847  bitmap &= ~BIT(AliRDHFCuts::kTooFewVtxContrib);
3848  return label.Data();
3849  }
3850  if (bitmap & BIT(AliRDHFCuts::kZVtxOutFid)) {
3851  label = "ZVtxOutFid";
3852  bitmap &= ~BIT(AliRDHFCuts::kZVtxOutFid);
3853  return label.Data();
3854  }
3855  if (bitmap & BIT(AliRDHFCuts::kPileup)) {
3856  label = "Pileup";
3857  bitmap &= ~BIT(AliRDHFCuts::kPileup);
3858  return label.Data();
3859  }
3860  if (bitmap & BIT(AliRDHFCuts::kOutsideCentrality)) {
3861  label = "OutsideCentrality";
3862  bitmap &= ~BIT(AliRDHFCuts::kOutsideCentrality);
3863  return label.Data();
3864  }
3865  if (bitmap & BIT(AliRDHFCuts::kPhysicsSelection)) {
3866  label = "PhysicsSelection";
3867  bitmap &= ~BIT(AliRDHFCuts::kPhysicsSelection);
3868  return label.Data();
3869  }
3870  if (bitmap & BIT(AliRDHFCuts::kBadSPDVertex)) {
3871  label = "BadSPDVertex";
3872  bitmap &= ~BIT(AliRDHFCuts::kBadSPDVertex);
3873  return label.Data();
3874  }
3875  if (bitmap & BIT(AliRDHFCuts::kZVtxSPDOutFid)) {
3876  label = "ZVtxSPDOutFid";
3877  bitmap &= ~BIT(AliRDHFCuts::kZVtxSPDOutFid);
3878  return label.Data();
3879  }
3880  if (bitmap & BIT(AliRDHFCuts::kCentralityFlattening)) {
3881  label = "CentralityFlattening";
3882  bitmap &= ~BIT(AliRDHFCuts::kCentralityFlattening);
3883  return label.Data();
3884  }
3885  if (bitmap & BIT(AliRDHFCuts::kBadTrackV0Correl)) {
3886  label = "BadTrackV0Correl";
3887  bitmap &= ~BIT(AliRDHFCuts::kBadTrackV0Correl);
3888  return label.Data();
3889  }
3890 
3891  return label.Data();
3892 }
3893 
3900 {
3901  if (handler->GetDataSlotNumber() >= 0 && handler->GetOutputObject()) {
3902  PostData(handler->GetDataSlotNumber(), handler->GetOutputObject());
3903  return handler->GetDataSlotNumber();
3904  }
3905  else {
3906  return -1;
3907  }
3908 }
3909 
3919 {
3920  // Get the pointer to the existing analysis manager via the static access method
3921  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
3922  if (!mgr) {
3923  ::Error("AddTaskDmesonJetsSub", "No analysis manager to connect to.");
3924  return NULL;
3925  }
3926 
3927  // Check the analysis type using the event handlers connected to the analysis manager
3928  AliVEventHandler* handler = mgr->GetInputEventHandler();
3929  if (!handler) {
3930  ::Error("AddTaskDmesonJetsSub", "This task requires an input event handler");
3931  return NULL;
3932  }
3933 
3935 
3936  if (handler->InheritsFrom("AliESDInputHandler")) {
3937  dataType = kESD;
3938  }
3939  else if (handler->InheritsFrom("AliAODInputHandler")) {
3940  dataType = kAOD;
3941  }
3942 
3943  // Init the task and do settings
3944  if (ntracks == "usedefault") {
3945  if (dataType == kESD) {
3946  ntracks = "Tracks";
3947  }
3948  else if (dataType == kAOD) {
3949  ntracks = "tracks";
3950  }
3951  else {
3952  ntracks = "";
3953  }
3954  }
3955 
3956  if (nclusters == "usedefault") {
3957  if (dataType == kESD) {
3958  nclusters = "CaloClusters";
3959  }
3960  else if (dataType == kAOD) {
3961  nclusters = "caloClusters";
3962  }
3963  else {
3964  nclusters = "";
3965  }
3966  }
3967 
3968  if (nMCpart == "usedefault") {
3969  nMCpart = "mcparticles"; // Always needs AliAODMCParticle objects
3970  }
3971 
3972  TString name("AliAnalysisTaskDmesonJetsSub");
3973  if (strcmp(suffix, "") != 0) {
3974  name += TString::Format("_%s", suffix.Data());
3975  }
3976 
3977  AliAnalysisTaskDmesonJetsSub* jetTask = new AliAnalysisTaskDmesonJetsSub(name, nMaxTrees);
3978 
3979  if (!ntracks.IsNull()) {
3980  AliHFTrackContainer* trackCont = new AliHFTrackContainer(ntracks);
3981  jetTask->AdoptParticleContainer(trackCont);
3982  }
3983 
3984  if (!nMCpart.IsNull()) {
3985  AliMCParticleContainer* partCont = new AliHFAODMCParticleContainer(nMCpart);
3986  partCont->SetEtaLimits(-1.5, 1.5);
3987  partCont->SetPtLimits(0, 1000);
3988  jetTask->AdoptParticleContainer(partCont);
3989  }
3990 
3991  jetTask->AddClusterContainer(nclusters.Data());
3992 
3993  // Final settings, pass to manager and set the containers
3994  mgr->AddTask(jetTask);
3995 
3996  // Create containers for input/output
3997  AliAnalysisDataContainer* cinput1 = mgr->GetCommonInputContainer();
3998  TString contname1(name);
3999  contname1 += "_histos";
4000  AliAnalysisDataContainer* coutput1 = mgr->CreateContainer(contname1.Data(),
4001  TList::Class(), AliAnalysisManager::kOutputContainer,
4002  Form("%s", AliAnalysisManager::GetCommonFileName()));
4003 
4004  mgr->ConnectInput(jetTask, 0, cinput1);
4005  mgr->ConnectOutput(jetTask, 1, coutput1);
4006 
4007  for (Int_t i = 0; i < nMaxTrees; i++) {
4008  TString contname = TString::Format("%s_tree_%d", name.Data(), i);
4009  AliAnalysisDataContainer *coutput = mgr->CreateContainer(contname.Data(),
4010  TTree::Class(),AliAnalysisManager::kOutputContainer,
4011  Form("%s", AliAnalysisManager::GetCommonFileName()));
4012  mgr->ConnectOutput(jetTask, 2+i, coutput);
4013  }
4014  return jetTask;
4015 }
Double_t fSoftPionPt
! Transverse momentum of the soft pion of the D* candidate
TString fBranchName
AOD branch where the D meson candidate are found.
Class that encapsulates event properties in a very compact structure (useful for pp simulation analys...
Int_t fDataSlotNumber
! Data slot where the tree output is posted
Int_t pdg
static std::pair< AliAnalysisTaskDmesonJetsSub::EMesonOrigin_t, AliAODMCParticle * > IsPromptCharm(const AliAODMCParticle *part, TClonesArray *mcArray)
Double32_t fDeltaInvMass
< Difference between the Kpipi and the Kpi invariant masses in GeV/c2
void Reset()
Reset all fields to their default values.
friend bool operator==(const AnalysisEngine &lhs, const AnalysisEngine &rhs)
AliEMCALGeometry * fGeom
!emcal geometry
void SetDMesonCandidate(AliAODRecoDecay *c)
virtual void Set(const AliDmesonJetInfo &source, std::string n)
virtual void Set(const AliDmesonJetInfo &source)
Bool_t fInhibit
!inhibit execution of the task
Double32_t fCorrZ
Z of the D meson after subtracting average background.
virtual Int_t Run()
double Double_t
Definition: External.C:58
EventInfo * fEventInfo
! Object conatining the event information (centrality, pt hard, weight, etc.)
Bool_t ExtractRecoDecayAttributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
EDataType_t
Switch for the data type.
UInt_t fEnabledAxis
Use bit defined in EAxis_t to enable axis in the THnSparse.
virtual ~AliAnalysisTaskDmesonJetsSub()
This is the standard destructor.
virtual AliAODMCParticle * GetMCParticle(Int_t i=-1) const
void Getd0MeasMinusExpProng(Int_t ip, Double_t magf, Double_t &d0diff, Double_t &errd0diff) const
const char * title
Definition: MakeQAPdf.C:27
std::map< int, AliDmesonJetInfo > * fDmesonJets
! Array containing the D meson jets
Double_t fTrackEfficiency
! Artificial tracking inefficiency (0...1) -> set automatically at ExecOnce by AliAnalysisTaskDmesonJ...
std::map< int, AliDmesonJetInfo > fDmesonJets
! Array containing the D meson jets
AliHFJetDefinition * AddJetDefinition(EJetType_t type, Double_t r, EJetAlgo_t algo, ERecoScheme_t reco)
Double_t fEPV0
!event plane V0
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
const TObjArray & GetDaughterList() const
AliTLorentzVector fD
! 4-momentum of the D meson candidate
bidirectional stl iterator over the EMCAL iterable container
std::list< AnalysisEngine > fAnalysisEngines
Array of analysis parameters.
static AliAODMCParticle * FindParticleOrigin(const AliAODMCParticle *part, TClonesArray *mcArray, EFindParticleOriginMode_t mode, const std::set< UInt_t > &pdgSet)
Double_t fMinMass
! Min mass in histogram axis
Bool_t ExtractDstarAttributes(const AliAODRecoCascadeHF *DstarCand, AliDmesonJetInfo &DmesonJet, UInt_t i)
virtual void Set(const AliDmesonJetInfo &source)
Invariant mass of the D0 meson candidate in GeV/c2.
Double_t fJetGhostArea
Area of the ghost particles.
Double_t mass
Container with name, TClonesArray and cuts for particles.
Bool_t ExtractD0Attributes(const AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, UInt_t i)
TArrayI fPDGdaughters
List of the PDG code of the daughters.
EventInfo fEventInfo
! Event info (centrality, weight, pt hard etc.)
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
AliAODMCParticle * fParton
! pointer to the parton in the shower tree of the D meson (only for particle level D mesons) ...
static Int_t CheckMatchingAODdeltaAODevents()
AnalysisEngine()
This is the default constructor, used for ROOT I/O purposes.
virtual void UserCreateOutputObjects()
Creates the output containers.
Lightweight class that encapsulates D meson jets.
Double_t fTrackEfficiency
Artificial tracking inefficiency (0...1)
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
Int_t fDataSlotNumber
! Data slot where the tree output is posted
static UShort_t GetRejectionReasonBitPosition(UInt_t rejectionReason)
Returns the highest bit in the rejection map as reason why the object was rejected.
Double_t fMinNeutralPt
Minimum pt of the leading neutral particle (or cluster)
AliAnalysisTaskDmesonJetsSub()
This is the default constructor, used for ROOT I/O purposes.
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:122
Int_t fNOutputTrees
Maximum number of output trees.
Bool_t FindJet(AliAODRecoDecayHF2Prong *Dcand, AliDmesonJetInfo &DmesonJet, AliHFJetDefinition &jetDef, Int_t numcand)
static AliAnalysisTaskDmesonJetsSub * AddTaskDmesonJetsSub(TString ntracks="usedefault", TString nclusters="usedefault", TString nMCpart="", Int_t nMaxTrees=2, TString suffix="")
Bool_t fRejectISR
Reject initial state radiation.
void SetEtaLimits(Double_t min, Double_t max)
const AliVEvent * fEvent
! pointer to the ESD/AOD event
friend bool operator<(const AnalysisEngine &lhs, const AnalysisEngine &rhs)
Double_t fInvMass2Prong
! 2-prong mass of the D* candidate (w/o the soft pion)
static FJRecoScheme ConvertToFJRecoScheme(ERecoScheme_t reco)
Bool_t fD0D0bar
! kTRUE if selected both as D0 and D0bar
static TString GenerateJetName(EJetType_t jetType, EJetAlgo_t jetAlgo, ERecoScheme_t recoScheme, Double_t radius, AliParticleContainer *partCont, AliClusterContainer *clusCont, TString tag)
static void CalculateMassLimits(Double_t range, Int_t pdg, Int_t nbins, Double_t &minMass, Double_t &maxMass)
Double32_t fSelectionType
Selection type: D0, D0bar, both.
AliEventNotFound(const std::string &class_name, const std::string &method_name)
AliRDHFCuts * fRDHFCuts
D meson candidates cuts.
void SetR(Double_t r)
Definition: AliFJWrapper.h:127
Select tracks based on specific prescriptions of HF analysis.
ERecoScheme_t fRecoScheme
Jet recombination scheme (pt scheme, E scheme, ...)
AliHFJetDefinition & operator=(const AliHFJetDefinition &source)
Output handler for D meson jet analysis.
Lightweight class that encapsulates D meson jets.
EBeamType_t fForceBeamType
forced beam type
UInt_t fAcceptedDecay
Bit mask with D meson decays that are accepted (only used for particle-level analysis) ...
const Int_t nPtBins
AliHFAODMCParticleContainer * fMCContainer
! MC particle container
virtual Bool_t AcceptObject(Int_t i, UInt_t &rejectionReason) const =0
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:34
std::vector< AliTrackContainer * > fTrackContainers
! Track containers
AliDmesonInfoSummary * fCurrentDmesonJetInfo
! Current D meson jet info
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
Create a new TH2 within the container.
TObject * FindObject(const char *name) const
Find an object inside the container.
int Int_t
Definition: External.C:63
Double_t GetDistance(const AliJetInfo &jet, Double_t &deta, Double_t &dphi) const
void RunParticleLevelAnalysis()
Run a particle level analysis.
void Init(const AliEMCALGeometry *const geom, Int_t runNumber)
Initialize the analysis engine.
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
Double_t CosThetaStarD0bar() const
angle of K
AliAODMCParticle * fAncestor
! pointer to the ancestor particle in the shower tree of the D meson (only for particle level D meson...
AliDmesonJetInfo & operator=(const AliDmesonJetInfo &source)
UInt_t fRejectedOrigin
Bit mask with D meson origins that are rejected (used for MC analysis, i.e. signal-only, background-only and particle-level)
std::vector< DMESONTYPE > fCurrentDmesonInfo
! Current D meson jet info
Int_t fMCLabel
! MC label, i.e. index of the generator level D meson (only for detector level D meson candidates) ...
Double32_t fPt
Transverse momentum of the jet in GeV/c.
AnalysisEngine * AddAnalysisEngine(ECandidateType_t type, TString cutfname, TString cuttype, EMCMode_t bkgMode, EJetType_t jettype, Double_t jetradius, TString rhoName="")
void IterativeDeclustering(Int_t jetnum, Double_t type, AliHFJetDefinition &jetDef, Double_t invm)
Base task in the EMCAL framework (lighter version of AliAnalysisTaskEmcal)
friend bool operator<(const AliHFJetDefinition &lhs, const AliHFJetDefinition &rhs)
Double_t Phi_0_2pi() const
UInt_t fEnabledAxis
! Use bit defined in EAxis_t to enable axis in the THnSparse
Bool_t fApplyKinematicCuts
Apply jet kinematic cuts.
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:121
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
Bool_t IsAnyJetInAcceptance(const AliDmesonJetInfo &dMesonJet) const
AliHFJetDefinition()
This is the default constructor, used for ROOT I/O purposes.
Bool_t fRejectISR
! Reject initial state radiation
AliRDHFCuts * LoadDMesonCutsFromFile(TString cutfname, TString cutsname)
Base class for container structures within the EMCAL framework.
AliAODTrack * GetBachelor() const
TClonesArray * fCandidateArray
! D meson candidate array
Bool_t fD0Extended
Store extended information in the tree (only for D0 mesons)
std::vector< AliClusterContainer * > fClusterContainers
! Cluster containers
static OutputHandlerTTreeExtendedBase * GenerateOutputHandler(AnalysisEngine *eng)
EJetType_t fJetType
Jet type (charged, full, neutral)
Bool_t fD0Extended
! Store extended information in the tree (only for D0 mesons)
TString fRhoName
Name of the object that holds the average background value.
AliFJWrapper * fFastJetWrapper
! Fastjet wrapper
AnalysisEngine & operator=(const AnalysisEngine &source)
TClonesArray * GetArray() const
void FillPartonLevelHistograms()
Fill histograms with parton-level information.
std::string fClassName
Class name where the event was not found.
virtual void Clear(const Option_t *="")
Double_t fMaxChargedPt
Transverse momentum of the leading charged particle (or track)
void RunDetectorLevelAnalysis()
Run a detector level analysis.
Int_t mode
Definition: anaM.C:41
std::vector< AliHFJetDefinition > fJetDefinitions
Jet definitions.
void SetPtLimits(Double_t min, Double_t max)
AliJetInfoSummary ** fCurrentJetInfo
! Current jet info
const char * GetName() const
Generate a name for this jet definition.
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
std::map< std::string, std::vector< JETTYPE > > fCurrentJetInfo
! Current jet info
void AddInputVectors(AliEmcalContainer *cont, Int_t offset, TH2 *rejectHist=0, Double_t eff=0.)
Double32_t fR
Distance between D meson and jet axis.
const AliMCParticleIterableMomentumContainer all_momentum() const
Select MC particles based on specific prescriptions of HF analysis.
virtual void Set(const AliDmesonJetInfo &source, std::string n)
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:125
virtual void BuildOutputObject(const char *taskName)
Int_t nbinsx
void Print() const
Prints the content of this object in the standard output.
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
static const char * GetHFEventRejectionReasonLabel(UInt_t &bitmap)
Double_t fMaxChargedPt
Maximum pt of the leading charged particle (or track)
static FJJetAlgo ConvertToFJAlgo(EJetAlgo_t algo)
std::map< std::string, AliJetInfo > fJets
! list of jets
decay
Definition: HFPtSpectrum.C: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:293
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:314
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.
const AliEmcalIterableMomentumContainer all_momentum() const
Create an iterable container interface over all objects in the EMCAL container.
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