AliPhysics  9b6b435 (9b6b435)
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==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 
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();
527  fShapesVar_Splittings_Zg.clear();
537  }
538  }
539  //delete D_Candidate;
540  }
541  if(N_DMesons==0) fhEvent->Fill(16);
542  if(N_DMesons==1) fhEvent->Fill(17);
543  if(N_DMesons==2) fhEvent->Fill(18);
544  if(N_DMesons==3) fhEvent->Fill(19);
545  if(N_DMesons==4) fhEvent->Fill(20);
546  if(N_DMesons==5) fhEvent->Fill(21);
547  if(N_DMesons==6) fhEvent->Fill(22);
548 
549  if (fIncludeInclusive){
551 
552  for (UInt_t i_D_Found=0; i_D_Found<D_Candidates_Vector.size(); i_D_Found++){
553  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);
554  //Track_Container->SetDMesonCandidate(D_Candidates_Vector[i_D_Found]);
555  }
556 
557  AliAODTrack *Track = NULL;
558  Bool_t DMeson_Daughter_Track=kFALSE;
559  for (Int_t i_Track=0; i_Track<Track_Container->GetNTracks(); i_Track++){
560  for (UInt_t i_D_Found=0; i_D_Found<D_Candidates_Vector.size(); i_D_Found++){
561  Track_Container->SetDMesonCandidate(D_Candidates_Vector[i_D_Found]);
562  Track = static_cast<AliAODTrack*>(Track_Container->GetAcceptParticle(i_Track));
563  if(!Track) DMeson_Daughter_Track=kTRUE;
564  }
565  if (DMeson_Daughter_Track) continue;
566  Track = static_cast<AliAODTrack*>(Track_Container->GetAcceptParticle(i_Track));
567  fFastJetWrapper->AddInputVector(Track->Px(), Track->Py(), Track->Pz(), Track->E(),i_Track+100);
568  }
569  //delete Track;
570  fFastJetWrapper->Run();
571  std::vector<fastjet::PseudoJet> Inclusive_Jets = fFastJetWrapper->GetInclusiveJets();
572  for (UInt_t i_Jet=0; i_Jet < Inclusive_Jets.size(); i_Jet++){
573  if (Inclusive_Jets[i_Jet].perp()<fJetMinPt) continue;
574  if (TMath::Abs(Inclusive_Jets[i_Jet].pseudorapidity()) > 0.9-fJetRadius) continue;
575  Bool_t Is_D_Jet = kFALSE;
576  std::vector<fastjet::PseudoJet> Constituents(fFastJetWrapper->GetJetConstituents(i_Jet));
577  for (UInt_t i_Constituents = 0; i_Constituents < Constituents.size(); i_Constituents++) {
578  for (UInt_t i_D_Found=0; i_D_Found<D_Candidates_Vector.size(); i_D_Found++){
579  if (Constituents[i_Constituents].user_index() == i_D_Found) {
580  Is_D_Jet = kTRUE;
581  }
582  }
583  }
584  if (Is_D_Jet) continue;
585  fhEvent->Fill(23);
586  std::vector<Double_t> Splittings_Zg;
587  std::vector<Double_t> Splittings_DeltaR;
588  std::vector<Double_t> Splittings_LeadingSubJetpT;
589  std::vector<Double_t> Splittings_HardestSubJetD0;
590  std::vector<Double_t> Splittings_RadiatorE;
591  std::vector<Double_t> Splittings_RadiatorpT;
592 
593 
594  fastjet::JetDefinition Jet_Definition(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
595 
596  try{
597  std::vector<fastjet::PseudoJet> Reclustered_Particles(fFastJetWrapper->GetJetConstituents(i_Jet));
598  fastjet::ClusterSequence Cluster_Sequence_CA(Reclustered_Particles, Jet_Definition);
599  std::vector<fastjet::PseudoJet> Reclustered_Jet = Cluster_Sequence_CA.inclusive_jets(0.0);
600  Reclustered_Jet = sorted_by_pt(Reclustered_Jet);
601 
602 
603  fastjet::PseudoJet Daughter_Jet = Reclustered_Jet[0];
604  fastjet::PseudoJet Parent_SubJet_1;
605  fastjet::PseudoJet Parent_SubJet_2;
606 
607  while(Daughter_Jet.has_parents(Parent_SubJet_1,Parent_SubJet_2)){
608  if(Parent_SubJet_1.perp() < Parent_SubJet_2.perp()) std::swap(Parent_SubJet_1,Parent_SubJet_2);
609  Splittings_LeadingSubJetpT.push_back(Parent_SubJet_1.perp());
610  Splittings_HardestSubJetD0.push_back(0.0);
611  Splittings_DeltaR.push_back(Parent_SubJet_1.delta_R(Parent_SubJet_2));
612  Splittings_Zg.push_back(Parent_SubJet_2.perp()/(Parent_SubJet_1.perp()+Parent_SubJet_2.perp()));
613  Splittings_RadiatorE.push_back(Daughter_Jet.E());
614  Splittings_RadiatorpT.push_back(Daughter_Jet.perp());
615  Daughter_Jet=Parent_SubJet_1;
616  }
617 
618 
619  } catch (fastjet::Error) { /*return -1;*/ }
620 
621 
622 
623  fShapesVar[0] = Inclusive_Jets[i_Jet].perp();
624  fShapesVar[1] = 0.0;
625  fShapesVar[2] = 0.0;
626  fShapesVar[3] = 0.0;
627  fShapesVar[4] = 0.0;
628  fShapesVar[5] = 0.0;
629  fShapesVar[6] = 0.0;
630  fShapesVar[7] = 0.0;
631  fShapesVar[8] = 0.0;
632  fShapesVar[9] = 0.0;
633 
634 
635  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR);
636  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR);
637  fShapesVar_Splittings_Zg.push_back(Splittings_Zg);
638  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg);
639  fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT);
640  fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT);
641  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0);
642  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0);
643  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE);
644  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE);
645  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT);
646  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT);
647 
648 
649  fTreeResponseMatrixAxis->Fill();
650  fTreeSplittings->Fill();
651 
654  fShapesVar_Splittings_Zg.clear();
664 
665  }
666  }
667  // delete Track_Container;
668  // if (fJetShapeType != kData) delete Particle_Container;
669  }
670 
671 
672 
673  if (fJetShapeType == kTrueDet){
674 
675 
677  Particle_Container->SetSpecialPDG(fCandidatePDG);
678  Particle_Container->SetRejectedOriginMap(fRejectedOrigin);
679  Particle_Container->SetAcceptedDecayMap(fAcceptedDecay);
680  Particle_Container->SetRejectISR(fRejectISR);
681 
682  std::vector<fastjet::PseudoJet> Inclusive_Jets_Truth;
683  std::vector<std::pair<Int_t, Int_t>> Inclusive_Jets_Truth_Labels;
684  std::vector<Int_t> Unmatched_Truth_Level_D;
685  Int_t NMatched_DMeson_Jets=0;
686 
687  if (Particle_Container->IsSpecialPDGFound()){
688 
689  fhEvent->Fill(2);
690 
691  fFastJetWrapper_Truth->SetAreaType(fastjet::active_area);
694  fFastJetWrapper_Truth->SetAlgorithm(fastjet::antikt_algorithm);
695  fFastJetWrapper_Truth->SetRecombScheme(static_cast<fastjet::RecombinationScheme>(0));
697 
698  AliAODMCParticle* Truth_Particle=NULL;
699  Int_t NTruthD=0;
700  for (Int_t i_Particle=0; i_Particle<Particle_Container->GetNParticles(); i_Particle++){
701  Truth_Particle = static_cast<AliAODMCParticle*>(Particle_Container->GetAcceptParticle(i_Particle));
702  if (!Truth_Particle) continue;
703  if (TMath::Abs(Truth_Particle->PdgCode())==fCandidatePDG){
704  std::pair<Int_t, Int_t> Inclusive_Jet_Truth_Labels;
705  Inclusive_Jet_Truth_Labels.first=Truth_Particle->GetLabel();
706  Inclusive_Jet_Truth_Labels.second=NTruthD;
707  Inclusive_Jets_Truth_Labels.push_back(Inclusive_Jet_Truth_Labels);
708  Unmatched_Truth_Level_D.push_back(NTruthD);
709  fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),NTruthD);
710  NTruthD++;
711  }
712  else fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),i_Particle+100);
713  }
714  // delete Truth_Particle;
716  Inclusive_Jets_Truth = fFastJetWrapper_Truth->GetInclusiveJets();
717  if (NTruthD==0) fhEvent->Fill(3);
718  if (NTruthD==1) fhEvent->Fill(4);
719  if (NTruthD==2) fhEvent->Fill(5);
720  if (NTruthD==3) fhEvent->Fill(6);
721  if (NTruthD==4) fhEvent->Fill(7);
722  if (NTruthD==5) fhEvent->Fill(8);
723  if (NTruthD==6) fhEvent->Fill(9);
724  }
725 
726 
727 
728 
729 
730  fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject("D0toKpi"));
732  if (!Track_Container) return kTRUE;
733  Track_Container->SetDMesonCandidate(NULL);
734 
735  for (Int_t i_D = 0; i_D < fCandidateArray->GetEntriesFast(); i_D++) {
736  AliAODRecoDecayHF2Prong* D_Candidate = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(i_D));
737  if (!D_Candidate) continue;
738  if (!fRDHFCuts->IsInFiducialAcceptance(D_Candidate->Pt(), D_Candidate->Y(fCandidatePDG))) continue;
739 
740  Int_t Mass_Hypo_Type=fRDHFCuts->IsSelected(D_Candidate, AliRDHFCuts::kAll, fAodEvent);
741  if (Mass_Hypo_Type <= 0 || Mass_Hypo_Type>3) continue;
742  fhEvent->Fill(10);
743 
744  Int_t Matched_Truth_Particle_PDG=0;
745 
746 
747  const Int_t D_Candidtae_N_Daughters=2;
748  Int_t D_Candidtae_Daughters_PDG[D_Candidtae_N_Daughters] = {211,321};
749  Int_t D_Candidate_MatchedTruth_Label = D_Candidate->MatchToMC(fCandidatePDG, Particle_Container->GetArray(), D_Candidtae_N_Daughters, D_Candidtae_Daughters_PDG);
750  Bool_t Is_Prompt_Correct_Quark=kFALSE;
751  Int_t Is_Prompt_Correct_Quark_PDG=-1;
752  AliAODMCParticle* Matched_Truth_Particle;
753  if (D_Candidate_MatchedTruth_Label >= 0) {
754  Matched_Truth_Particle = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(D_Candidate_MatchedTruth_Label));
755  if (Matched_Truth_Particle) {
756  fhEvent->Fill(11);
757 
758  Int_t Matched_Truth_Particle_Mother_Label=Matched_Truth_Particle->GetMother();
759  while (Matched_Truth_Particle_Mother_Label >= 0) {
760  AliAODMCParticle* Matched_Truth_Particle_Mother = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Matched_Truth_Particle_Mother_Label));
761  if (Matched_Truth_Particle_Mother){
762  Int_t Original_Quark_PDG=4;
763  if (fIsBDecay) Original_Quark_PDG=5;
764  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode()) == 4){
765  fhEvent->Fill(12);
766  Is_Prompt_Correct_Quark_PDG=4;
767  }
768  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode()) == 5){
769  fhEvent->Fill(13);
770  Is_Prompt_Correct_Quark_PDG=5;
771  }
772  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode())==Original_Quark_PDG) Is_Prompt_Correct_Quark=kTRUE;
773  if (Matched_Truth_Particle_Mother_Label==Matched_Truth_Particle_Mother->GetMother()) break;
774  Matched_Truth_Particle_Mother_Label=Matched_Truth_Particle_Mother->GetMother();
775  }
776  else break;
777  // delete Matched_Truth_Particle_Mother;
778  }
779  Matched_Truth_Particle_PDG = Matched_Truth_Particle->PdgCode();
780  }
781  }
782  else continue;
783  if (fPromptReject && !Is_Prompt_Correct_Quark) continue;
784 
785 
786  if (Mass_Hypo_Type==1 && Matched_Truth_Particle_PDG!=fCandidatePDG){
787  fhEvent->Fill(14);
788  continue;
789  }
790  else fhEvent->Fill(15);
791  if (Mass_Hypo_Type==2 && Matched_Truth_Particle_PDG!=-fCandidatePDG){
792  fhEvent->Fill(16);
793  continue;
794  }
795  else fhEvent->Fill(17);
796  if (Mass_Hypo_Type==3 && TMath::Abs(Matched_Truth_Particle_PDG)!=fCandidatePDG){
797  fhEvent->Fill(18);
798  continue;
799  }
800  else{
801  fhEvent->Fill(19);
802  if (Matched_Truth_Particle_PDG==fCandidatePDG) fhEvent->Fill(20);
803  if (Matched_Truth_Particle_PDG==-fCandidatePDG) fhEvent->Fill(21);
804  }
805 
806 
807  Double_t Inv_Mass_D=0.0;
808  Double_t Inv_Mass_D_Truth=0.0;
809 
810  if (Mass_Hypo_Type==1 || (Mass_Hypo_Type==3 && Matched_Truth_Particle_PDG==fCandidatePDG)){
811  Inv_Mass_D=D_Candidate->InvMassD0();
812  //Inv_Mass_D_Truth=Matched_Truth_Particle->InvMassD0();
813  Inv_Mass_D_Truth=0.0;
814  }
815  if (Mass_Hypo_Type==2 || (Mass_Hypo_Type==3 && Matched_Truth_Particle_PDG==-fCandidatePDG)){
816  Inv_Mass_D=D_Candidate->InvMassD0bar();
817  //Inv_Mass_D_Truth=Matched_Truth_Particle->InvMassD0bar();
818  Inv_Mass_D_Truth=0.0;
819  }
820 
821 
822 
823 
825  AliTLorentzVector D_Candidate_LorentzVector(0,0,0,0);
826  D_Candidate_LorentzVector.SetPtEtaPhiM(D_Candidate->Pt(), D_Candidate->Eta(), D_Candidate->Phi(), Inv_Mass_D);
827  fFastJetWrapper->AddInputVector(D_Candidate_LorentzVector.Px(), D_Candidate_LorentzVector.Py(), D_Candidate_LorentzVector.Pz(), D_Candidate_LorentzVector.E(), 0);
828 
829 
830  if (!Track_Container) continue;
831  Track_Container->SetDMesonCandidate(D_Candidate);
832  AliAODTrack *Track = NULL;
833  for (Int_t i_Track=0; i_Track<Track_Container->GetNTracks(); i_Track++){
834  Track = static_cast<AliAODTrack*>(Track_Container->GetAcceptParticle(i_Track));
835  if(!Track) continue;
836  fFastJetWrapper->AddInputVector(Track->Px(), Track->Py(), Track->Pz(), Track->E(),i_Track+100);
837  }
838  // delete Track;
839  fFastJetWrapper->Run();
840 
841 
842  std::vector<fastjet::PseudoJet> Inclusive_Jets = fFastJetWrapper->GetInclusiveJets();
843  for (UInt_t i_Jet=0; i_Jet < Inclusive_Jets.size(); i_Jet++){
844  Bool_t Is_D_Jet=kFALSE;
845  if (Inclusive_Jets[i_Jet].perp()<fJetMinPt) continue;
846  if (TMath::Abs(Inclusive_Jets[i_Jet].pseudorapidity()) > 0.9-fJetRadius) continue;
847  std::vector<fastjet::PseudoJet> Constituents(fFastJetWrapper->GetJetConstituents(i_Jet));
848  for (UInt_t i_Constituents = 0; i_Constituents < Constituents.size(); i_Constituents++) {
849  if (Constituents[i_Constituents].user_index() == 0) {
850  Is_D_Jet = kTRUE;
851  }
852  }
853 
854 
855  if (!Is_D_Jet) continue;
856  fhEvent->Fill(22);
857 
858 
859  Int_t i_Matched_D_Jet_Truth=-1;
860  for (UInt_t k=0; k< Inclusive_Jets_Truth_Labels.size(); k++){
861  if(Inclusive_Jets_Truth_Labels[k].first==D_Candidate_MatchedTruth_Label) i_Matched_D_Jet_Truth=Inclusive_Jets_Truth_Labels[k].second;
862  }
863 
864  for (UInt_t i_Jet_Truth=0; i_Jet_Truth < Inclusive_Jets_Truth.size(); i_Jet_Truth++){
865  Bool_t Is_Jet_Truth_Matched=kFALSE;
866  if (TMath::Abs(Inclusive_Jets_Truth[i_Jet_Truth].pseudorapidity()) > 0.9-fJetRadius) continue;
867  std::vector<fastjet::PseudoJet> Constituents_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
868  for (UInt_t i_Constituents_Truth = 0; i_Constituents_Truth < Constituents_Truth.size(); i_Constituents_Truth++) {
869  if (Constituents_Truth[i_Constituents_Truth].user_index() == i_Matched_D_Jet_Truth) {
870  Is_Jet_Truth_Matched=kTRUE;
871  for (UInt_t i_Unmacthed_D=0; i_Unmacthed_D<Unmatched_Truth_Level_D.size(); i_Unmacthed_D++){
872  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);
873  }
874  }
875  }
876  if (!Is_Jet_Truth_Matched) continue;
877  fhEvent->Fill(23);
878  NMatched_DMeson_Jets++;
879 
880 
881 
882 
883  std::vector<Double_t> Splittings_Zg;
884  std::vector<Double_t> Splittings_DeltaR;
885  std::vector<Double_t> Splittings_LeadingSubJetpT;
886  std::vector<Double_t> Splittings_HardestSubJetD0;
887  std::vector<Double_t> Splittings_RadiatorE;
888  std::vector<Double_t> Splittings_RadiatorpT;
889 
890  Bool_t Is_D_SubJet=kFALSE;
891  fastjet::JetDefinition Jet_Definition(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
892 
893  try{
894  std::vector<fastjet::PseudoJet> Reclustered_Particles(fFastJetWrapper->GetJetConstituents(i_Jet));
895  fastjet::ClusterSequence Cluster_Sequence_CA(Reclustered_Particles, Jet_Definition);
896  std::vector<fastjet::PseudoJet> Reclustered_Jet = Cluster_Sequence_CA.inclusive_jets(0.0);
897  Reclustered_Jet = sorted_by_pt(Reclustered_Jet);
898 
899 
900  fastjet::PseudoJet Daughter_Jet = Reclustered_Jet[0];
901  fastjet::PseudoJet Parent_SubJet_1;
902  fastjet::PseudoJet Parent_SubJet_2;
903 
904  while(Daughter_Jet.has_parents(Parent_SubJet_1,Parent_SubJet_2)){
905  if(Parent_SubJet_1.perp() < Parent_SubJet_2.perp()) std::swap(Parent_SubJet_1,Parent_SubJet_2);
906  Splittings_LeadingSubJetpT.push_back(Parent_SubJet_1.perp());
907  vector < fastjet::PseudoJet > Hard_SubJet_Constituents = sorted_by_pt(Parent_SubJet_1.constituents());
908  Is_D_SubJet=kFALSE;
909  for(UInt_t j=0;j<Hard_SubJet_Constituents.size();j++){
910  if(Hard_SubJet_Constituents[j].user_index()==0) Is_D_SubJet=kTRUE;
911  }
912 
913 
914  if (!Is_D_SubJet) Splittings_HardestSubJetD0.push_back(1.0);
915  else Splittings_HardestSubJetD0.push_back(2.0);
916 
917  // if(!Is_D_SubJet) std::swap(Parent_SubJet_1,Parent_SubJet_2);
918  Splittings_DeltaR.push_back(Parent_SubJet_1.delta_R(Parent_SubJet_2));
919  Splittings_Zg.push_back(Parent_SubJet_2.perp()/(Parent_SubJet_1.perp()+Parent_SubJet_2.perp()));
920  Splittings_RadiatorE.push_back(Daughter_Jet.E());
921  Splittings_RadiatorpT.push_back(Daughter_Jet.perp());
922  Daughter_Jet=Parent_SubJet_1;
923  }
924 
925 
926  } catch (fastjet::Error) { /*return -1;*/ }
927 
928 
929  std::vector<Double_t> Splittings_Zg_Truth;
930  std::vector<Double_t> Splittings_DeltaR_Truth;
931  std::vector<Double_t> Splittings_LeadingSubJetpT_Truth;
932  std::vector<Double_t> Splittings_HardestSubJetD0_Truth;
933  std::vector<Double_t> Splittings_RadiatorE_Truth;
934  std::vector<Double_t> Splittings_RadiatorpT_Truth;
935 
936 
937  Bool_t Is_D_SubJet_Truth=kFALSE;
938  fastjet::JetDefinition Jet_Definition_Truth(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
939 
940  try{
941  std::vector<fastjet::PseudoJet> Reclustered_Particles_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
942  fastjet::ClusterSequence Cluster_Sequence_CA_Truth(Reclustered_Particles_Truth, Jet_Definition_Truth);
943  std::vector<fastjet::PseudoJet> Reclustered_Jet_Truth = Cluster_Sequence_CA_Truth.inclusive_jets(0.0);
944  Reclustered_Jet_Truth = sorted_by_pt(Reclustered_Jet_Truth);
945 
946 
947  fastjet::PseudoJet Daughter_Jet_Truth = Reclustered_Jet_Truth[0];
948  fastjet::PseudoJet Parent_SubJet_1_Truth;
949  fastjet::PseudoJet Parent_SubJet_2_Truth;
950 
951  while(Daughter_Jet_Truth.has_parents(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth)){
952  if(Parent_SubJet_1_Truth.perp() < Parent_SubJet_2_Truth.perp()) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
953  Splittings_LeadingSubJetpT_Truth.push_back(Parent_SubJet_1_Truth.perp());
954  vector < fastjet::PseudoJet > Hard_SubJet_Constituents_Truth = sorted_by_pt(Parent_SubJet_1_Truth.constituents());
955  Is_D_SubJet_Truth=kFALSE;
956  for(UInt_t j=0;j<Hard_SubJet_Constituents_Truth.size();j++){
957  if(Hard_SubJet_Constituents_Truth[j].user_index()==i_Matched_D_Jet_Truth) Is_D_SubJet_Truth=kTRUE;
958  }
959 
960 
961  if (!Is_D_SubJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(1.0);
962  else Splittings_HardestSubJetD0_Truth.push_back(2.0);
963 
964  // if(!Is_D_SubJet_Truth) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
965  Splittings_DeltaR_Truth.push_back(Parent_SubJet_1_Truth.delta_R(Parent_SubJet_2_Truth));
966  Splittings_Zg_Truth.push_back(Parent_SubJet_2_Truth.perp()/(Parent_SubJet_1_Truth.perp()+Parent_SubJet_2_Truth.perp()));
967  Splittings_RadiatorE_Truth.push_back(Daughter_Jet_Truth.E());
968  Splittings_RadiatorpT_Truth.push_back(Daughter_Jet_Truth.perp());
969  Daughter_Jet_Truth=Parent_SubJet_1_Truth;
970  }
971 
972 
973  } catch (fastjet::Error) { /*return -1;*/ }
974 
975 
976 
977 
978  //set detector level flags
979  Double_t Flag_D=-1.0;
980  if (Mass_Hypo_Type ==1) Flag_D=1.0;
981  else if (Mass_Hypo_Type ==2) Flag_D=2.0;
982  else if (Mass_Hypo_Type ==3 && Matched_Truth_Particle->GetPdgCode()==fCandidatePDG) Flag_D=3.0;
983  else if (Mass_Hypo_Type ==3 && Matched_Truth_Particle->GetPdgCode()==-fCandidatePDG) Flag_D=4.0;
984 
985  Double_t Flag_D_Truth=-1.0;
986  if (Matched_Truth_Particle->GetPdgCode()==fCandidatePDG) Flag_D_Truth=1.0;
987  if (Matched_Truth_Particle->GetPdgCode()==-fCandidatePDG) Flag_D_Truth=2.0;
988 
989 
990  fShapesVar[0] = Inclusive_Jets[i_Jet].perp();
991  fShapesVar[1] = Inclusive_Jets_Truth[i_Jet_Truth].perp();
992  fShapesVar[2] = D_Candidate->Pt();
993  fShapesVar[3] = Matched_Truth_Particle->Pt();
994  fShapesVar[4] = Inv_Mass_D;
995  fShapesVar[5] = Inv_Mass_D_Truth;
996  fShapesVar[6] = Flag_D;
997  fShapesVar[7] = Flag_D_Truth;
998  fShapesVar[8] = Is_Prompt_Correct_Quark_PDG;
999  fShapesVar[9] = 0.0;
1000 
1001 
1002 
1003  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR);
1004  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR_Truth);
1005  fShapesVar_Splittings_Zg.push_back(Splittings_Zg);
1006  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg_Truth);
1007  fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT);
1008  fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT_Truth);
1009  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0);
1010  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0_Truth);
1011  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE);
1012  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE_Truth);
1013  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT);
1014  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT_Truth);
1015 
1016 
1017  fTreeResponseMatrixAxis->Fill();
1018  fTreeSplittings->Fill();
1019 
1022  fShapesVar_Splittings_Zg.clear();
1032  }
1033  }
1034  // delete Matched_Truth_Particle;
1035  // delete D_Candidate;
1036  }
1037  // delete Track_Container;
1038  if(NMatched_DMeson_Jets==0) fhEvent->Fill(24);
1039  if(NMatched_DMeson_Jets==1) fhEvent->Fill(25);
1040  if(NMatched_DMeson_Jets==2) fhEvent->Fill(26);
1041  if(NMatched_DMeson_Jets==3) fhEvent->Fill(27);
1042  if(NMatched_DMeson_Jets==4) fhEvent->Fill(28);
1043  if(NMatched_DMeson_Jets==5) fhEvent->Fill(29);
1044  if(NMatched_DMeson_Jets==6) fhEvent->Fill(30);
1045 
1046  if (fIncludeInclusive){
1047 
1048  if(Unmatched_Truth_Level_D.size()==0) fhEvent->Fill(31);
1049  if(Unmatched_Truth_Level_D.size()==1) fhEvent->Fill(32);
1050  if(Unmatched_Truth_Level_D.size()==2) fhEvent->Fill(33);
1051  if(Unmatched_Truth_Level_D.size()==3) fhEvent->Fill(34);
1052  if(Unmatched_Truth_Level_D.size()==4) fhEvent->Fill(35);
1053  if(Unmatched_Truth_Level_D.size()==5) fhEvent->Fill(36);
1054  if(Unmatched_Truth_Level_D.size()==6) fhEvent->Fill(37);
1055 
1056  for (UInt_t i_Jet_Truth=0; i_Jet_Truth < Inclusive_Jets_Truth.size(); i_Jet_Truth++){
1057  AliAODMCParticle* Truth_D_Particle = NULL;
1058  Bool_t Is_Unmatched_D=kFALSE;
1059  Int_t D_Meson_Matched_Index=-1;
1060  if (TMath::Abs(Inclusive_Jets_Truth[i_Jet_Truth].pseudorapidity()) > 0.9-fJetRadius) continue;
1061  if (Inclusive_Jets_Truth[i_Jet_Truth].perp()<fJetMinPt) continue;
1062  std::vector<fastjet::PseudoJet> Constituents_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1063  for (UInt_t i_Constituents_Truth = 0; i_Constituents_Truth < Constituents_Truth.size(); i_Constituents_Truth++) {
1064  for (UInt_t i_Unmacthed_D=0; i_Unmacthed_D<Unmatched_Truth_Level_D.size(); i_Unmacthed_D++){
1065  if(Constituents_Truth[i_Constituents_Truth].user_index() == Unmatched_Truth_Level_D[i_Unmacthed_D]){
1066  Is_Unmatched_D=kTRUE;
1067  D_Meson_Matched_Index=Constituents_Truth[i_Constituents_Truth].user_index();
1068  for(UInt_t i_MC_Label=0; i_MC_Label<Inclusive_Jets_Truth_Labels.size(); i_MC_Label++){
1069  if (Inclusive_Jets_Truth_Labels[i_MC_Label].second==Constituents_Truth[i_Constituents_Truth].user_index()){
1070  Truth_D_Particle=static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Inclusive_Jets_Truth_Labels[i_MC_Label].first));
1071  }
1072  }
1073  }
1074  }
1075  }
1076  if (!Is_Unmatched_D) continue;
1077  fhEvent->Fill(38);
1078 
1079 
1080  Bool_t Is_Prompt_Correct_Quark=kFALSE;
1081  Int_t Is_Prompt_Correct_Quark_PDG=-1;
1082  Int_t Truth_D_Particle_Mother_Label=Truth_D_Particle->GetMother();
1083  while (Truth_D_Particle_Mother_Label >= 0) {
1084  AliAODMCParticle* Truth_D_Particle_Mother = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Truth_D_Particle_Mother_Label));
1085  if (Truth_D_Particle_Mother){
1086  Int_t Original_Quark_PDG=4;
1087  if (fIsBDecay) Original_Quark_PDG=5;
1088  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) == 4){
1089  fhEvent->Fill(39);
1090  Is_Prompt_Correct_Quark_PDG=4;
1091  }
1092  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) == 5){
1093  fhEvent->Fill(40);
1094  Is_Prompt_Correct_Quark_PDG=5;
1095  }
1096  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode())==Original_Quark_PDG) Is_Prompt_Correct_Quark=kTRUE;
1097  if (Truth_D_Particle_Mother_Label==Truth_D_Particle_Mother->GetMother()) break;
1098  Truth_D_Particle_Mother_Label=Truth_D_Particle_Mother->GetMother();
1099  }
1100  else break;
1101  // delete Truth_D_Particle_Mother;
1102  }
1103 
1104  if (fPromptReject && !Is_Prompt_Correct_Quark) continue;
1105 
1106 
1107  Bool_t Is_D_SubJet_Truth=kFALSE;
1108  fastjet::JetDefinition Jet_Definition_Truth(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
1109 
1110  std::vector<Double_t> Splittings_Zg_Truth;
1111  std::vector<Double_t> Splittings_DeltaR_Truth;
1112  std::vector<Double_t> Splittings_LeadingSubJetpT_Truth;
1113  std::vector<Double_t> Splittings_HardestSubJetD0_Truth;
1114  std::vector<Double_t> Splittings_RadiatorE_Truth;
1115  std::vector<Double_t> Splittings_RadiatorpT_Truth;
1116 
1117 
1118  try{
1119  std::vector<fastjet::PseudoJet> Reclustered_Particles_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1120  fastjet::ClusterSequence Cluster_Sequence_CA_Truth(Reclustered_Particles_Truth, Jet_Definition_Truth);
1121  std::vector<fastjet::PseudoJet> Reclustered_Jet_Truth = Cluster_Sequence_CA_Truth.inclusive_jets(0.0);
1122  Reclustered_Jet_Truth = sorted_by_pt(Reclustered_Jet_Truth);
1123 
1124 
1125  fastjet::PseudoJet Daughter_Jet_Truth = Reclustered_Jet_Truth[0];
1126  fastjet::PseudoJet Parent_SubJet_1_Truth;
1127  fastjet::PseudoJet Parent_SubJet_2_Truth;
1128 
1129  while(Daughter_Jet_Truth.has_parents(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth)){
1130  if(Parent_SubJet_1_Truth.perp() < Parent_SubJet_2_Truth.perp()) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1131  Splittings_LeadingSubJetpT_Truth.push_back(Parent_SubJet_1_Truth.perp());
1132  vector < fastjet::PseudoJet > Hard_SubJet_Constituents_Truth = sorted_by_pt(Parent_SubJet_1_Truth.constituents());
1133  Is_D_SubJet_Truth=kFALSE;
1134  for(UInt_t j=0;j<Hard_SubJet_Constituents_Truth.size();j++){
1135  if(Hard_SubJet_Constituents_Truth[j].user_index()==D_Meson_Matched_Index) Is_D_SubJet_Truth=kTRUE;
1136  }
1137 
1138 
1139  if (!Is_D_SubJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(1.0);
1140  else Splittings_HardestSubJetD0_Truth.push_back(2.0);
1141 
1142  // if(!Is_D_SubJet_Truth) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1143  Splittings_DeltaR_Truth.push_back(Parent_SubJet_1_Truth.delta_R(Parent_SubJet_2_Truth));
1144  Splittings_Zg_Truth.push_back(Parent_SubJet_2_Truth.perp()/(Parent_SubJet_1_Truth.perp()+Parent_SubJet_2_Truth.perp()));
1145  Splittings_RadiatorE_Truth.push_back(Daughter_Jet_Truth.E());
1146  Splittings_RadiatorpT_Truth.push_back(Daughter_Jet_Truth.perp());
1147  Daughter_Jet_Truth=Parent_SubJet_1_Truth;
1148  }
1149 
1150 
1151  } catch (fastjet::Error) { /*return -1;*/ }
1152 
1153  Double_t Inv_Mass_D_Truth=0.0;
1154  Double_t Flag_D_Truth=-1.0;
1155  if (Truth_D_Particle->GetPdgCode()==fCandidatePDG){
1156  //Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0();
1157  Inv_Mass_D_Truth=0.0;
1158  Flag_D_Truth=3.0;
1159  }
1160  if (Truth_D_Particle->GetPdgCode()==-fCandidatePDG){
1161  // Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0bar();
1162  Inv_Mass_D_Truth=0.0;
1163  Flag_D_Truth=4.0;
1164  }
1165 
1166 
1167  fShapesVar[0] = 0.0;
1168  fShapesVar[1] = Inclusive_Jets_Truth[i_Jet_Truth].perp();
1169  fShapesVar[2] = 0.0;
1170  fShapesVar[3] = Truth_D_Particle->Pt();
1171  fShapesVar[4] = 0.0;
1172  fShapesVar[5] = Inv_Mass_D_Truth;
1173  fShapesVar[6] = 0.0;
1174  fShapesVar[7] = Flag_D_Truth;
1175  fShapesVar[8] = Is_Prompt_Correct_Quark_PDG;
1176  fShapesVar[9] = 0.0;
1177 
1178 
1179  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR_Truth);
1180  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR_Truth);
1181  fShapesVar_Splittings_Zg.push_back(Splittings_Zg_Truth);
1182  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg_Truth);
1183  fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT_Truth);
1184  fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT_Truth);
1185  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0_Truth);
1186  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0_Truth);
1187  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE_Truth);
1188  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE_Truth);
1189  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT_Truth);
1190  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT_Truth);
1191 
1192 
1193  fTreeResponseMatrixAxis->Fill();
1194  fTreeSplittings->Fill();
1195 
1198  fShapesVar_Splittings_Zg.clear();
1208  // delete Truth_D_Particle;
1209  }
1210  }
1211  // delete Particle_Container;
1212  }
1213 
1214 
1215  if (fJetShapeType == kTrue){
1216 
1217 
1218  //truth level only jet finding
1219 
1220 
1222  Particle_Container->SetSpecialPDG(fCandidatePDG);
1223  Particle_Container->SetRejectedOriginMap(fRejectedOrigin);
1224  Particle_Container->SetAcceptedDecayMap(kDecayD0toKpi);
1225  Particle_Container->SetRejectISR(fRejectISR);
1226 
1227  std::vector<fastjet::PseudoJet> Inclusive_Jets_Truth;
1228  std::vector<std::pair<Int_t, Int_t>> Inclusive_Jets_Truth_Labels;
1229 
1230  if (Particle_Container->IsSpecialPDGFound()){
1231  fhEvent->Fill(2);
1232 
1233  fFastJetWrapper_Truth->SetAreaType(fastjet::active_area);
1236  fFastJetWrapper_Truth->SetAlgorithm(fastjet::antikt_algorithm);
1237  fFastJetWrapper_Truth->SetRecombScheme(static_cast<fastjet::RecombinationScheme>(0));
1239 
1240  AliAODMCParticle* Truth_Particle=NULL;
1241  Int_t NTruthD=0;
1242  for (Int_t i_Particle=0; i_Particle<Particle_Container->GetNParticles(); i_Particle++){
1243  Truth_Particle = static_cast<AliAODMCParticle*>(Particle_Container->GetAcceptParticle(i_Particle));
1244  if (!Truth_Particle) continue;
1245  if (TMath::Abs(Truth_Particle->GetPdgCode())==fCandidatePDG){
1246  fhEvent->Fill(3);
1247  std::pair<Int_t, Int_t> Inclusive_Jet_Truth_Labels;
1248  Inclusive_Jet_Truth_Labels.first=Truth_Particle->GetLabel();
1249  Inclusive_Jet_Truth_Labels.second=NTruthD;
1250  Inclusive_Jets_Truth_Labels.push_back(Inclusive_Jet_Truth_Labels);
1251  fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),NTruthD);
1252  NTruthD++;
1253  }
1254  else fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),i_Particle+100);
1255  }
1256  // delete Truth_Particle;
1258  Inclusive_Jets_Truth = fFastJetWrapper_Truth->GetInclusiveJets();
1259 
1260  for (UInt_t i_Jet_Truth=0; i_Jet_Truth < Inclusive_Jets_Truth.size(); i_Jet_Truth++){
1261  AliAODMCParticle* Truth_D_Particle = NULL;
1262  Bool_t Is_DJet_Truth=kFALSE;
1263  Int_t D_Meson_Matched_Index=-1;
1264  if (TMath::Abs(Inclusive_Jets_Truth[i_Jet_Truth].pseudorapidity()) > 0.9-fJetRadius) continue;
1265  std::vector<fastjet::PseudoJet> Constituents_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1266  for (UInt_t i_Constituents_Truth = 0; i_Constituents_Truth < Constituents_Truth.size(); i_Constituents_Truth++) {
1267 
1268  if(Constituents_Truth[i_Constituents_Truth].user_index() >=0 && Constituents_Truth[i_Constituents_Truth].user_index() < NTruthD){
1269  D_Meson_Matched_Index=Constituents_Truth[i_Constituents_Truth].user_index();
1270  Is_DJet_Truth=kTRUE;
1271  for(UInt_t i_MC_Label=0; i_MC_Label<Inclusive_Jets_Truth_Labels.size(); i_MC_Label++){
1272  if (Inclusive_Jets_Truth_Labels[i_MC_Label].second==Constituents_Truth[i_Constituents_Truth].user_index()){
1273  Truth_D_Particle=static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Inclusive_Jets_Truth_Labels[i_MC_Label].first));
1274  }
1275  }
1276 
1277  }
1278  }
1279 
1280  if (!Is_DJet_Truth) fhEvent->Fill(4);
1281  if (Is_DJet_Truth) fhEvent->Fill(5);
1282 
1283  if (!fIncludeInclusive && !Is_DJet_Truth) continue;
1284 
1285  Bool_t Is_Prompt_Correct_Quark=kFALSE;
1286  Int_t Is_Prompt_Correct_Quark_PDG=-1;
1287  Int_t Truth_D_Particle_Mother_Label=Truth_D_Particle->GetMother();
1288  while (Truth_D_Particle_Mother_Label >= 0) {
1289  AliAODMCParticle* Truth_D_Particle_Mother = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Truth_D_Particle_Mother_Label));
1290  if (Truth_D_Particle_Mother){
1291  Int_t Original_Quark_PDG=4;
1292  if (fIsBDecay) Original_Quark_PDG=5;
1293  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) ==Original_Quark_PDG) Is_Prompt_Correct_Quark=kTRUE;
1294  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) ==4){
1295  fhEvent->Fill(6);
1296  Is_Prompt_Correct_Quark_PDG=4;
1297  }
1298  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) ==5){
1299  fhEvent->Fill(7);
1300  Is_Prompt_Correct_Quark_PDG=5;
1301  }
1302  if (Truth_D_Particle_Mother_Label==Truth_D_Particle_Mother->GetMother()) break;
1303  Truth_D_Particle_Mother_Label=Truth_D_Particle_Mother->GetMother();
1304  }
1305  else break;
1306  // delete Truth_D_Particle_Mother;
1307  }
1308 
1309  if(fPromptReject && !Is_Prompt_Correct_Quark) continue;
1310  fhEvent->Fill(8);
1311 
1312  Bool_t Is_D_SubJet_Truth=kFALSE;
1313  fastjet::JetDefinition Jet_Definition_Truth(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
1314 
1315  std::vector<Double_t> Splittings_Zg_Truth;
1316  std::vector<Double_t> Splittings_DeltaR_Truth;
1317  std::vector<Double_t> Splittings_LeadingSubJetpT_Truth;
1318  std::vector<Double_t> Splittings_HardestSubJetD0_Truth;
1319  std::vector<Double_t> Splittings_RadiatorE_Truth;
1320  std::vector<Double_t> Splittings_RadiatorpT_Truth;
1321 
1322  try{
1323  std::vector<fastjet::PseudoJet> Reclustered_Particles_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1324  fastjet::ClusterSequence Cluster_Sequence_CA_Truth(Reclustered_Particles_Truth, Jet_Definition_Truth);
1325  std::vector<fastjet::PseudoJet> Reclustered_Jet_Truth = Cluster_Sequence_CA_Truth.inclusive_jets(0.0);
1326  Reclustered_Jet_Truth = sorted_by_pt(Reclustered_Jet_Truth);
1327 
1328 
1329  fastjet::PseudoJet Daughter_Jet_Truth = Reclustered_Jet_Truth[0];
1330  fastjet::PseudoJet Parent_SubJet_1_Truth;
1331  fastjet::PseudoJet Parent_SubJet_2_Truth;
1332 
1333  while(Daughter_Jet_Truth.has_parents(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth)){
1334  if(Parent_SubJet_1_Truth.perp() < Parent_SubJet_2_Truth.perp()) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1335  Splittings_LeadingSubJetpT_Truth.push_back(Parent_SubJet_1_Truth.perp());
1336  vector < fastjet::PseudoJet > Hard_SubJet_Constituents_Truth = sorted_by_pt(Parent_SubJet_1_Truth.constituents());
1337  Is_D_SubJet_Truth=kFALSE;
1338  for(UInt_t j=0;j<Hard_SubJet_Constituents_Truth.size();j++){
1339  if(Hard_SubJet_Constituents_Truth[j].user_index()==D_Meson_Matched_Index) Is_D_SubJet_Truth=kTRUE;
1340  }
1341 
1342  if (fIncludeInclusive && !Is_DJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(0.0);
1343  else if (Is_DJet_Truth && !Is_D_SubJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(1.0);
1344  else if (Is_DJet_Truth && Is_D_SubJet_Truth)Splittings_HardestSubJetD0_Truth.push_back(2.0);
1345 
1346  // if(!Is_D_SubJet_Truth) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1347  Splittings_DeltaR_Truth.push_back(Parent_SubJet_1_Truth.delta_R(Parent_SubJet_2_Truth));
1348  Splittings_Zg_Truth.push_back(Parent_SubJet_2_Truth.perp()/(Parent_SubJet_1_Truth.perp()+Parent_SubJet_2_Truth.perp()));
1349  Splittings_RadiatorE_Truth.push_back(Daughter_Jet_Truth.E());
1350  Splittings_RadiatorpT_Truth.push_back(Daughter_Jet_Truth.perp());
1351  Daughter_Jet_Truth=Parent_SubJet_1_Truth;
1352  }
1353 
1354 
1355  } catch (fastjet::Error) { /*return -1;*/ }
1356 
1357  Double_t Inv_Mass_D_Truth=0.0;
1358  Double_t Flag_D_Truth=-1.0;
1359  Double_t D_Pt=-1.0;
1360  if (fIncludeInclusive && !Is_DJet_Truth){
1361  Flag_D_Truth=0.0;
1362  }
1363  else if (Is_DJet_Truth){
1364  if (Truth_D_Particle->GetPdgCode()==fCandidatePDG){
1365  //Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0();
1366  Inv_Mass_D_Truth=0.0;
1367  Flag_D_Truth=1.0;
1368  D_Pt=Truth_D_Particle->Pt();
1369  fhEvent->Fill(9);
1370  }
1371  if (Truth_D_Particle->GetPdgCode()==-fCandidatePDG){
1372  // Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0bar();
1373  Inv_Mass_D_Truth=0.0;
1374  Flag_D_Truth=2.0;
1375  D_Pt=Truth_D_Particle->Pt();
1376  fhEvent->Fill(10);
1377  }
1378  }
1379 
1380 
1381  fShapesVar[0] = 0.0;
1382  fShapesVar[1] = Inclusive_Jets_Truth[i_Jet_Truth].perp();
1383  fShapesVar[2] = 0.0;
1384  fShapesVar[3] = D_Pt;
1385  fShapesVar[4] = 0.0;
1386  fShapesVar[5] = Inv_Mass_D_Truth;
1387  fShapesVar[6] = 0.0;
1388  fShapesVar[7] = Flag_D_Truth;
1389  fShapesVar[8] = 0.0;
1390  fShapesVar[9] = Is_Prompt_Correct_Quark_PDG;
1391 
1392 
1393 
1394  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR_Truth);
1395  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR_Truth);
1396  fShapesVar_Splittings_Zg.push_back(Splittings_Zg_Truth);
1397  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg_Truth);
1398  fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT_Truth);
1399  fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT_Truth);
1400  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0_Truth);
1401  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0_Truth);
1402  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE_Truth);
1403  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE_Truth);
1404  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT_Truth);
1405  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT_Truth);
1406 
1407 
1408  fTreeResponseMatrixAxis->Fill();
1409  fTreeSplittings->Fill();
1410 
1413  fShapesVar_Splittings_Zg.clear();
1423 
1424 
1425  // delete Truth_D_Particle;
1426  }
1427  if (NTruthD==0) fhEvent->Fill(11);
1428  if (NTruthD==1) fhEvent->Fill(12);
1429  if (NTruthD==2) fhEvent->Fill(13);
1430  if (NTruthD==3) fhEvent->Fill(14);
1431  if (NTruthD==4) fhEvent->Fill(15);
1432  if (NTruthD==5) fhEvent->Fill(16);
1433  if (NTruthD==6) fhEvent->Fill(17);
1434  }
1435  //delete Particle_Container;
1436  }
1437 
1438  delete fFastJetWrapper_Truth;
1439  delete fFastJetWrapper;
1440 
1441  return kTRUE;
1442 
1443 
1444 }
1445 
1446 
1447 
1448 
1449 
1450 
1451 //________________________________________________________________________
1453  //
1454  // retrieve event objects
1455  //
1457  return kFALSE;
1458 
1459  return kTRUE;
1460 }
1461 
1462 
1463 //_______________________________________________________________________
1465 {
1466  // Called once at the end of the analysis.
1467 
1468 }
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:293
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:314
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