AliPhysics  a60a912 (a60a912)
AliAnalysisTaskHFSubstructure.cxx
Go to the documentation of this file.
1 /*************************************************************************
2 * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15 
16 //C++
17 #include <sstream>
18 #include <array>
19 
20 // Root
21 #include <TClonesArray.h>
22 #include <TDatabasePDG.h>
23 #include <TParticlePDG.h>
24 #include <TVector3.h>
25 #include <THnSparse.h>
26 #include <TParticle.h>
27 #include <TMath.h>
28 #include <THashList.h>
29 #include <TFile.h>
30 #include <TRandom3.h>
31 #include <TGrid.h>
32 #include <TObjArray.h>
33 
34 // Aliroot general
35 #include "AliLog.h"
36 #include "AliEMCALGeometry.h"
37 #include "AliAnalysisManager.h"
38 #include "AliVEventHandler.h"
39 #include "AliAnalysisDataSlot.h"
40 #include "AliAnalysisDataContainer.h"
41 
42 // Aliroot HF
43 #include "AliAODRecoCascadeHF.h"
44 #include "AliAODRecoDecayHF2Prong.h"
45 #include "AliRDHFCutsD0toKpi.h"
46 #include "AliRDHFCutsDStartoKpipi.h"
47 #include "AliHFAODMCParticleContainer.h"
48 #include "AliHFTrackContainer.h"
49 #include "AliAnalysisVertexingHF.h"
50 
51 // Aliroot EMCal jet framework
52 #include "AliEmcalJetTask.h"
53 #include "AliEmcalJet.h"
54 #include "AliJetContainer.h"
55 #include "AliParticleContainer.h"
56 #include "AliEmcalParticle.h"
57 #include "AliFJWrapper.h"
58 #include "AliRhoParameter.h"
59 
60 #include "AliAnalysisTaskHFSubstructure.h"
61 
62 
63 #include "FJ_includes.h"
64 
65 
66 
67 
68 
69 using std::cout;
70 using std::endl;
71 
72 
73 ClassImp(AliAnalysisTaskHFSubstructure)
74 
75 //________________________________________________________________________
76 AliAnalysisTaskHFSubstructure::AliAnalysisTaskHFSubstructure() :
77 AliAnalysisTaskEmcal("AliAnalysisTaskHFSubstructure", kTRUE),
78  fJetShapeType(kData),
79  fECandidateType(kD0toKpi),
80  fIncludeInclusive(kFALSE),
81  fIsBDecay(kFALSE),
82  fRejectISR(kFALSE),
83  fPromptReject(kFALSE),
84  fAlienConnect(kFALSE),
85  fBranchName(0),
86  fCutsType(0),
87  fCandidatePDG(421),
88  fRejectedOrigin(0),
89  fAcceptedDecay(kDecayD0toKpi),
90  fJetRadius(0.4),
91  fJetMinPt(0.0),
92  fCandidateArray(0),
93  fAodEvent(0),
94  fRDHFCuts(0),
95  fFastJetWrapper(0),
96  fFastJetWrapper_Truth(0),
97  fShapesVar_Splittings_DeltaR(0),
98  fShapesVar_Splittings_DeltaR_Truth(0),
99  fShapesVar_Splittings_Zg(0),
100  fShapesVar_Splittings_Zg_Truth(0),
101  fShapesVar_Splittings_LeadingSubJetpT(0),
102  fShapesVar_Splittings_LeadingSubJetpT_Truth(0),
103  fShapesVar_Splittings_HardestSubJetD0(0),
104  fShapesVar_Splittings_HardestSubJetD0_Truth(0),
105  fShapesVar_Splittings_RadiatorE(0),
106  fShapesVar_Splittings_RadiatorE_Truth(0),
107  fShapesVar_Splittings_RadiatorpT(0),
108  fShapesVar_Splittings_RadiatorpT_Truth(0),
109 
110  fTreeResponseMatrixAxis(0),
111  fTreeSplittings(0),
112  fhEvent(0x0)
113 
114 {
115  for(Int_t i=0;i<nVar;i++){
116  fShapesVar[i]=0;
117  }
118  SetMakeGeneralHistograms(kTRUE);
119  DefineOutput(1, TList::Class());
120  DefineOutput(2, TTree::Class());
121  DefineOutput(3, TTree::Class());
122 }
123 
124 //________________________________________________________________________
125 AliAnalysisTaskHFSubstructure::AliAnalysisTaskHFSubstructure(const char *name) :
126  AliAnalysisTaskEmcal(name, kTRUE),
127  fJetShapeType(kData),
128  fECandidateType(kD0toKpi),
129  fIncludeInclusive(kFALSE),
130  fIsBDecay(kFALSE),
131  fRejectISR(kFALSE),
132  fPromptReject(kFALSE),
133  fAlienConnect(kFALSE),
134  fBranchName(0),
135  fCutsType(0),
136  fCandidatePDG(421),
137  fRejectedOrigin(0),
138  fAcceptedDecay(kDecayD0toKpi),
139  fJetRadius(0.4),
140  fJetMinPt(0.0),
141  fCandidateArray(0),
142  fAodEvent(0),
143  fRDHFCuts(0),
144  fFastJetWrapper(0),
145  fFastJetWrapper_Truth(0),
146  fShapesVar_Splittings_DeltaR(0),
147  fShapesVar_Splittings_DeltaR_Truth(0),
148  fShapesVar_Splittings_Zg(0),
149  fShapesVar_Splittings_Zg_Truth(0),
150  fShapesVar_Splittings_LeadingSubJetpT(0),
151  fShapesVar_Splittings_LeadingSubJetpT_Truth(0),
152  fShapesVar_Splittings_HardestSubJetD0(0),
153  fShapesVar_Splittings_HardestSubJetD0_Truth(0),
154  fShapesVar_Splittings_RadiatorE(0),
155  fShapesVar_Splittings_RadiatorE_Truth(0),
156  fShapesVar_Splittings_RadiatorpT(0),
157  fShapesVar_Splittings_RadiatorpT_Truth(0),
158 
159  fTreeResponseMatrixAxis(0),
160  fTreeSplittings(0),
161  fhEvent(0x0)
162 {
163  // Standard constructor.
164  for(Int_t i=0;i<nVar;i++){
165  fShapesVar[i]=0;
166  }
167  SetMakeGeneralHistograms(kTRUE);
168  DefineOutput(1, TList::Class());
169  DefineOutput(2, TTree::Class());
170  DefineOutput(3, TTree::Class());
171 }
172 
173 //________________________________________________________________________
174 AliAnalysisTaskHFSubstructure::~AliAnalysisTaskHFSubstructure()
175 {
176  // Destructor.
177 }
178 
179 //________________________________________________________________________
180 void AliAnalysisTaskHFSubstructure::UserCreateOutputObjects()
181 {
182  // Create user output.
183 
184  AliAnalysisTaskEmcal::UserCreateOutputObjects();
185 
186  Bool_t oldStatus = TH1::AddDirectoryStatus();
187  TH1::AddDirectory(kFALSE);
188  TH1::AddDirectory(oldStatus);
189  //create a tree used for the MC data and making a 4D response matrix
190 
191 
192 
193  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
194  fTreeResponseMatrixAxis = new TTree(nameoutput, nameoutput);
195  TString *fShapesVarNames = new TString [nVar];
196  fShapesVarNames[0] = "pT_Jet";
197  fShapesVarNames[1] = "pT_Jet_Truth";
198  fShapesVarNames[2] = "pT_D";
199  fShapesVarNames[3] = "pT_D_Truth";
200  fShapesVarNames[4] = "Inv_M_D";
201  fShapesVarNames[5] = "Inv_M_D_Truth";
202  fShapesVarNames[6] = "Flag_D";
203  fShapesVarNames[7] = "Flag_D_Truth";
204  fShapesVarNames[8] = "Prompt_PDG";
205  fShapesVarNames[9] = "Prompt_PDG_Truth";
206 
207  for(Int_t ivar=0; ivar < nVar; ivar++){
208  cout<<"looping over variables"<<endl;
209  fTreeResponseMatrixAxis->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/D", fShapesVarNames[ivar].Data()));
210  }
211 
212  const char* nameoutput_Splittings = GetOutputSlot(3)->GetContainer()->GetName();
213  fTreeSplittings = new TTree(nameoutput_Splittings, nameoutput_Splittings);
214  TString *fShapesVarNames_Splittings=new TString[nVar_Splittings];
215  fShapesVarNames_Splittings[0] = "DeltaR";
216  fShapesVarNames_Splittings[1] = "DeltaR_Truth";
217  fShapesVarNames_Splittings[2] = "Zg";
218  fShapesVarNames_Splittings[3] = "Zg_Truth";
219  fShapesVarNames_Splittings[4] = "LeadingSubJetpT";
220  fShapesVarNames_Splittings[5] = "LeadingSubJetpT_Truth";
221  fShapesVarNames_Splittings[6] = "HardestSubJetD0";
222  fShapesVarNames_Splittings[7] = "HardestSubJetD0_Truth";
223  fShapesVarNames_Splittings[8] = "RadiatorE";
224  fShapesVarNames_Splittings[9] = "RadiatorE_Truth";
225  fShapesVarNames_Splittings[10] = "RadiatorpT";
226  fShapesVarNames_Splittings[11] = "RadiatorpT_Truth";
227  fTreeSplittings->Branch(fShapesVarNames_Splittings[0].Data(), &fShapesVar_Splittings_DeltaR, 0,1);
228  fTreeSplittings->Branch(fShapesVarNames_Splittings[1].Data(), &fShapesVar_Splittings_DeltaR_Truth, 0,1);
229  fTreeSplittings->Branch(fShapesVarNames_Splittings[2].Data(), &fShapesVar_Splittings_Zg, 0,1);
230  fTreeSplittings->Branch(fShapesVarNames_Splittings[3].Data(), &fShapesVar_Splittings_Zg_Truth, 0,1);
231  fTreeSplittings->Branch(fShapesVarNames_Splittings[4].Data(), &fShapesVar_Splittings_LeadingSubJetpT, 0,1);
232  fTreeSplittings->Branch(fShapesVarNames_Splittings[5].Data(), &fShapesVar_Splittings_LeadingSubJetpT_Truth, 0,1);
233  fTreeSplittings->Branch(fShapesVarNames_Splittings[6].Data(), &fShapesVar_Splittings_HardestSubJetD0, 0,1);
234  fTreeSplittings->Branch(fShapesVarNames_Splittings[7].Data(), &fShapesVar_Splittings_HardestSubJetD0_Truth, 0,1);
235  fTreeSplittings->Branch(fShapesVarNames_Splittings[8].Data(), &fShapesVar_Splittings_RadiatorE, 0,1);
236  fTreeSplittings->Branch(fShapesVarNames_Splittings[9].Data(), &fShapesVar_Splittings_RadiatorE_Truth, 0,1);
237  fTreeSplittings->Branch(fShapesVarNames_Splittings[10].Data(), &fShapesVar_Splittings_RadiatorpT, 0,1);
238  fTreeSplittings->Branch(fShapesVarNames_Splittings[11].Data(), &fShapesVar_Splittings_RadiatorpT_Truth, 0,1);
239 
240  fhEvent=new TH1D("fhEvent","fhEvent",40,-0.5,39.5);
241  fOutput->Add(fhEvent);
242 
243  PostData(1,fOutput);
244  PostData(2,fTreeResponseMatrixAxis);
245  PostData(3,fTreeSplittings);
246  // delete [] fShapesVarNames;
247 }
248 
249 //________________________________________________________________________
250 Bool_t AliAnalysisTaskHFSubstructure::Run()
251 {
252  // Run analysis code here, if needed. It will be executed before FillHistograms().
253 
254 
255  return kTRUE;
256 }
257 //________________________________________________________________________
258 Bool_t AliAnalysisTaskHFSubstructure::FillHistograms()
259 {
260 
261  Int_t Matching_AOD_deltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
262  if (Matching_AOD_deltaAODlevel <= 0) return kTRUE;
263 
264  fhEvent->Fill(0);
265  //if(fCutsFileName.Contains("alien://") && fAlienConnect) TGrid::Connect("alien://");
266  //TFile* Cuts_File = TFile::Open(fCutsFileName);
267  //TString cutsname="D0toKpiCuts";
268  //if (fCutsType!="") cutsname += TString::Format("_%s", fCutsType.Data());
269  //fRDHFCuts = dynamic_cast<AliRDHFCuts*>(fCutsFile->Get(cutsname));
270 
271 
272  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
273 
274 
275  if(!fRDHFCuts->IsEventSelected(fAodEvent)) return kTRUE;
276  fhEvent->Fill(1);
277  fFastJetWrapper=new AliFJWrapper("fastjetwrapper","fastjetwrapper");
278  fFastJetWrapper->SetAreaType(fastjet::active_area);
279  fFastJetWrapper->SetGhostArea(0.005);
280  fFastJetWrapper->SetR(fJetRadius);
281  fFastJetWrapper->SetAlgorithm(fastjet::antikt_algorithm);
282  fFastJetWrapper->SetRecombScheme(static_cast<fastjet::RecombinationScheme>(0));
283 
284  fFastJetWrapper_Truth=new AliFJWrapper("fastjetwrapper_truth","fastjetwrapper_truth");
285 
286  if (fJetShapeType == kData || fJetShapeType == kDetSignal || fJetShapeType == kDetBackground || fJetShapeType == kDetReflection){
287 
288  fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject("D0toKpi"));
289 
290  AliHFAODMCParticleContainer *Particle_Container = NULL;
291  if (fJetShapeType != kData) Particle_Container = (AliHFAODMCParticleContainer *) GetParticleContainer(1);
292  // if (!Particle_Container) continue;
293 
294 
295  std::vector<AliAODRecoDecayHF2Prong*> D_Candidates_Vector;
296  D_Candidates_Vector.clear();
297 
298 
299  Int_t N_DMesons=0;
300  //AliHFTrackContainer* Track_Container=(AliHFTrackContainer *) GetTrackContainer(0);
301  AliHFTrackContainer* Track_Container=dynamic_cast<AliHFTrackContainer*>(GetTrackContainer(0));
302  if (!Track_Container) return kTRUE;
303  Track_Container->SetDMesonCandidate(NULL);
304  for (Int_t i_D = 0; i_D < fCandidateArray->GetEntriesFast(); i_D++) {
305  AliAODRecoDecayHF2Prong* D_Candidate = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(i_D));
306  if (!D_Candidate) continue;
307  if (!fRDHFCuts->IsInFiducialAcceptance(D_Candidate->Pt(), D_Candidate->Y(fCandidatePDG))) continue;
308 
309  Int_t Mass_Hypo_Type=fRDHFCuts->IsSelected(D_Candidate, AliRDHFCuts::kAll, fAodEvent);
310  Int_t N_Mass_Hypotheses=1;
311  if (Mass_Hypo_Type <= 0 || Mass_Hypo_Type>3) continue;
312  else if (Mass_Hypo_Type ==3) N_Mass_Hypotheses=2;
313  fhEvent->Fill(2);
314  Int_t Matched_Truth_Particle_PDG=0;
315  Int_t Is_Prompt_Correct_Quark_PDG=-1;
316  if (fJetShapeType != kData){
317 
318 
319  const Int_t D_Candidtae_N_Daughters=2;
320  Int_t D_Candidtae_Daughters_PDG[D_Candidtae_N_Daughters] = {211,321};
321  Int_t D_Candidate_MatchedTruth_Label = D_Candidate->MatchToMC(fCandidatePDG, Particle_Container->GetArray(), D_Candidtae_N_Daughters, D_Candidtae_Daughters_PDG);
322  Bool_t Is_Prompt_Correct_Quark=kFALSE;
323  if (D_Candidate_MatchedTruth_Label >= 0) {
324  AliAODMCParticle* Matched_Truth_Particle = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(D_Candidate_MatchedTruth_Label));
325  if (Matched_Truth_Particle) {
326  fhEvent->Fill(3);
327 
328  Int_t Matched_Truth_Particle_Mother_Label=Matched_Truth_Particle->GetMother();
329  while (Matched_Truth_Particle_Mother_Label >= 0) {
330  AliAODMCParticle* Matched_Truth_Particle_Mother = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Matched_Truth_Particle_Mother_Label));
331  if (Matched_Truth_Particle_Mother){
332  Int_t Original_Quark_PDG=4;
333  if (fIsBDecay) Original_Quark_PDG=5;
334  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode())==Original_Quark_PDG) Is_Prompt_Correct_Quark=kTRUE;
335  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode()) == 4){
336  Is_Prompt_Correct_Quark_PDG=4;
337  fhEvent->Fill(4);
338  }
339  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode()) == 5){
340  Is_Prompt_Correct_Quark_PDG=5;
341  fhEvent->Fill(5);
342  }
343  if (Matched_Truth_Particle_Mother_Label==Matched_Truth_Particle_Mother->GetMother()) break;
344  Matched_Truth_Particle_Mother_Label=Matched_Truth_Particle_Mother->GetMother();
345  }
346  else break;
347  //delete Matched_Truth_Particle_Mother;
348  }
349  Matched_Truth_Particle_PDG = Matched_Truth_Particle->PdgCode();
350  }
351  //delete Matched_Truth_Particle;
352  }
353  //else continue;
354  if (fPromptReject && !Is_Prompt_Correct_Quark) continue;
355 
356  //if (TMath::Abs(Matched_Truth_Particle_PDG)!=fCandidatePDG) continue;
357 
358  }
359 
360 
361  Double_t Inv_Mass_D=0.0;
362  if (Mass_Hypo_Type==1){
363  if (fJetShapeType==kData || (fJetShapeType==kDetSignal && Matched_Truth_Particle_PDG==fCandidatePDG) || (fJetShapeType==kDetBackground && Matched_Truth_Particle_PDG!=fCandidatePDG) || (fJetShapeType==kDetReflection && Matched_Truth_Particle_PDG==-fCandidatePDG)){
364  Inv_Mass_D=D_Candidate->InvMassD0();
365  fhEvent->Fill(6);
366  }
367  else{
368  fhEvent->Fill(7);
369  continue;
370  }
371  }
372 
373  if (Mass_Hypo_Type==2){
374  if (fJetShapeType==kData || (fJetShapeType==kDetSignal && Matched_Truth_Particle_PDG==-fCandidatePDG) || (fJetShapeType==kDetBackground && Matched_Truth_Particle_PDG!=-fCandidatePDG) || (fJetShapeType==kDetReflection && Matched_Truth_Particle_PDG==fCandidatePDG)){
375  Inv_Mass_D=D_Candidate->InvMassD0bar();
376  fhEvent->Fill(8);
377  }
378  else{
379  fhEvent->Fill(9);
380  continue;
381  }
382  }
383 
384 
385 
386  for (Int_t i_Mass_Hypotheses=0; i_Mass_Hypotheses<N_Mass_Hypotheses; i_Mass_Hypotheses++){
387 
388 
389  if (Mass_Hypo_Type==3){
390  if(i_Mass_Hypotheses==0){
391  if (fJetShapeType==kData || (fJetShapeType==kDetSignal && Matched_Truth_Particle_PDG==fCandidatePDG) || (fJetShapeType==kDetBackground && Matched_Truth_Particle_PDG!=fCandidatePDG) || (fJetShapeType==kDetReflection && Matched_Truth_Particle_PDG==-fCandidatePDG)){
392  Inv_Mass_D=D_Candidate->InvMassD0();
393  fhEvent->Fill(11);
394  }
395  else{
396  fhEvent->Fill(12);
397  continue;
398  }
399  }
400  if(i_Mass_Hypotheses==1){
401  if (fJetShapeType==kData || (fJetShapeType==kDetSignal && Matched_Truth_Particle_PDG==-fCandidatePDG) || (fJetShapeType==kDetBackground && Matched_Truth_Particle_PDG!=-fCandidatePDG) || (fJetShapeType==kDetReflection && Matched_Truth_Particle_PDG==fCandidatePDG)){
402  Inv_Mass_D=D_Candidate->InvMassD0bar();
403  fhEvent->Fill(13);
404  }
405  else{
406  fhEvent->Fill(14);
407  continue;
408  }
409  }
410  }
411 
412  fFastJetWrapper->Clear();
413  AliTLorentzVector D_Candidate_LorentzVector(0,0,0,0);
414  D_Candidate_LorentzVector.SetPtEtaPhiM(D_Candidate->Pt(), D_Candidate->Eta(), D_Candidate->Phi(), Inv_Mass_D);
415  fFastJetWrapper->AddInputVector(D_Candidate_LorentzVector.Px(), D_Candidate_LorentzVector.Py(), D_Candidate_LorentzVector.Pz(), D_Candidate_LorentzVector.E(), 0);
416  Track_Container->SetDMesonCandidate(D_Candidate);
417  AliAODTrack *Track = NULL;
418  for (Int_t i_Track=0; i_Track<Track_Container->GetNTracks(); i_Track++){
419  Track = static_cast<AliAODTrack*>(Track_Container->GetAcceptParticle(i_Track));
420  if(!Track) continue;
421  fFastJetWrapper->AddInputVector(Track->Px(), Track->Py(), Track->Pz(), Track->E(),i_Track+100);
422  }
423  // delete Track;
424 
425  fFastJetWrapper->Run();
426  std::vector<fastjet::PseudoJet> Inclusive_Jets = fFastJetWrapper->GetInclusiveJets();
427  for (UInt_t i_Jet=0; i_Jet < Inclusive_Jets.size(); i_Jet++){
428  Bool_t Is_D_Jet=kFALSE;
429  if (Inclusive_Jets[i_Jet].perp()<fJetMinPt) continue;
430  if (TMath::Abs(Inclusive_Jets[i_Jet].pseudorapidity()) > 0.9-fJetRadius) continue;
431  std::vector<fastjet::PseudoJet> Constituents(fFastJetWrapper->GetJetConstituents(i_Jet));
432  for (UInt_t i_Constituents = 0; i_Constituents < Constituents.size(); i_Constituents++) {
433  if (Constituents[i_Constituents].user_index() == 0) {
434  Is_D_Jet = kTRUE;
435  }
436  }
437  if (!Is_D_Jet) continue;
438  fhEvent->Fill(10);
439  std::vector<Double_t> Splittings_Zg;
440  std::vector<Double_t> Splittings_DeltaR;
441  std::vector<Double_t> Splittings_LeadingSubJetpT;
442  std::vector<Double_t> Splittings_HardestSubJetD0;
443  std::vector<Double_t> Splittings_RadiatorE;
444  std::vector<Double_t> Splittings_RadiatorpT;
445 
446  Bool_t Is_D_SubJet=kFALSE;
447  fastjet::JetDefinition Jet_Definition(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
448  try{
449  std::vector<fastjet::PseudoJet> Reclustered_Particles(fFastJetWrapper->GetJetConstituents(i_Jet));
450  fastjet::ClusterSequence Cluster_Sequence_CA(Reclustered_Particles, Jet_Definition);
451  std::vector<fastjet::PseudoJet> Reclustered_Jet = Cluster_Sequence_CA.inclusive_jets(0.0);
452  Reclustered_Jet = sorted_by_pt(Reclustered_Jet);
453 
454 
455  fastjet::PseudoJet Daughter_Jet = Reclustered_Jet[0];
456  fastjet::PseudoJet Parent_SubJet_1;
457  fastjet::PseudoJet Parent_SubJet_2;
458 
459  while(Daughter_Jet.has_parents(Parent_SubJet_1,Parent_SubJet_2)){
460  if(Parent_SubJet_1.perp() < Parent_SubJet_2.perp()) std::swap(Parent_SubJet_1,Parent_SubJet_2);
461  Splittings_LeadingSubJetpT.push_back(Parent_SubJet_1.perp());
462  vector < fastjet::PseudoJet > Hard_SubJet_Constituents = sorted_by_pt(Parent_SubJet_1.constituents());
463  Is_D_SubJet=kFALSE;
464  for(UInt_t j=0;j<Hard_SubJet_Constituents.size();j++){
465  if(Hard_SubJet_Constituents[j].user_index()==0) Is_D_SubJet=kTRUE;
466  }
467 
468  if (!Is_D_SubJet) Splittings_HardestSubJetD0.push_back(1.0);
469  else Splittings_HardestSubJetD0.push_back(2.0);
470 
471  // if(!Is_D_SubJet) std::swap(Parent_SubJet_1,Parent_SubJet_2);
472  Splittings_DeltaR.push_back(Parent_SubJet_1.delta_R(Parent_SubJet_2));
473  Splittings_Zg.push_back(Parent_SubJet_2.perp()/(Parent_SubJet_1.perp()+Parent_SubJet_2.perp()));
474  Splittings_RadiatorE.push_back(Daughter_Jet.E());
475  Splittings_RadiatorpT.push_back(Daughter_Jet.perp());
476  Daughter_Jet=Parent_SubJet_1;
477  }
478 
479 
480  } catch (fastjet::Error) { /*return -1;*/ }
481 
482 
483 
484  // for (Int_t i_Mass_Hypotheses=0; i_Mass_Hypotheses<N_Mass_Hypotheses; i_Mass_Hypotheses++){
485 
486 
487  D_Candidates_Vector.push_back(D_Candidate);
488  if (i_Mass_Hypotheses==0) fhEvent->Fill(15);
489  N_DMesons++;
490 
491  Double_t Flag_D=-1.0;
492  if (Mass_Hypo_Type ==1) Flag_D=1.0;
493  else if (Mass_Hypo_Type ==2) Flag_D=2.0;
494  else if (Mass_Hypo_Type ==3 && i_Mass_Hypotheses==0) Flag_D=3.0;
495  else if (Mass_Hypo_Type ==3 && i_Mass_Hypotheses==1) Flag_D=4.0;
496 
497 
498  fShapesVar[0] = Inclusive_Jets[i_Jet].perp();
499  fShapesVar[1] = 0.0;
500  fShapesVar[2] = D_Candidate->Pt();
501  fShapesVar[3] = 0.0;
502  fShapesVar[4] = Inv_Mass_D;
503  fShapesVar[5] = 0.0;
504  fShapesVar[6] = Flag_D;
505  fShapesVar[7] = 0.0;
506  fShapesVar[8] = Is_Prompt_Correct_Quark_PDG;
507  fShapesVar[9] = 0.0;
508 
509 
510  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR);
511  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR);
512  fShapesVar_Splittings_Zg.push_back(Splittings_Zg);
513  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg);
514  fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT);
515  fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT);
516  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0);
517  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0);
518  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE);
519  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE);
520  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT);
521  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT);
522 
523  fTreeResponseMatrixAxis->Fill();
524  fTreeSplittings->Fill();
525  }
526  }
527  //delete D_Candidate;
528  }
529  if(N_DMesons==0) fhEvent->Fill(16);
530  if(N_DMesons==1) fhEvent->Fill(17);
531  if(N_DMesons==2) fhEvent->Fill(18);
532  if(N_DMesons==3) fhEvent->Fill(19);
533  if(N_DMesons==4) fhEvent->Fill(20);
534  if(N_DMesons==5) fhEvent->Fill(21);
535  if(N_DMesons==6) fhEvent->Fill(22);
536 
537  if (fIncludeInclusive){
538  fFastJetWrapper->Clear();
539 
540  for (UInt_t i_D_Found=0; i_D_Found<D_Candidates_Vector.size(); i_D_Found++){
541  fFastJetWrapper->AddInputVector(D_Candidates_Vector[i_D_Found]->Px(), D_Candidates_Vector[i_D_Found]->Py(), D_Candidates_Vector[i_D_Found]->Pz(), D_Candidates_Vector[i_D_Found]->AliAODRecoDecay::E(fCandidatePDG), i_D_Found);
542  //Track_Container->SetDMesonCandidate(D_Candidates_Vector[i_D_Found]);
543  }
544 
545  AliAODTrack *Track = NULL;
546  Bool_t DMeson_Daughter_Track=kFALSE;
547  for (Int_t i_Track=0; i_Track<Track_Container->GetNTracks(); i_Track++){
548  for (UInt_t i_D_Found=0; i_D_Found<D_Candidates_Vector.size(); i_D_Found++){
549  Track_Container->SetDMesonCandidate(D_Candidates_Vector[i_D_Found]);
550  Track = static_cast<AliAODTrack*>(Track_Container->GetAcceptParticle(i_Track));
551  if(!Track) DMeson_Daughter_Track=kTRUE;
552  }
553  if (DMeson_Daughter_Track) continue;
554  Track = static_cast<AliAODTrack*>(Track_Container->GetAcceptParticle(i_Track));
555  fFastJetWrapper->AddInputVector(Track->Px(), Track->Py(), Track->Pz(), Track->E(),i_Track+100);
556  }
557  //delete Track;
558  fFastJetWrapper->Run();
559  std::vector<fastjet::PseudoJet> Inclusive_Jets = fFastJetWrapper->GetInclusiveJets();
560  for (UInt_t i_Jet=0; i_Jet < Inclusive_Jets.size(); i_Jet++){
561  if (Inclusive_Jets[i_Jet].perp()<fJetMinPt) continue;
562  if (TMath::Abs(Inclusive_Jets[i_Jet].pseudorapidity()) > 0.9-fJetRadius) continue;
563  Bool_t Is_D_Jet = kFALSE;
564  std::vector<fastjet::PseudoJet> Constituents(fFastJetWrapper->GetJetConstituents(i_Jet));
565  for (UInt_t i_Constituents = 0; i_Constituents < Constituents.size(); i_Constituents++) {
566  for (UInt_t i_D_Found=0; i_D_Found<D_Candidates_Vector.size(); i_D_Found++){
567  if (Constituents[i_Constituents].user_index() == i_D_Found) {
568  Is_D_Jet = kTRUE;
569  }
570  }
571  }
572  if (Is_D_Jet) continue;
573  fhEvent->Fill(23);
574  std::vector<Double_t> Splittings_Zg;
575  std::vector<Double_t> Splittings_DeltaR;
576  std::vector<Double_t> Splittings_LeadingSubJetpT;
577  std::vector<Double_t> Splittings_HardestSubJetD0;
578  std::vector<Double_t> Splittings_RadiatorE;
579  std::vector<Double_t> Splittings_RadiatorpT;
580 
581 
582  fastjet::JetDefinition Jet_Definition(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
583 
584  try{
585  std::vector<fastjet::PseudoJet> Reclustered_Particles(fFastJetWrapper->GetJetConstituents(i_Jet));
586  fastjet::ClusterSequence Cluster_Sequence_CA(Reclustered_Particles, Jet_Definition);
587  std::vector<fastjet::PseudoJet> Reclustered_Jet = Cluster_Sequence_CA.inclusive_jets(0.0);
588  Reclustered_Jet = sorted_by_pt(Reclustered_Jet);
589 
590 
591  fastjet::PseudoJet Daughter_Jet = Reclustered_Jet[0];
592  fastjet::PseudoJet Parent_SubJet_1;
593  fastjet::PseudoJet Parent_SubJet_2;
594 
595  while(Daughter_Jet.has_parents(Parent_SubJet_1,Parent_SubJet_2)){
596  if(Parent_SubJet_1.perp() < Parent_SubJet_2.perp()) std::swap(Parent_SubJet_1,Parent_SubJet_2);
597  Splittings_LeadingSubJetpT.push_back(Parent_SubJet_1.perp());
598  Splittings_HardestSubJetD0.push_back(0.0);
599  Splittings_DeltaR.push_back(Parent_SubJet_1.delta_R(Parent_SubJet_2));
600  Splittings_Zg.push_back(Parent_SubJet_2.perp()/(Parent_SubJet_1.perp()+Parent_SubJet_2.perp()));
601  Splittings_RadiatorE.push_back(Daughter_Jet.E());
602  Splittings_RadiatorpT.push_back(Daughter_Jet.perp());
603  Daughter_Jet=Parent_SubJet_1;
604  }
605 
606 
607  } catch (fastjet::Error) { /*return -1;*/ }
608 
609 
610 
611  fShapesVar[0] = Inclusive_Jets[i_Jet].perp();
612  fShapesVar[1] = 0.0;
613  fShapesVar[2] = 0.0;
614  fShapesVar[3] = 0.0;
615  fShapesVar[4] = 0.0;
616  fShapesVar[5] = 0.0;
617  fShapesVar[6] = 0.0;
618  fShapesVar[7] = 0.0;
619  fShapesVar[8] = 0.0;
620  fShapesVar[9] = 0.0;
621 
622 
623  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR);
624  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR);
625  fShapesVar_Splittings_Zg.push_back(Splittings_Zg);
626  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg);
627  fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT);
628  fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT);
629  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0);
630  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0);
631  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE);
632  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE);
633  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT);
634  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT);
635 
636 
637  fTreeResponseMatrixAxis->Fill();
638  fTreeSplittings->Fill();
639 
640  }
641  }
642  // delete Track_Container;
643  // if (fJetShapeType != kData) delete Particle_Container;
644  }
645 
646 
647 
648  if (fJetShapeType == kTrueDet){
649 
650 
651  AliHFAODMCParticleContainer *Particle_Container = (AliHFAODMCParticleContainer*) GetParticleContainer(1);
652  Particle_Container->SetSpecialPDG(fCandidatePDG);
653  Particle_Container->SetRejectedOriginMap(fRejectedOrigin);
654  Particle_Container->SetAcceptedDecayMap(fAcceptedDecay);
655  Particle_Container->SetRejectISR(fRejectISR);
656 
657  std::vector<fastjet::PseudoJet> Inclusive_Jets_Truth;
658  std::vector<std::pair<Int_t, Int_t>> Inclusive_Jets_Truth_Labels;
659  std::vector<Int_t> Unmatched_Truth_Level_D;
660  Int_t NMatched_DMeson_Jets=0;
661 
662  if (Particle_Container->IsSpecialPDGFound()){
663 
664  fhEvent->Fill(2);
665 
666  fFastJetWrapper_Truth->SetAreaType(fastjet::active_area);
667  fFastJetWrapper_Truth->SetGhostArea(0.005);
668  fFastJetWrapper_Truth->SetR(fJetRadius);
669  fFastJetWrapper_Truth->SetAlgorithm(fastjet::antikt_algorithm);
670  fFastJetWrapper_Truth->SetRecombScheme(static_cast<fastjet::RecombinationScheme>(0));
671  fFastJetWrapper_Truth->Clear();
672 
673  AliAODMCParticle* Truth_Particle=NULL;
674  Int_t NTruthD=0;
675  for (Int_t i_Particle=0; i_Particle<Particle_Container->GetNParticles(); i_Particle++){
676  Truth_Particle = static_cast<AliAODMCParticle*>(Particle_Container->GetAcceptParticle(i_Particle));
677  if (!Truth_Particle) continue;
678  if (TMath::Abs(Truth_Particle->PdgCode())==fCandidatePDG){
679  std::pair<Int_t, Int_t> Inclusive_Jet_Truth_Labels;
680  Inclusive_Jet_Truth_Labels.first=Truth_Particle->GetLabel();
681  Inclusive_Jet_Truth_Labels.second=NTruthD;
682  Inclusive_Jets_Truth_Labels.push_back(Inclusive_Jet_Truth_Labels);
683  Unmatched_Truth_Level_D.push_back(NTruthD);
684  fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),NTruthD);
685  NTruthD++;
686  }
687  else fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),i_Particle+100);
688  }
689  // delete Truth_Particle;
690  fFastJetWrapper_Truth->Run();
691  Inclusive_Jets_Truth = fFastJetWrapper_Truth->GetInclusiveJets();
692  if (NTruthD==0) fhEvent->Fill(3);
693  if (NTruthD==1) fhEvent->Fill(4);
694  if (NTruthD==2) fhEvent->Fill(5);
695  if (NTruthD==3) fhEvent->Fill(6);
696  if (NTruthD==4) fhEvent->Fill(7);
697  if (NTruthD==5) fhEvent->Fill(8);
698  if (NTruthD==6) fhEvent->Fill(9);
699  }
700 
701 
702 
703 
704 
705  fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject("D0toKpi"));
706  AliHFTrackContainer* Track_Container=(AliHFTrackContainer *) GetTrackContainer(0);
707  if (!Track_Container) return kTRUE;
708  Track_Container->SetDMesonCandidate(NULL);
709 
710  for (Int_t i_D = 0; i_D < fCandidateArray->GetEntriesFast(); i_D++) {
711  AliAODRecoDecayHF2Prong* D_Candidate = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(i_D));
712  if (!D_Candidate) continue;
713  if (!fRDHFCuts->IsInFiducialAcceptance(D_Candidate->Pt(), D_Candidate->Y(fCandidatePDG))) continue;
714 
715  Int_t Mass_Hypo_Type=fRDHFCuts->IsSelected(D_Candidate, AliRDHFCuts::kAll, fAodEvent);
716  if (Mass_Hypo_Type <= 0 || Mass_Hypo_Type>3) continue;
717  fhEvent->Fill(10);
718 
719  Int_t Matched_Truth_Particle_PDG=0;
720 
721 
722  const Int_t D_Candidtae_N_Daughters=2;
723  Int_t D_Candidtae_Daughters_PDG[D_Candidtae_N_Daughters] = {211,321};
724  Int_t D_Candidate_MatchedTruth_Label = D_Candidate->MatchToMC(fCandidatePDG, Particle_Container->GetArray(), D_Candidtae_N_Daughters, D_Candidtae_Daughters_PDG);
725  Bool_t Is_Prompt_Correct_Quark=kFALSE;
726  Int_t Is_Prompt_Correct_Quark_PDG=-1;
727  AliAODMCParticle* Matched_Truth_Particle;
728  if (D_Candidate_MatchedTruth_Label >= 0) {
729  Matched_Truth_Particle = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(D_Candidate_MatchedTruth_Label));
730  if (Matched_Truth_Particle) {
731  fhEvent->Fill(11);
732 
733  Int_t Matched_Truth_Particle_Mother_Label=Matched_Truth_Particle->GetMother();
734  while (Matched_Truth_Particle_Mother_Label >= 0) {
735  AliAODMCParticle* Matched_Truth_Particle_Mother = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Matched_Truth_Particle_Mother_Label));
736  if (Matched_Truth_Particle_Mother){
737  Int_t Original_Quark_PDG=4;
738  if (fIsBDecay) Original_Quark_PDG=5;
739  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode()) == 4){
740  fhEvent->Fill(12);
741  Is_Prompt_Correct_Quark_PDG=4;
742  }
743  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode()) == 5){
744  fhEvent->Fill(13);
745  Is_Prompt_Correct_Quark_PDG=5;
746  }
747  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode())==Original_Quark_PDG) Is_Prompt_Correct_Quark=kTRUE;
748  if (Matched_Truth_Particle_Mother_Label==Matched_Truth_Particle_Mother->GetMother()) break;
749  Matched_Truth_Particle_Mother_Label=Matched_Truth_Particle_Mother->GetMother();
750  }
751  else break;
752  // delete Matched_Truth_Particle_Mother;
753  }
754  Matched_Truth_Particle_PDG = Matched_Truth_Particle->PdgCode();
755  }
756  }
757  else continue;
758  if (fPromptReject && !Is_Prompt_Correct_Quark) continue;
759 
760 
761  if (Mass_Hypo_Type==1 && Matched_Truth_Particle_PDG!=fCandidatePDG){
762  fhEvent->Fill(14);
763  continue;
764  }
765  else fhEvent->Fill(15);
766  if (Mass_Hypo_Type==2 && Matched_Truth_Particle_PDG!=-fCandidatePDG){
767  fhEvent->Fill(16);
768  continue;
769  }
770  else fhEvent->Fill(17);
771  if (Mass_Hypo_Type==3 && TMath::Abs(Matched_Truth_Particle_PDG)!=fCandidatePDG){
772  fhEvent->Fill(18);
773  continue;
774  }
775  else{
776  fhEvent->Fill(19);
777  if (Matched_Truth_Particle_PDG==fCandidatePDG) fhEvent->Fill(20);
778  if (Matched_Truth_Particle_PDG==-fCandidatePDG) fhEvent->Fill(21);
779  }
780 
781 
782  Double_t Inv_Mass_D=0.0;
783  Double_t Inv_Mass_D_Truth=0.0;
784 
785  if (Mass_Hypo_Type==1 || (Mass_Hypo_Type==3 && Matched_Truth_Particle_PDG==fCandidatePDG)){
786  Inv_Mass_D=D_Candidate->InvMassD0();
787  //Inv_Mass_D_Truth=Matched_Truth_Particle->InvMassD0();
788  Inv_Mass_D_Truth=0.0;
789  }
790  if (Mass_Hypo_Type==2 || (Mass_Hypo_Type==3 && Matched_Truth_Particle_PDG==-fCandidatePDG)){
791  Inv_Mass_D=D_Candidate->InvMassD0bar();
792  //Inv_Mass_D_Truth=Matched_Truth_Particle->InvMassD0bar();
793  Inv_Mass_D_Truth=0.0;
794  }
795 
796 
797 
798 
799  fFastJetWrapper->Clear();
800  AliTLorentzVector D_Candidate_LorentzVector(0,0,0,0);
801  D_Candidate_LorentzVector.SetPtEtaPhiM(D_Candidate->Pt(), D_Candidate->Eta(), D_Candidate->Phi(), Inv_Mass_D);
802  fFastJetWrapper->AddInputVector(D_Candidate_LorentzVector.Px(), D_Candidate_LorentzVector.Py(), D_Candidate_LorentzVector.Pz(), D_Candidate_LorentzVector.E(), 0);
803 
804 
805  if (!Track_Container) continue;
806  Track_Container->SetDMesonCandidate(D_Candidate);
807  AliAODTrack *Track = NULL;
808  for (Int_t i_Track=0; i_Track<Track_Container->GetNTracks(); i_Track++){
809  Track = static_cast<AliAODTrack*>(Track_Container->GetAcceptParticle(i_Track));
810  if(!Track) continue;
811  fFastJetWrapper->AddInputVector(Track->Px(), Track->Py(), Track->Pz(), Track->E(),i_Track+100);
812  }
813  // delete Track;
814  fFastJetWrapper->Run();
815 
816 
817  std::vector<fastjet::PseudoJet> Inclusive_Jets = fFastJetWrapper->GetInclusiveJets();
818  for (UInt_t i_Jet=0; i_Jet < Inclusive_Jets.size(); i_Jet++){
819  Bool_t Is_D_Jet=kFALSE;
820  if (Inclusive_Jets[i_Jet].perp()<fJetMinPt) continue;
821  if (TMath::Abs(Inclusive_Jets[i_Jet].pseudorapidity()) > 0.9-fJetRadius) continue;
822  std::vector<fastjet::PseudoJet> Constituents(fFastJetWrapper->GetJetConstituents(i_Jet));
823  for (UInt_t i_Constituents = 0; i_Constituents < Constituents.size(); i_Constituents++) {
824  if (Constituents[i_Constituents].user_index() == 0) {
825  Is_D_Jet = kTRUE;
826  }
827  }
828 
829 
830  if (!Is_D_Jet) continue;
831  fhEvent->Fill(22);
832 
833 
834  Int_t i_Matched_D_Jet_Truth=-1;
835  for (UInt_t k=0; k< Inclusive_Jets_Truth_Labels.size(); k++){
836  if(Inclusive_Jets_Truth_Labels[k].first==D_Candidate_MatchedTruth_Label) i_Matched_D_Jet_Truth=Inclusive_Jets_Truth_Labels[k].second;
837  }
838 
839  for (UInt_t i_Jet_Truth=0; i_Jet_Truth < Inclusive_Jets_Truth.size(); i_Jet_Truth++){
840  Bool_t Is_Jet_Truth_Matched=kFALSE;
841  if (TMath::Abs(Inclusive_Jets_Truth[i_Jet_Truth].pseudorapidity()) > 0.9-fJetRadius) continue;
842  std::vector<fastjet::PseudoJet> Constituents_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
843  for (UInt_t i_Constituents_Truth = 0; i_Constituents_Truth < Constituents_Truth.size(); i_Constituents_Truth++) {
844  if (Constituents_Truth[i_Constituents_Truth].user_index() == i_Matched_D_Jet_Truth) {
845  Is_Jet_Truth_Matched=kTRUE;
846  for (UInt_t i_Unmacthed_D=0; i_Unmacthed_D<Unmatched_Truth_Level_D.size(); i_Unmacthed_D++){
847  if (Unmatched_Truth_Level_D[i_Unmacthed_D]==i_Matched_D_Jet_Truth) Unmatched_Truth_Level_D.erase(Unmatched_Truth_Level_D.begin()+i_Unmacthed_D);
848  }
849  }
850  }
851  if (!Is_Jet_Truth_Matched) continue;
852  fhEvent->Fill(23);
853  NMatched_DMeson_Jets++;
854 
855 
856 
857 
858  std::vector<Double_t> Splittings_Zg;
859  std::vector<Double_t> Splittings_DeltaR;
860  std::vector<Double_t> Splittings_LeadingSubJetpT;
861  std::vector<Double_t> Splittings_HardestSubJetD0;
862  std::vector<Double_t> Splittings_RadiatorE;
863  std::vector<Double_t> Splittings_RadiatorpT;
864 
865  Bool_t Is_D_SubJet=kFALSE;
866  fastjet::JetDefinition Jet_Definition(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
867 
868  try{
869  std::vector<fastjet::PseudoJet> Reclustered_Particles(fFastJetWrapper->GetJetConstituents(i_Jet));
870  fastjet::ClusterSequence Cluster_Sequence_CA(Reclustered_Particles, Jet_Definition);
871  std::vector<fastjet::PseudoJet> Reclustered_Jet = Cluster_Sequence_CA.inclusive_jets(0.0);
872  Reclustered_Jet = sorted_by_pt(Reclustered_Jet);
873 
874 
875  fastjet::PseudoJet Daughter_Jet = Reclustered_Jet[0];
876  fastjet::PseudoJet Parent_SubJet_1;
877  fastjet::PseudoJet Parent_SubJet_2;
878 
879  while(Daughter_Jet.has_parents(Parent_SubJet_1,Parent_SubJet_2)){
880  if(Parent_SubJet_1.perp() < Parent_SubJet_2.perp()) std::swap(Parent_SubJet_1,Parent_SubJet_2);
881  Splittings_LeadingSubJetpT.push_back(Parent_SubJet_1.perp());
882  vector < fastjet::PseudoJet > Hard_SubJet_Constituents = sorted_by_pt(Parent_SubJet_1.constituents());
883  Is_D_SubJet=kFALSE;
884  for(UInt_t j=0;j<Hard_SubJet_Constituents.size();j++){
885  if(Hard_SubJet_Constituents[j].user_index()==0) Is_D_SubJet=kTRUE;
886  }
887 
888 
889  if (!Is_D_SubJet) Splittings_HardestSubJetD0.push_back(1.0);
890  else Splittings_HardestSubJetD0.push_back(2.0);
891 
892  // if(!Is_D_SubJet) std::swap(Parent_SubJet_1,Parent_SubJet_2);
893  Splittings_DeltaR.push_back(Parent_SubJet_1.delta_R(Parent_SubJet_2));
894  Splittings_Zg.push_back(Parent_SubJet_2.perp()/(Parent_SubJet_1.perp()+Parent_SubJet_2.perp()));
895  Splittings_RadiatorE.push_back(Daughter_Jet.E());
896  Splittings_RadiatorpT.push_back(Daughter_Jet.perp());
897  Daughter_Jet=Parent_SubJet_1;
898  }
899 
900 
901  } catch (fastjet::Error) { /*return -1;*/ }
902 
903 
904  std::vector<Double_t> Splittings_Zg_Truth;
905  std::vector<Double_t> Splittings_DeltaR_Truth;
906  std::vector<Double_t> Splittings_LeadingSubJetpT_Truth;
907  std::vector<Double_t> Splittings_HardestSubJetD0_Truth;
908  std::vector<Double_t> Splittings_RadiatorE_Truth;
909  std::vector<Double_t> Splittings_RadiatorpT_Truth;
910 
911 
912  Bool_t Is_D_SubJet_Truth=kFALSE;
913  fastjet::JetDefinition Jet_Definition_Truth(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
914 
915  try{
916  std::vector<fastjet::PseudoJet> Reclustered_Particles_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
917  fastjet::ClusterSequence Cluster_Sequence_CA_Truth(Reclustered_Particles_Truth, Jet_Definition_Truth);
918  std::vector<fastjet::PseudoJet> Reclustered_Jet_Truth = Cluster_Sequence_CA_Truth.inclusive_jets(0.0);
919  Reclustered_Jet_Truth = sorted_by_pt(Reclustered_Jet_Truth);
920 
921 
922  fastjet::PseudoJet Daughter_Jet_Truth = Reclustered_Jet_Truth[0];
923  fastjet::PseudoJet Parent_SubJet_1_Truth;
924  fastjet::PseudoJet Parent_SubJet_2_Truth;
925 
926  while(Daughter_Jet_Truth.has_parents(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth)){
927  if(Parent_SubJet_1_Truth.perp() < Parent_SubJet_2_Truth.perp()) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
928  Splittings_LeadingSubJetpT_Truth.push_back(Parent_SubJet_1_Truth.perp());
929  vector < fastjet::PseudoJet > Hard_SubJet_Constituents_Truth = sorted_by_pt(Parent_SubJet_1_Truth.constituents());
930  Is_D_SubJet_Truth=kFALSE;
931  for(UInt_t j=0;j<Hard_SubJet_Constituents_Truth.size();j++){
932  if(Hard_SubJet_Constituents_Truth[j].user_index()==i_Matched_D_Jet_Truth) Is_D_SubJet_Truth=kTRUE;
933  }
934 
935 
936  if (!Is_D_SubJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(1.0);
937  else Splittings_HardestSubJetD0_Truth.push_back(2.0);
938 
939  // if(!Is_D_SubJet_Truth) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
940  Splittings_DeltaR_Truth.push_back(Parent_SubJet_1_Truth.delta_R(Parent_SubJet_2_Truth));
941  Splittings_Zg_Truth.push_back(Parent_SubJet_2_Truth.perp()/(Parent_SubJet_1_Truth.perp()+Parent_SubJet_2_Truth.perp()));
942  Splittings_RadiatorE_Truth.push_back(Daughter_Jet_Truth.E());
943  Splittings_RadiatorpT_Truth.push_back(Daughter_Jet_Truth.perp());
944  Daughter_Jet_Truth=Parent_SubJet_1_Truth;
945  }
946 
947 
948  } catch (fastjet::Error) { /*return -1;*/ }
949 
950 
951 
952 
953  //set detector level flags
954  Double_t Flag_D=-1.0;
955  if (Mass_Hypo_Type ==1) Flag_D=1.0;
956  else if (Mass_Hypo_Type ==2) Flag_D=2.0;
957  else if (Mass_Hypo_Type ==3 && Matched_Truth_Particle->GetPdgCode()==fCandidatePDG) Flag_D=3.0;
958  else if (Mass_Hypo_Type ==3 && Matched_Truth_Particle->GetPdgCode()==-fCandidatePDG) Flag_D=4.0;
959 
960  Double_t Flag_D_Truth=-1.0;
961  if (Matched_Truth_Particle->GetPdgCode()==fCandidatePDG) Flag_D_Truth=1.0;
962  if (Matched_Truth_Particle->GetPdgCode()==-fCandidatePDG) Flag_D_Truth=2.0;
963 
964 
965  fShapesVar[0] = Inclusive_Jets[i_Jet].perp();
966  fShapesVar[1] = Inclusive_Jets_Truth[i_Jet_Truth].perp();
967  fShapesVar[2] = D_Candidate->Pt();
968  fShapesVar[3] = Matched_Truth_Particle->Pt();
969  fShapesVar[4] = Inv_Mass_D;
970  fShapesVar[5] = Inv_Mass_D_Truth;
971  fShapesVar[6] = Flag_D;
972  fShapesVar[7] = Flag_D_Truth;
973  fShapesVar[8] = Is_Prompt_Correct_Quark_PDG;
974  fShapesVar[9] = 0.0;
975 
976 
977 
978  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR);
979  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR_Truth);
980  fShapesVar_Splittings_Zg.push_back(Splittings_Zg);
981  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg_Truth);
982  fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT);
983  fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT_Truth);
984  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0);
985  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0_Truth);
986  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE);
987  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE_Truth);
988  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT);
989  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT_Truth);
990 
991 
992  fTreeResponseMatrixAxis->Fill();
993  fTreeSplittings->Fill();
994 
995  }
996  }
997  // delete Matched_Truth_Particle;
998  // delete D_Candidate;
999  }
1000  // delete Track_Container;
1001  if(NMatched_DMeson_Jets==0) fhEvent->Fill(24);
1002  if(NMatched_DMeson_Jets==1) fhEvent->Fill(25);
1003  if(NMatched_DMeson_Jets==2) fhEvent->Fill(26);
1004  if(NMatched_DMeson_Jets==3) fhEvent->Fill(27);
1005  if(NMatched_DMeson_Jets==4) fhEvent->Fill(28);
1006  if(NMatched_DMeson_Jets==5) fhEvent->Fill(29);
1007  if(NMatched_DMeson_Jets==6) fhEvent->Fill(30);
1008 
1009  if (fIncludeInclusive){
1010 
1011  if(Unmatched_Truth_Level_D.size()==0) fhEvent->Fill(31);
1012  if(Unmatched_Truth_Level_D.size()==1) fhEvent->Fill(32);
1013  if(Unmatched_Truth_Level_D.size()==2) fhEvent->Fill(33);
1014  if(Unmatched_Truth_Level_D.size()==3) fhEvent->Fill(34);
1015  if(Unmatched_Truth_Level_D.size()==4) fhEvent->Fill(35);
1016  if(Unmatched_Truth_Level_D.size()==5) fhEvent->Fill(36);
1017  if(Unmatched_Truth_Level_D.size()==6) fhEvent->Fill(37);
1018 
1019  for (UInt_t i_Jet_Truth=0; i_Jet_Truth < Inclusive_Jets_Truth.size(); i_Jet_Truth++){
1020  AliAODMCParticle* Truth_D_Particle = NULL;
1021  Bool_t Is_Unmatched_D=kFALSE;
1022  Int_t D_Meson_Matched_Index=-1;
1023  if (TMath::Abs(Inclusive_Jets_Truth[i_Jet_Truth].pseudorapidity()) > 0.9-fJetRadius) continue;
1024  std::vector<fastjet::PseudoJet> Constituents_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1025  for (UInt_t i_Constituents_Truth = 0; i_Constituents_Truth < Constituents_Truth.size(); i_Constituents_Truth++) {
1026  for (UInt_t i_Unmacthed_D=0; i_Unmacthed_D<Unmatched_Truth_Level_D.size(); i_Unmacthed_D++){
1027  if(Constituents_Truth[i_Constituents_Truth].user_index() == Unmatched_Truth_Level_D[i_Unmacthed_D]){
1028  Is_Unmatched_D=kTRUE;
1029  D_Meson_Matched_Index=Constituents_Truth[i_Constituents_Truth].user_index();
1030  for(UInt_t i_MC_Label=0; i_MC_Label<Inclusive_Jets_Truth_Labels.size(); i_MC_Label++){
1031  if (Inclusive_Jets_Truth_Labels[i_MC_Label].second==Constituents_Truth[i_Constituents_Truth].user_index()){
1032  Truth_D_Particle=static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Inclusive_Jets_Truth_Labels[i_MC_Label].first));
1033  }
1034  }
1035  }
1036  }
1037  }
1038  if (!Is_Unmatched_D) continue;
1039  fhEvent->Fill(38);
1040 
1041 
1042  Bool_t Is_Prompt_Correct_Quark=kFALSE;
1043  Int_t Is_Prompt_Correct_Quark_PDG=-1;
1044  Int_t Truth_D_Particle_Mother_Label=Truth_D_Particle->GetMother();
1045  while (Truth_D_Particle_Mother_Label >= 0) {
1046  AliAODMCParticle* Truth_D_Particle_Mother = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Truth_D_Particle_Mother_Label));
1047  if (Truth_D_Particle_Mother){
1048  Int_t Original_Quark_PDG=4;
1049  if (fIsBDecay) Original_Quark_PDG=5;
1050  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) == 4){
1051  fhEvent->Fill(39);
1052  Is_Prompt_Correct_Quark_PDG=4;
1053  }
1054  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) == 5){
1055  fhEvent->Fill(40);
1056  Is_Prompt_Correct_Quark_PDG=5;
1057  }
1058  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode())==Original_Quark_PDG) Is_Prompt_Correct_Quark=kTRUE;
1059  if (Truth_D_Particle_Mother_Label==Truth_D_Particle_Mother->GetMother()) break;
1060  Truth_D_Particle_Mother_Label=Truth_D_Particle_Mother->GetMother();
1061  }
1062  else break;
1063  // delete Truth_D_Particle_Mother;
1064  }
1065 
1066  if (fPromptReject && !Is_Prompt_Correct_Quark) continue;
1067 
1068 
1069  Bool_t Is_D_SubJet_Truth=kFALSE;
1070  fastjet::JetDefinition Jet_Definition_Truth(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
1071 
1072  std::vector<Double_t> Splittings_Zg_Truth;
1073  std::vector<Double_t> Splittings_DeltaR_Truth;
1074  std::vector<Double_t> Splittings_LeadingSubJetpT_Truth;
1075  std::vector<Double_t> Splittings_HardestSubJetD0_Truth;
1076  std::vector<Double_t> Splittings_RadiatorE_Truth;
1077  std::vector<Double_t> Splittings_RadiatorpT_Truth;
1078 
1079 
1080  try{
1081  std::vector<fastjet::PseudoJet> Reclustered_Particles_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1082  fastjet::ClusterSequence Cluster_Sequence_CA_Truth(Reclustered_Particles_Truth, Jet_Definition_Truth);
1083  std::vector<fastjet::PseudoJet> Reclustered_Jet_Truth = Cluster_Sequence_CA_Truth.inclusive_jets(0.0);
1084  Reclustered_Jet_Truth = sorted_by_pt(Reclustered_Jet_Truth);
1085 
1086 
1087  fastjet::PseudoJet Daughter_Jet_Truth = Reclustered_Jet_Truth[0];
1088  fastjet::PseudoJet Parent_SubJet_1_Truth;
1089  fastjet::PseudoJet Parent_SubJet_2_Truth;
1090 
1091  while(Daughter_Jet_Truth.has_parents(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth)){
1092  if(Parent_SubJet_1_Truth.perp() < Parent_SubJet_2_Truth.perp()) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1093  Splittings_LeadingSubJetpT_Truth.push_back(Parent_SubJet_1_Truth.perp());
1094  vector < fastjet::PseudoJet > Hard_SubJet_Constituents_Truth = sorted_by_pt(Parent_SubJet_1_Truth.constituents());
1095  Is_D_SubJet_Truth=kFALSE;
1096  for(UInt_t j=0;j<Hard_SubJet_Constituents_Truth.size();j++){
1097  if(Hard_SubJet_Constituents_Truth[j].user_index()==D_Meson_Matched_Index) Is_D_SubJet_Truth=kTRUE;
1098  }
1099 
1100 
1101  if (!Is_D_SubJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(1.0);
1102  else Splittings_HardestSubJetD0_Truth.push_back(2.0);
1103 
1104  // if(!Is_D_SubJet_Truth) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1105  Splittings_DeltaR_Truth.push_back(Parent_SubJet_1_Truth.delta_R(Parent_SubJet_2_Truth));
1106  Splittings_Zg_Truth.push_back(Parent_SubJet_2_Truth.perp()/(Parent_SubJet_1_Truth.perp()+Parent_SubJet_2_Truth.perp()));
1107  Splittings_RadiatorE_Truth.push_back(Daughter_Jet_Truth.E());
1108  Splittings_RadiatorpT_Truth.push_back(Daughter_Jet_Truth.perp());
1109  Daughter_Jet_Truth=Parent_SubJet_1_Truth;
1110  }
1111 
1112 
1113  } catch (fastjet::Error) { /*return -1;*/ }
1114 
1115  Double_t Inv_Mass_D_Truth=0.0;
1116  Double_t Flag_D_Truth=-1.0;
1117  if (Truth_D_Particle->GetPdgCode()==fCandidatePDG){
1118  //Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0();
1119  Inv_Mass_D_Truth=0.0;
1120  Flag_D_Truth=3.0;
1121  }
1122  if (Truth_D_Particle->GetPdgCode()==-fCandidatePDG){
1123  // Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0bar();
1124  Inv_Mass_D_Truth=0.0;
1125  Flag_D_Truth=4.0;
1126  }
1127 
1128 
1129  fShapesVar[0] = 0.0;
1130  fShapesVar[1] = Inclusive_Jets_Truth[i_Jet_Truth].perp();
1131  fShapesVar[2] = 0.0;
1132  fShapesVar[3] = Truth_D_Particle->Pt();
1133  fShapesVar[4] = 0.0;
1134  fShapesVar[5] = Inv_Mass_D_Truth;
1135  fShapesVar[6] = 0.0;
1136  fShapesVar[7] = Flag_D_Truth;
1137  fShapesVar[8] = Is_Prompt_Correct_Quark_PDG;
1138  fShapesVar[9] = 0.0;
1139 
1140 
1141  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR_Truth);
1142  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR_Truth);
1143  fShapesVar_Splittings_Zg.push_back(Splittings_Zg_Truth);
1144  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg_Truth);
1145  fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT_Truth);
1146  fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT_Truth);
1147  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0_Truth);
1148  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0_Truth);
1149  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE_Truth);
1150  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE_Truth);
1151  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT_Truth);
1152  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT_Truth);
1153 
1154 
1155  fTreeResponseMatrixAxis->Fill();
1156  fTreeSplittings->Fill();
1157  // delete Truth_D_Particle;
1158  }
1159  }
1160  // delete Particle_Container;
1161  }
1162 
1163 
1164  if (fJetShapeType == kTrue){
1165 
1166 
1167  //truth level only jet finding
1168 
1169 
1170  AliHFAODMCParticleContainer *Particle_Container = (AliHFAODMCParticleContainer*) GetParticleContainer(0);
1171  Particle_Container->SetSpecialPDG(fCandidatePDG);
1172  Particle_Container->SetRejectedOriginMap(fRejectedOrigin);
1173  Particle_Container->SetAcceptedDecayMap(kDecayD0toKpi);
1174  Particle_Container->SetRejectISR(fRejectISR);
1175 
1176  std::vector<fastjet::PseudoJet> Inclusive_Jets_Truth;
1177  std::vector<std::pair<Int_t, Int_t>> Inclusive_Jets_Truth_Labels;
1178 
1179  if (Particle_Container->IsSpecialPDGFound()){
1180  fhEvent->Fill(2);
1181 
1182  fFastJetWrapper_Truth->SetAreaType(fastjet::active_area);
1183  fFastJetWrapper_Truth->SetGhostArea(0.005);
1184  fFastJetWrapper_Truth->SetR(fJetRadius);
1185  fFastJetWrapper_Truth->SetAlgorithm(fastjet::antikt_algorithm);
1186  fFastJetWrapper_Truth->SetRecombScheme(static_cast<fastjet::RecombinationScheme>(0));
1187  fFastJetWrapper_Truth->Clear();
1188 
1189  AliAODMCParticle* Truth_Particle=NULL;
1190  Int_t NTruthD=0;
1191  for (Int_t i_Particle=0; i_Particle<Particle_Container->GetNParticles(); i_Particle++){
1192  Truth_Particle = static_cast<AliAODMCParticle*>(Particle_Container->GetAcceptParticle(i_Particle));
1193  if (!Truth_Particle) continue;
1194  if (TMath::Abs(Truth_Particle->GetPdgCode())==fCandidatePDG){
1195  fhEvent->Fill(3);
1196  std::pair<Int_t, Int_t> Inclusive_Jet_Truth_Labels;
1197  Inclusive_Jet_Truth_Labels.first=Truth_Particle->GetLabel();
1198  Inclusive_Jet_Truth_Labels.second=NTruthD;
1199  Inclusive_Jets_Truth_Labels.push_back(Inclusive_Jet_Truth_Labels);
1200  fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),NTruthD);
1201  NTruthD++;
1202  }
1203  else fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),i_Particle+100);
1204  }
1205  // delete Truth_Particle;
1206  fFastJetWrapper_Truth->Run();
1207  Inclusive_Jets_Truth = fFastJetWrapper_Truth->GetInclusiveJets();
1208 
1209  for (UInt_t i_Jet_Truth=0; i_Jet_Truth < Inclusive_Jets_Truth.size(); i_Jet_Truth++){
1210  AliAODMCParticle* Truth_D_Particle = NULL;
1211  Bool_t Is_DJet_Truth=kFALSE;
1212  Int_t D_Meson_Matched_Index=-1;
1213  if (TMath::Abs(Inclusive_Jets_Truth[i_Jet_Truth].pseudorapidity()) > 0.9-fJetRadius) continue;
1214  std::vector<fastjet::PseudoJet> Constituents_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1215  for (UInt_t i_Constituents_Truth = 0; i_Constituents_Truth < Constituents_Truth.size(); i_Constituents_Truth++) {
1216 
1217  if(Constituents_Truth[i_Constituents_Truth].user_index() >=0 && Constituents_Truth[i_Constituents_Truth].user_index() < NTruthD){
1218  D_Meson_Matched_Index=Constituents_Truth[i_Constituents_Truth].user_index();
1219  Is_DJet_Truth=kTRUE;
1220  for(UInt_t i_MC_Label=0; i_MC_Label<Inclusive_Jets_Truth_Labels.size(); i_MC_Label++){
1221  if (Inclusive_Jets_Truth_Labels[i_MC_Label].second==Constituents_Truth[i_Constituents_Truth].user_index()){
1222  Truth_D_Particle=static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Inclusive_Jets_Truth_Labels[i_MC_Label].first));
1223  }
1224  }
1225 
1226  }
1227  }
1228 
1229  if (!Is_DJet_Truth) fhEvent->Fill(4);
1230  if (Is_DJet_Truth) fhEvent->Fill(5);
1231 
1232  if (!fIncludeInclusive && !Is_DJet_Truth) continue;
1233 
1234  Bool_t Is_Prompt_Correct_Quark=kFALSE;
1235  Int_t Is_Prompt_Correct_Quark_PDG=-1;
1236  Int_t Truth_D_Particle_Mother_Label=Truth_D_Particle->GetMother();
1237  while (Truth_D_Particle_Mother_Label >= 0) {
1238  AliAODMCParticle* Truth_D_Particle_Mother = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Truth_D_Particle_Mother_Label));
1239  if (Truth_D_Particle_Mother){
1240  Int_t Original_Quark_PDG=4;
1241  if (fIsBDecay) Original_Quark_PDG=5;
1242  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) ==Original_Quark_PDG) Is_Prompt_Correct_Quark=kTRUE;
1243  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) ==4){
1244  fhEvent->Fill(6);
1245  Is_Prompt_Correct_Quark_PDG=4;
1246  }
1247  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) ==5){
1248  fhEvent->Fill(7);
1249  Is_Prompt_Correct_Quark_PDG=5;
1250  }
1251  if (Truth_D_Particle_Mother_Label==Truth_D_Particle_Mother->GetMother()) break;
1252  Truth_D_Particle_Mother_Label=Truth_D_Particle_Mother->GetMother();
1253  }
1254  else break;
1255  // delete Truth_D_Particle_Mother;
1256  }
1257 
1258  if(fPromptReject && !Is_Prompt_Correct_Quark) continue;
1259  fhEvent->Fill(8);
1260 
1261  Bool_t Is_D_SubJet_Truth=kFALSE;
1262  fastjet::JetDefinition Jet_Definition_Truth(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
1263 
1264  std::vector<Double_t> Splittings_Zg_Truth;
1265  std::vector<Double_t> Splittings_DeltaR_Truth;
1266  std::vector<Double_t> Splittings_LeadingSubJetpT_Truth;
1267  std::vector<Double_t> Splittings_HardestSubJetD0_Truth;
1268  std::vector<Double_t> Splittings_RadiatorE_Truth;
1269  std::vector<Double_t> Splittings_RadiatorpT_Truth;
1270 
1271  try{
1272  std::vector<fastjet::PseudoJet> Reclustered_Particles_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1273  fastjet::ClusterSequence Cluster_Sequence_CA_Truth(Reclustered_Particles_Truth, Jet_Definition_Truth);
1274  std::vector<fastjet::PseudoJet> Reclustered_Jet_Truth = Cluster_Sequence_CA_Truth.inclusive_jets(0.0);
1275  Reclustered_Jet_Truth = sorted_by_pt(Reclustered_Jet_Truth);
1276 
1277 
1278  fastjet::PseudoJet Daughter_Jet_Truth = Reclustered_Jet_Truth[0];
1279  fastjet::PseudoJet Parent_SubJet_1_Truth;
1280  fastjet::PseudoJet Parent_SubJet_2_Truth;
1281 
1282  while(Daughter_Jet_Truth.has_parents(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth)){
1283  if(Parent_SubJet_1_Truth.perp() < Parent_SubJet_2_Truth.perp()) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1284  Splittings_LeadingSubJetpT_Truth.push_back(Parent_SubJet_1_Truth.perp());
1285  vector < fastjet::PseudoJet > Hard_SubJet_Constituents_Truth = sorted_by_pt(Parent_SubJet_1_Truth.constituents());
1286  Is_D_SubJet_Truth=kFALSE;
1287  for(UInt_t j=0;j<Hard_SubJet_Constituents_Truth.size();j++){
1288  if(Hard_SubJet_Constituents_Truth[j].user_index()==D_Meson_Matched_Index) Is_D_SubJet_Truth=kTRUE;
1289  }
1290 
1291  if (fIncludeInclusive && !Is_DJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(0.0);
1292  else if (Is_DJet_Truth && !Is_D_SubJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(1.0);
1293  else if (Is_DJet_Truth && Is_D_SubJet_Truth)Splittings_HardestSubJetD0_Truth.push_back(2.0);
1294 
1295  // if(!Is_D_SubJet_Truth) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1296  Splittings_DeltaR_Truth.push_back(Parent_SubJet_1_Truth.delta_R(Parent_SubJet_2_Truth));
1297  Splittings_Zg_Truth.push_back(Parent_SubJet_2_Truth.perp()/(Parent_SubJet_1_Truth.perp()+Parent_SubJet_2_Truth.perp()));
1298  Splittings_RadiatorE_Truth.push_back(Daughter_Jet_Truth.E());
1299  Splittings_RadiatorpT_Truth.push_back(Daughter_Jet_Truth.perp());
1300  Daughter_Jet_Truth=Parent_SubJet_1_Truth;
1301  }
1302 
1303 
1304  } catch (fastjet::Error) { /*return -1;*/ }
1305 
1306  Double_t Inv_Mass_D_Truth=0.0;
1307  Double_t Flag_D_Truth=-1.0;
1308  Double_t D_Pt=-1.0;
1309  if (fIncludeInclusive && !Is_DJet_Truth){
1310  Flag_D_Truth=0.0;
1311  }
1312  else if (Is_DJet_Truth){
1313  if (Truth_D_Particle->GetPdgCode()==fCandidatePDG){
1314  //Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0();
1315  Inv_Mass_D_Truth=0.0;
1316  Flag_D_Truth=1.0;
1317  D_Pt=Truth_D_Particle->Pt();
1318  fhEvent->Fill(9);
1319  }
1320  if (Truth_D_Particle->GetPdgCode()==-fCandidatePDG){
1321  // Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0bar();
1322  Inv_Mass_D_Truth=0.0;
1323  Flag_D_Truth=2.0;
1324  D_Pt=Truth_D_Particle->Pt();
1325  fhEvent->Fill(10);
1326  }
1327  }
1328 
1329 
1330  fShapesVar[0] = 0.0;
1331  fShapesVar[1] = Inclusive_Jets_Truth[i_Jet_Truth].perp();
1332  fShapesVar[2] = 0.0;
1333  fShapesVar[3] = D_Pt;
1334  fShapesVar[4] = 0.0;
1335  fShapesVar[5] = Inv_Mass_D_Truth;
1336  fShapesVar[6] = 0.0;
1337  fShapesVar[7] = Flag_D_Truth;
1338  fShapesVar[8] = 0.0;
1339  fShapesVar[9] = Is_Prompt_Correct_Quark_PDG;
1340 
1341 
1342 
1343  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR_Truth);
1344  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR_Truth);
1345  fShapesVar_Splittings_Zg.push_back(Splittings_Zg_Truth);
1346  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg_Truth);
1347  fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT_Truth);
1348  fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT_Truth);
1349  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0_Truth);
1350  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0_Truth);
1351  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE_Truth);
1352  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE_Truth);
1353  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT_Truth);
1354  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT_Truth);
1355 
1356 
1357  fTreeResponseMatrixAxis->Fill();
1358  fTreeSplittings->Fill();
1359 
1360 
1361  // delete Truth_D_Particle;
1362  }
1363  if (NTruthD==0) fhEvent->Fill(11);
1364  if (NTruthD==1) fhEvent->Fill(12);
1365  if (NTruthD==2) fhEvent->Fill(13);
1366  if (NTruthD==3) fhEvent->Fill(14);
1367  if (NTruthD==4) fhEvent->Fill(15);
1368  if (NTruthD==5) fhEvent->Fill(16);
1369  if (NTruthD==6) fhEvent->Fill(17);
1370  }
1371  //delete Particle_Container;
1372  }
1373 
1374  delete fFastJetWrapper_Truth;
1375  delete fFastJetWrapper;
1376 
1377  return kTRUE;
1378 
1379 
1380 }
1381 
1382 
1383 
1384 
1385 
1386 
1387 //________________________________________________________________________
1388 Bool_t AliAnalysisTaskHFSubstructure::RetrieveEventObjects() {
1389  //
1390  // retrieve event objects
1391  //
1392  if (!AliAnalysisTaskEmcal::RetrieveEventObjects())
1393  return kFALSE;
1394 
1395  return kTRUE;
1396 }
1397 
1398 
1399 //_______________________________________________________________________
1400 void AliAnalysisTaskHFSubstructure::Terminate(Option_t *)
1401 {
1402  // Called once at the end of the analysis.
1403 
1404 }
AliAnalysisTaskHFSubstructure::fShapesVar_Splittings_HardestSubJetD0_Truth
std::vector< std::vector< Double_t > > fShapesVar_Splittings_HardestSubJetD0_Truth
Definition: AliAnalysisTaskHFSubstructure.h:221
AliRDHFCutsD0toKpi.h
AliHFTrackContainer::SetDMesonCandidate
void SetDMesonCandidate(AliAODRecoDecay *c)
Definition: AliHFTrackContainer.cxx:92
AliFJWrapper::Run
virtual Int_t Run()
Double_t
double Double_t
Definition: External.C:58
kD0toKpi
Definition: FitMassSpectra.C:38
AliTrackContainer::GetNTracks
Int_t GetNTracks() const
Definition: AliTrackContainer.h:140
AliAnalysisTaskHFSubstructure.h
AliAnalysisTaskHFSubstructure::fShapesVar_Splittings_RadiatorE_Truth
std::vector< std::vector< Double_t > > fShapesVar_Splittings_RadiatorE_Truth
Definition: AliAnalysisTaskHFSubstructure.h:223
AliAODRecoCascadeHF.h
AliAnalysisTaskEmcal
Base task in the EMCAL framework.
Definition: AliAnalysisTaskEmcal.h:94
AliAnalysisVertexingHF.h
AliAnalysisTaskHFSubstructure::fShapesVar_Splittings_HardestSubJetD0
std::vector< std::vector< Double_t > > fShapesVar_Splittings_HardestSubJetD0
Definition: AliAnalysisTaskHFSubstructure.h:220
AliParticleContainer::GetNParticles
Int_t GetNParticles() const
Definition: AliParticleContainer.h:72
AliRDHFCuts::CheckMatchingAODdeltaAODevents
static Int_t CheckMatchingAODdeltaAODevents()
Definition: AliRDHFCuts.cxx:764
AliAnalysisTaskHFSubstructure::fShapesVar_Splittings_LeadingSubJetpT
std::vector< std::vector< Double_t > > fShapesVar_Splittings_LeadingSubJetpT
Definition: AliAnalysisTaskHFSubstructure.h:218
AliAnalysisTaskHFSubstructure::fShapesVar_Splittings_RadiatorpT_Truth
std::vector< std::vector< Double_t > > fShapesVar_Splittings_RadiatorpT_Truth
Definition: AliAnalysisTaskHFSubstructure.h:225
AliAnalysisTaskHFSubstructure::fShapesVar_Splittings_Zg
std::vector< std::vector< Double_t > > fShapesVar_Splittings_Zg
Definition: AliAnalysisTaskHFSubstructure.h:216
AliEmcalParticle.h
AliAnalysisTaskHFSubstructure::FillHistograms
Bool_t FillHistograms()
Function filling histograms.
Definition: AliAnalysisTaskHFSubstructure.cxx:258
AliEmcalJetTask.h
AliFJWrapper::SetRecombScheme
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:122
AliHFTrackContainer.h
AliFJWrapper.h
FJ_includes.h
AliFJWrapper::SetR
void SetR(Double_t r)
Definition: AliFJWrapper.h:127
AliAnalysisTaskHFSubstructure::fShapesVar_Splittings_DeltaR_Truth
std::vector< std::vector< Double_t > > fShapesVar_Splittings_DeltaR_Truth
Definition: AliAnalysisTaskHFSubstructure.h:215
kData
TString kData
Declare data MC or deltaAOD.
Definition: AddTaskCaloTrackCorr.C:63
AliHFTrackContainer
Select tracks based on specific prescriptions of HF analysis.
Definition: AliHFTrackContainer.h:36
AliAnalysisTaskEmcal::GetParticleContainer
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Definition: AliAnalysisTaskEmcal.cxx:1671
AliFJWrapper::GetInclusiveJets
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:34
AliTrackContainer::GetAcceptParticle
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
Definition: AliTrackContainer.h:123
Int_t
int Int_t
Definition: External.C:63
UInt_t
unsigned int UInt_t
Definition: External.C:33
AliMCParticleContainer::GetAcceptParticle
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
Definition: AliMCParticleContainer.h:46
AliAnalysisTaskHFSubstructure::Run
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
Definition: AliAnalysisTaskHFSubstructure.cxx:250
AliParticleContainer.h
AliFJWrapper::SetAlgorithm
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:121
AliAnalysisTaskHFSubstructure::RetrieveEventObjects
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Definition: AliAnalysisTaskHFSubstructure.cxx:1388
TH1D
Definition: External.C:212
AliJetContainer.h
AliRhoParameter.h
AliHFAODMCParticleContainer.h
AliFJWrapper::Clear
virtual void Clear(const Option_t *="")
AliAODRecoDecayHF2Prong.h
AliAnalysisTaskEmcal::RetrieveEventObjects
virtual Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Definition: AliAnalysisTaskEmcal.cxx:1414
Data
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Definition: charmCutsOptimization.C:371
AliHFAODMCParticleContainer
Select MC particles based on specific prescriptions of HF analysis.
Definition: AliHFAODMCParticleContainer.h:34
AliAnalysisTaskHFSubstructure::fShapesVar_Splittings_Zg_Truth
std::vector< std::vector< Double_t > > fShapesVar_Splittings_Zg_Truth
Definition: AliAnalysisTaskHFSubstructure.h:217
AliFJWrapper::SetGhostArea
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:125
AliFJWrapper::GetJetConstituents
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
Track
void Track()
Definition: DrawProductionComparison.C:1494
AliRDHFCuts::IsEventSelected
Bool_t IsEventSelected(AliVEvent *event)
Definition: AliRDHFCuts.cxx:503
AliAnalysisTaskEmcal::fOutput
AliEmcalList * fOutput
!output list
Definition: AliAnalysisTaskEmcal.h:990
AliAnalysisTaskEmcal::GetTrackContainer
AliTrackContainer * GetTrackContainer(Int_t i=0) const
Definition: AliAnalysisTaskEmcal.h:248
AliAnalysisTaskEmcal::SetMakeGeneralHistograms
void SetMakeGeneralHistograms(Bool_t g)
Definition: AliAnalysisTaskEmcal.h:283
AliRDHFCuts::IsSelected
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:293
AliAnalysisTaskHFSubstructure::fShapesVar_Splittings_DeltaR
std::vector< std::vector< Double_t > > fShapesVar_Splittings_DeltaR
Definition: AliAnalysisTaskHFSubstructure.h:214
AliRDHFCutsDStartoKpipi.h
Option_t
const char Option_t
Definition: External.C:48
AliFJWrapper::AddInputVector
virtual void AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
AliAnalysisTaskEmcal::UserCreateOutputObjects
void UserCreateOutputObjects()
Main initialization function on the worker.
Definition: AliAnalysisTaskEmcal.cxx:362
Bool_t
bool Bool_t
Definition: External.C:53
AliFJWrapper::SetAreaType
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:123
AliRDHFCuts::IsInFiducialAcceptance
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:314
AliAnalysisTaskHFSubstructure::fShapesVar_Splittings_LeadingSubJetpT_Truth
std::vector< std::vector< Double_t > > fShapesVar_Splittings_LeadingSubJetpT_Truth
Definition: AliAnalysisTaskHFSubstructure.h:219
swap
void swap(PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &first, PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &second)
Definition: AliAnalysisTaskEmcalJetHPerformance.cxx:756
AliEmcalJet.h
AliAnalysisTaskHFSubstructure::fShapesVar_Splittings_RadiatorpT
std::vector< std::vector< Double_t > > fShapesVar_Splittings_RadiatorpT
Definition: AliAnalysisTaskHFSubstructure.h:224
AliAnalysisTaskHFSubstructure::fShapesVar_Splittings_RadiatorE
std::vector< std::vector< Double_t > > fShapesVar_Splittings_RadiatorE
Definition: AliAnalysisTaskHFSubstructure.h:222