AliPhysics  608b256 (608b256)
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"
45 #include "AliRDHFCutsD0toKpi.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 
61 
62 
63 #include "FJ_includes.h"
64 
65 
66 
67 
68 
69 using std::cout;
70 using std::endl;
71 
72 
74 
75 //________________________________________________________________________
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 //________________________________________________________________________
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  }
168  DefineOutput(1, TList::Class());
169  DefineOutput(2, TTree::Class());
170  DefineOutput(3, TTree::Class());
171 }
172 
173 //________________________________________________________________________
175 {
176  // Destructor.
177 }
178 
179 //________________________________________________________________________
181 {
182  // Create user output.
183 
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 //________________________________________________________________________
251 {
252  // Run analysis code here, if needed. It will be executed before FillHistograms().
253 
254 
255  return kTRUE;
256 }
257 //________________________________________________________________________
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);
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 
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 == kDet || (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 == kDet || (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 == kDet || (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 == kDet || (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 
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  if (fJetShapeType == kData) fShapesVar[7] = 0.0;
506  else fShapesVar[7] = Matched_Truth_Particle_PDG;
507  fShapesVar[8] = Is_Prompt_Correct_Quark_PDG;
508  fShapesVar[9] = 0.0;
509 
510 
511  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR);
512  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR);
513  fShapesVar_Splittings_Zg.push_back(Splittings_Zg);
514  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg);
515  fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT);
516  fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT);
517  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0);
518  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0);
519  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE);
520  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE);
521  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT);
522  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT);
523 
524  fTreeResponseMatrixAxis->Fill();
525  fTreeSplittings->Fill();
528  fShapesVar_Splittings_Zg.clear();
538  }
539  }
540  //delete D_Candidate;
541  }
542  if(N_DMesons==0) fhEvent->Fill(16);
543  if(N_DMesons==1) fhEvent->Fill(17);
544  if(N_DMesons==2) fhEvent->Fill(18);
545  if(N_DMesons==3) fhEvent->Fill(19);
546  if(N_DMesons==4) fhEvent->Fill(20);
547  if(N_DMesons==5) fhEvent->Fill(21);
548  if(N_DMesons==6) fhEvent->Fill(22);
549 
550  if (fIncludeInclusive){
552  /*
553  for (UInt_t i_D_Found=0; i_D_Found<D_Candidates_Vector.size(); i_D_Found++){
554  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);
555  //Track_Container->SetDMesonCandidate(D_Candidates_Vector[i_D_Found]);
556  }
557  */
558  AliAODTrack *Track = NULL;
559  Bool_t DMeson_Daughter_Track=kFALSE;
560  for (Int_t i_Track=0; i_Track<Track_Container->GetNTracks(); i_Track++){
561  /* for (UInt_t i_D_Found=0; i_D_Found<D_Candidates_Vector.size(); i_D_Found++){
562  Track_Container->SetDMesonCandidate(D_Candidates_Vector[i_D_Found]);
563  Track = static_cast<AliAODTrack*>(Track_Container->GetAcceptParticle(i_Track));
564  if(!Track) DMeson_Daughter_Track=kTRUE;
565  }
566  if (DMeson_Daughter_Track) continue;*/
567  Track = static_cast<AliAODTrack*>(Track_Container->GetAcceptParticle(i_Track));
568  if(!Track) continue;
569  fFastJetWrapper->AddInputVector(Track->Px(), Track->Py(), Track->Pz(), Track->E(),i_Track+100);
570  }
571  //delete Track;
572  fFastJetWrapper->Run();
573  std::vector<fastjet::PseudoJet> Inclusive_Jets = fFastJetWrapper->GetInclusiveJets();
574  for (UInt_t i_Jet=0; i_Jet < Inclusive_Jets.size(); i_Jet++){
575  if (Inclusive_Jets[i_Jet].perp()<fJetMinPt) continue;
576  if (TMath::Abs(Inclusive_Jets[i_Jet].pseudorapidity()) > 0.9-fJetRadius) continue;
577  /* Bool_t Is_D_Jet = kFALSE;
578  std::vector<fastjet::PseudoJet> Constituents(fFastJetWrapper->GetJetConstituents(i_Jet));
579  for (UInt_t i_Constituents = 0; i_Constituents < Constituents.size(); i_Constituents++) {
580  for (UInt_t i_D_Found=0; i_D_Found<D_Candidates_Vector.size(); i_D_Found++){
581  if (Constituents[i_Constituents].user_index() == i_D_Found) {
582  Is_D_Jet = kTRUE;
583  }
584  }
585  }
586  if (Is_D_Jet) continue; */
587  fhEvent->Fill(23);
588  std::vector<Double_t> Splittings_Zg;
589  std::vector<Double_t> Splittings_DeltaR;
590  std::vector<Double_t> Splittings_LeadingSubJetpT;
591  std::vector<Double_t> Splittings_HardestSubJetD0;
592  std::vector<Double_t> Splittings_RadiatorE;
593  std::vector<Double_t> Splittings_RadiatorpT;
594 
595  fastjet::JetDefinition Jet_Definition(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
596 
597  try{
598  std::vector<fastjet::PseudoJet> Reclustered_Particles(fFastJetWrapper->GetJetConstituents(i_Jet));
599  fastjet::ClusterSequence Cluster_Sequence_CA(Reclustered_Particles, Jet_Definition);
600  std::vector<fastjet::PseudoJet> Reclustered_Jet = Cluster_Sequence_CA.inclusive_jets(0.0);
601  Reclustered_Jet = sorted_by_pt(Reclustered_Jet);
602 
603 
604  fastjet::PseudoJet Daughter_Jet = Reclustered_Jet[0];
605  fastjet::PseudoJet Parent_SubJet_1;
606  fastjet::PseudoJet Parent_SubJet_2;
607 
608  while(Daughter_Jet.has_parents(Parent_SubJet_1,Parent_SubJet_2)){
609  if(Parent_SubJet_1.perp() < Parent_SubJet_2.perp()) std::swap(Parent_SubJet_1,Parent_SubJet_2);
610  Splittings_LeadingSubJetpT.push_back(Parent_SubJet_1.perp());
611  Splittings_HardestSubJetD0.push_back(0.0);
612  Splittings_DeltaR.push_back(Parent_SubJet_1.delta_R(Parent_SubJet_2));
613  Splittings_Zg.push_back(Parent_SubJet_2.perp()/(Parent_SubJet_1.perp()+Parent_SubJet_2.perp()));
614  Splittings_RadiatorE.push_back(Daughter_Jet.E());
615  Splittings_RadiatorpT.push_back(Daughter_Jet.perp());
616  Daughter_Jet=Parent_SubJet_1;
617  }
618 
619 
620  } catch (fastjet::Error) { /*return -1;*/ }
621 
622 
623 
624  fShapesVar[0] = Inclusive_Jets[i_Jet].perp();
625  fShapesVar[1] = 0.0;
626  fShapesVar[2] = 0.0;
627  fShapesVar[3] = 0.0;
628  fShapesVar[4] = 0.0;
629  fShapesVar[5] = 0.0;
630  fShapesVar[6] = 0.0;
631  fShapesVar[7] = 0.0;
632  fShapesVar[8] = 0.0;
633  fShapesVar[9] = 0.0;
634 
635 
636  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR);
637  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR);
638  fShapesVar_Splittings_Zg.push_back(Splittings_Zg);
639  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg);
640  fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT);
641  fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT);
642  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0);
643  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0);
644  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE);
645  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE);
646  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT);
647  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT);
648 
649 
650  fTreeResponseMatrixAxis->Fill();
651  fTreeSplittings->Fill();
652 
655  fShapesVar_Splittings_Zg.clear();
665 
666  }
667  }
668  // delete Track_Container;
669  // if (fJetShapeType != kData) delete Particle_Container;
670 
671  }
672 
673 
674 
675  if (fJetShapeType == kTrueDet){
676 
677 
679  Particle_Container->SetSpecialPDG(fCandidatePDG);
680  Particle_Container->SetRejectedOriginMap(fRejectedOrigin);
681  Particle_Container->SetAcceptedDecayMap(fAcceptedDecay);
682  Particle_Container->SetRejectISR(fRejectISR);
683 
684  std::vector<fastjet::PseudoJet> Inclusive_Jets_Truth;
685  std::vector<std::pair<Int_t, Int_t>> Inclusive_Jets_Truth_Labels;
686  std::vector<Int_t> Unmatched_Truth_Level_D;
687  Int_t NMatched_DMeson_Jets=0;
688 
689  if (Particle_Container->IsSpecialPDGFound()){
690 
691  fhEvent->Fill(2);
692 
693  fFastJetWrapper_Truth->SetAreaType(fastjet::active_area);
696  fFastJetWrapper_Truth->SetAlgorithm(fastjet::antikt_algorithm);
697  fFastJetWrapper_Truth->SetRecombScheme(static_cast<fastjet::RecombinationScheme>(0));
699 
700  AliAODMCParticle* Truth_Particle=NULL;
701  Int_t NTruthD=0;
702  for (Int_t i_Particle=0; i_Particle<Particle_Container->GetNParticles(); i_Particle++){
703  Truth_Particle = static_cast<AliAODMCParticle*>(Particle_Container->GetAcceptParticle(i_Particle));
704  if (!Truth_Particle) continue;
705  if (TMath::Abs(Truth_Particle->PdgCode())==fCandidatePDG){
706  std::pair<Int_t, Int_t> Inclusive_Jet_Truth_Labels;
707  Inclusive_Jet_Truth_Labels.first=Truth_Particle->GetLabel();
708  Inclusive_Jet_Truth_Labels.second=NTruthD;
709  Inclusive_Jets_Truth_Labels.push_back(Inclusive_Jet_Truth_Labels);
710  Unmatched_Truth_Level_D.push_back(NTruthD);
711  fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),NTruthD);
712  NTruthD++;
713  }
714  else fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),i_Particle+100);
715  }
716  // delete Truth_Particle;
718  Inclusive_Jets_Truth = fFastJetWrapper_Truth->GetInclusiveJets();
719  if (NTruthD==0) fhEvent->Fill(3);
720  if (NTruthD==1) fhEvent->Fill(4);
721  if (NTruthD==2) fhEvent->Fill(5);
722  if (NTruthD==3) fhEvent->Fill(6);
723  if (NTruthD==4) fhEvent->Fill(7);
724  if (NTruthD==5) fhEvent->Fill(8);
725  if (NTruthD==6) fhEvent->Fill(9);
726  }
727 
728 
729 
730 
731 
732  fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject("D0toKpi"));
734  if (!Track_Container) return kTRUE;
735  Track_Container->SetDMesonCandidate(NULL);
736 
737  for (Int_t i_D = 0; i_D < fCandidateArray->GetEntriesFast(); i_D++) {
738  AliAODRecoDecayHF2Prong* D_Candidate = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(i_D));
739  if (!D_Candidate) continue;
740  if (!fRDHFCuts->IsInFiducialAcceptance(D_Candidate->Pt(), D_Candidate->Y(fCandidatePDG))) continue;
741 
742  Int_t Mass_Hypo_Type=fRDHFCuts->IsSelected(D_Candidate, AliRDHFCuts::kAll, fAodEvent);
743  if (Mass_Hypo_Type <= 0 || Mass_Hypo_Type>3) continue;
744  fhEvent->Fill(10);
745 
746  Int_t Matched_Truth_Particle_PDG=0;
747 
748 
749  const Int_t D_Candidtae_N_Daughters=2;
750  Int_t D_Candidtae_Daughters_PDG[D_Candidtae_N_Daughters] = {211,321};
751  Int_t D_Candidate_MatchedTruth_Label = D_Candidate->MatchToMC(fCandidatePDG, Particle_Container->GetArray(), D_Candidtae_N_Daughters, D_Candidtae_Daughters_PDG);
752  Bool_t Is_Prompt_Correct_Quark=kFALSE;
753  Int_t Is_Prompt_Correct_Quark_PDG=-1;
754  AliAODMCParticle* Matched_Truth_Particle;
755  if (D_Candidate_MatchedTruth_Label >= 0) {
756  Matched_Truth_Particle = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(D_Candidate_MatchedTruth_Label));
757  if (Matched_Truth_Particle) {
758  fhEvent->Fill(11);
759 
760  Int_t Matched_Truth_Particle_Mother_Label=Matched_Truth_Particle->GetMother();
761  while (Matched_Truth_Particle_Mother_Label >= 0) {
762  AliAODMCParticle* Matched_Truth_Particle_Mother = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Matched_Truth_Particle_Mother_Label));
763  if (Matched_Truth_Particle_Mother){
764  Int_t Original_Quark_PDG=4;
765  if (fIsBDecay) Original_Quark_PDG=5;
766  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode()) == 4){
767  fhEvent->Fill(12);
768  Is_Prompt_Correct_Quark_PDG=4;
769  }
770  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode()) == 5){
771  fhEvent->Fill(13);
772  Is_Prompt_Correct_Quark_PDG=5;
773  }
774  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode())==Original_Quark_PDG) Is_Prompt_Correct_Quark=kTRUE;
775  if (Matched_Truth_Particle_Mother_Label==Matched_Truth_Particle_Mother->GetMother()) break;
776  Matched_Truth_Particle_Mother_Label=Matched_Truth_Particle_Mother->GetMother();
777  }
778  else break;
779  // delete Matched_Truth_Particle_Mother;
780  }
781  Matched_Truth_Particle_PDG = Matched_Truth_Particle->PdgCode();
782  }
783  }
784  else continue;
785  if (fPromptReject && !Is_Prompt_Correct_Quark) continue;
786 
787 
788  if (Mass_Hypo_Type==1 && Matched_Truth_Particle_PDG!=fCandidatePDG){
789  fhEvent->Fill(14);
790  continue;
791  }
792  else fhEvent->Fill(15);
793  if (Mass_Hypo_Type==2 && Matched_Truth_Particle_PDG!=-fCandidatePDG){
794  fhEvent->Fill(16);
795  continue;
796  }
797  else fhEvent->Fill(17);
798  if (Mass_Hypo_Type==3 && TMath::Abs(Matched_Truth_Particle_PDG)!=fCandidatePDG){
799  fhEvent->Fill(18);
800  continue;
801  }
802  else{
803  fhEvent->Fill(19);
804  if (Matched_Truth_Particle_PDG==fCandidatePDG) fhEvent->Fill(20);
805  if (Matched_Truth_Particle_PDG==-fCandidatePDG) fhEvent->Fill(21);
806  }
807 
808 
809  Double_t Inv_Mass_D=0.0;
810  Double_t Inv_Mass_D_Truth=0.0;
811 
812  if (Mass_Hypo_Type==1 || (Mass_Hypo_Type==3 && Matched_Truth_Particle_PDG==fCandidatePDG)){
813  Inv_Mass_D=D_Candidate->InvMassD0();
814  //Inv_Mass_D_Truth=Matched_Truth_Particle->InvMassD0();
815  Inv_Mass_D_Truth=0.0;
816  }
817  if (Mass_Hypo_Type==2 || (Mass_Hypo_Type==3 && Matched_Truth_Particle_PDG==-fCandidatePDG)){
818  Inv_Mass_D=D_Candidate->InvMassD0bar();
819  //Inv_Mass_D_Truth=Matched_Truth_Particle->InvMassD0bar();
820  Inv_Mass_D_Truth=0.0;
821  }
822 
823 
824 
825 
827  AliTLorentzVector D_Candidate_LorentzVector(0,0,0,0);
828  D_Candidate_LorentzVector.SetPtEtaPhiM(D_Candidate->Pt(), D_Candidate->Eta(), D_Candidate->Phi(), Inv_Mass_D);
829  fFastJetWrapper->AddInputVector(D_Candidate_LorentzVector.Px(), D_Candidate_LorentzVector.Py(), D_Candidate_LorentzVector.Pz(), D_Candidate_LorentzVector.E(), 0);
830 
831 
832  if (!Track_Container) continue;
833  Track_Container->SetDMesonCandidate(D_Candidate);
834  AliAODTrack *Track = NULL;
835  for (Int_t i_Track=0; i_Track<Track_Container->GetNTracks(); i_Track++){
836  Track = static_cast<AliAODTrack*>(Track_Container->GetAcceptParticle(i_Track));
837  if(!Track) continue;
838  fFastJetWrapper->AddInputVector(Track->Px(), Track->Py(), Track->Pz(), Track->E(),i_Track+100);
839  }
840  // delete Track;
841  fFastJetWrapper->Run();
842 
843 
844  std::vector<fastjet::PseudoJet> Inclusive_Jets = fFastJetWrapper->GetInclusiveJets();
845  for (UInt_t i_Jet=0; i_Jet < Inclusive_Jets.size(); i_Jet++){
846  Bool_t Is_D_Jet=kFALSE;
847  if (Inclusive_Jets[i_Jet].perp()<fJetMinPt) continue;
848  if (TMath::Abs(Inclusive_Jets[i_Jet].pseudorapidity()) > 0.9-fJetRadius) continue;
849  std::vector<fastjet::PseudoJet> Constituents(fFastJetWrapper->GetJetConstituents(i_Jet));
850  for (UInt_t i_Constituents = 0; i_Constituents < Constituents.size(); i_Constituents++) {
851  if (Constituents[i_Constituents].user_index() == 0) {
852  Is_D_Jet = kTRUE;
853  }
854  }
855 
856 
857  if (!Is_D_Jet) continue;
858  fhEvent->Fill(22);
859 
860 
861  Int_t i_Matched_D_Jet_Truth=-1;
862  for (UInt_t k=0; k< Inclusive_Jets_Truth_Labels.size(); k++){
863  if(Inclusive_Jets_Truth_Labels[k].first==D_Candidate_MatchedTruth_Label) i_Matched_D_Jet_Truth=Inclusive_Jets_Truth_Labels[k].second;
864  }
865 
866  for (UInt_t i_Jet_Truth=0; i_Jet_Truth < Inclusive_Jets_Truth.size(); i_Jet_Truth++){
867  Bool_t Is_Jet_Truth_Matched=kFALSE;
868  if (TMath::Abs(Inclusive_Jets_Truth[i_Jet_Truth].pseudorapidity()) > 0.9-fJetRadius) continue;
869  std::vector<fastjet::PseudoJet> Constituents_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
870  for (UInt_t i_Constituents_Truth = 0; i_Constituents_Truth < Constituents_Truth.size(); i_Constituents_Truth++) {
871  if (Constituents_Truth[i_Constituents_Truth].user_index() == i_Matched_D_Jet_Truth) {
872  Is_Jet_Truth_Matched=kTRUE;
873  for (UInt_t i_Unmacthed_D=0; i_Unmacthed_D<Unmatched_Truth_Level_D.size(); i_Unmacthed_D++){
874  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);
875  }
876  }
877  }
878  if (!Is_Jet_Truth_Matched) continue;
879  fhEvent->Fill(23);
880  NMatched_DMeson_Jets++;
881 
882 
883 
884 
885  std::vector<Double_t> Splittings_Zg;
886  std::vector<Double_t> Splittings_DeltaR;
887  std::vector<Double_t> Splittings_LeadingSubJetpT;
888  std::vector<Double_t> Splittings_HardestSubJetD0;
889  std::vector<Double_t> Splittings_RadiatorE;
890  std::vector<Double_t> Splittings_RadiatorpT;
891 
892  Bool_t Is_D_SubJet=kFALSE;
893  fastjet::JetDefinition Jet_Definition(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
894 
895  try{
896  std::vector<fastjet::PseudoJet> Reclustered_Particles(fFastJetWrapper->GetJetConstituents(i_Jet));
897  fastjet::ClusterSequence Cluster_Sequence_CA(Reclustered_Particles, Jet_Definition);
898  std::vector<fastjet::PseudoJet> Reclustered_Jet = Cluster_Sequence_CA.inclusive_jets(0.0);
899  Reclustered_Jet = sorted_by_pt(Reclustered_Jet);
900 
901 
902  fastjet::PseudoJet Daughter_Jet = Reclustered_Jet[0];
903  fastjet::PseudoJet Parent_SubJet_1;
904  fastjet::PseudoJet Parent_SubJet_2;
905 
906  while(Daughter_Jet.has_parents(Parent_SubJet_1,Parent_SubJet_2)){
907  if(Parent_SubJet_1.perp() < Parent_SubJet_2.perp()) std::swap(Parent_SubJet_1,Parent_SubJet_2);
908  Splittings_LeadingSubJetpT.push_back(Parent_SubJet_1.perp());
909  vector < fastjet::PseudoJet > Hard_SubJet_Constituents = sorted_by_pt(Parent_SubJet_1.constituents());
910  Is_D_SubJet=kFALSE;
911  for(UInt_t j=0;j<Hard_SubJet_Constituents.size();j++){
912  if(Hard_SubJet_Constituents[j].user_index()==0) Is_D_SubJet=kTRUE;
913  }
914 
915 
916  if (!Is_D_SubJet) Splittings_HardestSubJetD0.push_back(1.0);
917  else Splittings_HardestSubJetD0.push_back(2.0);
918 
919  // if(!Is_D_SubJet) std::swap(Parent_SubJet_1,Parent_SubJet_2);
920  Splittings_DeltaR.push_back(Parent_SubJet_1.delta_R(Parent_SubJet_2));
921  Splittings_Zg.push_back(Parent_SubJet_2.perp()/(Parent_SubJet_1.perp()+Parent_SubJet_2.perp()));
922  Splittings_RadiatorE.push_back(Daughter_Jet.E());
923  Splittings_RadiatorpT.push_back(Daughter_Jet.perp());
924  Daughter_Jet=Parent_SubJet_1;
925  }
926 
927 
928  } catch (fastjet::Error) { /*return -1;*/ }
929 
930 
931  std::vector<Double_t> Splittings_Zg_Truth;
932  std::vector<Double_t> Splittings_DeltaR_Truth;
933  std::vector<Double_t> Splittings_LeadingSubJetpT_Truth;
934  std::vector<Double_t> Splittings_HardestSubJetD0_Truth;
935  std::vector<Double_t> Splittings_RadiatorE_Truth;
936  std::vector<Double_t> Splittings_RadiatorpT_Truth;
937 
938 
939  Bool_t Is_D_SubJet_Truth=kFALSE;
940  fastjet::JetDefinition Jet_Definition_Truth(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
941 
942  try{
943  std::vector<fastjet::PseudoJet> Reclustered_Particles_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
944  fastjet::ClusterSequence Cluster_Sequence_CA_Truth(Reclustered_Particles_Truth, Jet_Definition_Truth);
945  std::vector<fastjet::PseudoJet> Reclustered_Jet_Truth = Cluster_Sequence_CA_Truth.inclusive_jets(0.0);
946  Reclustered_Jet_Truth = sorted_by_pt(Reclustered_Jet_Truth);
947 
948 
949  fastjet::PseudoJet Daughter_Jet_Truth = Reclustered_Jet_Truth[0];
950  fastjet::PseudoJet Parent_SubJet_1_Truth;
951  fastjet::PseudoJet Parent_SubJet_2_Truth;
952 
953  while(Daughter_Jet_Truth.has_parents(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth)){
954  if(Parent_SubJet_1_Truth.perp() < Parent_SubJet_2_Truth.perp()) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
955  Splittings_LeadingSubJetpT_Truth.push_back(Parent_SubJet_1_Truth.perp());
956  vector < fastjet::PseudoJet > Hard_SubJet_Constituents_Truth = sorted_by_pt(Parent_SubJet_1_Truth.constituents());
957  Is_D_SubJet_Truth=kFALSE;
958  for(UInt_t j=0;j<Hard_SubJet_Constituents_Truth.size();j++){
959  if(Hard_SubJet_Constituents_Truth[j].user_index()==i_Matched_D_Jet_Truth) Is_D_SubJet_Truth=kTRUE;
960  }
961 
962 
963  if (!Is_D_SubJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(1.0);
964  else Splittings_HardestSubJetD0_Truth.push_back(2.0);
965 
966  // if(!Is_D_SubJet_Truth) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
967  Splittings_DeltaR_Truth.push_back(Parent_SubJet_1_Truth.delta_R(Parent_SubJet_2_Truth));
968  Splittings_Zg_Truth.push_back(Parent_SubJet_2_Truth.perp()/(Parent_SubJet_1_Truth.perp()+Parent_SubJet_2_Truth.perp()));
969  Splittings_RadiatorE_Truth.push_back(Daughter_Jet_Truth.E());
970  Splittings_RadiatorpT_Truth.push_back(Daughter_Jet_Truth.perp());
971  Daughter_Jet_Truth=Parent_SubJet_1_Truth;
972  }
973 
974 
975  } catch (fastjet::Error) { /*return -1;*/ }
976 
977 
978 
979 
980  //set detector level flags
981  Double_t Flag_D=-1.0;
982  if (Mass_Hypo_Type ==1) Flag_D=1.0;
983  else if (Mass_Hypo_Type ==2) Flag_D=2.0;
984  else if (Mass_Hypo_Type ==3 && Matched_Truth_Particle->GetPdgCode()==fCandidatePDG) Flag_D=3.0;
985  else if (Mass_Hypo_Type ==3 && Matched_Truth_Particle->GetPdgCode()==-fCandidatePDG) Flag_D=4.0;
986 
987  Double_t Flag_D_Truth=-1.0;
988  if (Matched_Truth_Particle->GetPdgCode()==fCandidatePDG) Flag_D_Truth=1.0;
989  if (Matched_Truth_Particle->GetPdgCode()==-fCandidatePDG) Flag_D_Truth=2.0;
990 
991 
992  fShapesVar[0] = Inclusive_Jets[i_Jet].perp();
993  fShapesVar[1] = Inclusive_Jets_Truth[i_Jet_Truth].perp();
994  fShapesVar[2] = D_Candidate->Pt();
995  fShapesVar[3] = Matched_Truth_Particle->Pt();
996  fShapesVar[4] = Inv_Mass_D;
997  fShapesVar[5] = Inv_Mass_D_Truth;
998  fShapesVar[6] = Flag_D;
999  fShapesVar[7] = Flag_D_Truth;
1000  fShapesVar[8] = Is_Prompt_Correct_Quark_PDG;
1001  fShapesVar[9] = 0.0;
1002 
1003 
1004 
1005  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR);
1006  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR_Truth);
1007  fShapesVar_Splittings_Zg.push_back(Splittings_Zg);
1008  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg_Truth);
1009  fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT);
1010  fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT_Truth);
1011  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0);
1012  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0_Truth);
1013  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE);
1014  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE_Truth);
1015  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT);
1016  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT_Truth);
1017 
1018 
1019  fTreeResponseMatrixAxis->Fill();
1020  fTreeSplittings->Fill();
1021 
1024  fShapesVar_Splittings_Zg.clear();
1034  }
1035  }
1036  // delete Matched_Truth_Particle;
1037  // delete D_Candidate;
1038  }
1039  // delete Track_Container;
1040  if(NMatched_DMeson_Jets==0) fhEvent->Fill(24);
1041  if(NMatched_DMeson_Jets==1) fhEvent->Fill(25);
1042  if(NMatched_DMeson_Jets==2) fhEvent->Fill(26);
1043  if(NMatched_DMeson_Jets==3) fhEvent->Fill(27);
1044  if(NMatched_DMeson_Jets==4) fhEvent->Fill(28);
1045  if(NMatched_DMeson_Jets==5) fhEvent->Fill(29);
1046  if(NMatched_DMeson_Jets==6) fhEvent->Fill(30);
1047 
1048  if (fIncludeInclusive){
1049 
1050  if(Unmatched_Truth_Level_D.size()==0) fhEvent->Fill(31);
1051  if(Unmatched_Truth_Level_D.size()==1) fhEvent->Fill(32);
1052  if(Unmatched_Truth_Level_D.size()==2) fhEvent->Fill(33);
1053  if(Unmatched_Truth_Level_D.size()==3) fhEvent->Fill(34);
1054  if(Unmatched_Truth_Level_D.size()==4) fhEvent->Fill(35);
1055  if(Unmatched_Truth_Level_D.size()==5) fhEvent->Fill(36);
1056  if(Unmatched_Truth_Level_D.size()==6) fhEvent->Fill(37);
1057 
1058  for (UInt_t i_Jet_Truth=0; i_Jet_Truth < Inclusive_Jets_Truth.size(); i_Jet_Truth++){
1059  AliAODMCParticle* Truth_D_Particle = NULL;
1060  Bool_t Is_Unmatched_D=kFALSE;
1061  Int_t D_Meson_Matched_Index=-1;
1062  if (TMath::Abs(Inclusive_Jets_Truth[i_Jet_Truth].pseudorapidity()) > 0.9-fJetRadius) continue;
1063  if (Inclusive_Jets_Truth[i_Jet_Truth].perp()<fJetMinPt) continue;
1064  std::vector<fastjet::PseudoJet> Constituents_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1065  for (UInt_t i_Constituents_Truth = 0; i_Constituents_Truth < Constituents_Truth.size(); i_Constituents_Truth++) {
1066  for (UInt_t i_Unmacthed_D=0; i_Unmacthed_D<Unmatched_Truth_Level_D.size(); i_Unmacthed_D++){
1067  if(Constituents_Truth[i_Constituents_Truth].user_index() == Unmatched_Truth_Level_D[i_Unmacthed_D]){
1068  Is_Unmatched_D=kTRUE;
1069  D_Meson_Matched_Index=Constituents_Truth[i_Constituents_Truth].user_index();
1070  for(UInt_t i_MC_Label=0; i_MC_Label<Inclusive_Jets_Truth_Labels.size(); i_MC_Label++){
1071  if (Inclusive_Jets_Truth_Labels[i_MC_Label].second==Constituents_Truth[i_Constituents_Truth].user_index()){
1072  Truth_D_Particle=static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Inclusive_Jets_Truth_Labels[i_MC_Label].first));
1073  }
1074  }
1075  }
1076  }
1077  }
1078  if (!Is_Unmatched_D) continue;
1079  fhEvent->Fill(38);
1080 
1081 
1082  Bool_t Is_Prompt_Correct_Quark=kFALSE;
1083  Int_t Is_Prompt_Correct_Quark_PDG=-1;
1084  Int_t Truth_D_Particle_Mother_Label=Truth_D_Particle->GetMother();
1085  while (Truth_D_Particle_Mother_Label >= 0) {
1086  AliAODMCParticle* Truth_D_Particle_Mother = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Truth_D_Particle_Mother_Label));
1087  if (Truth_D_Particle_Mother){
1088  Int_t Original_Quark_PDG=4;
1089  if (fIsBDecay) Original_Quark_PDG=5;
1090  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) == 4){
1091  fhEvent->Fill(39);
1092  Is_Prompt_Correct_Quark_PDG=4;
1093  }
1094  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) == 5){
1095  fhEvent->Fill(40);
1096  Is_Prompt_Correct_Quark_PDG=5;
1097  }
1098  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode())==Original_Quark_PDG) Is_Prompt_Correct_Quark=kTRUE;
1099  if (Truth_D_Particle_Mother_Label==Truth_D_Particle_Mother->GetMother()) break;
1100  Truth_D_Particle_Mother_Label=Truth_D_Particle_Mother->GetMother();
1101  }
1102  else break;
1103  // delete Truth_D_Particle_Mother;
1104  }
1105 
1106  if (fPromptReject && !Is_Prompt_Correct_Quark) continue;
1107 
1108 
1109  Bool_t Is_D_SubJet_Truth=kFALSE;
1110  fastjet::JetDefinition Jet_Definition_Truth(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
1111 
1112  std::vector<Double_t> Splittings_Zg_Truth;
1113  std::vector<Double_t> Splittings_DeltaR_Truth;
1114  std::vector<Double_t> Splittings_LeadingSubJetpT_Truth;
1115  std::vector<Double_t> Splittings_HardestSubJetD0_Truth;
1116  std::vector<Double_t> Splittings_RadiatorE_Truth;
1117  std::vector<Double_t> Splittings_RadiatorpT_Truth;
1118 
1119 
1120  try{
1121  std::vector<fastjet::PseudoJet> Reclustered_Particles_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1122  fastjet::ClusterSequence Cluster_Sequence_CA_Truth(Reclustered_Particles_Truth, Jet_Definition_Truth);
1123  std::vector<fastjet::PseudoJet> Reclustered_Jet_Truth = Cluster_Sequence_CA_Truth.inclusive_jets(0.0);
1124  Reclustered_Jet_Truth = sorted_by_pt(Reclustered_Jet_Truth);
1125 
1126 
1127  fastjet::PseudoJet Daughter_Jet_Truth = Reclustered_Jet_Truth[0];
1128  fastjet::PseudoJet Parent_SubJet_1_Truth;
1129  fastjet::PseudoJet Parent_SubJet_2_Truth;
1130 
1131  while(Daughter_Jet_Truth.has_parents(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth)){
1132  if(Parent_SubJet_1_Truth.perp() < Parent_SubJet_2_Truth.perp()) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1133  Splittings_LeadingSubJetpT_Truth.push_back(Parent_SubJet_1_Truth.perp());
1134  vector < fastjet::PseudoJet > Hard_SubJet_Constituents_Truth = sorted_by_pt(Parent_SubJet_1_Truth.constituents());
1135  Is_D_SubJet_Truth=kFALSE;
1136  for(UInt_t j=0;j<Hard_SubJet_Constituents_Truth.size();j++){
1137  if(Hard_SubJet_Constituents_Truth[j].user_index()==D_Meson_Matched_Index) Is_D_SubJet_Truth=kTRUE;
1138  }
1139 
1140 
1141  if (!Is_D_SubJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(1.0);
1142  else Splittings_HardestSubJetD0_Truth.push_back(2.0);
1143 
1144  // if(!Is_D_SubJet_Truth) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1145  Splittings_DeltaR_Truth.push_back(Parent_SubJet_1_Truth.delta_R(Parent_SubJet_2_Truth));
1146  Splittings_Zg_Truth.push_back(Parent_SubJet_2_Truth.perp()/(Parent_SubJet_1_Truth.perp()+Parent_SubJet_2_Truth.perp()));
1147  Splittings_RadiatorE_Truth.push_back(Daughter_Jet_Truth.E());
1148  Splittings_RadiatorpT_Truth.push_back(Daughter_Jet_Truth.perp());
1149  Daughter_Jet_Truth=Parent_SubJet_1_Truth;
1150  }
1151 
1152 
1153  } catch (fastjet::Error) { /*return -1;*/ }
1154 
1155  Double_t Inv_Mass_D_Truth=0.0;
1156  Double_t Flag_D_Truth=-1.0;
1157  if (Truth_D_Particle->GetPdgCode()==fCandidatePDG){
1158  //Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0();
1159  Inv_Mass_D_Truth=0.0;
1160  Flag_D_Truth=3.0;
1161  }
1162  if (Truth_D_Particle->GetPdgCode()==-fCandidatePDG){
1163  // Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0bar();
1164  Inv_Mass_D_Truth=0.0;
1165  Flag_D_Truth=4.0;
1166  }
1167 
1168 
1169  fShapesVar[0] = 0.0;
1170  fShapesVar[1] = Inclusive_Jets_Truth[i_Jet_Truth].perp();
1171  fShapesVar[2] = 0.0;
1172  fShapesVar[3] = Truth_D_Particle->Pt();
1173  fShapesVar[4] = 0.0;
1174  fShapesVar[5] = Inv_Mass_D_Truth;
1175  fShapesVar[6] = 0.0;
1176  fShapesVar[7] = Flag_D_Truth;
1177  fShapesVar[8] = Is_Prompt_Correct_Quark_PDG;
1178  fShapesVar[9] = 0.0;
1179 
1180 
1181  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR_Truth);
1182  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR_Truth);
1183  fShapesVar_Splittings_Zg.push_back(Splittings_Zg_Truth);
1184  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg_Truth);
1185  fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT_Truth);
1186  fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT_Truth);
1187  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0_Truth);
1188  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0_Truth);
1189  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE_Truth);
1190  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE_Truth);
1191  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT_Truth);
1192  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT_Truth);
1193 
1194 
1195  fTreeResponseMatrixAxis->Fill();
1196  fTreeSplittings->Fill();
1197 
1200  fShapesVar_Splittings_Zg.clear();
1210  // delete Truth_D_Particle;
1211  }
1212  }
1213  // delete Particle_Container;
1214  }
1215 
1216 
1217  if (fJetShapeType == kTrue){
1218 
1219 
1220  //truth level only jet finding
1221 
1222 
1224  Particle_Container->SetSpecialPDG(fCandidatePDG);
1225  Particle_Container->SetRejectedOriginMap(fRejectedOrigin);
1226  Particle_Container->SetAcceptedDecayMap(kDecayD0toKpi);
1227  Particle_Container->SetRejectISR(fRejectISR);
1228 
1229  std::vector<fastjet::PseudoJet> Inclusive_Jets_Truth;
1230  std::vector<std::pair<Int_t, Int_t>> Inclusive_Jets_Truth_Labels;
1231 
1232  if (Particle_Container->IsSpecialPDGFound()){
1233  fhEvent->Fill(2);
1234 
1235  fFastJetWrapper_Truth->SetAreaType(fastjet::active_area);
1238  fFastJetWrapper_Truth->SetAlgorithm(fastjet::antikt_algorithm);
1239  fFastJetWrapper_Truth->SetRecombScheme(static_cast<fastjet::RecombinationScheme>(0));
1241 
1242  AliAODMCParticle* Truth_Particle=NULL;
1243  Int_t NTruthD=0;
1244  for (Int_t i_Particle=0; i_Particle<Particle_Container->GetNParticles(); i_Particle++){
1245  Truth_Particle = static_cast<AliAODMCParticle*>(Particle_Container->GetAcceptParticle(i_Particle));
1246  if (!Truth_Particle) continue;
1247  if (TMath::Abs(Truth_Particle->GetPdgCode())==fCandidatePDG){
1248  fhEvent->Fill(3);
1249  std::pair<Int_t, Int_t> Inclusive_Jet_Truth_Labels;
1250  Inclusive_Jet_Truth_Labels.first=Truth_Particle->GetLabel();
1251  Inclusive_Jet_Truth_Labels.second=NTruthD;
1252  Inclusive_Jets_Truth_Labels.push_back(Inclusive_Jet_Truth_Labels);
1253  fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),NTruthD);
1254  NTruthD++;
1255  }
1256  else fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),i_Particle+100);
1257  }
1258  // delete Truth_Particle;
1260  Inclusive_Jets_Truth = fFastJetWrapper_Truth->GetInclusiveJets();
1261 
1262  for (UInt_t i_Jet_Truth=0; i_Jet_Truth < Inclusive_Jets_Truth.size(); i_Jet_Truth++){
1263  AliAODMCParticle* Truth_D_Particle = NULL;
1264  Bool_t Is_DJet_Truth=kFALSE;
1265  Int_t D_Meson_Matched_Index=-1;
1266  if (TMath::Abs(Inclusive_Jets_Truth[i_Jet_Truth].pseudorapidity()) > 0.9-fJetRadius) continue;
1267  std::vector<fastjet::PseudoJet> Constituents_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1268  for (UInt_t i_Constituents_Truth = 0; i_Constituents_Truth < Constituents_Truth.size(); i_Constituents_Truth++) {
1269 
1270  if(Constituents_Truth[i_Constituents_Truth].user_index() >=0 && Constituents_Truth[i_Constituents_Truth].user_index() < NTruthD){
1271  D_Meson_Matched_Index=Constituents_Truth[i_Constituents_Truth].user_index();
1272  Is_DJet_Truth=kTRUE;
1273  for(UInt_t i_MC_Label=0; i_MC_Label<Inclusive_Jets_Truth_Labels.size(); i_MC_Label++){
1274  if (Inclusive_Jets_Truth_Labels[i_MC_Label].second==Constituents_Truth[i_Constituents_Truth].user_index()){
1275  Truth_D_Particle=static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Inclusive_Jets_Truth_Labels[i_MC_Label].first));
1276  }
1277  }
1278 
1279  }
1280  }
1281 
1282  if (!Is_DJet_Truth) fhEvent->Fill(4);
1283  if (Is_DJet_Truth) fhEvent->Fill(5);
1284 
1285  if (!fIncludeInclusive && !Is_DJet_Truth) continue;
1286 
1287  Bool_t Is_Prompt_Correct_Quark=kFALSE;
1288  Int_t Is_Prompt_Correct_Quark_PDG=-1;
1289  Int_t Truth_D_Particle_Mother_Label=Truth_D_Particle->GetMother();
1290  while (Truth_D_Particle_Mother_Label >= 0) {
1291  AliAODMCParticle* Truth_D_Particle_Mother = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Truth_D_Particle_Mother_Label));
1292  if (Truth_D_Particle_Mother){
1293  Int_t Original_Quark_PDG=4;
1294  if (fIsBDecay) Original_Quark_PDG=5;
1295  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) ==Original_Quark_PDG) Is_Prompt_Correct_Quark=kTRUE;
1296  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) ==4){
1297  fhEvent->Fill(6);
1298  Is_Prompt_Correct_Quark_PDG=4;
1299  }
1300  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) ==5){
1301  fhEvent->Fill(7);
1302  Is_Prompt_Correct_Quark_PDG=5;
1303  }
1304  if (Truth_D_Particle_Mother_Label==Truth_D_Particle_Mother->GetMother()) break;
1305  Truth_D_Particle_Mother_Label=Truth_D_Particle_Mother->GetMother();
1306  }
1307  else break;
1308  // delete Truth_D_Particle_Mother;
1309  }
1310 
1311  if(fPromptReject && !Is_Prompt_Correct_Quark) continue;
1312  fhEvent->Fill(8);
1313 
1314  Bool_t Is_D_SubJet_Truth=kFALSE;
1315  fastjet::JetDefinition Jet_Definition_Truth(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
1316 
1317  std::vector<Double_t> Splittings_Zg_Truth;
1318  std::vector<Double_t> Splittings_DeltaR_Truth;
1319  std::vector<Double_t> Splittings_LeadingSubJetpT_Truth;
1320  std::vector<Double_t> Splittings_HardestSubJetD0_Truth;
1321  std::vector<Double_t> Splittings_RadiatorE_Truth;
1322  std::vector<Double_t> Splittings_RadiatorpT_Truth;
1323 
1324  try{
1325  std::vector<fastjet::PseudoJet> Reclustered_Particles_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1326  fastjet::ClusterSequence Cluster_Sequence_CA_Truth(Reclustered_Particles_Truth, Jet_Definition_Truth);
1327  std::vector<fastjet::PseudoJet> Reclustered_Jet_Truth = Cluster_Sequence_CA_Truth.inclusive_jets(0.0);
1328  Reclustered_Jet_Truth = sorted_by_pt(Reclustered_Jet_Truth);
1329 
1330 
1331  fastjet::PseudoJet Daughter_Jet_Truth = Reclustered_Jet_Truth[0];
1332  fastjet::PseudoJet Parent_SubJet_1_Truth;
1333  fastjet::PseudoJet Parent_SubJet_2_Truth;
1334 
1335  while(Daughter_Jet_Truth.has_parents(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth)){
1336  if(Parent_SubJet_1_Truth.perp() < Parent_SubJet_2_Truth.perp()) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1337  Splittings_LeadingSubJetpT_Truth.push_back(Parent_SubJet_1_Truth.perp());
1338  vector < fastjet::PseudoJet > Hard_SubJet_Constituents_Truth = sorted_by_pt(Parent_SubJet_1_Truth.constituents());
1339  Is_D_SubJet_Truth=kFALSE;
1340  for(UInt_t j=0;j<Hard_SubJet_Constituents_Truth.size();j++){
1341  if(Hard_SubJet_Constituents_Truth[j].user_index()==D_Meson_Matched_Index) Is_D_SubJet_Truth=kTRUE;
1342  }
1343 
1344  if (fIncludeInclusive && !Is_DJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(0.0);
1345  else if (Is_DJet_Truth && !Is_D_SubJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(1.0);
1346  else if (Is_DJet_Truth && Is_D_SubJet_Truth)Splittings_HardestSubJetD0_Truth.push_back(2.0);
1347 
1348  // if(!Is_D_SubJet_Truth) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1349  Splittings_DeltaR_Truth.push_back(Parent_SubJet_1_Truth.delta_R(Parent_SubJet_2_Truth));
1350  Splittings_Zg_Truth.push_back(Parent_SubJet_2_Truth.perp()/(Parent_SubJet_1_Truth.perp()+Parent_SubJet_2_Truth.perp()));
1351  Splittings_RadiatorE_Truth.push_back(Daughter_Jet_Truth.E());
1352  Splittings_RadiatorpT_Truth.push_back(Daughter_Jet_Truth.perp());
1353  Daughter_Jet_Truth=Parent_SubJet_1_Truth;
1354  }
1355 
1356 
1357  } catch (fastjet::Error) { /*return -1;*/ }
1358 
1359  Double_t Inv_Mass_D_Truth=0.0;
1360  Double_t Flag_D_Truth=-1.0;
1361  Double_t D_Pt=-1.0;
1362  if (fIncludeInclusive && !Is_DJet_Truth){
1363  Flag_D_Truth=0.0;
1364  }
1365  else if (Is_DJet_Truth){
1366  if (Truth_D_Particle->GetPdgCode()==fCandidatePDG){
1367  //Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0();
1368  Inv_Mass_D_Truth=0.0;
1369  Flag_D_Truth=1.0;
1370  D_Pt=Truth_D_Particle->Pt();
1371  fhEvent->Fill(9);
1372  }
1373  if (Truth_D_Particle->GetPdgCode()==-fCandidatePDG){
1374  // Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0bar();
1375  Inv_Mass_D_Truth=0.0;
1376  Flag_D_Truth=2.0;
1377  D_Pt=Truth_D_Particle->Pt();
1378  fhEvent->Fill(10);
1379  }
1380  }
1381 
1382 
1383  fShapesVar[0] = 0.0;
1384  fShapesVar[1] = Inclusive_Jets_Truth[i_Jet_Truth].perp();
1385  fShapesVar[2] = 0.0;
1386  fShapesVar[3] = D_Pt;
1387  fShapesVar[4] = 0.0;
1388  fShapesVar[5] = Inv_Mass_D_Truth;
1389  fShapesVar[6] = 0.0;
1390  fShapesVar[7] = Flag_D_Truth;
1391  fShapesVar[8] = 0.0;
1392  fShapesVar[9] = Is_Prompt_Correct_Quark_PDG;
1393 
1394 
1395 
1396  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR_Truth);
1397  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR_Truth);
1398  fShapesVar_Splittings_Zg.push_back(Splittings_Zg_Truth);
1399  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg_Truth);
1400  fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT_Truth);
1401  fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT_Truth);
1402  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0_Truth);
1403  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0_Truth);
1404  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE_Truth);
1405  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE_Truth);
1406  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT_Truth);
1407  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT_Truth);
1408 
1409 
1410  fTreeResponseMatrixAxis->Fill();
1411  fTreeSplittings->Fill();
1412 
1415  fShapesVar_Splittings_Zg.clear();
1425 
1426 
1427  // delete Truth_D_Particle;
1428  }
1429  if (NTruthD==0) fhEvent->Fill(11);
1430  if (NTruthD==1) fhEvent->Fill(12);
1431  if (NTruthD==2) fhEvent->Fill(13);
1432  if (NTruthD==3) fhEvent->Fill(14);
1433  if (NTruthD==4) fhEvent->Fill(15);
1434  if (NTruthD==5) fhEvent->Fill(16);
1435  if (NTruthD==6) fhEvent->Fill(17);
1436  }
1437  //delete Particle_Container;
1438  }
1439 
1440 
1441 
1442 
1443  if (fJetShapeType == kDataInclusive){
1444 
1445  AliHFTrackContainer* Track_Container=dynamic_cast<AliHFTrackContainer*>(GetTrackContainer(0));
1446  if (!Track_Container) return kTRUE;
1447  //Track_Container->SetDMesonCandidate(NULL);
1449  AliAODTrack *Track = NULL;
1450  for (Int_t i_Track=0; i_Track<Track_Container->GetNTracks(); i_Track++){
1451  Track = static_cast<AliAODTrack*>(Track_Container->GetAcceptParticle(i_Track));
1452  if(!Track) continue;
1453  fFastJetWrapper->AddInputVector(Track->Px(), Track->Py(), Track->Pz(), Track->E(),i_Track+100);
1454  }
1455  fFastJetWrapper->Run();
1456  std::vector<fastjet::PseudoJet> Inclusive_Jets = fFastJetWrapper->GetInclusiveJets();
1457  for (UInt_t i_Jet=0; i_Jet < Inclusive_Jets.size(); i_Jet++){
1458  if (Inclusive_Jets[i_Jet].perp()<fJetMinPt) continue;
1459  if (TMath::Abs(Inclusive_Jets[i_Jet].pseudorapidity()) > 0.9-fJetRadius) continue;
1460  fhEvent->Fill(23);
1461  std::vector<Double_t> Splittings_Zg;
1462  std::vector<Double_t> Splittings_DeltaR;
1463  std::vector<Double_t> Splittings_LeadingSubJetpT;
1464  std::vector<Double_t> Splittings_HardestSubJetD0;
1465  std::vector<Double_t> Splittings_RadiatorE;
1466  std::vector<Double_t> Splittings_RadiatorpT;
1467 
1468  fastjet::JetDefinition Jet_Definition(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
1469 
1470  try{
1471  std::vector<fastjet::PseudoJet> Reclustered_Particles(fFastJetWrapper->GetJetConstituents(i_Jet));
1472  fastjet::ClusterSequence Cluster_Sequence_CA(Reclustered_Particles, Jet_Definition);
1473  std::vector<fastjet::PseudoJet> Reclustered_Jet = Cluster_Sequence_CA.inclusive_jets(0.0);
1474  Reclustered_Jet = sorted_by_pt(Reclustered_Jet);
1475 
1476 
1477  fastjet::PseudoJet Daughter_Jet = Reclustered_Jet[0];
1478  fastjet::PseudoJet Parent_SubJet_1;
1479  fastjet::PseudoJet Parent_SubJet_2;
1480 
1481  while(Daughter_Jet.has_parents(Parent_SubJet_1,Parent_SubJet_2)){
1482  if(Parent_SubJet_1.perp() < Parent_SubJet_2.perp()) std::swap(Parent_SubJet_1,Parent_SubJet_2);
1483  Double_t Leading_Track_pT=0.0;
1484  vector <fastjet::PseudoJet> Hard_SubJet_Constituents = sorted_by_pt(Parent_SubJet_1.constituents());
1485  for (Int_t i_Leading_Track=0; i_Leading_Track < Hard_SubJet_Constituents.size(); i_Leading_Track++){
1486  if (Hard_SubJet_Constituents[i_Leading_Track].perp() > Leading_Track_pT) Leading_Track_pT=Hard_SubJet_Constituents[i_Leading_Track].perp();
1487  }
1488  Splittings_LeadingSubJetpT.push_back(Parent_SubJet_1.perp());
1489  Splittings_HardestSubJetD0.push_back(Leading_Track_pT);
1490  Splittings_DeltaR.push_back(Parent_SubJet_1.delta_R(Parent_SubJet_2));
1491  Splittings_Zg.push_back(Parent_SubJet_2.perp()/(Parent_SubJet_1.perp()+Parent_SubJet_2.perp()));
1492  Splittings_RadiatorE.push_back(Daughter_Jet.E());
1493  Splittings_RadiatorpT.push_back(Daughter_Jet.perp());
1494  Daughter_Jet=Parent_SubJet_1;
1495  }
1496 
1497 
1498  } catch (fastjet::Error) { /*return -1;*/ }
1499 
1500 
1501 
1502  fShapesVar[0] = Inclusive_Jets[i_Jet].perp();
1503  fShapesVar[1] = 0.0;
1504  fShapesVar[2] = 0.0;
1505  fShapesVar[3] = 0.0;
1506  fShapesVar[4] = 0.0;
1507  fShapesVar[5] = 0.0;
1508  fShapesVar[6] = 0.0;
1509  fShapesVar[7] = 0.0;
1510  fShapesVar[8] = 0.0;
1511  fShapesVar[9] = 0.0;
1512 
1513 
1514  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR);
1515  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR);
1516  fShapesVar_Splittings_Zg.push_back(Splittings_Zg);
1517  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg);
1518  fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT);
1519  fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT);
1520  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0);
1521  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0);
1522  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE);
1523  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE);
1524  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT);
1525  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT);
1526 
1527 
1528  fTreeResponseMatrixAxis->Fill();
1529  fTreeSplittings->Fill();
1530 
1533  fShapesVar_Splittings_Zg.clear();
1543 
1544  }
1545  }
1546 
1547 
1548 
1549 
1550 
1551 
1552 
1553 
1554 
1555 
1556 
1557 
1558 
1559  delete fFastJetWrapper_Truth;
1560  delete fFastJetWrapper;
1561 
1562  return kTRUE;
1563 
1564 
1565 }
1566 
1567 
1568 
1569 
1570 
1571 
1572 //________________________________________________________________________
1574  //
1575  // retrieve event objects
1576  //
1578  return kFALSE;
1579 
1580  return kTRUE;
1581 }
1582 
1583 
1584 //_______________________________________________________________________
1586 {
1587  // Called once at the end of the analysis.
1588 
1589 }
std::vector< std::vector< Double_t > > fShapesVar_Splittings_HardestSubJetD0_Truth
void SetDMesonCandidate(AliAODRecoDecay *c)
virtual Int_t Run()
double Double_t
Definition: External.C:58
Int_t GetNTracks() const
std::vector< std::vector< Double_t > > fShapesVar_Splittings_RadiatorE_Truth
Base task in the EMCAL framework.
std::vector< std::vector< Double_t > > fShapesVar_Splittings_HardestSubJetD0
Int_t GetNParticles() const
static Int_t CheckMatchingAODdeltaAODevents()
std::vector< std::vector< Double_t > > fShapesVar_Splittings_LeadingSubJetpT
std::vector< std::vector< Double_t > > fShapesVar_Splittings_RadiatorpT_Truth
std::vector< std::vector< Double_t > > fShapesVar_Splittings_Zg
Bool_t FillHistograms()
Function filling histograms.
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:122
void SetR(Double_t r)
Definition: AliFJWrapper.h:127
std::vector< std::vector< Double_t > > fShapesVar_Splittings_DeltaR_Truth
TString kData
Declare data MC or deltaAOD.
Select tracks based on specific prescriptions of HF analysis.
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:34
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
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.
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:121
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Definition: External.C:212
TClonesArray * GetArray() const
virtual void Clear(const Option_t *="")
virtual Bool_t RetrieveEventObjects()
Retrieve common objects from event.
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)
Select MC particles based on specific prescriptions of HF analysis.
std::vector< std::vector< Double_t > > fShapesVar_Splittings_Zg_Truth
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:125
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
void Track()
Bool_t IsEventSelected(AliVEvent *event)
AliEmcalList * fOutput
!output list
AliTrackContainer * GetTrackContainer(Int_t i=0) const
void SetMakeGeneralHistograms(Bool_t g)
Enable general histograms.
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:306
std::vector< std::vector< Double_t > > fShapesVar_Splittings_DeltaR
const char Option_t
Definition: External.C:48
virtual void AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:123
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:328
std::vector< std::vector< Double_t > > fShapesVar_Splittings_LeadingSubJetpT_Truth
void swap(PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &first, PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &second)
std::vector< std::vector< Double_t > > fShapesVar_Splittings_RadiatorpT
std::vector< std::vector< Double_t > > fShapesVar_Splittings_RadiatorE