AliPhysics  master (3d17d9d)
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  fTrackingEfficiency(1.0),
93  fCandidateArray(0),
94  fAodEvent(0),
95  fRDHFCuts(0),
96  fFastJetWrapper(0),
97  fFastJetWrapper_Truth(0),
98  fShapesVar_Splittings_DeltaR(0),
99  fShapesVar_Splittings_DeltaR_Truth(0),
100  fShapesVar_Splittings_Zg(0),
101  fShapesVar_Splittings_Zg_Truth(0),
102  fShapesVar_Splittings_LeadingSubJetpT(0),
103  fShapesVar_Splittings_LeadingSubJetpT_Truth(0),
104  fShapesVar_Splittings_HardestSubJetD0(0),
105  fShapesVar_Splittings_HardestSubJetD0_Truth(0),
106  fShapesVar_Splittings_RadiatorE(0),
107  fShapesVar_Splittings_RadiatorE_Truth(0),
108  fShapesVar_Splittings_RadiatorpT(0),
109  fShapesVar_Splittings_RadiatorpT_Truth(0),
110 
111  fTreeResponseMatrixAxis(0),
112  fTreeSplittings(0),
113  fhEvent(0x0)
114 
115 {
116  for(Int_t i=0;i<nVar;i++){
117  fShapesVar[i]=0;
118  }
119  SetMakeGeneralHistograms(kTRUE);
120  DefineOutput(1, TList::Class());
121  DefineOutput(2, TTree::Class());
122  DefineOutput(3, TTree::Class());
123 }
124 
125 //________________________________________________________________________
127  AliAnalysisTaskEmcal(name, kTRUE),
128  fJetShapeType(kData),
129  fECandidateType(kD0toKpi),
130  fIncludeInclusive(kFALSE),
131  fIsBDecay(kFALSE),
132  fRejectISR(kFALSE),
133  fPromptReject(kFALSE),
134  fAlienConnect(kFALSE),
135  fBranchName(0),
136  fCutsType(0),
137  fCandidatePDG(421),
138  fRejectedOrigin(0),
139  fAcceptedDecay(kDecayD0toKpi),
140  fJetRadius(0.4),
141  fJetMinPt(0.0),
142  fTrackingEfficiency(1.0),
143  fCandidateArray(0),
144  fAodEvent(0),
145  fRDHFCuts(0),
146  fFastJetWrapper(0),
147  fFastJetWrapper_Truth(0),
148  fShapesVar_Splittings_DeltaR(0),
149  fShapesVar_Splittings_DeltaR_Truth(0),
150  fShapesVar_Splittings_Zg(0),
151  fShapesVar_Splittings_Zg_Truth(0),
152  fShapesVar_Splittings_LeadingSubJetpT(0),
153  fShapesVar_Splittings_LeadingSubJetpT_Truth(0),
154  fShapesVar_Splittings_HardestSubJetD0(0),
155  fShapesVar_Splittings_HardestSubJetD0_Truth(0),
156  fShapesVar_Splittings_RadiatorE(0),
157  fShapesVar_Splittings_RadiatorE_Truth(0),
158  fShapesVar_Splittings_RadiatorpT(0),
159  fShapesVar_Splittings_RadiatorpT_Truth(0),
160 
161  fTreeResponseMatrixAxis(0),
162  fTreeSplittings(0),
163  fhEvent(0x0)
164 {
165  // Standard constructor.
166  for(Int_t i=0;i<nVar;i++){
167  fShapesVar[i]=0;
168  }
170  DefineOutput(1, TList::Class());
171  DefineOutput(2, TTree::Class());
172  DefineOutput(3, TTree::Class());
173 }
174 
175 //________________________________________________________________________
177 {
178  // Destructor.
179 }
180 
181 //________________________________________________________________________
183 {
184  // Create user output.
185 
187 
188  Bool_t oldStatus = TH1::AddDirectoryStatus();
189  TH1::AddDirectory(kFALSE);
190  TH1::AddDirectory(oldStatus);
191  //create a tree used for the MC data and making a 4D response matrix
192 
193 
194 
195  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
196  fTreeResponseMatrixAxis = new TTree(nameoutput, nameoutput);
197  TString *fShapesVarNames = new TString [nVar];
198  fShapesVarNames[0] = "pT_Jet";
199  fShapesVarNames[1] = "pT_Jet_Truth";
200  fShapesVarNames[2] = "pT_D";
201  fShapesVarNames[3] = "pT_D_Truth";
202  fShapesVarNames[4] = "Inv_M_D";
203  fShapesVarNames[5] = "Inv_M_D_Truth";
204  fShapesVarNames[6] = "Flag_D";
205  fShapesVarNames[7] = "Flag_D_Truth";
206  fShapesVarNames[8] = "Prompt_PDG";
207  fShapesVarNames[9] = "Prompt_PDG_Truth";
208  // fShapesVarNames[10] = "NTracks";
209  //fShapesVarNames[11] = "NTracks_Truth";
210  // fShapesVarNames[12] = "Eta_Jet";
211  // fShapesVarNames[13] = "Eta_Jet_Truth";
212  // fShapesVarNames[14] = "Eta_D";
213  // fShapesVarNames[15] = "Eta_D_Truth";
214  // fShapesVarNames[16] = "Y_D";
215  // fShapesVarNames[17] = "Y_D_Truth";
216 
217  for(Int_t ivar=0; ivar < nVar; ivar++){
218  cout<<"looping over variables"<<endl;
219  fTreeResponseMatrixAxis->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/D", fShapesVarNames[ivar].Data()));
220  }
221 
222  const char* nameoutput_Splittings = GetOutputSlot(3)->GetContainer()->GetName();
223  fTreeSplittings = new TTree(nameoutput_Splittings, nameoutput_Splittings);
224  TString *fShapesVarNames_Splittings=new TString[nVar_Splittings];
225  fShapesVarNames_Splittings[0] = "DeltaR";
226  fShapesVarNames_Splittings[1] = "DeltaR_Truth";
227  fShapesVarNames_Splittings[2] = "Zg";
228  fShapesVarNames_Splittings[3] = "Zg_Truth";
229  // fShapesVarNames_Splittings[4] = "LeadingSubJetpT";
230  //fShapesVarNames_Splittings[5] = "LeadingSubJetpT_Truth";
231  fShapesVarNames_Splittings[4] = "HardestSubJetD0";
232  fShapesVarNames_Splittings[5] = "HardestSubJetD0_Truth";
233  fShapesVarNames_Splittings[6] = "RadiatorE";
234  fShapesVarNames_Splittings[7] = "RadiatorE_Truth";
235  fShapesVarNames_Splittings[8] = "RadiatorpT";
236  fShapesVarNames_Splittings[9] = "RadiatorpT_Truth";
237  fTreeSplittings->Branch(fShapesVarNames_Splittings[0].Data(), &fShapesVar_Splittings_DeltaR, 0,1);
238  fTreeSplittings->Branch(fShapesVarNames_Splittings[1].Data(), &fShapesVar_Splittings_DeltaR_Truth, 0,1);
239  fTreeSplittings->Branch(fShapesVarNames_Splittings[2].Data(), &fShapesVar_Splittings_Zg, 0,1);
240  fTreeSplittings->Branch(fShapesVarNames_Splittings[3].Data(), &fShapesVar_Splittings_Zg_Truth, 0,1);
241  // fTreeSplittings->Branch(fShapesVarNames_Splittings[4].Data(), &fShapesVar_Splittings_LeadingSubJetpT, 0,1);
242  // fTreeSplittings->Branch(fShapesVarNames_Splittings[5].Data(), &fShapesVar_Splittings_LeadingSubJetpT_Truth, 0,1);
243  fTreeSplittings->Branch(fShapesVarNames_Splittings[4].Data(), &fShapesVar_Splittings_HardestSubJetD0, 0,1);
244  fTreeSplittings->Branch(fShapesVarNames_Splittings[5].Data(), &fShapesVar_Splittings_HardestSubJetD0_Truth, 0,1);
245  fTreeSplittings->Branch(fShapesVarNames_Splittings[6].Data(), &fShapesVar_Splittings_RadiatorE, 0,1);
246  fTreeSplittings->Branch(fShapesVarNames_Splittings[7].Data(), &fShapesVar_Splittings_RadiatorE_Truth, 0,1);
247  fTreeSplittings->Branch(fShapesVarNames_Splittings[8].Data(), &fShapesVar_Splittings_RadiatorpT, 0,1);
248  fTreeSplittings->Branch(fShapesVarNames_Splittings[9].Data(), &fShapesVar_Splittings_RadiatorpT_Truth, 0,1);
249 
250  fhEvent=new TH1D("fhEvent","fhEvent",40,-0.5,39.5);
251  fOutput->Add(fhEvent);
252 
253  PostData(1,fOutput);
254  PostData(2,fTreeResponseMatrixAxis);
255  PostData(3,fTreeSplittings);
256  // delete [] fShapesVarNames;
257 }
258 
259 //________________________________________________________________________
261 {
262  // Run analysis code here, if needed. It will be executed before FillHistograms().
263 
264 
265  return kTRUE;
266 }
267 //________________________________________________________________________
269 {
270 
271  Int_t Matching_AOD_deltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
272  if (Matching_AOD_deltaAODlevel <= 0) return kTRUE;
273 
274  fhEvent->Fill(0);
275  //if(fCutsFileName.Contains("alien://") && fAlienConnect) TGrid::Connect("alien://");
276  //TFile* Cuts_File = TFile::Open(fCutsFileName);
277  //TString cutsname="D0toKpiCuts";
278  //if (fCutsType!="") cutsname += TString::Format("_%s", fCutsType.Data());
279  //fRDHFCuts = dynamic_cast<AliRDHFCuts*>(fCutsFile->Get(cutsname));
280 
281 
282  fAodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
283 
284 
285  if(!fRDHFCuts->IsEventSelected(fAodEvent)) return kTRUE;
286  fhEvent->Fill(1);
287  fFastJetWrapper=new AliFJWrapper("fastjetwrapper","fastjetwrapper");
288  fFastJetWrapper->SetAreaType(fastjet::active_area);
289  fFastJetWrapper->SetGhostArea(0.005);
291  fFastJetWrapper->SetAlgorithm(fastjet::antikt_algorithm);
292  fFastJetWrapper->SetRecombScheme(static_cast<fastjet::RecombinationScheme>(0));
293 
294  fFastJetWrapper_Truth=new AliFJWrapper("fastjetwrapper_truth","fastjetwrapper_truth");
295 
297 
298  TRandom3 Random;
299  Random.SetSeed(0);
300  Double_t Random_Number;
301 
302  fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject("D0toKpi"));
303 
304  AliHFAODMCParticleContainer *Particle_Container = NULL;
305  if (fJetShapeType != kData) Particle_Container = (AliHFAODMCParticleContainer *) GetParticleContainer(1);
306  // if (!Particle_Container) continue;
307 
308 
309  std::vector<AliAODRecoDecayHF2Prong*> D_Candidates_Vector;
310  D_Candidates_Vector.clear();
311 
312 
313  Int_t N_DMesons=0;
314  //AliHFTrackContainer* Track_Container=(AliHFTrackContainer *) GetTrackContainer(0);
315  AliHFTrackContainer* Track_Container=dynamic_cast<AliHFTrackContainer*>(GetTrackContainer(0));
316  if (!Track_Container) return kTRUE;
317  Track_Container->SetDMesonCandidate(NULL);
318  for (Int_t i_D = 0; i_D < fCandidateArray->GetEntriesFast(); i_D++) {
319  AliAODRecoDecayHF2Prong* D_Candidate = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(i_D));
320  if (!D_Candidate) continue;
321  if (!fRDHFCuts->IsInFiducialAcceptance(D_Candidate->Pt(), D_Candidate->Y(fCandidatePDG))) continue;
322 
323 
324 
325  Int_t Mass_Hypo_Type=fRDHFCuts->IsSelected(D_Candidate, AliRDHFCuts::kAll, fAodEvent);
326  Int_t N_Mass_Hypotheses=1;
327  if (Mass_Hypo_Type <= 0 || Mass_Hypo_Type>3) continue;
328  else if (Mass_Hypo_Type ==3) N_Mass_Hypotheses=2;
329  fhEvent->Fill(2);
330  Int_t Matched_Truth_Particle_PDG=0;
331  Int_t Is_Prompt_Correct_Quark_PDG=-1;
332  if (fJetShapeType != kData){
333 
334 
335  const Int_t D_Candidtae_N_Daughters=2;
336  Int_t D_Candidtae_Daughters_PDG[D_Candidtae_N_Daughters] = {211,321};
337  Int_t D_Candidate_MatchedTruth_Label = D_Candidate->MatchToMC(fCandidatePDG, Particle_Container->GetArray(), D_Candidtae_N_Daughters, D_Candidtae_Daughters_PDG);
338  Bool_t Is_Prompt_Correct_Quark=kFALSE;
339  if (D_Candidate_MatchedTruth_Label >= 0) {
340  AliAODMCParticle* Matched_Truth_Particle = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(D_Candidate_MatchedTruth_Label));
341  if (Matched_Truth_Particle) {
342  fhEvent->Fill(3);
343 
344  Int_t Matched_Truth_Particle_Mother_Label=Matched_Truth_Particle->GetMother();
345  while (Matched_Truth_Particle_Mother_Label >= 0) {
346  AliAODMCParticle* Matched_Truth_Particle_Mother = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Matched_Truth_Particle_Mother_Label));
347  if (Matched_Truth_Particle_Mother){
348  Int_t Original_Quark_PDG=4;
349  if (fIsBDecay) Original_Quark_PDG=5;
350  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode())==Original_Quark_PDG) Is_Prompt_Correct_Quark=kTRUE;
351  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode()) == 4){
352  Is_Prompt_Correct_Quark_PDG=4;
353  fhEvent->Fill(4);
354  break;
355  }
356  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode()) == 5){
357  Is_Prompt_Correct_Quark_PDG=5;
358  fhEvent->Fill(5);
359  break;
360  }
361  if (Matched_Truth_Particle_Mother_Label==Matched_Truth_Particle_Mother->GetMother()) break;
362  Matched_Truth_Particle_Mother_Label=Matched_Truth_Particle_Mother->GetMother();
363  }
364  else break;
365  //delete Matched_Truth_Particle_Mother;
366  }
367  Matched_Truth_Particle_PDG = Matched_Truth_Particle->PdgCode();
368  }
369  //delete Matched_Truth_Particle;
370  }
371  //else continue;
372  if (fPromptReject && !Is_Prompt_Correct_Quark) continue;
373 
374  //if (TMath::Abs(Matched_Truth_Particle_PDG)!=fCandidatePDG) continue;
375 
376  }
377 
378 
379  Double_t Inv_Mass_D=0.0;
380  if (Mass_Hypo_Type==1){
381  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)){
382  Inv_Mass_D=D_Candidate->InvMassD0();
383  fhEvent->Fill(6);
384  }
385  else{
386  fhEvent->Fill(7);
387  continue;
388  }
389  }
390 
391  if (Mass_Hypo_Type==2){
392  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)){
393  Inv_Mass_D=D_Candidate->InvMassD0bar();
394  fhEvent->Fill(8);
395  }
396  else{
397  fhEvent->Fill(9);
398  continue;
399  }
400  }
401 
402 
403 
404  for (Int_t i_Mass_Hypotheses=0; i_Mass_Hypotheses<N_Mass_Hypotheses; i_Mass_Hypotheses++){
405 
406 
407  if (Mass_Hypo_Type==3){
408  if(i_Mass_Hypotheses==0){
409  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)){
410  Inv_Mass_D=D_Candidate->InvMassD0();
411  fhEvent->Fill(11);
412  }
413  else{
414  fhEvent->Fill(12);
415  continue;
416  }
417  }
418  if(i_Mass_Hypotheses==1){
419  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)){
420  Inv_Mass_D=D_Candidate->InvMassD0bar();
421  fhEvent->Fill(13);
422  }
423  else{
424  fhEvent->Fill(14);
425  continue;
426  }
427  }
428  }
429  //Random_Number=Random.Rndm();
430  //if(Random_Number > fTrackingEfficiency*fTrackingEfficiency) continue; // here it shows that the D did not get reconstructed cause one of the daughters was missing...however should we do this before incase the same daughter is involved multiple times?
432  AliTLorentzVector D_Candidate_LorentzVector(0,0,0,0);
433  D_Candidate_LorentzVector.SetPtEtaPhiM(D_Candidate->Pt(), D_Candidate->Eta(), D_Candidate->Phi(), Inv_Mass_D);
434  fFastJetWrapper->AddInputVector(D_Candidate_LorentzVector.Px(), D_Candidate_LorentzVector.Py(), D_Candidate_LorentzVector.Pz(), D_Candidate_LorentzVector.E(), 0);
435  Track_Container->SetDMesonCandidate(D_Candidate);
436  AliAODTrack *Track = NULL;
437  for (Int_t i_Track=0; i_Track<Track_Container->GetNTracks(); i_Track++){
438  Track = static_cast<AliAODTrack*>(Track_Container->GetAcceptParticle(i_Track));
439  if(!Track) continue;
440  if (Track->Pt() > 100.0 || TMath::Abs(Track->Eta()) > 0.9) continue;
441  Random_Number=Random.Rndm();
442  if(Random_Number > fTrackingEfficiency) continue;
443  fFastJetWrapper->AddInputVector(Track->Px(), Track->Py(), Track->Pz(), Track->E(),i_Track+100);
444  }
445  // delete Track;
446 
447  fFastJetWrapper->Run();
448  std::vector<fastjet::PseudoJet> Inclusive_Jets = fFastJetWrapper->GetInclusiveJets();
449  for (UInt_t i_Jet=0; i_Jet < Inclusive_Jets.size(); i_Jet++){
450  Bool_t Is_D_Jet=kFALSE;
451  if (Inclusive_Jets[i_Jet].perp()<fJetMinPt) continue;
452  if (TMath::Abs(Inclusive_Jets[i_Jet].pseudorapidity()) > 0.9-fJetRadius) continue;
453  std::vector<fastjet::PseudoJet> Constituents(fFastJetWrapper->GetJetConstituents(i_Jet));
454  for (UInt_t i_Constituents = 0; i_Constituents < Constituents.size(); i_Constituents++) {
455  if (Constituents[i_Constituents].user_index() == 0) Is_D_Jet = kTRUE;
456  }
457  if (!Is_D_Jet) continue;
458  fhEvent->Fill(10);
459  std::vector<Double_t> Splittings_Zg;
460  std::vector<Double_t> Splittings_DeltaR;
461  std::vector<Double_t> Splittings_LeadingSubJetpT;
462  std::vector<Double_t> Splittings_HardestSubJetD0;
463  std::vector<Double_t> Splittings_RadiatorE;
464  std::vector<Double_t> Splittings_RadiatorpT;
465 
466  Bool_t Is_D_SubJet=kFALSE;
467  fastjet::JetDefinition Jet_Definition(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
468  try{
469  std::vector<fastjet::PseudoJet> Reclustered_Particles(fFastJetWrapper->GetJetConstituents(i_Jet));
470  fastjet::ClusterSequence Cluster_Sequence_CA(Reclustered_Particles, Jet_Definition);
471  std::vector<fastjet::PseudoJet> Reclustered_Jet = Cluster_Sequence_CA.inclusive_jets(0.0);
472  Reclustered_Jet = sorted_by_pt(Reclustered_Jet);
473 
474 
475  fastjet::PseudoJet Daughter_Jet = Reclustered_Jet[0];
476  fastjet::PseudoJet Parent_SubJet_1;
477  fastjet::PseudoJet Parent_SubJet_2;
478 
479  while(Daughter_Jet.has_parents(Parent_SubJet_1,Parent_SubJet_2)){
480  if(Parent_SubJet_1.perp() < Parent_SubJet_2.perp()) std::swap(Parent_SubJet_1,Parent_SubJet_2);
481  // Splittings_LeadingSubJetpT.push_back(Parent_SubJet_1.perp());
482  vector < fastjet::PseudoJet > Hard_SubJet_Constituents = sorted_by_pt(Parent_SubJet_1.constituents());
483  Is_D_SubJet=kFALSE;
484  for(UInt_t j=0;j<Hard_SubJet_Constituents.size();j++){
485  if(Hard_SubJet_Constituents[j].user_index()==0) Is_D_SubJet=kTRUE;
486  }
487 
488  if (!Is_D_SubJet) Splittings_HardestSubJetD0.push_back(1.0);
489  else Splittings_HardestSubJetD0.push_back(2.0);
490 
491  // if(!Is_D_SubJet) std::swap(Parent_SubJet_1,Parent_SubJet_2);
492  Splittings_DeltaR.push_back(Parent_SubJet_1.delta_R(Parent_SubJet_2));
493  Splittings_Zg.push_back(Parent_SubJet_2.perp()/(Parent_SubJet_1.perp()+Parent_SubJet_2.perp()));
494  Splittings_RadiatorE.push_back(Daughter_Jet.E());
495  Splittings_RadiatorpT.push_back(Daughter_Jet.perp());
496  Daughter_Jet=Parent_SubJet_1;
497  }
498 
499 
500  } catch (fastjet::Error) { /*return -1;*/ }
501 
502 
503 
504  // for (Int_t i_Mass_Hypotheses=0; i_Mass_Hypotheses<N_Mass_Hypotheses; i_Mass_Hypotheses++){
505 
506 
507  D_Candidates_Vector.push_back(D_Candidate);
508  if (i_Mass_Hypotheses==0) fhEvent->Fill(15);
509  N_DMesons++;
510 
511  Double_t Flag_D=-1.0;
512  if (Mass_Hypo_Type ==1) Flag_D=1.0;
513  else if (Mass_Hypo_Type ==2) Flag_D=2.0;
514  else if (Mass_Hypo_Type ==3 && i_Mass_Hypotheses==0) Flag_D=3.0;
515  else if (Mass_Hypo_Type ==3 && i_Mass_Hypotheses==1) Flag_D=4.0;
516 
517 
518  fShapesVar[0] = Inclusive_Jets[i_Jet].perp();
519  fShapesVar[1] = 0.0;
520  fShapesVar[2] = D_Candidate->Pt();
521  fShapesVar[3] = 0.0;
522  fShapesVar[4] = Inv_Mass_D;
523  fShapesVar[5] = 0.0;
524  fShapesVar[6] = Flag_D;
525  if (fJetShapeType == kData) fShapesVar[7] = 0.0;
526  else fShapesVar[7] = Matched_Truth_Particle_PDG;
527  fShapesVar[8] = Is_Prompt_Correct_Quark_PDG;
528  fShapesVar[9] = 0.0;
529  // fShapesVar[10] = NTracks;
530  //fShapesVar[11] = 0.0;
531  // fShapesVar[12] = TMath::Abs(Inclusive_Jets[i_Jet].pseudorapidity());
532  // fShapesVar[13] = 0.0;
533  // fShapesVar[14] = Dmeson_Eta;
534  //fShapesVar[15] = 0.0;
535  //fShapesVar[16] = Dmeson_Y;
536  //fShapesVar[17] = 0.0;
537 
538 
539 
540  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR);
541  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR);
542  fShapesVar_Splittings_Zg.push_back(Splittings_Zg);
543  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg);
544  // fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT);
545  //fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT);
546  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0);
547  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0);
548  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE);
549  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE);
550  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT);
551  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT);
552 
553  fTreeResponseMatrixAxis->Fill();
554  fTreeSplittings->Fill();
557  fShapesVar_Splittings_Zg.clear();
559  // fShapesVar_Splittings_LeadingSubJetpT.clear();
560  //fShapesVar_Splittings_LeadingSubJetpT_Truth.clear();
567  }
568  }
569  //delete D_Candidate;
570  }
571  if(N_DMesons==0) fhEvent->Fill(16);
572  if(N_DMesons==1) fhEvent->Fill(17);
573  if(N_DMesons==2) fhEvent->Fill(18);
574  if(N_DMesons==3) fhEvent->Fill(19);
575  if(N_DMesons==4) fhEvent->Fill(20);
576  if(N_DMesons==5) fhEvent->Fill(21);
577  if(N_DMesons==6) fhEvent->Fill(22);
578 
579  // delete Track_Container;
580  // if (fJetShapeType != kData) delete Particle_Container;
581 
582  }
583 
584 
585 
586  if (fJetShapeType == kTrueDet){
587 
588  TRandom3 Random;
589  Random.SetSeed(0);
590  Double_t Random_Number;
591 
593  Particle_Container->SetSpecialPDG(fCandidatePDG);
594  Particle_Container->SetRejectedOriginMap(fRejectedOrigin);
595  Particle_Container->SetAcceptedDecayMap(fAcceptedDecay);
596  Particle_Container->SetRejectISR(fRejectISR);
597  Particle_Container->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
598 
599  std::vector<fastjet::PseudoJet> Inclusive_Jets_Truth;
600  std::vector<std::pair<Int_t, Int_t>> Inclusive_Jets_Truth_Labels;
601  std::vector<Int_t> Unmatched_Truth_Level_D;
602  Int_t NMatched_DMeson_Jets=0;
603 
604  // Double_t NTracks=0;
605  // Double_t NTracks_Truth=0;
606  //Double_t Jet_Eta=0;
607  //Double_t Jet_Eta_Truth=0;
608  //Double_t Dmeson_Eta=0;
609  //Double_t Dmeson_Eta_Truth=0;
610  //Double_t Dmeson_Y=0;
611  //Double_t Dmeson_Y_Truth=0;
612 
613 
614 
615  if (Particle_Container->IsSpecialPDGFound()){
616 
617  fhEvent->Fill(2);
618 
619  fFastJetWrapper_Truth->SetAreaType(fastjet::active_area);
622  fFastJetWrapper_Truth->SetAlgorithm(fastjet::antikt_algorithm);
623  fFastJetWrapper_Truth->SetRecombScheme(static_cast<fastjet::RecombinationScheme>(0));
625 
626  AliAODMCParticle* Truth_Particle=NULL;
627  Int_t NTruthD=0;
628  for (Int_t i_Particle=0; i_Particle<Particle_Container->GetNParticles(); i_Particle++){
629  Truth_Particle = static_cast<AliAODMCParticle*>(Particle_Container->GetAcceptMCParticle(i_Particle));
630  if (!Truth_Particle) continue;
631  // if (TMath::Abs(Truth_Particle->Eta())>0.9) continue;
632  if (TMath::Abs(Truth_Particle->PdgCode())==fCandidatePDG){
633  if (Truth_Particle->Pt() > 5.0){
634  if (TMath::Abs(Truth_Particle->Y()) > 0.8) continue;
635  }
636  else{
637  if(Truth_Particle->Y() < 0.2/15*Truth_Particle->Pt()*Truth_Particle->Pt()-1.9/15*Truth_Particle->Pt()-0.5 || Truth_Particle->Y() > -0.2/15*Truth_Particle->Pt()*Truth_Particle->Pt()+1.9/15*Truth_Particle->Pt()+0.5) continue;
638  }
639  std::pair<Int_t, Int_t> Inclusive_Jet_Truth_Labels;
640  Inclusive_Jet_Truth_Labels.first=Truth_Particle->GetLabel();
641  Inclusive_Jet_Truth_Labels.second=NTruthD;
642  Inclusive_Jets_Truth_Labels.push_back(Inclusive_Jet_Truth_Labels);
643  Unmatched_Truth_Level_D.push_back(NTruthD);
644  fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),NTruthD);
645  NTruthD++;
646  }
647  else{
648  if (TMath::Abs(Truth_Particle->Eta()) > 0.9) continue;
649  fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),i_Particle+100);
650  }
651  }
652  // delete Truth_Particle;
654  Inclusive_Jets_Truth = fFastJetWrapper_Truth->GetInclusiveJets();
655  if (NTruthD==0) fhEvent->Fill(3);
656  if (NTruthD==1) fhEvent->Fill(4);
657  if (NTruthD==2) fhEvent->Fill(5);
658  if (NTruthD==3) fhEvent->Fill(6);
659  if (NTruthD==4) fhEvent->Fill(7);
660  if (NTruthD==5) fhEvent->Fill(8);
661  if (NTruthD==6) fhEvent->Fill(9);
662  }
663 
664 
665 
666 
667 
668  fCandidateArray = dynamic_cast<TClonesArray*>(fAodEvent->GetList()->FindObject("D0toKpi"));
670  if (!Track_Container) return kTRUE;
671  Track_Container->SetDMesonCandidate(NULL);
672 
673  for (Int_t i_D = 0; i_D < fCandidateArray->GetEntriesFast(); i_D++) {
674  AliAODRecoDecayHF2Prong* D_Candidate = static_cast<AliAODRecoDecayHF2Prong*>(fCandidateArray->At(i_D));
675  if (!D_Candidate) continue;
676  if (!fRDHFCuts->IsInFiducialAcceptance(D_Candidate->Pt(), D_Candidate->Y(fCandidatePDG))) continue;
677 
678  Int_t Mass_Hypo_Type=fRDHFCuts->IsSelected(D_Candidate, AliRDHFCuts::kAll, fAodEvent);
679  if (Mass_Hypo_Type <= 0 || Mass_Hypo_Type>3) continue;
680  fhEvent->Fill(10);
681 
682  Int_t Matched_Truth_Particle_PDG=0;
683 
684 
685  const Int_t D_Candidtae_N_Daughters=2;
686  Int_t D_Candidtae_Daughters_PDG[D_Candidtae_N_Daughters] = {211,321};
687  Int_t D_Candidate_MatchedTruth_Label = D_Candidate->MatchToMC(fCandidatePDG, Particle_Container->GetArray(), D_Candidtae_N_Daughters, D_Candidtae_Daughters_PDG);
688  Bool_t Is_Prompt_Correct_Quark=kFALSE;
689  Int_t Is_Prompt_Correct_Quark_PDG=-1;
690  AliAODMCParticle* Matched_Truth_Particle;
691  if (D_Candidate_MatchedTruth_Label >= 0) {
692  Matched_Truth_Particle = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(D_Candidate_MatchedTruth_Label));
693  if (Matched_Truth_Particle) {
694  fhEvent->Fill(11);
695 
696  Int_t Matched_Truth_Particle_Mother_Label=Matched_Truth_Particle->GetMother();
697  while (Matched_Truth_Particle_Mother_Label >= 0) {
698  AliAODMCParticle* Matched_Truth_Particle_Mother = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Matched_Truth_Particle_Mother_Label));
699  if (Matched_Truth_Particle_Mother){
700  Int_t Original_Quark_PDG=4;
701  if (fIsBDecay) Original_Quark_PDG=5;
702  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode()) == 4){
703  fhEvent->Fill(12);
704  Is_Prompt_Correct_Quark_PDG=4;
705  break;
706  }
707  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode()) == 5){
708  fhEvent->Fill(13);
709  Is_Prompt_Correct_Quark_PDG=5;
710  break;
711  }
712  if (TMath::Abs(Matched_Truth_Particle_Mother->GetPdgCode())==Original_Quark_PDG) Is_Prompt_Correct_Quark=kTRUE;
713  if (Matched_Truth_Particle_Mother_Label==Matched_Truth_Particle_Mother->GetMother()) break;
714  Matched_Truth_Particle_Mother_Label=Matched_Truth_Particle_Mother->GetMother();
715  }
716  else break;
717  // delete Matched_Truth_Particle_Mother;
718  }
719  Matched_Truth_Particle_PDG = Matched_Truth_Particle->PdgCode();
720  }
721  }
722  else continue;
723  if (fPromptReject && !Is_Prompt_Correct_Quark) continue;
724 
725 
726  if (Mass_Hypo_Type==1 && Matched_Truth_Particle_PDG!=fCandidatePDG){
727  fhEvent->Fill(14);
728  continue;
729  }
730  else fhEvent->Fill(15);
731  if (Mass_Hypo_Type==2 && Matched_Truth_Particle_PDG!=-fCandidatePDG){
732  fhEvent->Fill(16);
733  continue;
734  }
735  else fhEvent->Fill(17);
736  if (Mass_Hypo_Type==3 && TMath::Abs(Matched_Truth_Particle_PDG)!=fCandidatePDG){
737  fhEvent->Fill(18);
738  continue;
739  }
740  else{
741  fhEvent->Fill(19);
742  if (Matched_Truth_Particle_PDG==fCandidatePDG) fhEvent->Fill(20);
743  if (Matched_Truth_Particle_PDG==-fCandidatePDG) fhEvent->Fill(21);
744  }
745 
746 
747  Double_t Inv_Mass_D=0.0;
748  Double_t Inv_Mass_D_Truth=0.0;
749 
750  if (Mass_Hypo_Type==1 || (Mass_Hypo_Type==3 && Matched_Truth_Particle_PDG==fCandidatePDG)){
751  Inv_Mass_D=D_Candidate->InvMassD0();
752  //Inv_Mass_D_Truth=Matched_Truth_Particle->InvMassD0();
753  Inv_Mass_D_Truth=0.0;
754  }
755  if (Mass_Hypo_Type==2 || (Mass_Hypo_Type==3 && Matched_Truth_Particle_PDG==-fCandidatePDG)){
756  Inv_Mass_D=D_Candidate->InvMassD0bar();
757  //Inv_Mass_D_Truth=Matched_Truth_Particle->InvMassD0bar();
758  Inv_Mass_D_Truth=0.0;
759  }
760 
761  //Random_Number=Random.Rndm();
762  //if(Random_Number > fTrackingEfficiency*fTrackingEfficiency) continue;
763 
765  AliTLorentzVector D_Candidate_LorentzVector(0,0,0,0);
766  D_Candidate_LorentzVector.SetPtEtaPhiM(D_Candidate->Pt(), D_Candidate->Eta(), D_Candidate->Phi(), Inv_Mass_D);
767  //if (TMath::Abs(D_Candidate->Eta())>0.9) continue;
768  // Dmeson_Eta=TMath::Abs(D_Candidate->Eta());
769  //Dmeson_Y=TMath::Abs(D_Candidate->Y(fCandidatePDG));
770  fFastJetWrapper->AddInputVector(D_Candidate_LorentzVector.Px(), D_Candidate_LorentzVector.Py(), D_Candidate_LorentzVector.Pz(), D_Candidate_LorentzVector.E(), 0);
771 
772 
773  if (!Track_Container) continue;
774  Track_Container->SetDMesonCandidate(D_Candidate);
775  AliAODTrack *Track = NULL;
776  for (Int_t i_Track=0; i_Track<Track_Container->GetNTracks(); i_Track++){
777  Track = static_cast<AliAODTrack*>(Track_Container->GetAcceptParticle(i_Track));
778  if(!Track) continue;
779  if (Track->Pt() > 100.0 || TMath::Abs(Track->Eta())>0.9) continue;
780  Random_Number=Random.Rndm();
781  if(Random_Number > fTrackingEfficiency) continue;
782  fFastJetWrapper->AddInputVector(Track->Px(), Track->Py(), Track->Pz(), Track->E(),i_Track+100);
783  }
784  // delete Track;
785  fFastJetWrapper->Run();
786 
787 
788  std::vector<fastjet::PseudoJet> Inclusive_Jets = fFastJetWrapper->GetInclusiveJets();
789  for (UInt_t i_Jet=0; i_Jet < Inclusive_Jets.size(); i_Jet++){
790  Bool_t Is_D_Jet=kFALSE;
791  if (Inclusive_Jets[i_Jet].perp()<fJetMinPt) continue;
792  //Jet_Eta=TMath::Abs(Inclusive_Jets[i_Jet].pseudorapidity());
793  if (TMath::Abs(Inclusive_Jets[i_Jet].pseudorapidity()) > 0.9-fJetRadius) continue;
794  std::vector<fastjet::PseudoJet> Constituents(fFastJetWrapper->GetJetConstituents(i_Jet));
795  // NTracks=Constituents.size();
796  for (UInt_t i_Constituents = 0; i_Constituents < Constituents.size(); i_Constituents++) {
797  if (Constituents[i_Constituents].user_index() == 0) {
798  Is_D_Jet = kTRUE;
799  }
800  }
801 
802 
803  if (!Is_D_Jet) continue;
804  fhEvent->Fill(22);
805 
806 
807  Int_t i_Matched_D_Jet_Truth=-1;
808  for (UInt_t k=0; k< Inclusive_Jets_Truth_Labels.size(); k++){
809  if(Inclusive_Jets_Truth_Labels[k].first==D_Candidate_MatchedTruth_Label) i_Matched_D_Jet_Truth=Inclusive_Jets_Truth_Labels[k].second;
810  }
811 
812  for (UInt_t i_Jet_Truth=0; i_Jet_Truth < Inclusive_Jets_Truth.size(); i_Jet_Truth++){
813  Bool_t Is_Jet_Truth_Matched=kFALSE;
814  if (TMath::Abs(Inclusive_Jets_Truth[i_Jet_Truth].pseudorapidity()) > 0.9-fJetRadius) continue;
815  //Jet_Eta_Truth=TMath::Abs(Inclusive_Jets_Truth[i_Jet_Truth].pseudorapidity());
816  std::vector<fastjet::PseudoJet> Constituents_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
817  // NTracks_Truth=Constituents_Truth.size();
818  for (UInt_t i_Constituents_Truth = 0; i_Constituents_Truth < Constituents_Truth.size(); i_Constituents_Truth++) {
819  if (Constituents_Truth[i_Constituents_Truth].user_index() == i_Matched_D_Jet_Truth) {
820  Is_Jet_Truth_Matched=kTRUE;
821  // Dmeson_Eta_Truth=TMath::Abs(Constituents_Truth[i_Constituents_Truth].pseudorapidity());
822  //Dmeson_Y_Truth=TMath::Abs(Constituents_Truth[i_Constituents_Truth].rapidity());
823  for (UInt_t i_Unmacthed_D=0; i_Unmacthed_D<Unmatched_Truth_Level_D.size(); i_Unmacthed_D++){
824  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);
825  }
826  }
827  }
828  if (!Is_Jet_Truth_Matched) continue;
829  fhEvent->Fill(23);
830  NMatched_DMeson_Jets++;
831 
832 
833 
834 
835  std::vector<Double_t> Splittings_Zg;
836  std::vector<Double_t> Splittings_DeltaR;
837  // std::vector<Double_t> Splittings_LeadingSubJetpT;
838  std::vector<Double_t> Splittings_HardestSubJetD0;
839  std::vector<Double_t> Splittings_RadiatorE;
840  std::vector<Double_t> Splittings_RadiatorpT;
841 
842  Bool_t Is_D_SubJet=kFALSE;
843  fastjet::JetDefinition Jet_Definition(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
844 
845  try{
846  std::vector<fastjet::PseudoJet> Reclustered_Particles(fFastJetWrapper->GetJetConstituents(i_Jet));
847  fastjet::ClusterSequence Cluster_Sequence_CA(Reclustered_Particles, Jet_Definition);
848  std::vector<fastjet::PseudoJet> Reclustered_Jet = Cluster_Sequence_CA.inclusive_jets(0.0);
849  Reclustered_Jet = sorted_by_pt(Reclustered_Jet);
850 
851 
852  fastjet::PseudoJet Daughter_Jet = Reclustered_Jet[0];
853  fastjet::PseudoJet Parent_SubJet_1;
854  fastjet::PseudoJet Parent_SubJet_2;
855 
856  while(Daughter_Jet.has_parents(Parent_SubJet_1,Parent_SubJet_2)){
857  if(Parent_SubJet_1.perp() < Parent_SubJet_2.perp()) std::swap(Parent_SubJet_1,Parent_SubJet_2);
858  // Splittings_LeadingSubJetpT.push_back(Parent_SubJet_1.perp());
859  vector < fastjet::PseudoJet > Hard_SubJet_Constituents = sorted_by_pt(Parent_SubJet_1.constituents());
860  Is_D_SubJet=kFALSE;
861  for(UInt_t j=0;j<Hard_SubJet_Constituents.size();j++){
862  if(Hard_SubJet_Constituents[j].user_index()==0) Is_D_SubJet=kTRUE;
863  }
864 
865 
866  if (!Is_D_SubJet) Splittings_HardestSubJetD0.push_back(1.0);
867  else Splittings_HardestSubJetD0.push_back(2.0);
868 
869  // if(!Is_D_SubJet) std::swap(Parent_SubJet_1,Parent_SubJet_2);
870  Splittings_DeltaR.push_back(Parent_SubJet_1.delta_R(Parent_SubJet_2));
871  Splittings_Zg.push_back(Parent_SubJet_2.perp()/(Parent_SubJet_1.perp()+Parent_SubJet_2.perp()));
872  Splittings_RadiatorE.push_back(Daughter_Jet.E());
873  Splittings_RadiatorpT.push_back(Daughter_Jet.perp());
874  Daughter_Jet=Parent_SubJet_1;
875  }
876 
877 
878  } catch (fastjet::Error) { /*return -1;*/ }
879 
880  std::vector<Double_t> Splittings_Zg_Truth;
881  std::vector<Double_t> Splittings_DeltaR_Truth;
882  // std::vector<Double_t> Splittings_LeadingSubJetpT_Truth;
883  std::vector<Double_t> Splittings_HardestSubJetD0_Truth;
884  std::vector<Double_t> Splittings_RadiatorE_Truth;
885  std::vector<Double_t> Splittings_RadiatorpT_Truth;
886 
887 
888  Bool_t Is_D_SubJet_Truth=kFALSE;
889  fastjet::JetDefinition Jet_Definition_Truth(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
890 
891  try{
892  std::vector<fastjet::PseudoJet> Reclustered_Particles_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
893  fastjet::ClusterSequence Cluster_Sequence_CA_Truth(Reclustered_Particles_Truth, Jet_Definition_Truth);
894  std::vector<fastjet::PseudoJet> Reclustered_Jet_Truth = Cluster_Sequence_CA_Truth.inclusive_jets(0.0);
895  Reclustered_Jet_Truth = sorted_by_pt(Reclustered_Jet_Truth);
896 
897 
898  fastjet::PseudoJet Daughter_Jet_Truth = Reclustered_Jet_Truth[0];
899  fastjet::PseudoJet Parent_SubJet_1_Truth;
900  fastjet::PseudoJet Parent_SubJet_2_Truth;
901 
902  while(Daughter_Jet_Truth.has_parents(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth)){
903  if(Parent_SubJet_1_Truth.perp() < Parent_SubJet_2_Truth.perp()) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
904  // Splittings_LeadingSubJetpT_Truth.push_back(Parent_SubJet_1_Truth.perp());
905  vector < fastjet::PseudoJet > Hard_SubJet_Constituents_Truth = sorted_by_pt(Parent_SubJet_1_Truth.constituents());
906  Is_D_SubJet_Truth=kFALSE;
907  for(UInt_t j=0;j<Hard_SubJet_Constituents_Truth.size();j++){
908  if(Hard_SubJet_Constituents_Truth[j].user_index()==i_Matched_D_Jet_Truth) Is_D_SubJet_Truth=kTRUE;
909  }
910 
911 
912  if (!Is_D_SubJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(1.0);
913  else Splittings_HardestSubJetD0_Truth.push_back(2.0);
914 
915  // if(!Is_D_SubJet_Truth) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
916  Splittings_DeltaR_Truth.push_back(Parent_SubJet_1_Truth.delta_R(Parent_SubJet_2_Truth));
917  Splittings_Zg_Truth.push_back(Parent_SubJet_2_Truth.perp()/(Parent_SubJet_1_Truth.perp()+Parent_SubJet_2_Truth.perp()));
918  Splittings_RadiatorE_Truth.push_back(Daughter_Jet_Truth.E());
919  Splittings_RadiatorpT_Truth.push_back(Daughter_Jet_Truth.perp());
920  Daughter_Jet_Truth=Parent_SubJet_1_Truth;
921  }
922 
923 
924  } catch (fastjet::Error) { /*return -1;*/ }
925 
926 
927  //set detector level flags
928  Double_t Flag_D=-1.0;
929  if (Mass_Hypo_Type ==1) Flag_D=1.0;
930  else if (Mass_Hypo_Type ==2) Flag_D=2.0;
931  else if (Mass_Hypo_Type ==3 && Matched_Truth_Particle->GetPdgCode()==fCandidatePDG) Flag_D=3.0;
932  else if (Mass_Hypo_Type ==3 && Matched_Truth_Particle->GetPdgCode()==-fCandidatePDG) Flag_D=4.0;
933 
934  Double_t Flag_D_Truth=-1.0;
935  if (Matched_Truth_Particle->GetPdgCode()==fCandidatePDG) Flag_D_Truth=1.0;
936  if (Matched_Truth_Particle->GetPdgCode()==-fCandidatePDG) Flag_D_Truth=2.0;
937 
938 
939  fShapesVar[0] = Inclusive_Jets[i_Jet].perp();
940  fShapesVar[1] = Inclusive_Jets_Truth[i_Jet_Truth].perp();
941  fShapesVar[2] = D_Candidate->Pt();
942  fShapesVar[3] = Matched_Truth_Particle->Pt();
943  fShapesVar[4] = Inv_Mass_D;
944  fShapesVar[5] = Inv_Mass_D_Truth;
945  fShapesVar[6] = Flag_D;
946  fShapesVar[7] = Flag_D_Truth;
947  fShapesVar[8] = Is_Prompt_Correct_Quark_PDG;
948  fShapesVar[9] = 0.0;
949  // fShapesVar[10] = NTracks;
950  //fShapesVar[11] = NTracks_Truth;
951  //fShapesVar[12] = Jet_Eta;
952  //fShapesVar[13] = Jet_Eta_Truth;
953  //fShapesVar[14] = Dmeson_Eta;
954  //fShapesVar[15] = Dmeson_Eta_Truth;
955  //fShapesVar[16] = Dmeson_Y;
956  //fShapesVar[17] = Dmeson_Y_Truth;
957 
958 
959 
960  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR);
961  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR_Truth);
962  fShapesVar_Splittings_Zg.push_back(Splittings_Zg);
963  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg_Truth);
964  // fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT);
965  //fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT_Truth);
966  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0);
967  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0_Truth);
968  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE);
969  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE_Truth);
970  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT);
971  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT_Truth);
972 
973 
974  fTreeResponseMatrixAxis->Fill();
975  fTreeSplittings->Fill();
976 
979  fShapesVar_Splittings_Zg.clear();
981  // fShapesVar_Splittings_LeadingSubJetpT.clear();
982  //fShapesVar_Splittings_LeadingSubJetpT_Truth.clear();
989  }
990  }
991  // delete Matched_Truth_Particle;
992  // delete D_Candidate;
993  }
994  // delete Track_Container;
995  if(NMatched_DMeson_Jets==0) fhEvent->Fill(24);
996  if(NMatched_DMeson_Jets==1) fhEvent->Fill(25);
997  if(NMatched_DMeson_Jets==2) fhEvent->Fill(26);
998  if(NMatched_DMeson_Jets==3) fhEvent->Fill(27);
999  if(NMatched_DMeson_Jets==4) fhEvent->Fill(28);
1000  if(NMatched_DMeson_Jets==5) fhEvent->Fill(29);
1001  if(NMatched_DMeson_Jets==6) fhEvent->Fill(30);
1002 
1003  if (fIncludeInclusive){
1004 
1005  if(Unmatched_Truth_Level_D.size()==0) fhEvent->Fill(31);
1006  if(Unmatched_Truth_Level_D.size()==1) fhEvent->Fill(32);
1007  if(Unmatched_Truth_Level_D.size()==2) fhEvent->Fill(33);
1008  if(Unmatched_Truth_Level_D.size()==3) fhEvent->Fill(34);
1009  if(Unmatched_Truth_Level_D.size()==4) fhEvent->Fill(35);
1010  if(Unmatched_Truth_Level_D.size()==5) fhEvent->Fill(36);
1011  if(Unmatched_Truth_Level_D.size()==6) fhEvent->Fill(37);
1012 
1013  for (UInt_t i_Jet_Truth=0; i_Jet_Truth < Inclusive_Jets_Truth.size(); i_Jet_Truth++){
1014  AliAODMCParticle* Truth_D_Particle = NULL;
1015  Bool_t Is_Unmatched_D=kFALSE;
1016  Int_t D_Meson_Matched_Index=-1;
1017  if (TMath::Abs(Inclusive_Jets_Truth[i_Jet_Truth].pseudorapidity()) > 0.9-fJetRadius) continue;
1018  //Jet_Eta_Truth=TMath::Abs(Inclusive_Jets_Truth[i_Jet_Truth].pseudorapidity());
1019  if (Inclusive_Jets_Truth[i_Jet_Truth].perp()<fJetMinPt) continue;
1020  std::vector<fastjet::PseudoJet> Constituents_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1021  //NTracks_Truth=Constituents_Truth.size();
1022  for (UInt_t i_Constituents_Truth = 0; i_Constituents_Truth < Constituents_Truth.size(); i_Constituents_Truth++) {
1023  for (UInt_t i_Unmacthed_D=0; i_Unmacthed_D<Unmatched_Truth_Level_D.size(); i_Unmacthed_D++){
1024  if(Constituents_Truth[i_Constituents_Truth].user_index() == Unmatched_Truth_Level_D[i_Unmacthed_D]){
1025  Is_Unmatched_D=kTRUE;
1026  D_Meson_Matched_Index=Constituents_Truth[i_Constituents_Truth].user_index();
1027  for(UInt_t i_MC_Label=0; i_MC_Label<Inclusive_Jets_Truth_Labels.size(); i_MC_Label++){
1028  if (Inclusive_Jets_Truth_Labels[i_MC_Label].second==Constituents_Truth[i_Constituents_Truth].user_index()){
1029  Truth_D_Particle=static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Inclusive_Jets_Truth_Labels[i_MC_Label].first));
1030  // Dmeson_Eta_Truth=TMath::Abs(Truth_D_Particle->Eta());
1031  //Dmeson_Y_Truth=TMath::Abs(Truth_D_Particle->Y());
1032  }
1033  }
1034  }
1035  }
1036  }
1037  if (!Is_Unmatched_D) continue;
1038  fhEvent->Fill(38);
1039 
1040 
1041  Bool_t Is_Prompt_Correct_Quark=kFALSE;
1042  Int_t Is_Prompt_Correct_Quark_PDG=-1;
1043  Int_t Truth_D_Particle_Mother_Label=Truth_D_Particle->GetMother();
1044  while (Truth_D_Particle_Mother_Label >= 0) {
1045  AliAODMCParticle* Truth_D_Particle_Mother = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Truth_D_Particle_Mother_Label));
1046  if (Truth_D_Particle_Mother){
1047  Int_t Original_Quark_PDG=4;
1048  if (fIsBDecay) Original_Quark_PDG=5;
1049  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) == 4){
1050  fhEvent->Fill(39);
1051  Is_Prompt_Correct_Quark_PDG=4;
1052  break;
1053  }
1054  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) == 5){
1055  fhEvent->Fill(40);
1056  Is_Prompt_Correct_Quark_PDG=5;
1057  break;
1058  }
1059  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode())==Original_Quark_PDG) Is_Prompt_Correct_Quark=kTRUE;
1060  if (Truth_D_Particle_Mother_Label==Truth_D_Particle_Mother->GetMother()) break;
1061  Truth_D_Particle_Mother_Label=Truth_D_Particle_Mother->GetMother();
1062  }
1063  else break;
1064  // delete Truth_D_Particle_Mother;
1065  }
1066 
1067  if (fPromptReject && !Is_Prompt_Correct_Quark) continue;
1068 
1069 
1070  Bool_t Is_D_SubJet_Truth=kFALSE;
1071  fastjet::JetDefinition Jet_Definition_Truth(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
1072 
1073  std::vector<Double_t> Splittings_Zg_Truth;
1074  std::vector<Double_t> Splittings_DeltaR_Truth;
1075  // std::vector<Double_t> Splittings_LeadingSubJetpT_Truth;
1076  std::vector<Double_t> Splittings_HardestSubJetD0_Truth;
1077  std::vector<Double_t> Splittings_RadiatorE_Truth;
1078  std::vector<Double_t> Splittings_RadiatorpT_Truth;
1079 
1080 
1081  try{
1082  std::vector<fastjet::PseudoJet> Reclustered_Particles_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1083  fastjet::ClusterSequence Cluster_Sequence_CA_Truth(Reclustered_Particles_Truth, Jet_Definition_Truth);
1084  std::vector<fastjet::PseudoJet> Reclustered_Jet_Truth = Cluster_Sequence_CA_Truth.inclusive_jets(0.0);
1085  Reclustered_Jet_Truth = sorted_by_pt(Reclustered_Jet_Truth);
1086 
1087 
1088  fastjet::PseudoJet Daughter_Jet_Truth = Reclustered_Jet_Truth[0];
1089  fastjet::PseudoJet Parent_SubJet_1_Truth;
1090  fastjet::PseudoJet Parent_SubJet_2_Truth;
1091 
1092  while(Daughter_Jet_Truth.has_parents(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth)){
1093  if(Parent_SubJet_1_Truth.perp() < Parent_SubJet_2_Truth.perp()) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1094  // Splittings_LeadingSubJetpT_Truth.push_back(Parent_SubJet_1_Truth.perp());
1095  vector < fastjet::PseudoJet > Hard_SubJet_Constituents_Truth = sorted_by_pt(Parent_SubJet_1_Truth.constituents());
1096  Is_D_SubJet_Truth=kFALSE;
1097  for(UInt_t j=0;j<Hard_SubJet_Constituents_Truth.size();j++){
1098  if(Hard_SubJet_Constituents_Truth[j].user_index()==D_Meson_Matched_Index) Is_D_SubJet_Truth=kTRUE;
1099  }
1100 
1101 
1102  if (!Is_D_SubJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(1.0);
1103  else Splittings_HardestSubJetD0_Truth.push_back(2.0);
1104 
1105  // if(!Is_D_SubJet_Truth) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1106  Splittings_DeltaR_Truth.push_back(Parent_SubJet_1_Truth.delta_R(Parent_SubJet_2_Truth));
1107  Splittings_Zg_Truth.push_back(Parent_SubJet_2_Truth.perp()/(Parent_SubJet_1_Truth.perp()+Parent_SubJet_2_Truth.perp()));
1108  Splittings_RadiatorE_Truth.push_back(Daughter_Jet_Truth.E());
1109  Splittings_RadiatorpT_Truth.push_back(Daughter_Jet_Truth.perp());
1110  Daughter_Jet_Truth=Parent_SubJet_1_Truth;
1111  }
1112 
1113 
1114  } catch (fastjet::Error) { /*return -1;*/ }
1115 
1116  Double_t Inv_Mass_D_Truth=0.0;
1117  Double_t Flag_D_Truth=-1.0;
1118  if (Truth_D_Particle->GetPdgCode()==fCandidatePDG){
1119  //Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0();
1120  Inv_Mass_D_Truth=0.0;
1121  Flag_D_Truth=3.0;
1122  }
1123  if (Truth_D_Particle->GetPdgCode()==-fCandidatePDG){
1124  // Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0bar();
1125  Inv_Mass_D_Truth=0.0;
1126  Flag_D_Truth=4.0;
1127  }
1128 
1129 
1130  fShapesVar[0] = 0.0;
1131  fShapesVar[1] = Inclusive_Jets_Truth[i_Jet_Truth].perp();
1132  fShapesVar[2] = 0.0;
1133  fShapesVar[3] = Truth_D_Particle->Pt();
1134  fShapesVar[4] = 0.0;
1135  fShapesVar[5] = Inv_Mass_D_Truth;
1136  fShapesVar[6] = 0.0;
1137  fShapesVar[7] = Flag_D_Truth;
1138  fShapesVar[8] = Is_Prompt_Correct_Quark_PDG;
1139  fShapesVar[9] = 0.0;
1140  // fShapesVar[10] = 0.0;
1141  // fShapesVar[11] = NTracks_Truth;
1142  // fShapesVar[12] = 0.0;
1143  // fShapesVar[13] = Jet_Eta_Truth;
1144  // fShapesVar[14] = 0.0;
1145  // fShapesVar[15] = Dmeson_Eta_Truth;
1146  // fShapesVar[16] = 0.0;
1147  // fShapesVar[17] = Dmeson_Y_Truth;
1148 
1149  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR_Truth);
1150  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR_Truth);
1151  fShapesVar_Splittings_Zg.push_back(Splittings_Zg_Truth);
1152  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg_Truth);
1153  // fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT_Truth);
1154  // fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT_Truth);
1155  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0_Truth);
1156  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0_Truth);
1157  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE_Truth);
1158  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE_Truth);
1159  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT_Truth);
1160  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT_Truth);
1161 
1162 
1163  fTreeResponseMatrixAxis->Fill();
1164  fTreeSplittings->Fill();
1165 
1168  fShapesVar_Splittings_Zg.clear();
1170  // fShapesVar_Splittings_LeadingSubJetpT.clear();
1171  // fShapesVar_Splittings_LeadingSubJetpT_Truth.clear();
1178  // delete Truth_D_Particle;
1179  }
1180  }
1181  // delete Particle_Container;
1182  }
1183 
1184 
1185  if (fJetShapeType == kTrue){
1186 
1187 
1188  //truth level only jet finding
1189 
1190  Double_t NTracks_Truth=0;
1191  Double_t Jet_Eta_Truth=-5.0;
1192  Double_t Dmeson_Eta_Truth=-5.0;
1193  Double_t Dmeson_Y_Truth=-5.0;
1194 
1196  Particle_Container->SetSpecialPDG(fCandidatePDG);
1197  Particle_Container->SetRejectedOriginMap(fRejectedOrigin);
1198  Particle_Container->SetAcceptedDecayMap(fAcceptedDecay);
1199  Particle_Container->SetRejectISR(fRejectISR);
1200  Particle_Container->SetCharge(AliParticleContainer::EChargeCut_t::kCharged);
1201 
1202  std::vector<fastjet::PseudoJet> Inclusive_Jets_Truth;
1203  std::vector<std::pair<Int_t, Int_t>> Inclusive_Jets_Truth_Labels;
1204 
1205  if (Particle_Container->IsSpecialPDGFound()){
1206  fhEvent->Fill(2);
1207 
1208  fFastJetWrapper_Truth->SetAreaType(fastjet::active_area);
1211  fFastJetWrapper_Truth->SetAlgorithm(fastjet::antikt_algorithm);
1212  fFastJetWrapper_Truth->SetRecombScheme(static_cast<fastjet::RecombinationScheme>(1));
1214 
1215  AliAODMCParticle* Truth_Particle=NULL;
1216  Int_t NTruthD=0;
1217  for (Int_t i_Particle=0; i_Particle<Particle_Container->GetNParticles(); i_Particle++){
1218  Truth_Particle = static_cast<AliAODMCParticle*>(Particle_Container->GetAcceptMCParticle(i_Particle));
1219  if (!Truth_Particle) continue;
1220  if (TMath::Abs(Truth_Particle->GetPdgCode())==fCandidatePDG){
1221  fhEvent->Fill(3);
1222  std::pair<Int_t, Int_t> Inclusive_Jet_Truth_Labels;
1223  Inclusive_Jet_Truth_Labels.first=Truth_Particle->GetLabel();
1224  Inclusive_Jet_Truth_Labels.second=NTruthD;
1225  Inclusive_Jets_Truth_Labels.push_back(Inclusive_Jet_Truth_Labels);
1226  fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),NTruthD);
1227  NTruthD++;
1228  }
1229  else fFastJetWrapper_Truth->AddInputVector(Truth_Particle->Px(), Truth_Particle->Py(), Truth_Particle->Pz(), Truth_Particle->E(),i_Particle+100);
1230  }
1231  // delete Truth_Particle;
1233  Inclusive_Jets_Truth = fFastJetWrapper_Truth->GetInclusiveJets();
1234 
1235  for (UInt_t i_Jet_Truth=0; i_Jet_Truth < Inclusive_Jets_Truth.size(); i_Jet_Truth++){
1236  AliAODMCParticle* Truth_D_Particle = NULL;
1237  Bool_t Is_DJet_Truth=kFALSE;
1238  Int_t D_Meson_Matched_Index=-1;
1239  //if (TMath::Abs(Inclusive_Jets_Truth[i_Jet_Truth].pseudorapidity()) > 0.9-fJetRadius) continue;
1240  Jet_Eta_Truth=TMath::Abs(Inclusive_Jets_Truth[i_Jet_Truth].pseudorapidity());
1241  std::vector<fastjet::PseudoJet> Constituents_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1242  NTracks_Truth=Constituents_Truth.size();
1243  for (UInt_t i_Constituents_Truth = 0; i_Constituents_Truth < Constituents_Truth.size(); i_Constituents_Truth++) {
1244 
1245  if(Constituents_Truth[i_Constituents_Truth].user_index() >=0 && Constituents_Truth[i_Constituents_Truth].user_index() < NTruthD){
1246  D_Meson_Matched_Index=Constituents_Truth[i_Constituents_Truth].user_index();
1247  Is_DJet_Truth=kTRUE;
1248  for(UInt_t i_MC_Label=0; i_MC_Label<Inclusive_Jets_Truth_Labels.size(); i_MC_Label++){
1249  if (Inclusive_Jets_Truth_Labels[i_MC_Label].second==Constituents_Truth[i_Constituents_Truth].user_index()){
1250  Truth_D_Particle=static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Inclusive_Jets_Truth_Labels[i_MC_Label].first));
1251  Dmeson_Eta_Truth=TMath::Abs(Truth_D_Particle->Eta());
1252  Dmeson_Y_Truth=TMath::Abs(Truth_D_Particle->Y());
1253  }
1254  }
1255 
1256  }
1257  }
1258 
1259  if (!Is_DJet_Truth) fhEvent->Fill(4);
1260  if (Is_DJet_Truth) fhEvent->Fill(5);
1261 
1262  if (!fIncludeInclusive && !Is_DJet_Truth) continue;
1263 
1264  Bool_t Is_Prompt_Correct_Quark=kFALSE;
1265  Int_t Is_Prompt_Correct_Quark_PDG=-1;
1266  Int_t Truth_D_Particle_Mother_Label=Truth_D_Particle->GetMother();
1267  while (Truth_D_Particle_Mother_Label >= 0) {
1268  AliAODMCParticle* Truth_D_Particle_Mother = static_cast<AliAODMCParticle*>(Particle_Container->GetArray()->At(Truth_D_Particle_Mother_Label));
1269  if (Truth_D_Particle_Mother){
1270  Int_t Original_Quark_PDG=4;
1271  if (fIsBDecay) Original_Quark_PDG=5;
1272  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) ==Original_Quark_PDG) Is_Prompt_Correct_Quark=kTRUE;
1273  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) ==4){
1274  fhEvent->Fill(6);
1275  Is_Prompt_Correct_Quark_PDG=4;
1276  break;
1277  }
1278  if (TMath::Abs(Truth_D_Particle_Mother->GetPdgCode()) ==5){
1279  fhEvent->Fill(7);
1280  Is_Prompt_Correct_Quark_PDG=5;
1281  break;
1282  }
1283  if (Truth_D_Particle_Mother_Label==Truth_D_Particle_Mother->GetMother()) break;
1284  Truth_D_Particle_Mother_Label=Truth_D_Particle_Mother->GetMother();
1285  }
1286  else break;
1287  // delete Truth_D_Particle_Mother;
1288  }
1289 
1290  if(fPromptReject && !Is_Prompt_Correct_Quark) continue;
1291  fhEvent->Fill(8);
1292 
1293  Bool_t Is_D_SubJet_Truth=kFALSE;
1294  fastjet::JetDefinition Jet_Definition_Truth(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(1), fastjet::Best);
1295 
1296  std::vector<Double_t> Splittings_Zg_Truth;
1297  std::vector<Double_t> Splittings_DeltaR_Truth;
1298  // std::vector<Double_t> Splittings_LeadingSubJetpT_Truth;
1299  std::vector<Double_t> Splittings_HardestSubJetD0_Truth;
1300  std::vector<Double_t> Splittings_RadiatorE_Truth;
1301  std::vector<Double_t> Splittings_RadiatorpT_Truth;
1302 
1303  try{
1304  std::vector<fastjet::PseudoJet> Reclustered_Particles_Truth(fFastJetWrapper_Truth->GetJetConstituents(i_Jet_Truth));
1305  fastjet::ClusterSequence Cluster_Sequence_CA_Truth(Reclustered_Particles_Truth, Jet_Definition_Truth);
1306  std::vector<fastjet::PseudoJet> Reclustered_Jet_Truth = Cluster_Sequence_CA_Truth.inclusive_jets(0.0);
1307  Reclustered_Jet_Truth = sorted_by_pt(Reclustered_Jet_Truth);
1308 
1309 
1310  fastjet::PseudoJet Daughter_Jet_Truth = Reclustered_Jet_Truth[0];
1311  fastjet::PseudoJet Parent_SubJet_1_Truth;
1312  fastjet::PseudoJet Parent_SubJet_2_Truth;
1313 
1314  while(Daughter_Jet_Truth.has_parents(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth)){
1315  if(Parent_SubJet_1_Truth.perp() < Parent_SubJet_2_Truth.perp()) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1316  // Splittings_LeadingSubJetpT_Truth.push_back(Parent_SubJet_1_Truth.perp());
1317  vector < fastjet::PseudoJet > Hard_SubJet_Constituents_Truth = sorted_by_pt(Parent_SubJet_1_Truth.constituents());
1318  Is_D_SubJet_Truth=kFALSE;
1319  for(UInt_t j=0;j<Hard_SubJet_Constituents_Truth.size();j++){
1320  if(Hard_SubJet_Constituents_Truth[j].user_index()==D_Meson_Matched_Index) Is_D_SubJet_Truth=kTRUE;
1321  }
1322 
1323  if (fIncludeInclusive && !Is_DJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(0.0);
1324  else if (Is_DJet_Truth && !Is_D_SubJet_Truth) Splittings_HardestSubJetD0_Truth.push_back(1.0);
1325  else if (Is_DJet_Truth && Is_D_SubJet_Truth)Splittings_HardestSubJetD0_Truth.push_back(2.0);
1326 
1327  // if(!Is_D_SubJet_Truth) std::swap(Parent_SubJet_1_Truth,Parent_SubJet_2_Truth);
1328  Splittings_DeltaR_Truth.push_back(Parent_SubJet_1_Truth.delta_R(Parent_SubJet_2_Truth));
1329  Splittings_Zg_Truth.push_back(Parent_SubJet_2_Truth.perp()/(Parent_SubJet_1_Truth.perp()+Parent_SubJet_2_Truth.perp()));
1330  Splittings_RadiatorE_Truth.push_back(Daughter_Jet_Truth.E());
1331  Splittings_RadiatorpT_Truth.push_back(Daughter_Jet_Truth.perp());
1332  Daughter_Jet_Truth=Parent_SubJet_1_Truth;
1333  }
1334 
1335 
1336  } catch (fastjet::Error) { /*return -1;*/ }
1337 
1338  Double_t Inv_Mass_D_Truth=0.0;
1339  Double_t Flag_D_Truth=-1.0;
1340  Double_t D_Pt=-1.0;
1341  if (fIncludeInclusive && !Is_DJet_Truth){
1342  Flag_D_Truth=0.0;
1343  }
1344  else if (Is_DJet_Truth){
1345  if (Truth_D_Particle->GetPdgCode()==fCandidatePDG){
1346  //Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0();
1347  Inv_Mass_D_Truth=0.0;
1348  Flag_D_Truth=1.0;
1349  D_Pt=Truth_D_Particle->Pt();
1350  fhEvent->Fill(9);
1351  }
1352  if (Truth_D_Particle->GetPdgCode()==-fCandidatePDG){
1353  // Inv_Mass_D_Truth=Truth_D_Particle->InvMassD0bar();
1354  Inv_Mass_D_Truth=0.0;
1355  Flag_D_Truth=2.0;
1356  D_Pt=Truth_D_Particle->Pt();
1357  fhEvent->Fill(10);
1358  }
1359  }
1360 
1361 
1362  fShapesVar[0] = 0.0;
1363  fShapesVar[1] = Inclusive_Jets_Truth[i_Jet_Truth].perp();
1364  fShapesVar[2] = 0.0;
1365  fShapesVar[3] = D_Pt;
1366  fShapesVar[4] = 0.0;
1367  fShapesVar[5] = Inv_Mass_D_Truth;
1368  fShapesVar[6] = 0.0;
1369  fShapesVar[7] = Flag_D_Truth;
1370  fShapesVar[8] = 0.0;
1371  fShapesVar[9] = Is_Prompt_Correct_Quark_PDG;
1372  // fShapesVar[10] = 0.0;
1373  //fShapesVar[11] = NTracks_Truth;
1374  //fShapesVar[12] = 0.0;
1375  //fShapesVar[13] = Jet_Eta_Truth;
1376  //fShapesVar[14] = 0.0;
1377  // fShapesVar[15] = Dmeson_Eta_Truth;
1378  // fShapesVar[16] = 0.0;
1379  // fShapesVar[17] = Dmeson_Y_Truth;
1380 
1381 
1382 
1383  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR_Truth);
1384  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR_Truth);
1385  fShapesVar_Splittings_Zg.push_back(Splittings_Zg_Truth);
1386  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg_Truth);
1387  // fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT_Truth);
1388  //fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT_Truth);
1389  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0_Truth);
1390  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0_Truth);
1391  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE_Truth);
1392  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE_Truth);
1393  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT_Truth);
1394  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT_Truth);
1395 
1396 
1397  fTreeResponseMatrixAxis->Fill();
1398  fTreeSplittings->Fill();
1399 
1402  fShapesVar_Splittings_Zg.clear();
1404  // fShapesVar_Splittings_LeadingSubJetpT.clear();
1405  //fShapesVar_Splittings_LeadingSubJetpT_Truth.clear();
1412 
1413 
1414  // delete Truth_D_Particle;
1415  }
1416  if (NTruthD==0) fhEvent->Fill(11);
1417  if (NTruthD==1) fhEvent->Fill(12);
1418  if (NTruthD==2) fhEvent->Fill(13);
1419  if (NTruthD==3) fhEvent->Fill(14);
1420  if (NTruthD==4) fhEvent->Fill(15);
1421  if (NTruthD==5) fhEvent->Fill(16);
1422  if (NTruthD==6) fhEvent->Fill(17);
1423  }
1424  //delete Particle_Container;
1425  }
1426 
1427 
1428 
1429 
1430  if (fJetShapeType == kDataInclusive){
1431 
1432  TRandom3 Random;
1433  Random.SetSeed(0);
1434  Double_t Random_Number;
1435 
1436  AliHFTrackContainer* Track_Container=dynamic_cast<AliHFTrackContainer*>(GetTrackContainer(0));
1437  if (!Track_Container) return kTRUE;
1438  //Track_Container->SetDMesonCandidate(NULL);
1440  // Double_t NTracks=0;
1441  //Double_t Jet_Eta=-5.0;
1442  //Double_t HardestTrack_Eta=-5.0;
1443  //Double_t HardestTrack_Y=-5.0;
1444  Double_t HardestTrack_Pt=-5.0;
1445  AliAODTrack *Track = NULL;
1446  for (Int_t i_Track=0; i_Track<Track_Container->GetNTracks(); i_Track++){
1447  Track = static_cast<AliAODTrack*>(Track_Container->GetAcceptParticle(i_Track));
1448  if(!Track) continue;
1449  if(Track->Pt() > 100.0 || TMath::Abs(Track->Eta()) > 0.9) continue;
1450  Random_Number=Random.Rndm();
1451  if(Random_Number > fTrackingEfficiency) continue;
1452  fFastJetWrapper->AddInputVector(Track->Px(), Track->Py(), Track->Pz(), Track->E(),i_Track+100);
1453  }
1454  fFastJetWrapper->Run();
1455  std::vector<fastjet::PseudoJet> Inclusive_Jets = fFastJetWrapper->GetInclusiveJets();
1456  for (UInt_t i_Jet=0; i_Jet < Inclusive_Jets.size(); i_Jet++){
1457  if (Inclusive_Jets[i_Jet].perp()<fJetMinPt) continue;
1458  if (TMath::Abs(Inclusive_Jets[i_Jet].pseudorapidity()) > 0.9-fJetRadius) continue;
1459  // Jet_Eta=TMath::Abs(Inclusive_Jets[i_Jet].pseudorapidity());
1460  fhEvent->Fill(23);
1461 
1462  HardestTrack_Pt=-5.0;
1463  std::vector<fastjet::PseudoJet> Constituents(fFastJetWrapper->GetJetConstituents(i_Jet));
1464  // NTracks=Constituents.size();
1465  for (UInt_t i_Constituents = 0; i_Constituents < Constituents.size(); i_Constituents++) {
1466  if (Constituents[i_Constituents].perp() > HardestTrack_Pt){
1467  HardestTrack_Pt=Constituents[i_Constituents].perp();
1468  //HardestTrack_Eta=TMath::Abs(Constituents[i_Constituents].pseudorapidity());
1469  //HardestTrack_Y=TMath::Abs(Constituents[i_Constituents].rapidity());
1470  }
1471  }
1472 
1473 
1474  std::vector<Double_t> Splittings_Zg;
1475  std::vector<Double_t> Splittings_DeltaR;
1476  // std::vector<Double_t> Splittings_LeadingSubJetpT;
1477  std::vector<Double_t> Splittings_HardestSubJetD0;
1478  std::vector<Double_t> Splittings_RadiatorE;
1479  std::vector<Double_t> Splittings_RadiatorpT;
1480 
1481  fastjet::JetDefinition Jet_Definition(fastjet::cambridge_algorithm, fJetRadius*2.5,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
1482 
1483  try{
1484  std::vector<fastjet::PseudoJet> Reclustered_Particles(fFastJetWrapper->GetJetConstituents(i_Jet));
1485  fastjet::ClusterSequence Cluster_Sequence_CA(Reclustered_Particles, Jet_Definition);
1486  std::vector<fastjet::PseudoJet> Reclustered_Jet = Cluster_Sequence_CA.inclusive_jets(0.0);
1487  Reclustered_Jet = sorted_by_pt(Reclustered_Jet);
1488 
1489 
1490  fastjet::PseudoJet Daughter_Jet = Reclustered_Jet[0];
1491  fastjet::PseudoJet Parent_SubJet_1;
1492  fastjet::PseudoJet Parent_SubJet_2;
1493 
1494  while(Daughter_Jet.has_parents(Parent_SubJet_1,Parent_SubJet_2)){
1495  if(Parent_SubJet_1.perp() < Parent_SubJet_2.perp()) std::swap(Parent_SubJet_1,Parent_SubJet_2);
1496  Double_t Leading_Track_pT=0.0;
1497  vector <fastjet::PseudoJet> Hard_SubJet_Constituents = sorted_by_pt(Parent_SubJet_1.constituents());
1498  for (Int_t i_Leading_Track=0; i_Leading_Track < Hard_SubJet_Constituents.size(); i_Leading_Track++){
1499  if (Hard_SubJet_Constituents[i_Leading_Track].perp() > Leading_Track_pT) Leading_Track_pT=Hard_SubJet_Constituents[i_Leading_Track].perp();
1500  }
1501  // Splittings_LeadingSubJetpT.push_back(Parent_SubJet_1.perp());
1502  Splittings_HardestSubJetD0.push_back(Leading_Track_pT);
1503  Splittings_DeltaR.push_back(Parent_SubJet_1.delta_R(Parent_SubJet_2));
1504  Splittings_Zg.push_back(Parent_SubJet_2.perp()/(Parent_SubJet_1.perp()+Parent_SubJet_2.perp()));
1505  Splittings_RadiatorE.push_back(Daughter_Jet.E());
1506  Splittings_RadiatorpT.push_back(Daughter_Jet.perp());
1507  Daughter_Jet=Parent_SubJet_1;
1508  }
1509 
1510 
1511  } catch (fastjet::Error) { /*return -1;*/ }
1512 
1513 
1514 
1515  fShapesVar[0] = Inclusive_Jets[i_Jet].perp();
1516  fShapesVar[1] = 0.0;
1517  fShapesVar[2] = HardestTrack_Pt;
1518  fShapesVar[3] = 0.0;
1519  fShapesVar[4] = 0.0;
1520  fShapesVar[5] = 0.0;
1521  fShapesVar[6] = 0.0;
1522  fShapesVar[7] = 0.0;
1523  fShapesVar[8] = 0.0;
1524  fShapesVar[9] = 0.0;
1525  // fShapesVar[10] = NTracks;
1526  //fShapesVar[11] = 0.0;
1527  //fShapesVar[12] = Jet_Eta;
1528  //fShapesVar[13] = 0.0;
1529  //fShapesVar[14] = HardestTrack_Eta;
1530  //fShapesVar[15] = 0.0;
1531  //fShapesVar[16] = HardestTrack_Y;
1532  //fShapesVar[17] = 0.0;
1533 
1534 
1535  fShapesVar_Splittings_DeltaR.push_back(Splittings_DeltaR);
1536  fShapesVar_Splittings_DeltaR_Truth.push_back(Splittings_DeltaR);
1537  fShapesVar_Splittings_Zg.push_back(Splittings_Zg);
1538  fShapesVar_Splittings_Zg_Truth.push_back(Splittings_Zg);
1539  // fShapesVar_Splittings_LeadingSubJetpT.push_back(Splittings_LeadingSubJetpT);
1540  //fShapesVar_Splittings_LeadingSubJetpT_Truth.push_back(Splittings_LeadingSubJetpT);
1541  fShapesVar_Splittings_HardestSubJetD0.push_back(Splittings_HardestSubJetD0);
1542  fShapesVar_Splittings_HardestSubJetD0_Truth.push_back(Splittings_HardestSubJetD0);
1543  fShapesVar_Splittings_RadiatorE.push_back(Splittings_RadiatorE);
1544  fShapesVar_Splittings_RadiatorE_Truth.push_back(Splittings_RadiatorE);
1545  fShapesVar_Splittings_RadiatorpT.push_back(Splittings_RadiatorpT);
1546  fShapesVar_Splittings_RadiatorpT_Truth.push_back(Splittings_RadiatorpT);
1547 
1548 
1549  fTreeResponseMatrixAxis->Fill();
1550  fTreeSplittings->Fill();
1551 
1554  fShapesVar_Splittings_Zg.clear();
1556  // fShapesVar_Splittings_LeadingSubJetpT.clear();
1557  //fShapesVar_Splittings_LeadingSubJetpT_Truth.clear();
1564 
1565  }
1566  }
1567 
1568 
1569 
1570 
1571 
1572 
1573 
1574 
1575 
1576 
1577 
1578 
1579 
1580  delete fFastJetWrapper_Truth;
1581  delete fFastJetWrapper;
1582 
1583  return kTRUE;
1584 
1585 
1586 }
1587 
1588 
1589 
1590 
1591 
1592 
1593 //________________________________________________________________________
1595  //
1596  // retrieve event objects
1597  //
1599  return kFALSE;
1600 
1601  return kTRUE;
1602 }
1603 
1604 
1605 //_______________________________________________________________________
1607 {
1608  // Called once at the end of the analysis.
1609 
1610 }
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()
virtual AliAODMCParticle * GetAcceptMCParticle(Int_t i=-1) const
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
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:312
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:338
void swap(PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &first, PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &second)
void SetCharge(EChargeCut_t c)
std::vector< std::vector< Double_t > > fShapesVar_Splittings_RadiatorpT
std::vector< std::vector< Double_t > > fShapesVar_Splittings_RadiatorE