AliPhysics  master (3d17d9d)
AliAnalysisTaskJetExtractor.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2018, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: R. Haake. *
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 #include <iostream>
17 #include <cstdlib>
18 #include <algorithm>
19 #include <vector>
20 #include <TClonesArray.h>
21 #include <TH1F.h>
22 #include <TH2F.h>
23 #include <TH3F.h>
24 #include <TProfile.h>
25 #include <TTree.h>
26 #include <TList.h>
27 #include <TLorentzVector.h>
28 #include <TNamed.h>
29 #include <TGrid.h>
30 #include <TFile.h>
31 #include <TSystem.h>
32 //#include <TTimeStamp.h>
33 
34 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
35  #include <TPython.h>
36 #endif
37 
38 #include "AliMCEvent.h"
39 #include "AliESDVertex.h"
40 #include "AliAODVertex.h"
41 #include "AliAODPid.h"
42 
43 #include "AliRhoParameter.h"
44 
45 #include "AliVTrack.h"
46 #include "AliVHeader.h"
47 #include "AliEmcalJet.h"
48 #include "AliLog.h"
49 #include "AliJetContainer.h"
50 #include "AliParticleContainer.h"
51 #include "AliClusterContainer.h"
52 #include "AliAODTrack.h"
53 #include "AliVParticle.h"
54 #include "TRandom3.h"
55 #include "AliEmcalPythiaInfo.h"
57 #include "AliAnalysisManager.h"
58 #include "AliEmcalContainerUtils.h"
59 #include "AliFJWrapper.h"
60 
61 
62 #include "AliGenHepMCEventHeader.h"
63 #include "AliGenHijingEventHeader.h"
64 #include "AliHFJetsTaggingVertex.h"
65 #include "AliRDHFJetsCutsVertex.h"
66 
68 
70 ClassImp(AliEmcalJetTree)
72 
76 
77 //________________________________________________________________________
78 AliEmcalJetTree::AliEmcalJetTree() : TNamed("CustomTree", "CustomTree"), fJetTree(0), fInitialized(0), fExtractionPercentages(), fExtractionPercentagePtBins(), fExtractionJetTypes_HM(), fExtractionJetTypes_PM()
79 {
80  // For these arrays, we need to reserve memory
81  fBuffer_Track_Pt = new Float_t[kMaxNumConstituents];
82  fBuffer_Track_Eta = new Float_t[kMaxNumConstituents];
83  fBuffer_Track_Phi = new Float_t[kMaxNumConstituents];
84  fBuffer_Track_Charge = new Float_t[kMaxNumConstituents];
85  fBuffer_Track_Label = new Int_t [kMaxNumConstituents];
86  fBuffer_Track_ProdVtx_X = new Float_t[kMaxNumConstituents];
87  fBuffer_Track_ProdVtx_Y = new Float_t[kMaxNumConstituents];
88  fBuffer_Track_ProdVtx_Z = new Float_t[kMaxNumConstituents];
89 
90  fBuffer_Cluster_Pt = new Float_t[kMaxNumConstituents];
91  fBuffer_Cluster_E = new Float_t[kMaxNumConstituents];
92  fBuffer_Cluster_Eta = new Float_t[kMaxNumConstituents];
93  fBuffer_Cluster_Phi = new Float_t[kMaxNumConstituents];
94  fBuffer_Cluster_M02 = new Float_t[kMaxNumConstituents];
95  fBuffer_Cluster_Time = new Float_t[kMaxNumConstituents];
96  fBuffer_Cluster_Label = new Int_t[kMaxNumConstituents];
97 }
98 
99 //________________________________________________________________________
100 AliEmcalJetTree::AliEmcalJetTree(const char* name) : TNamed(name, name), fJetTree(0), fInitialized(0), fExtractionPercentages(), fExtractionPercentagePtBins(), fExtractionJetTypes_HM(), fExtractionJetTypes_PM()
101 {
102  // For these arrays, we need to reserve memory
111 
119 }
120 
121 //________________________________________________________________________
122 Bool_t AliEmcalJetTree::AddJetToTree(AliEmcalJet* jet, Bool_t saveConstituents, Bool_t saveConstituentsIP, Bool_t saveCaloClusters, Double_t* vertex, Float_t rho, Float_t rhoMass, Float_t centrality, Int_t multiplicity, Long64_t eventID, Float_t magField)
123 {
124  if(!fInitialized)
125  AliFatal("Tree is not initialized.");
126 
127  fBuffer_JetPt = jet->Pt() - rho*jet->Area();
128 
129  // Check if jet type is contained in extraction list
130  if( (fExtractionJetTypes_PM.size() || fExtractionJetTypes_HM.size()) &&
133  return kFALSE;
134 
135  // Check acceptance percentage for the given jet and discard statistically on demand
136  Bool_t inPtRange = kFALSE;
137  for(size_t i = 0; i<fExtractionPercentages.size(); i++)
138  {
140  {
141  inPtRange = kTRUE;
142  if(fRandomGenerator->Rndm() >= fExtractionPercentages[i])
143  return kFALSE;
144  break;
145  }
146  }
147  if(!inPtRange && fExtractionPercentagePtBins.size()) // only discard jet if fExtractionPercentagePtBins was given
148  return kFALSE;
149 
150  // Set basic properties
151  fBuffer_JetEta = jet->Eta();
152  fBuffer_JetPhi = jet->Phi();
153  fBuffer_JetArea = jet->Area();
154 
155  // Set event properties
158  fBuffer_Event_Vertex_X = vertex ? vertex[0] : 0;
159  fBuffer_Event_Vertex_Y = vertex ? vertex[1] : 0;
160  fBuffer_Event_Vertex_Z = vertex ? vertex[2] : 0;
162  fBuffer_Event_Multiplicity = multiplicity;
163  fBuffer_Event_ID = eventID;
164  fBuffer_Event_MagneticField = magField;
165 
166  // Extract basic constituent track properties directly from AliEmcalJet object
167  fBuffer_NumTracks = 0;
168  if(saveConstituents || saveConstituentsIP)
169  for(Int_t i = 0; i < jet->GetNumberOfParticleConstituents(); i++)
170  {
171  const AliVParticle* particle = jet->GetParticleConstituents()[i].GetParticle();
172  if(!particle) continue;
173 
174  if(saveConstituents)
175  {
176  fBuffer_Track_Pt[fBuffer_NumTracks] = particle->Pt();
177  fBuffer_Track_Eta[fBuffer_NumTracks] = particle->Eta();
178  fBuffer_Track_Phi[fBuffer_NumTracks] = particle->Phi();
179  fBuffer_Track_Charge[fBuffer_NumTracks] = particle->Charge();
180  fBuffer_Track_Label[fBuffer_NumTracks] = particle->GetLabel();
181  }
182  if(saveConstituentsIP) // track production vertices
183  {
184  fBuffer_Track_ProdVtx_X[fBuffer_NumTracks] = particle->Xv();
185  fBuffer_Track_ProdVtx_Y[fBuffer_NumTracks] = particle->Yv();
186  fBuffer_Track_ProdVtx_Z[fBuffer_NumTracks] = particle->Zv();
187  }
189  }
190 
191  // Extract basic constituent cluster properties directly from AliEmcalJet object
193  if(saveCaloClusters)
194  for(Int_t i = 0; i < jet->GetNumberOfClusterConstituents(); i++)
195  {
196  const AliVCluster* cluster = jet->GetClusterConstituents()[i].GetCluster();
197  if(!cluster) continue;
198 
199  // #### Retrieve cluster pT
200  TLorentzVector clusterMomentum;
201  cluster->GetMomentum(clusterMomentum, vertex);
202  // ####
203 
204  fBuffer_Cluster_Pt[fBuffer_NumClusters] = clusterMomentum.Perp();
205  fBuffer_Cluster_E[fBuffer_NumClusters] = cluster->E();
206  fBuffer_Cluster_Eta[fBuffer_NumClusters] = clusterMomentum.Eta();
207  fBuffer_Cluster_Phi[fBuffer_NumClusters] = clusterMomentum.Phi();
208  fBuffer_Cluster_M02[fBuffer_NumClusters] = cluster->GetM02();
209  fBuffer_Cluster_Time[fBuffer_NumClusters] = cluster->GetTOF();
210  fBuffer_Cluster_Label[fBuffer_NumClusters] = cluster->GetLabel();
212  }
213 
214 
215  // Add all buffers to tree
216  fJetTree->Fill();
217 
218  return kTRUE;
219 }
220 
221 //________________________________________________________________________
222 void AliEmcalJetTree::FillBuffer_TriggerTracks(std::vector<Float_t>& triggerTrackPt, std::vector<Float_t>& triggerTrackDeltaEta, std::vector<Float_t>& triggerTrackDeltaPhi)
223 {
224  fBuffer_NumTriggerTracks = triggerTrackPt.size();
225  fJetTree->SetBranchAddress("Jet_TriggerTrack_Pt", triggerTrackPt.data());
226  fJetTree->SetBranchAddress("Jet_TriggerTrack_dEta", triggerTrackDeltaEta.data());
227  fJetTree->SetBranchAddress("Jet_TriggerTrack_dPhi", triggerTrackDeltaPhi.data());
228 }
229 
230 //________________________________________________________________________
231 void AliEmcalJetTree::FillBuffer_ImpactParameters(std::vector<Float_t>& trackIP_d0, std::vector<Float_t>& trackIP_z0, std::vector<Float_t>& trackIP_d0cov, std::vector<Float_t>& trackIP_z0cov)
232 {
233  fJetTree->SetBranchAddress("Jet_Track_CovIPd", trackIP_d0cov.data());
234  fJetTree->SetBranchAddress("Jet_Track_CovIPz", trackIP_z0cov.data());
235  fJetTree->SetBranchAddress("Jet_Track_IPd", trackIP_d0.data());
236  fJetTree->SetBranchAddress("Jet_Track_IPz", trackIP_z0.data());
237 }
238 
239 //________________________________________________________________________
240 void AliEmcalJetTree::FillBuffer_MonteCarlo(Int_t motherParton, Int_t motherHadron, Int_t partonInitialCollision, Float_t matchedJetDistance_Det, Float_t matchedJetPt_Det, Float_t matchedJetMass_Det, Float_t matchedJetAngularity_Det, Float_t matchedJetpTD_Det, Float_t matchedJetDistance_Part, Float_t matchedJetPt_Part, Float_t matchedJetMass_Part, Float_t matchedJetAngularity_Part, Float_t matchedJetpTD_Part, Float_t truePtFraction, Float_t truePtFraction_PartLevel, Float_t ptHard, Float_t eventWeight, Float_t impactParameter)
241 {
242  fBuffer_Jet_MC_MotherParton = motherParton;
243  fBuffer_Jet_MC_MotherHadron = motherHadron;
244  fBuffer_Jet_MC_MotherIC = partonInitialCollision;
245  fBuffer_Jet_MC_MatchedDetLevelJet_Distance = matchedJetDistance_Det;
246  fBuffer_Jet_MC_MatchedDetLevelJet_Pt = matchedJetPt_Det;
247  fBuffer_Jet_MC_MatchedDetLevelJet_Mass = matchedJetMass_Det;
248  fBuffer_Jet_MC_MatchedDetLevelJet_Angularity = matchedJetAngularity_Det;
249  fBuffer_Jet_MC_MatchedDetLevelJet_pTD = matchedJetpTD_Det;
250 
251  fBuffer_Jet_MC_MatchedPartLevelJet_Distance = matchedJetDistance_Part;
252  fBuffer_Jet_MC_MatchedPartLevelJet_Pt = matchedJetPt_Part;
253  fBuffer_Jet_MC_MatchedPartLevelJet_Mass = matchedJetMass_Part;
254  fBuffer_Jet_MC_MatchedPartLevelJet_Angularity = matchedJetAngularity_Part;
255  fBuffer_Jet_MC_MatchedPartLevelJet_pTD = matchedJetpTD_Part;
256 
257  fBuffer_Jet_MC_TruePtFraction = truePtFraction;
258  fBuffer_Jet_MC_TruePtFraction_PartLevel = truePtFraction_PartLevel;
259 
260 
261  fBuffer_Event_PtHard = ptHard;
262  fBuffer_Event_Weight = eventWeight;
263  fBuffer_Event_ImpactParameter = impactParameter;
264 }
265 //________________________________________________________________________
266 void AliEmcalJetTree::FillBuffer_PID(std::vector<Float_t>& trackPID_ITS, std::vector<Float_t>& trackPID_TPC, std::vector<Float_t>& trackPID_TOF, std::vector<Float_t>& trackPID_TRD, std::vector<Short_t>& trackPID_Reco, std::vector<Int_t>& trackPID_Truth)
267 {
268  fJetTree->SetBranchAddress("Jet_Track_PID_ITS", trackPID_ITS.data());
269  fJetTree->SetBranchAddress("Jet_Track_PID_TPC", trackPID_TPC.data());
270  fJetTree->SetBranchAddress("Jet_Track_PID_TOF", trackPID_TOF.data());
271  fJetTree->SetBranchAddress("Jet_Track_PID_TRD", trackPID_TRD.data());
272  fJetTree->SetBranchAddress("Jet_Track_PID_Reconstructed", trackPID_Reco.data());
273  if(trackPID_Truth.data())
274  fJetTree->SetBranchAddress("Jet_Track_PID_Truth", trackPID_Truth.data());
275 }
276 
277 //________________________________________________________________________
278 void AliEmcalJetTree::FillBuffer_JetShapes(AliEmcalJet* jet, Double_t leSub_noCorr, Double_t angularity, Double_t momentumDispersion, Double_t trackPtMean, Double_t trackPtMedian)
279 {
280  fBuffer_Shape_Mass_NoCorr = jet->M();
285  fBuffer_Shape_Angularity_NoCorr = angularity;
294  fBuffer_Shape_LeSub_NoCorr = leSub_noCorr;
295  fBuffer_Shape_MomentumDispersion = momentumDispersion;
296  fBuffer_Shape_TrackPtMean = trackPtMean;
297  fBuffer_Shape_TrackPtMedian = trackPtMedian;
298 }
299 
300 //________________________________________________________________________
301 void AliEmcalJetTree::FillBuffer_Splittings(std::vector<Float_t>& splittings_radiatorE, std::vector<Float_t>& splittings_kT, std::vector<Float_t>& splittings_theta, Bool_t saveSecondaryVertices, std::vector<Int_t>& splittings_secVtx_rank, std::vector<Int_t>& splittings_secVtx_index)
302 {
303  fBuffer_NumSplittings = splittings_radiatorE.size();
304  fJetTree->SetBranchAddress("Jet_Splitting_RadiatorE", splittings_radiatorE.data());
305  fJetTree->SetBranchAddress("Jet_Splitting_kT", splittings_kT.data());
306  fJetTree->SetBranchAddress("Jet_Splitting_Theta", splittings_theta.data());
307  if(saveSecondaryVertices)
308  {
309  fJetTree->SetBranchAddress("Jet_Splitting_SecVtx_Rank", splittings_secVtx_rank.data());
310  fJetTree->SetBranchAddress("Jet_Splitting_SecVtx_Index", splittings_secVtx_index.data());
311  }
312 }
313 
314 //________________________________________________________________________
315 void AliEmcalJetTree::FillBuffer_SecVertices(std::vector<Float_t>& secVtx_X, std::vector<Float_t>& secVtx_Y, std::vector<Float_t>& secVtx_Z, std::vector<Float_t>& secVtx_Mass, std::vector<Float_t>& secVtx_Lxy, std::vector<Float_t>& secVtx_SigmaLxy, std::vector<Float_t>& secVtx_Chi2, std::vector<Float_t>& secVtx_Dispersion)
316 {
317  fBuffer_NumSecVertices = secVtx_X.size();
318  fJetTree->SetBranchAddress("Jet_SecVtx_X", secVtx_X.data());
319  fJetTree->SetBranchAddress("Jet_SecVtx_Y", secVtx_Y.data());
320  fJetTree->SetBranchAddress("Jet_SecVtx_Z", secVtx_Z.data());
321  fJetTree->SetBranchAddress("Jet_SecVtx_Mass", secVtx_Mass.data());
322  fJetTree->SetBranchAddress("Jet_SecVtx_Lxy", secVtx_Lxy.data());
323  fJetTree->SetBranchAddress("Jet_SecVtx_SigmaLxy", secVtx_SigmaLxy.data());
324  fJetTree->SetBranchAddress("Jet_SecVtx_Chi2", secVtx_Chi2.data());
325  fJetTree->SetBranchAddress("Jet_SecVtx_Dispersion", secVtx_Dispersion.data());
326 }
327 
328 //________________________________________________________________________
329 void AliEmcalJetTree::InitializeTree(Bool_t saveCaloClusters, Bool_t saveMCInformation, Bool_t saveMatchedJets_Det, Bool_t saveMatchedJets_Part, Bool_t saveConstituents, Bool_t saveConstituentsIP, Bool_t saveConstituentPID, Bool_t saveJetShapes, Bool_t saveSplittings, Bool_t saveSecondaryVertices, Bool_t saveTriggerTracks)
330 {
331  // Create the tree with active branches
332 
333  void* dummy = 0; // for some branches, do not need a buffer pointer yet
334 
335 
336  if(fInitialized)
337  AliFatal("Tree is already initialized.");
338 
339  // ### Prepare the jet tree
340  fJetTree = new TTree(Form("JetTree_%s", GetName()), "");
341 
342  fJetTree->Branch("Jet_Pt",&fBuffer_JetPt,"Jet_Pt/F");
343  fJetTree->Branch("Jet_Phi",&fBuffer_JetPhi,"Jet_Phi/F");
344  fJetTree->Branch("Jet_Eta",&fBuffer_JetEta,"Jet_Eta/F");
345  fJetTree->Branch("Jet_Area",&fBuffer_JetArea,"Jet_Area/F");
346  fJetTree->Branch("Jet_NumTracks",&fBuffer_NumTracks,"Jet_NumTracks/I");
347  if(saveCaloClusters)
348  fJetTree->Branch("Jet_NumClusters",&fBuffer_NumClusters,"Jet_NumClusters/I");
349 
350  fJetTree->Branch("Event_BackgroundDensity",&fBuffer_Event_BackgroundDensity,"Event_BackgroundDensity/F");
351  fJetTree->Branch("Event_BackgroundDensityMass",&fBuffer_Event_BackgroundDensityMass,"Event_BackgroundDensityMass/F");
352  fJetTree->Branch("Event_Vertex_X",&fBuffer_Event_Vertex_X,"Event_Vertex_X/F");
353  fJetTree->Branch("Event_Vertex_Y",&fBuffer_Event_Vertex_Y,"Event_Vertex_Y/F");
354  fJetTree->Branch("Event_Vertex_Z",&fBuffer_Event_Vertex_Z,"Event_Vertex_Z/F");
355  fJetTree->Branch("Event_Centrality",&fBuffer_Event_Centrality,"Event_Centrality/F");
356  fJetTree->Branch("Event_Multiplicity",&fBuffer_Event_Multiplicity,"Event_Multiplicity/I");
357  fJetTree->Branch("Event_ID",&fBuffer_Event_ID,"Event_ID/L");
358  fJetTree->Branch("Event_MagneticField",&fBuffer_Event_MagneticField,"Event_MagneticField/F");
359 
360  if(saveMCInformation)
361  {
362  fJetTree->Branch("Event_PtHard",&fBuffer_Event_PtHard,"Event_PtHard/F");
363  fJetTree->Branch("Event_Weight",&fBuffer_Event_Weight,"Event_Weight/F");
364  fJetTree->Branch("Event_ImpactParameter",&fBuffer_Event_ImpactParameter,"Event_ImpactParameter/F");
365  }
366 
367  if(saveConstituents)
368  {
369  fJetTree->Branch("Jet_Track_Pt",fBuffer_Track_Pt,"Jet_Track_Pt[Jet_NumTracks]/F");
370  fJetTree->Branch("Jet_Track_Phi",fBuffer_Track_Phi,"Jet_Track_Phi[Jet_NumTracks]/F");
371  fJetTree->Branch("Jet_Track_Eta",fBuffer_Track_Eta,"Jet_Track_Eta[Jet_NumTracks]/F");
372  fJetTree->Branch("Jet_Track_Charge",fBuffer_Track_Charge,"Jet_Track_Charge[Jet_NumTracks]/F");
373  if(saveMCInformation)
374  fJetTree->Branch("Jet_Track_Label",fBuffer_Track_Label,"Jet_Track_Label[Jet_NumTracks]/I");
375  }
376 
377  if(saveCaloClusters)
378  {
379  fJetTree->Branch("Jet_Cluster_Pt",fBuffer_Cluster_Pt,"Jet_Cluster_Pt[Jet_NumClusters]/F");
380  fJetTree->Branch("Jet_Cluster_E",fBuffer_Cluster_E,"Jet_Cluster_Pt[Jet_NumClusters]/F");
381  fJetTree->Branch("Jet_Cluster_Phi",fBuffer_Cluster_Phi,"Jet_Cluster_Phi[Jet_NumClusters]/F");
382  fJetTree->Branch("Jet_Cluster_Eta",fBuffer_Cluster_Eta,"Jet_Cluster_Eta[Jet_NumClusters]/F");
383  fJetTree->Branch("Jet_Cluster_M02",fBuffer_Cluster_M02,"Jet_Cluster_M02[Jet_NumClusters]/F");
384  fJetTree->Branch("Jet_Cluster_Time",fBuffer_Cluster_Time,"Jet_Cluster_Time[Jet_NumClusters]/F");
385  if(saveMCInformation)
386  fJetTree->Branch("Jet_Cluster_Label",fBuffer_Cluster_Label,"Jet_Cluster_Label[Jet_NumClusters]/I");
387  }
388 
389  if(saveConstituentsIP)
390  {
391  fJetTree->Branch("Jet_Track_IPd",&dummy,"Jet_Track_IPd[Jet_NumTracks]/F");
392  fJetTree->Branch("Jet_Track_IPz",&dummy,"Jet_Track_IPz[Jet_NumTracks]/F");
393  fJetTree->Branch("Jet_Track_CovIPd",&dummy,"Jet_Track_CovIPd[Jet_NumTracks]/F");
394  fJetTree->Branch("Jet_Track_CovIPz",&dummy,"Jet_Track_CovIPz[Jet_NumTracks]/F");
395 
396  fJetTree->Branch("Jet_Track_ProdVtx_X",fBuffer_Track_ProdVtx_X,"Jet_Track_ProdVtx_X[Jet_NumTracks]/F");
397  fJetTree->Branch("Jet_Track_ProdVtx_Y",fBuffer_Track_ProdVtx_Y,"Jet_Track_ProdVtx_Y[Jet_NumTracks]/F");
398  fJetTree->Branch("Jet_Track_ProdVtx_Z",fBuffer_Track_ProdVtx_Z,"Jet_Track_ProdVtx_Z[Jet_NumTracks]/F");
399  }
400 
401  if(saveConstituentPID)
402  {
403  fJetTree->Branch("Jet_Track_PID_ITS",&dummy,"Jet_Track_PID_ITS[Jet_NumTracks]/F");
404  fJetTree->Branch("Jet_Track_PID_TPC",&dummy,"Jet_Track_PID_TPC[Jet_NumTracks]/F");
405  fJetTree->Branch("Jet_Track_PID_TOF",&dummy,"Jet_Track_PID_TOF[Jet_NumTracks]/F");
406  fJetTree->Branch("Jet_Track_PID_TRD",&dummy,"Jet_Track_PID_TRD[Jet_NumTracks]/F");
407 
408  fJetTree->Branch("Jet_Track_PID_Reconstructed",&dummy,"Jet_Track_PID_Reconstructed[Jet_NumTracks]/S");
409  if(saveMCInformation)
410  fJetTree->Branch("Jet_Track_PID_Truth",&dummy,"Jet_Track_PID_Truth[Jet_NumTracks]/I");
411  }
412 
413  if(saveJetShapes)
414  {
415  fJetTree->Branch("Jet_Shape_Mass_NoCorr",&fBuffer_Shape_Mass_NoCorr,"Jet_Shape_Mass_NoCorr/F");
416  fJetTree->Branch("Jet_Shape_Mass_DerivCorr_1",&fBuffer_Shape_Mass_DerivCorr_1,"Jet_Shape_Mass_DerivCorr_1/F");
417  fJetTree->Branch("Jet_Shape_Mass_DerivCorr_2",&fBuffer_Shape_Mass_DerivCorr_2,"Jet_Shape_Mass_DerivCorr_2/F");
418  fJetTree->Branch("Jet_Shape_pTD_DerivCorr_1",&fBuffer_Shape_pTD_DerivCorr_1,"Jet_Shape_pTD_DerivCorr_1/F");
419  fJetTree->Branch("Jet_Shape_pTD_DerivCorr_2",&fBuffer_Shape_pTD_DerivCorr_2,"Jet_Shape_pTD_DerivCorr_2/F");
420  fJetTree->Branch("Jet_Shape_LeSub_NoCorr",&fBuffer_Shape_LeSub_NoCorr,"Jet_Shape_LeSub_NoCorr/F");
421  fJetTree->Branch("Jet_Shape_LeSub_DerivCorr",&fBuffer_Shape_LeSub_DerivCorr,"Jet_Shape_LeSub_DerivCorr/F");
422  fJetTree->Branch("Jet_Shape_Angularity",&fBuffer_Shape_Angularity_NoCorr,"Jet_Shape_Angularity/F");
423  fJetTree->Branch("Jet_Shape_Angularity_DerivCorr_1",&fBuffer_Shape_Angularity_DerivCorr_1,"Jet_Shape_Angularity_DerivCorr_1/F");
424  fJetTree->Branch("Jet_Shape_Angularity_DerivCorr_2",&fBuffer_Shape_Angularity_DerivCorr_2,"Jet_Shape_Angularity_DerivCorr_2/F");
425  fJetTree->Branch("Jet_Shape_Circularity_DerivCorr_1",&fBuffer_Shape_Circularity_DerivCorr_1,"Jet_Shape_Circularity_DerivCorr_1/F");
426  fJetTree->Branch("Jet_Shape_Circularity_DerivCorr_2",&fBuffer_Shape_Circularity_DerivCorr_2,"Jet_Shape_Circularity_DerivCorr_2/F");
427  fJetTree->Branch("Jet_Shape_Sigma2_DerivCorr_1",&fBuffer_Shape_Sigma2_DerivCorr_1,"Jet_Shape_Sigma2_DerivCorr_1/F");
428  fJetTree->Branch("Jet_Shape_Sigma2_DerivCorr_2",&fBuffer_Shape_Sigma2_DerivCorr_2,"Jet_Shape_Sigma2_DerivCorr_2/F");
429  fJetTree->Branch("Jet_Shape_NumTracks_DerivCorr",&fBuffer_Shape_NumTracks_DerivCorr,"Jet_Shape_NumTracks_DerivCorr/F");
430  fJetTree->Branch("Jet_Shape_MomentumDispersion",&fBuffer_Shape_MomentumDispersion,"Jet_Shape_MomentumDispersion/F");
431  fJetTree->Branch("Jet_Shape_TrackPtMean",&fBuffer_Shape_TrackPtMean,"Jet_Shape_TrackPtMean/F");
432  fJetTree->Branch("Jet_Shape_TrackPtMedian",&fBuffer_Shape_TrackPtMedian,"Jet_Shape_TrackPtMedian/F");
433  }
434 
435  if(saveSplittings)
436  {
437  fJetTree->Branch("Jet_NumSplittings",&fBuffer_NumSplittings,"Jet_NumSplittings/I");
438  fJetTree->Branch("Jet_Splitting_Theta",&dummy,"Jet_Splitting_Theta[Jet_NumSplittings]/F");
439  fJetTree->Branch("Jet_Splitting_RadiatorE",&dummy,"Jet_Splitting_RadiatorE[Jet_NumSplittings]/F");
440  fJetTree->Branch("Jet_Splitting_kT",&dummy,"Jet_Splitting_kT[Jet_NumSplittings]/F");
441  if(saveSecondaryVertices)
442  {
443  fJetTree->Branch("Jet_Splitting_SecVtx_Rank",&dummy,"Jet_Splitting_SecVtx_Rank[Jet_NumSplittings]/I");
444  fJetTree->Branch("Jet_Splitting_SecVtx_Index",&dummy,"Jet_Splitting_SecVtx_Index[Jet_NumSplittings]/I");
445  }
446  }
447 
448  if(saveMCInformation)
449  {
450  fJetTree->Branch("Jet_MC_MotherParton",&fBuffer_Jet_MC_MotherParton,"Jet_MC_MotherParton/I");
451  fJetTree->Branch("Jet_MC_MotherHadron",&fBuffer_Jet_MC_MotherHadron,"Jet_MC_MotherHadron/I");
452  fJetTree->Branch("Jet_MC_MotherIC",&fBuffer_Jet_MC_MotherIC,"Jet_MC_MotherIC/I");
453 
454  if(saveMatchedJets_Det)
455  {
456  fJetTree->Branch("Jet_MC_MatchedDetLevelJet_Distance",&fBuffer_Jet_MC_MatchedDetLevelJet_Distance,"Jet_MC_MatchedDetLevelJet_Distance/F");
457  fJetTree->Branch("Jet_MC_MatchedDetLevelJet_Pt",&fBuffer_Jet_MC_MatchedDetLevelJet_Pt,"Jet_MC_MatchedDetLevelJet_Pt/F");
458  fJetTree->Branch("Jet_MC_MatchedDetLevelJet_Mass",&fBuffer_Jet_MC_MatchedDetLevelJet_Mass,"Jet_MC_MatchedDetLevelJet_Mass/F");
459  fJetTree->Branch("Jet_MC_MatchedDetLevelJet_Angularity",&fBuffer_Jet_MC_MatchedDetLevelJet_Angularity,"Jet_MC_MatchedDetLevelJet_Angularity/F");
460  fJetTree->Branch("Jet_MC_MatchedDetLevelJet_pTD",&fBuffer_Jet_MC_MatchedDetLevelJet_pTD,"Jet_MC_MatchedDetLevelJet_pTD/F");
461  }
462  if(saveMatchedJets_Part)
463  {
464  fJetTree->Branch("Jet_MC_MatchedPartLevelJet_Distance",&fBuffer_Jet_MC_MatchedPartLevelJet_Distance,"Jet_MC_MatchedPartLevelJet_Distance/F");
465  fJetTree->Branch("Jet_MC_MatchedPartLevelJet_Pt",&fBuffer_Jet_MC_MatchedPartLevelJet_Pt,"Jet_MC_MatchedPartLevelJet_Pt/F");
466  fJetTree->Branch("Jet_MC_MatchedPartLevelJet_Mass",&fBuffer_Jet_MC_MatchedPartLevelJet_Mass,"Jet_MC_MatchedPartLevelJet_Mass/F");
467  fJetTree->Branch("Jet_MC_MatchedPartLevelJet_Angularity",&fBuffer_Jet_MC_MatchedPartLevelJet_Angularity,"Jet_MC_MatchedPartLevelJet_Angularity/F");
468  fJetTree->Branch("Jet_MC_MatchedPartLevelJet_pTD",&fBuffer_Jet_MC_MatchedPartLevelJet_pTD,"Jet_MC_MatchedPartLevelJet_pTD/F");
469  }
470 
471  fJetTree->Branch("Jet_MC_TruePtFraction",&fBuffer_Jet_MC_TruePtFraction,"Jet_MC_TruePtFraction/F");
472  fJetTree->Branch("Jet_MC_TruePtFraction_PartLevel",&fBuffer_Jet_MC_TruePtFraction_PartLevel,"Jet_MC_TruePtFraction_PartLevel/F");
473  }
474 
475 
476  if(saveSecondaryVertices)
477  {
478  fJetTree->Branch("Jet_NumSecVertices",&fBuffer_NumSecVertices,"Jet_NumSecVertices/I");
479 
480  fJetTree->Branch("Jet_SecVtx_X",&dummy,"Jet_SecVtx_X[Jet_NumSecVertices]/F");
481  fJetTree->Branch("Jet_SecVtx_Y",&dummy,"Jet_SecVtx_Y[Jet_NumSecVertices]/F");
482  fJetTree->Branch("Jet_SecVtx_Z",&dummy,"Jet_SecVtx_Z[Jet_NumSecVertices]/F");
483  fJetTree->Branch("Jet_SecVtx_Mass",&dummy,"Jet_SecVtx_Mass[Jet_NumSecVertices]/F");
484  fJetTree->Branch("Jet_SecVtx_Lxy",&dummy,"Jet_SecVtx_Lxy[Jet_NumSecVertices]/F");
485  fJetTree->Branch("Jet_SecVtx_SigmaLxy",&dummy,"Jet_SecVtx_SigmaLxy[Jet_NumSecVertices]/F");
486  fJetTree->Branch("Jet_SecVtx_Chi2",&dummy,"Jet_SecVtx_Chi2[Jet_NumSecVertices]/F");
487  fJetTree->Branch("Jet_SecVtx_Dispersion",&dummy,"Jet_SecVtx_Dispersion[Jet_NumSecVertices]/F");
488  }
489 
490  // Trigger track requirement active -> Save trigger track
491  if(saveTriggerTracks)
492  {
493  fJetTree->Branch("Jet_NumTriggerTracks",&fBuffer_NumTriggerTracks,"Jet_NumTriggerTracks/I");
494  fJetTree->Branch("Jet_TriggerTrack_Pt",&dummy,"Jet_TriggerTrack_Pt[Jet_NumTriggerTracks]/F");
495  fJetTree->Branch("Jet_TriggerTrack_dEta",&dummy,"Jet_TriggerTrack_dEta[Jet_NumTriggerTracks]/F");
496  fJetTree->Branch("Jet_TriggerTrack_dPhi",&dummy,"Jet_TriggerTrack_dPhi[Jet_NumTriggerTracks]/F");
497  }
498 
499  fInitialized = kTRUE;
500 }
501 
502 //________________________________________________________________________
504  AliAnalysisTaskEmcalJet("AliAnalysisTaskJetExtractor", kTRUE),
505  fSaveConstituents(0), fSaveConstituentsIP(0), fSaveConstituentPID(0), fSaveJetShapes(0), fSaveJetSplittings(0), fSaveMCInformation(0), fSaveSecondaryVertices(0), fSaveTriggerTracks(0), fSaveCaloClusters(0),
506  fEventPercentage(1.0),
507  fEventCut_TriggerTrackMinPt(0),
508  fEventCut_TriggerTrackMaxPt(0),
509  fEventCut_TriggerTrackMinLabel(-9999999),
510  fEventCut_TriggerTrackMaxLabel(+9999999),
511  fEventCut_TriggerTrackOrigin(-1),
512  fTruthMinLabel(0),
513  fTruthMaxLabel(100000),
514  fHadronMatchingRadius(0.4),
515  fJetMatchingRadius(0.3),
516  fJetMatchingSharedPtFraction(0.5),
517  fMCParticleArrayName("mcparticles"),
518  fNeedEmbedClusterContainer(0),
519  fRandomSeed(0),
520  fRandomSeedCones(0),
521  fVertexerCuts(0),
522  fSecondaryVertexMaxChi2(1e10),
523  fSecondaryVertexMaxDispersion(0.05),
524  fCustomStartupScript(),
525  fSetEmcalJetFlavour(0),
526  fSaveTrackPDGCode(kTRUE),
527  fJetTree(0),
528  fEventWeight(0.),
529  fImpactParameter(0.),
530  fMultiplicity(0),
531  fTriggerTracks_Pt(),
532  fTriggerTracks_Eta(),
533  fTriggerTracks_Phi(),
534  fMCParticleArray(0),
535  fRandomGenerator(0),
536  fRandomGeneratorCones(0),
537  fVtxTagger(0),
538  fIsEmbeddedEvent(kFALSE),
539  fDoDetLevelMatching(kFALSE),
540  fDoPartLevelMatching(kFALSE),
541  fSimpleSecVertices()
542 {
543  fRandomGenerator = new TRandom3();
544  fRandomGeneratorCones = new TRandom3();
545 
547  fJetTree = new AliEmcalJetTree(GetName());
548  DefineOutput(2, TTree::Class());
549 }
550 
551 //________________________________________________________________________
553  AliAnalysisTaskEmcalJet(name, kTRUE),
555  fEventPercentage(1.0),
561  fTruthMinLabel(0),
562  fTruthMaxLabel(100000),
564  fJetMatchingRadius(0.3),
566  fMCParticleArrayName("mcparticles"),
568  fRandomSeed(0),
569  fRandomSeedCones(0),
570  fVertexerCuts(0),
575  fSaveTrackPDGCode(kTRUE),
576  fJetTree(0),
577  fEventWeight(0.),
578  fImpactParameter(0.),
579  fMultiplicity(0),
583  fMCParticleArray(0),
584  fRandomGenerator(0),
586  fVtxTagger(0),
587  fIsEmbeddedEvent(kFALSE),
588  fDoDetLevelMatching(kFALSE),
589  fDoPartLevelMatching(kFALSE),
591 {
592  fRandomGenerator = new TRandom3();
593  fRandomGeneratorCones = new TRandom3();
594 
596  fJetTree = new AliEmcalJetTree(GetName());
597  DefineOutput(2, TTree::Class());
598 }
599 
600 //________________________________________________________________________
602 {
603  // Destructor.
604 }
605 
606 //________________________________________________________________________
608 {
610 
611  // ### Basic container settings
612  if(!GetJetContainer(0))
613  AliFatal("Jet input container not found!");
615  if(!GetParticleContainer(0))
616  AliFatal("At least one particle input container needs to be added");
618  AliFatal("Cluster input container not found although cluster extraction demanded!");
619 
620  fRandomGenerator->SetSeed(fRandomSeed);
622 
623 
624  // Activate saving of trigger tracks if this is demanded
626  fSaveTriggerTracks = kTRUE;
627 
628  // Use the jet containers from the task to activate matching
629  if (GetJetContainer(1)) {
630  fDoDetLevelMatching = kTRUE;
631  }
632  if (GetJetContainer(2)) {
633  fDoPartLevelMatching = kTRUE;
634  }
635  // ### Initialize the jet tree (settings must all be given at this stage)
638  OpenFile(2);
639  PostData(2, fJetTree->GetTreePointer());
640 
641  // ### Add control histograms (already some created in base task)
642  AddHistogram2D<TH2D>("hTrackCount", "Number of tracks in acceptance vs. centrality", "COLZ", 500, 0., 5000., 100, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}");
643  AddHistogram2D<TH2D>("hBackgroundPt", "Background p_{T} distribution", "", 150, 0., 150., 100, 0, 100, "Background p_{T} (GeV/c)", "Centrality", "dN^{Events}/dp_{T}");
644 
645  AddHistogram2D<TH2D>("hJetPtRaw", "Jets p_{T} distribution (raw)", "COLZ", 300, 0., 300., 100, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
646  AddHistogram2D<TH2D>("hJetPt", "Jets p_{T} distribution (background subtracted)", "COLZ", 400, -100., 300., 100, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
647  AddHistogram2D<TH2D>("hJetPtExtracted", "Extracted jets p_{T} distribution (background subtracted)", "COLZ", 400, -100., 300., 100, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
648  AddHistogram2D<TH2D>("hJetPhiEta", "Jet angular distribution #phi/#eta", "COLZ", 180, 0., 2*TMath::Pi(), 100, -2.5, 2.5, "#phi", "#eta", "dN^{Jets}/d#phi d#eta");
649  AddHistogram2D<TH2D>("hJetArea", "Jet area", "COLZ", 200, 0., 2., 100, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
650  AddHistogram2D<TH2D>("hConstituentPt", "Jet constituent p_{T} distribution", "COLZ", 400, 0., 300., 100, 0, 100, "p_{T, const} (GeV/c)", "Centrality", "dN^{Const}/dp_{T}");
651  AddHistogram2D<TH2D>("hConstituentPhiEta", "Jet constituent relative #phi/#eta distribution", "COLZ", 120, -0.6, 0.6, 120, -0.6, 0.6, "#Delta#phi", "#Delta#eta", "dN^{Const}/d#phi d#eta");
652  AddHistogram1D<TH1D>("hExtractionPercentage", "Extracted jets p_{T} distribution (background subtracted)", "e1p", 400, -100., 600., "p_{T, jet} (GeV/c)", "Extracted percentage");
653 
654  // Track QA plots
655  AddHistogram2D<TH2D>("hTrackPt", "Tracks p_{T} distribution", "", 300, 0., 300., 100, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
656  AddHistogram2D<TH2D>("hTrackPhi", "Track angular distribution in #phi", "LEGO2", 180, 0., 2*TMath::Pi(), 100, 0, 100, "#phi", "Centrality", "dN^{Tracks}/(d#phi)");
657  AddHistogram2D<TH2D>("hTrackEta", "Track angular distribution in #eta", "LEGO2", 100, -2.5, 2.5, 100, 0, 100, "#eta", "Centrality", "dN^{Tracks}/(d#eta)");
658  AddHistogram2D<TH2D>("hTrackPhiEta", "Track angular distribution #phi/#eta", "COLZ", 180, 0., 2*TMath::Pi(), 100, -2.5, 2.5, "#phi", "#eta", "dN^{Tracks}/d#phi d#eta");
659 
660  AddHistogram2D<TH2D>("hTrackEtaPt", "Track angular distribution in #eta vs. p_{T}", "LEGO2", 100, -2.5, 2.5, 300, 0., 300., "#eta", "p_{T} (GeV/c)", "dN^{Tracks}/(d#eta dp_{T})");
661  AddHistogram2D<TH2D>("hTrackPhiPt", "Track angular distribution in #phi vs. p_{T}", "LEGO2", 180, 0, 2*TMath::Pi(), 300, 0., 300., "#phi", "p_{T} (GeV/c)", "dN^{Tracks}/(d#phi dp_{T})");
662 }
663 
664 
665 //________________________________________________________________________
667 {
669 
670 
671  // ### If there is an embedding track container, set flag that an embedded event is used
672  for(Int_t iCont=0; iCont<fParticleCollArray.GetEntriesFast(); iCont++)
673  if (GetParticleContainer(iCont)->GetIsEmbedding())
674  fIsEmbeddedEvent = kTRUE;
675 
676  // ### Need to explicitly tell jet container to load rho mass object
677  GetJetContainer(0)->LoadRhoMass(InputEvent());
678 
679  if (fMCParticleArrayName != "")
680  {
681  if(!fIsEmbeddedEvent)
682  fMCParticleArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCParticleArrayName.Data()));
683  else
684  {
685  // In case of embedding, the MC particle array needs to be fetched differently
686  AliVEvent* event = AliEmcalContainerUtils::GetEvent(InputEvent(), kTRUE);
687  fMCParticleArray = dynamic_cast<TClonesArray*>(event->FindListObject(fMCParticleArrayName.Data()));
688  }
689  }
690 
691  // ### Prepare vertexer
693  {
694  if(!fVertexerCuts)
695  AliFatal("VertexerCuts not given but secondary vertex calculation turned on.");
696  fVtxTagger = new AliHFJetsTaggingVertex();
697  fVtxTagger->SetCuts(fVertexerCuts);
698  }
699 
700 
701 
702  // ### Save tree extraction percentage to histogram
703  std::vector<Float_t> extractionPtBins = fJetTree->GetExtractionPercentagePtBins();
704  std::vector<Float_t> extractionPercentages = fJetTree->GetExtractionPercentages();
705 
706  for(Int_t i=0; i<static_cast<Int_t>(extractionPercentages.size()); i++)
707  {
708  Double_t percentage = extractionPercentages[i];
709  for(Int_t pt=static_cast<Int_t>(extractionPtBins[i*2]); pt<static_cast<Int_t>(extractionPtBins[i*2+1]); pt++)
710  {
711  FillHistogram("hExtractionPercentage", pt, percentage);
712  }
713  }
714 
715  // ### Add HIJING/JEWEL histograms (in case we need it)
716  if(MCEvent())
717  {
718  if(dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()))
719  AddHistogram1D<TH1D>("hHIJING_Ncoll", "N_{coll} from HIJING", "e1p", 1000, 0., 5000, "N_{coll}", "dN^{Events}/dN^{N_{coll}}");
720 
721  if(dynamic_cast<AliGenHepMCEventHeader*>(MCEvent()->GenEventHeader()))
722  {
723  AddHistogram1D<TH1D>("hJEWEL_NProduced", "NProduced from JEWEL", "e1p", 1000, 0., 1000, "Nproduced", "dN^{Events}/dN^{Nproduced}");
724  AddHistogram1D<TH1D>("hJEWEL_ImpactParameter", "Impact parameter from JEWEL", "e1p", 1000, 0., 1, "IP", "dN^{Events}/dN^{IP}");
725  TProfile* evweight = new TProfile("hJEWEL_EventWeight", "Event weight from JEWEL", 1, 0,1);
726  fOutput->Add(evweight);
727  }
728  }
729 
730  // ### Execute shell script at startup
731  if(!fCustomStartupScript.IsNull())
732  {
733  TGrid::Connect("alien://");
734  TFile::Cp(fCustomStartupScript.Data(), "./myScript.sh");
735  gSystem->Exec("bash ./myScript.sh");
736  }
737 
738  PrintConfig();
739 
740 }
741 
742 //________________________________________________________________________
744 {
745  // ################################### EVENT SELECTION
746  // For debugging
747  //auto startTime = std::chrono::high_resolution_clock::now();
748  //auto elapsedTime = std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::high_resolution_clock::now() - startTime).count();
749  //std::cout << Form("%s, Time (start): %E s", GetName(), elapsedTime/100000000.) << std::endl;
750 
751  if(!IsTriggerTrackInEvent())
752  return kFALSE;
753  if(fRandomGenerator->Rndm() >= fEventPercentage)
754  return kFALSE;
755 
756  // ################################### EVENT PROPERTIES
758 
759  // Load vertex if possible
760  Long64_t eventID = 0;
761  const AliVVertex* myVertex = InputEvent()->GetPrimaryVertex();
762  if(!myVertex && MCEvent())
763  myVertex = MCEvent()->GetPrimaryVertex();
764  Double_t vtx[3] = {0, 0, 0};
765  if(myVertex)
766  {
767  vtx[0] = myVertex->GetX(); vtx[1] = myVertex->GetY(); vtx[2] = myVertex->GetZ();
768  }
769 
770  // Get event ID from header
771  AliVHeader* eventIDHeader = InputEvent()->GetHeader();
772  if(eventIDHeader)
773  eventID = eventIDHeader->GetEventIdAsLong();
774 
775  // If available, get Ncoll from HIJING
776  if(MCEvent())
777  {
778  if(AliGenHijingEventHeader* mcHeader = dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()))
779  {
780  Double_t ncoll = mcHeader->NN() + mcHeader->NNw() + mcHeader->NwN() + mcHeader->NwNw();
781  FillHistogram("hHIJING_Ncoll", ncoll);
782  }
783  if(AliGenHepMCEventHeader* mcHeader = dynamic_cast<AliGenHepMCEventHeader*>(MCEvent()->GenEventHeader()))
784  {
785  fEventWeight = mcHeader->EventWeight();
786  fImpactParameter = mcHeader->impact_parameter();
787 
788  TProfile* tmpHist = static_cast<TProfile*>(fOutput->FindObject("hJEWEL_EventWeight"));
789  tmpHist->Fill(0., fEventWeight);
790  FillHistogram("hJEWEL_NProduced", mcHeader->NProduced());
791  FillHistogram("hJEWEL_ImpactParameter", fImpactParameter);
792  }
793  }
794 
795  // Perform the matching before the main jet loop
797 
798  // ################################### MAIN JET LOOP
800  Int_t jetCount = 0;
801  while(AliEmcalJet *jet = GetJetContainer(0)->GetNextAcceptJet())
802  {
803  // LOCAL BUFFER (Will be set for each jet, only added to tree if jet is accepted )
804  std::vector<Float_t> vecSigITS; std::vector<Float_t> vecSigTPC; std::vector<Float_t> vecSigTRD; std::vector<Float_t> vecSigTOF; std::vector<Short_t> vecRecoPID; std::vector<Int_t> vecTruePID;
805  std::vector<Float_t> vec_d0; std::vector<Float_t> vec_d0cov; std::vector<Float_t> vec_z0; std::vector<Float_t> vec_z0cov;
806  std::vector<Float_t> secVtx_X; std::vector<Float_t> secVtx_Y; std::vector<Float_t> secVtx_Z; std::vector<Float_t> secVtx_Mass; std::vector<Float_t> secVtx_Lxy; std::vector<Float_t> secVtx_SigmaLxy; std::vector<Float_t> secVtx_Chi2; std::vector<Float_t> secVtx_Dispersion;
807  std::vector<Float_t> triggerTracks_dEta(fTriggerTracks_Eta);
808  std::vector<Float_t> triggerTracks_dPhi(fTriggerTracks_Phi);
809  std::vector<Float_t> splittings_radiatorE; std::vector<Float_t> splittings_kT; std::vector<Float_t> splittings_theta; std::vector<Int_t> splittings_secVtx_rank; std::vector<Int_t> splittings_secVtx_index;
810 
813  {
814  Double_t matchedJetDistance_Det = 0;
815  Double_t matchedJetPt_Det = 0;
816  Double_t matchedJetMass_Det = 0;
817  Double_t matchedJetAngularity_Det = 0;
818  Double_t matchedJetpTD_Det = 0;
819  Double_t matchedJetDistance_Part = 0;
820  Double_t matchedJetPt_Part = 0;
821  Double_t matchedJetMass_Part = 0;
822  Double_t matchedJetAngularity_Part = 0;
823  Double_t matchedJetpTD_Part = 0;
824  Double_t truePtFraction = 0;
825  Double_t truePtFraction_PartLevel = 0;
826  Int_t currentJetType_HM = 0;
827  Int_t currentJetType_PM = 0;
828  Int_t currentJetType_IC = 0;
829 
830  // Get jet type from MC (hadron matching, parton matching definition - for HF jets)
831  GetJetType(jet, currentJetType_HM, currentJetType_PM, currentJetType_IC);
832  // Get true estimators: for pt, jet mass, ...
833  GetTrueJetPtFraction(jet, truePtFraction, truePtFraction_PartLevel);
834 
835  GetMatchedJetObservables(jet, matchedJetPt_Det, matchedJetPt_Part, matchedJetDistance_Det, matchedJetDistance_Part, matchedJetMass_Det, matchedJetMass_Part, matchedJetAngularity_Det, matchedJetAngularity_Part, matchedJetpTD_Det,matchedJetpTD_Part);
836 
837 
838  fJetTree->FillBuffer_MonteCarlo(currentJetType_PM,currentJetType_HM,currentJetType_IC,
839  matchedJetDistance_Det,matchedJetPt_Det,matchedJetMass_Det,matchedJetAngularity_Det,matchedJetpTD_Det,
840  matchedJetDistance_Part,matchedJetPt_Part,matchedJetMass_Part,matchedJetAngularity_Part,matchedJetpTD_Part,
841  truePtFraction,truePtFraction_PartLevel,fPtHard,fEventWeight,fImpactParameter);
842  }
843 
844  // ### CONSTITUENT LOOP: Retrieve PID values + impact parameters
846  {
847  for(Int_t i = 0; i < jet->GetNumberOfParticleConstituents(); i++)
848  {
849  const AliVParticle* particle = jet->GetParticleConstituents()[i].GetParticle();
850  if(!particle) continue;
852  {
853  Float_t sigITS = 0; Float_t sigTPC = 0; Float_t sigTOF = 0; Float_t sigTRD = 0; Short_t recoPID = 0; Int_t truePID = 0;
854  AddPIDInformation(const_cast<AliVParticle*>(particle), sigITS, sigTPC, sigTOF, sigTRD, recoPID, truePID);
855  vecSigITS.push_back(sigITS); vecSigTPC.push_back(sigTPC); vecSigTOF.push_back(sigTOF); vecSigTRD.push_back(sigTRD); vecRecoPID.push_back(recoPID); vecTruePID.push_back(truePID);
856  fJetTree->FillBuffer_PID(vecSigITS, vecSigTPC, vecSigTOF, vecSigTRD, vecRecoPID, vecTruePID);
857  }
859  {
860  Float_t d0 = 0; Float_t d0cov = 0; Float_t z0 = 0; Float_t z0cov = 0;
861  GetTrackImpactParameters(myVertex, dynamic_cast<const AliAODTrack*>(particle), d0, d0cov, z0, z0cov);
862  vec_d0.push_back(d0); vec_d0cov.push_back(d0cov); vec_z0.push_back(z0); vec_z0cov.push_back(z0cov);
863  fJetTree->FillBuffer_ImpactParameters(vec_d0, vec_z0, vec_d0cov, vec_z0cov);
864  }
865  }
866  }
867 
868  // Reconstruct secondary vertices
870  {
871  ReconstructSecondaryVertices(myVertex, jet, secVtx_X, secVtx_Y, secVtx_Z, secVtx_Mass, secVtx_Lxy, secVtx_SigmaLxy, secVtx_Chi2, secVtx_Dispersion);
872  fJetTree->FillBuffer_SecVertices(secVtx_X, secVtx_Y, secVtx_Z, secVtx_Mass, secVtx_Lxy, secVtx_SigmaLxy, secVtx_Chi2, secVtx_Dispersion);
873  }
874 
875  // Now change trigger tracks eta/phi to dEta/dPhi relative to jet and save them
877  {
878  for(UInt_t i=0; i<triggerTracks_dEta.size(); i++)
879  {
880  triggerTracks_dEta[i] = jet->Eta() - fTriggerTracks_Eta[i];
881  triggerTracks_dPhi[i] = TMath::Min(TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]), TMath::TwoPi() - TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]));
882  if( ((TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]) <= TMath::Pi()) && (jet->Phi() - fTriggerTracks_Phi[i] > 0))
883  || ((TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]) > TMath::Pi()) && (jet->Phi() - fTriggerTracks_Phi[i] < 0)) )
884  triggerTracks_dPhi[i] = -triggerTracks_dPhi[i];
885  }
886  fJetTree->FillBuffer_TriggerTracks(fTriggerTracks_Pt, triggerTracks_dEta, triggerTracks_dPhi);
887  }
888 
889  if(fSaveJetShapes)
890  {
891  // Calculate jet shapes and set them in the tree (some are retrieved in the tree itself)
892  Double_t leSub_noCorr = 0;
893  Double_t angularity = 0;
894  Double_t momentumDispersion = 0;
895  Double_t trackPtMean = 0;
896  Double_t trackPtMedian = 0;
897  CalculateJetShapes(jet, leSub_noCorr, angularity, momentumDispersion, trackPtMean, trackPtMedian);
898  fJetTree->FillBuffer_JetShapes(jet, leSub_noCorr, angularity, momentumDispersion, trackPtMean, trackPtMedian);
899  }
900 
902  {
903  GetJetSplittings(jet, splittings_radiatorE, splittings_kT, splittings_theta, splittings_secVtx_rank, splittings_secVtx_index);
904  fJetTree->FillBuffer_Splittings(splittings_radiatorE, splittings_kT, splittings_theta, fSaveSecondaryVertices, splittings_secVtx_rank, splittings_secVtx_index);
905  }
906 
907  // Fill jet to tree (here adding the minimum properties)
908  Bool_t accepted = fJetTree->AddJetToTree(jet, fSaveConstituents, fSaveConstituentsIP, fSaveCaloClusters, vtx, GetJetContainer(0)->GetRhoVal(), GetJetContainer(0)->GetRhoMassVal(), fCent, fMultiplicity, eventID, InputEvent()->GetMagneticField());
909  if(accepted)
910  FillHistogram("hJetPtExtracted", jet->Pt() - GetJetContainer(0)->GetRhoVal()*jet->Area(), fCent);
911  jetCount++;
912  }
913 
914  // ################################### PARTICLE PROPERTIES
915  for(Int_t iCont=0; iCont<fParticleCollArray.GetEntriesFast(); iCont++)
916  {
918  while(AliVTrack *track = static_cast<AliVTrack*>(GetParticleContainer(iCont)->GetNextAcceptParticle()))
920  }
921 
922  // For debugging
923  //elapsedTime = std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::high_resolution_clock::now() - startTime).count();
924  //std::cout << Form("Time (end): %E s", elapsedTime/100000000.) << std::endl;
925  return kTRUE;
926 }
927 
928 
929 //________________________________________________________________________
931 {
932  // Cut for trigger track requirement
934  {
935  // Clear vector of trigger tracks
936  fTriggerTracks_Pt.clear();
937  fTriggerTracks_Eta.clear();
938  fTriggerTracks_Phi.clear();
939 
940  // Go through all tracks, in all attached track containers, and check whether trigger tracks can be found
941  for(Int_t iCont=0; iCont<fParticleCollArray.GetEntriesFast(); iCont++)
942  {
943  AliParticleContainer* partCont = GetParticleContainer(iCont);
944  partCont->ResetCurrentID();
945  if((fEventCut_TriggerTrackOrigin == 1) && !partCont->GetIsEmbedding())
946  continue;
947  else if((fEventCut_TriggerTrackOrigin == 2) && partCont->GetIsEmbedding())
948  continue;
949 
950  while(AliVTrack *track = static_cast<AliVTrack*>(partCont->GetNextAcceptParticle()))
951  {
952  if( (track->GetLabel() >= fEventCut_TriggerTrackMinLabel) && (track->GetLabel() < fEventCut_TriggerTrackMaxLabel) )
953  if( (track->Pt() >= fEventCut_TriggerTrackMinPt) && (track->Pt() < fEventCut_TriggerTrackMaxPt) )
954  {
955  fTriggerTracks_Pt.push_back(track->Pt());
956  fTriggerTracks_Eta.push_back(track->Eta());
957  fTriggerTracks_Phi.push_back(track->Phi());
958  }
959  }
960  }
961  // No particle has been found that fulfills requirement -> Do not accept event
962  if(fTriggerTracks_Pt.size() == 0)
963  return kFALSE;
964  }
965  return kTRUE;
966 }
967 
968 //________________________________________________________________________
969 void AliAnalysisTaskJetExtractor::CalculateJetShapes(AliEmcalJet* jet, Double_t& leSub_noCorr, Double_t& angularity, Double_t& momentumDispersion, Double_t& trackPtMean, Double_t& trackPtMedian)
970 {
971  // #### Calculate mean, median of constituents, radial moment (angularity), momentum dispersion, leSub (no correction)
972  Double_t jetLeadingHadronPt = -999.;
973  Double_t jetSubleadingHadronPt = -999.;
974  Double_t jetSummedPt = 0;
975  Double_t jetSummedPt2 = 0;
976  trackPtMean = 0;
977  trackPtMedian = 0;
978  angularity = 0;
979  momentumDispersion = 0;
980  std::vector<PWG::JETFW::AliEmcalParticleJetConstituent> tracks_sorted = jet->GetParticleConstituents();
981  std::sort(tracks_sorted.rbegin(), tracks_sorted.rend());
982  Int_t numTracks = tracks_sorted.size();
983  if(!numTracks) return;
984  Double_t* constPts = new Double_t[numTracks];
985 
986  // Loop over all constituents and do jet shape calculations
987  for (Int_t i=0;i<numTracks;i++)
988  {
989  const AliVParticle* particle = tracks_sorted[i].GetParticle();
990  trackPtMean += particle->Pt();
991  constPts[i] = particle->Pt();
992  if(particle->Pt() > jetLeadingHadronPt)
993  {
994  jetSubleadingHadronPt = jetLeadingHadronPt;
995  jetLeadingHadronPt = particle->Pt();
996  }
997  else if(particle->Pt() > jetSubleadingHadronPt)
998  jetSubleadingHadronPt = particle->Pt();
999 
1000  Double_t deltaR = GetDistance(particle->Eta(), jet->Eta(), particle->Phi(), jet->Phi());
1001  jetSummedPt += particle->Pt();
1002  jetSummedPt2 += particle->Pt()*particle->Pt();
1003  angularity += particle->Pt() * deltaR;
1004  }
1005 
1006  if(numTracks)
1007  {
1008  trackPtMean /= numTracks;
1009  trackPtMedian = TMath::Median(numTracks, constPts);
1010  }
1011 
1012  if(numTracks > 1)
1013  leSub_noCorr = jetLeadingHadronPt - jetSubleadingHadronPt;
1014  else
1015  leSub_noCorr = jetLeadingHadronPt;
1016 
1017  if(jetSummedPt)
1018  {
1019  momentumDispersion = TMath::Sqrt(jetSummedPt2)/jetSummedPt;
1020  angularity /= jetSummedPt;
1021  }
1022 }
1023 
1024 //________________________________________________________________________
1025 void AliAnalysisTaskJetExtractor::GetTrueJetPtFraction(AliEmcalJet* jet, Double_t& truePtFraction, Double_t& truePtFraction_mcparticles)
1026 {
1027  // #################################################################################
1028  // ##### FRACTION OF TRUE PT IN JET: Defined as "not from toy"
1029  Double_t pt_truth = 0.;
1030  Double_t pt_truth_mcparticles = 0.;
1031  Double_t pt_all = 0.;
1032  truePtFraction = 0;
1033  truePtFraction_mcparticles = 0;
1034 
1035  // ### Loop over all tracks constituents
1036  for(Int_t iConst = 0; iConst < jet->GetNumberOfParticleConstituents(); iConst++)
1037  {
1038  const AliVParticle* particle = jet->GetParticleConstituents()[iConst].GetParticle();
1039  if(!particle) continue;
1040  if(particle->Pt() < 1e-6) continue;
1041 
1042  // Particles marked w/ labels within label range OR explicitly set as embedded tracks are considered to be from truth
1043  if ( (fIsEmbeddedEvent && jet->GetParticleConstituents()[iConst].IsFromEmbeddedEvent()) ||
1044  (!fIsEmbeddedEvent && ((particle->GetLabel() >= fTruthMinLabel) && (particle->GetLabel() < fTruthMaxLabel))) ){
1045  pt_truth += particle->Pt();
1046  }
1047  pt_all += particle->Pt();
1048  }
1049  // Get the Primary Vertex
1050  Double_t vtxX = 0;
1051  Double_t vtxY = 0;
1052  Double_t vtxZ = 0;
1053  const AliVVertex* myVertex = InputEvent()->GetPrimaryVertex();
1054  if(!myVertex && MCEvent())
1055  myVertex = MCEvent()->GetPrimaryVertex();
1056  if(myVertex)
1057  {
1058  vtxX = myVertex->GetX();
1059  vtxY = myVertex->GetY();
1060  vtxZ = myVertex->GetZ();
1061  }
1062 
1063  Double_t primVtx[3] = {vtxX, vtxY, vtxZ};
1064 
1065  // ### Loop over all cluster constituents
1066  for(Int_t iConst = 0; iConst < jet->GetNumberOfClusterConstituents(); iConst++)
1067  {
1068  const AliVCluster* cluster = jet->GetClusterConstituents()[iConst].GetCluster();
1069  if(!cluster) continue;
1070  if(cluster->E() < 1e-6) continue;
1071 
1072  // #### Retrieve cluster pT
1073  TLorentzVector clusterMomentum;
1074  cluster->GetMomentum(clusterMomentum, primVtx);
1075  // Get the total pT from clusters in the jet
1076  Double_t ClusterPt = clusterMomentum.Perp();
1077  pt_all += ClusterPt;
1078  } // end loop over clusters
1079 
1080  // Calculate the true pt fraction for the case where we have clusters
1081  if (GetClusterContainer(0)){
1082  // get the cluster container from the input event corresponding
1083  AliClusterContainer* clusterCont = 0;
1085  // if hybrid event, take external cluster container
1086  clusterCont = GetClusterContainer(1);
1087  }
1088  else{
1089  clusterCont = GetClusterContainer(0);
1090  }
1091  // loop over clusters in the input event
1092  // check to make sure that we are getting a cluster container
1093  if(clusterCont){
1094  for(int k=0; k< clusterCont->GetNClusters(); k++){
1095  const AliVCluster* cluster = clusterCont->GetAcceptCluster(k);
1096  if (!cluster)continue;
1097  TLorentzVector clusterMomentum;
1098  cluster->GetMomentum(clusterMomentum, primVtx);
1099  Double_t ClusterPt = clusterMomentum.Perp();
1100  if(IsClusterInCone(clusterMomentum, jet->Eta(), jet->Phi(), GetJetContainer(0)->GetJetRadius()) ){
1101  pt_truth += ClusterPt;
1102  }
1103  } // end loop over clusters
1104  }// end if clusterCont
1105  }
1106 
1107  // ### Loop over all primary (charged) MC particles and check if they have a corresponding track/cluster
1108  // Correspondence is checked geometrically, sum of matched particles pT is truth
1110  Double_t jetRadius = GetJetContainer(0)->GetJetRadius();
1111  if(fMCParticleArray)
1112  for(Int_t iPart=0; iPart<fMCParticleArray->GetEntriesFast();iPart++)
1113  {
1114  AliAODMCParticle* part = (AliAODMCParticle*)fMCParticleArray->At(iPart);
1115  if(!part) continue;
1116  if(!part->IsPhysicalPrimary()) continue;
1117  if(!fulljets && !part->Charge()) continue;
1118  if(part->Pt() < 1e-6) continue;
1119 
1120  if(IsTrackInCone(part, jet->Eta(), jet->Phi(), jetRadius))
1121  pt_truth_mcparticles += part->Pt();
1122  }
1123 
1124  if(pt_all)
1125  {
1126  truePtFraction = (pt_truth/pt_all);
1127  truePtFraction_mcparticles = (pt_truth_mcparticles/pt_all);
1128  }
1129 }
1130 //________________________________________________________________________
1132  // Perform the matching before the main jet loop (Only if we use the tagger for matching)
1133  AliJetContainer * jetsHybrid = GetJetContainer(0);
1134  AliJetContainer * jetsDetLevel = GetJetContainer(1);
1135  AliJetContainer * jetsPartLevel = GetJetContainer(2);
1136 
1137 
1138  // Now, begin the actual matching.
1139  // Hybrid <-> det first
1140  AliDebugStream(2) << "Matching hybrid to detector level jets.\n";
1141  // First, we reset the tagging
1142  for(auto j : jetsHybrid->all()){
1143  j->ResetMatching();
1144  }
1145  for(auto j : jetsDetLevel->all()){
1146  j->ResetMatching();
1147  }
1148  // Next, we perform the matching
1149  PerformGeometricalJetMatching(*jetsHybrid, *jetsDetLevel, fJetMatchingRadius);
1150  // Now, begin the next matching stage
1151  // det <-> particle
1152  AliDebugStream(2) << "Matching detector level to particle level jets.\n";
1153  // First, we reset the tagging. We need to reset the det matching again to ensure
1154  // that it doesn't accidentally keep some latent matches to the hybrid jets.
1155  for(auto j : jetsDetLevel->all()){
1156  j->ResetMatching();
1157  }
1158  // if we do not need to do particle level matching return
1159  if(!fDoPartLevelMatching) return;
1160 
1161  for(auto j : jetsPartLevel->all()){
1162  j->ResetMatching();
1163  }
1164  // Next, we perform the matching
1165  PerformGeometricalJetMatching(*jetsHybrid, *jetsDetLevel, fJetMatchingRadius);
1166  // Now, begin the next matching stage
1167  // det <-> particle
1168  AliDebugStream(2) << "Matching detector level to particle level jets.\n";
1169  // First, we reset the tagging. We need to reset the det matching again to ensure
1170  // that it doesn't accidentally keep some latent matches to the hybrid jets.
1171  for(auto j : jetsDetLevel->all()){
1172  j->ResetMatching();
1173  }
1174  for(auto j : jetsPartLevel->all()){
1175  j->ResetMatching();
1176  }
1177  // Next, we perform the matching
1178  PerformGeometricalJetMatching(*jetsDetLevel, *jetsPartLevel, fJetMatchingRadius);
1179 }
1180 
1181 //________________________________________________________________________
1183  AliJetContainer& contTag, double maxDist)
1184 {
1185  // Note that this function is also utilized in /PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetHPerformance.cxx. For more details, see this file.
1186  // Setup
1187  const Int_t kNacceptedBase = contBase.GetNAcceptedJets(), kNacceptedTag = contTag.GetNAcceptedJets();
1188  if (!(kNacceptedBase && kNacceptedTag)) {
1189  return false;
1190  }
1191 
1192  // Build up vectors of jet pointers to use when assigning the closest jets.
1193  // The storages are needed later for applying the tagging, in order to avoid multiple occurrence of jet selection
1194  std::vector<AliEmcalJet*> jetsBase(kNacceptedBase), jetsTag(kNacceptedTag);
1195 
1196  int countBase(0), countTag(0);
1197  for (auto jb : contBase.accepted()) {
1198  jetsBase[countBase] = jb;
1199  countBase++;
1200  }
1201  for (auto jt : contTag.accepted()) {
1202  jetsTag[countTag] = jt;
1203  countTag++;
1204  }
1205 
1206  TArrayI faMatchIndexTag(kNacceptedBase), faMatchIndexBase(kNacceptedTag);
1207  faMatchIndexBase.Reset(-1);
1208  faMatchIndexTag.Reset(-1);
1209 
1210  // find the closest distance to the base jet
1211  countBase = 0;
1212  for (auto jet1 : contBase.accepted()) {
1213  double distance = maxDist;
1214 
1215  // Loop over all accepted jets and brute force search for the closest jet.
1216  // NOTE: current_index() returns the jet index in the underlying array, not
1217  // the index within the accepted jets that are returned.
1218  int contTagAcceptedIndex = 0;
1219  for (auto jet2 : contTag.accepted()) {
1220  double dR = jet1->DeltaR(jet2);
1221  if (dR < distance && dR < maxDist) {
1222  faMatchIndexTag[countBase] = contTagAcceptedIndex;
1223  distance = dR;
1224  }
1225  contTagAcceptedIndex++;
1226  }
1227 
1228 
1229  countBase++;
1230  }
1231 
1232  // other way around
1233  countTag = 0;
1234  for (auto jet1 : contTag.accepted()) {
1235  double distance = maxDist;
1236 
1237  // Loop over all accepted jets and brute force search for the closest jet.
1238  // NOTE: current_index() returns the jet index in the underlying array, not
1239  // the index within the accepted jets that are returned.
1240  int contBaseAcceptedIndex = 0;
1241  for (auto jet2 : contBase.accepted()) {
1242  double dR = jet1->DeltaR(jet2);
1243  if (dR < distance && dR < maxDist) {
1244  faMatchIndexBase[countTag] = contBaseAcceptedIndex;
1245  distance = dR;
1246  }
1247  contBaseAcceptedIndex++;
1248  }
1249  countTag++;
1250  }
1251 
1252  // check for "true" correlations
1253  // these are pairs where the base jet is the closest to the tag jet and vice versa
1254  // As the lists are linear a loop over the outer base jet is sufficient.
1255  AliDebugStream(1) << "Starting true jet loop: nbase(" << kNacceptedBase << "), ntag(" << kNacceptedTag << ")\n";
1256  for (int ibase = 0; ibase < kNacceptedBase; ibase++) {
1257  AliDebugStream(2) << "base jet " << ibase << ": match index in tag jet container " << faMatchIndexTag[ibase]
1258  << "\n";
1259  if (faMatchIndexTag[ibase] > -1) {
1260  AliDebugStream(2) << "tag jet " << faMatchIndexTag[ibase] << ": matched base jet " << faMatchIndexBase[faMatchIndexTag[ibase]] << "\n";
1261  }
1262  // We have a true correlation where each jet points to the other.
1263  if (faMatchIndexTag[ibase] > -1 && faMatchIndexBase[faMatchIndexTag[ibase]] == ibase) {
1264  AliDebugStream(2) << "found a true match \n";
1265  AliEmcalJet *jetBase = jetsBase[ibase], *jetTag = jetsTag[faMatchIndexTag[ibase]];
1266  // We have a valid pair of matched jets, so set the closest jet properties.
1267  if (jetBase && jetTag) {
1268  Double_t dR = jetBase->DeltaR(jetTag);
1269  jetBase->SetClosestJet(jetTag, dR);
1270  jetTag->SetClosestJet(jetBase, dR);
1271  }
1272  }
1273  }
1274  return true;
1275 }
1276 //________________________________________________________________________
1277 void AliAnalysisTaskJetExtractor::GetMatchedJetObservables(AliEmcalJet* jet, Double_t& detJetPt, Double_t& partJetPt, Double_t& detJetDistance, Double_t& partJetDistance, Double_t& detJetMass, Double_t& partJetMass, Double_t& detJetAngularity, Double_t& partJetAngularity, Double_t& detJetpTD, Double_t& partJetpTD)
1278 {
1279 
1280 
1281  // Get the Matched Observables
1282  AliJetContainer * hybridJetCont = GetJetContainer(0);
1283  AliJetContainer * detJetCont = GetJetContainer(1);
1284  AliJetContainer * partJetCont = GetJetContainer(2);
1285 
1286  // hybrid to detector level matching
1287  AliEmcalJet * jet2 = jet->ClosestJet();
1288  if (!jet2) { return; }// if there is no match return.
1289  Double_t ptJet2 = jet2->Pt() - detJetCont->GetRhoVal() * jet2->Area();
1290  // This will retrieve the fraction of jet2's momentum in jet1.
1291  Double_t fraction = hybridJetCont->GetFractionSharedPt(jet);
1292  if (fraction < fJetMatchingSharedPtFraction) {
1293  return;
1294  }
1295 
1296  // Set temp jet shape parameters
1297  Double_t leSub_noCorr = 0;
1298  Double_t trackPtMean = 0;
1299  Double_t trackPtMedian = 0;
1300 
1301  // if we are not doing the particle level match, fill all observables here and return
1302  if(!fDoPartLevelMatching){
1303  detJetPt = ptJet2;
1304  detJetDistance = jet->DeltaR(jet2);
1305  detJetMass = jet2->M();
1306 
1307  // Get Jet Shape parameters
1308  detJetAngularity = 0;
1309  detJetpTD = 0;
1310  CalculateJetShapes(jet2, leSub_noCorr, detJetAngularity, detJetpTD, trackPtMean, trackPtMedian);
1311  return;
1312  }
1313 
1314  // detector to particle matching
1315  AliEmcalJet * jet3 = jet2->ClosestJet();
1316  if(!jet3){return;}
1317  Double_t ptJet3 = jet3->Pt() - partJetCont->GetRhoVal() * jet3->Area();
1318  // make no fraction cut on the detector to particle
1319 
1320  // if we are doing particle level matching, require that there is both a particle and detector level match
1321  detJetPt = ptJet2;
1322  partJetPt = ptJet3;
1323  detJetDistance = jet->DeltaR(jet2);
1324  partJetDistance = jet2->DeltaR(jet3);
1325  detJetMass = jet2->M();
1326  partJetMass = jet3->M();
1327 
1328  // Get Jet Shape Parameters
1329  detJetAngularity = 0;
1330  detJetpTD = 0;
1331  partJetAngularity = 0;
1332  partJetpTD = 0;
1333 
1334  CalculateJetShapes(jet2, leSub_noCorr, detJetAngularity, detJetpTD, trackPtMean, trackPtMedian);
1335 
1336  // reset the temp variables
1337  leSub_noCorr = 0;
1338  trackPtMean = 0;
1339  trackPtMedian = 0;
1340 
1341  CalculateJetShapes(jet3, leSub_noCorr, partJetAngularity, partJetpTD, trackPtMean, trackPtMedian);
1342 }
1343 
1344 
1345 
1346 
1347 //________________________________________________________________________
1349 {
1350  if(!fMCParticleArray)
1351  return;
1352 
1353  typeHM = 0;
1354  typePM = 0;
1355 
1356  AliAODMCParticle* parton[2];
1357  parton[0] = (AliAODMCParticle*) fVtxTagger->IsMCJetParton(fMCParticleArray, jet, fHadronMatchingRadius); // method 1 (parton matching)
1358  parton[1] = (AliAODMCParticle*) fVtxTagger->IsMCJetMeson(fMCParticleArray, jet, fHadronMatchingRadius); // method 2 (hadron matching)
1359 
1360  if (parton[0]) {
1361  Int_t pdg = TMath::Abs(parton[0]->PdgCode());
1362  typePM = pdg;
1363  }
1364 
1365  if (!parton[1])
1366  {
1367  // No HF jet (according to hadron matching) -- now also separate jets in udg (1) and s-jets (3)
1368  if(IsStrangeJet(jet))
1369  typeHM = 3;
1370  else
1371  typeHM = 1;
1372  }
1373  else {
1374  Int_t pdg = TMath::Abs(parton[1]->PdgCode());
1375  if(fVtxTagger->IsDMeson(pdg)) typeHM = 4;
1376  else if (fVtxTagger->IsBMeson(pdg)) typeHM = 5;
1377  }
1378 
1379  // Set flavour of AliEmcalJet object (set ith bit while i corresponds to type)
1381  jet->AddFlavourTag(static_cast<Int_t>(TMath::Power(2, typeHM)));
1382 
1383  const AliEmcalPythiaInfo* partonsInfo = GetPythiaInfo();
1384  typeIC = 0;
1385  if (partonsInfo)
1386  {
1387  // Get primary partons directions
1388  Double_t parton1phi = partonsInfo->GetPartonPhi6();
1389  Double_t parton1eta = partonsInfo->GetPartonEta6();
1390  Double_t parton2phi = partonsInfo->GetPartonPhi7();
1391  Double_t parton2eta = partonsInfo->GetPartonEta7();
1392 
1393  Double_t delta1R = GetDistance(parton1eta, jet->Eta(), parton1phi, jet->Phi());
1394  Double_t delta2R = GetDistance(parton2eta, jet->Eta(), parton2phi, jet->Phi());
1395 
1396  // Check if one of the partons if closer than matching criterion
1397  Bool_t matched = (delta1R < GetJetContainer(0)->GetJetRadius()/2.) || (delta2R < GetJetContainer(0)->GetJetRadius()/2.);
1398 
1399  // Matching criterion fulfilled -> Set flag to closest
1400  if(matched)
1401  {
1402  if(delta1R < delta2R)
1403  typeIC = partonsInfo->GetPartonFlag6();
1404  else
1405  typeIC = partonsInfo->GetPartonFlag7();
1406  }
1407  }
1408 }
1409 
1410 //________________________________________________________________________
1411 void AliAnalysisTaskJetExtractor::GetTrackImpactParameters(const AliVVertex* vtx, const AliAODTrack* track, Float_t& d0, Float_t& d0cov, Float_t& z0, Float_t& z0cov)
1412 {
1413  if (track)
1414  {
1415  AliAODTrack localTrack(*track);
1416  Double_t d0rphiz[2],covd0[3];
1417  Bool_t isDCA=localTrack.PropagateToDCA(vtx,InputEvent()->GetMagneticField(),3.0,d0rphiz,covd0);
1418  if(isDCA)
1419  {
1420  d0 = d0rphiz[0];
1421  z0 = d0rphiz[1];
1422  d0cov = covd0[0];
1423  z0cov = covd0[2];
1424  }
1425  }
1426 }
1427 
1428 //________________________________________________________________________
1429 void AliAnalysisTaskJetExtractor::ReconstructSecondaryVertices(const AliVVertex* primVtx, const AliEmcalJet* jet, std::vector<Float_t>& secVtx_X, std::vector<Float_t>& secVtx_Y, std::vector<Float_t>& secVtx_Z, std::vector<Float_t>& secVtx_Mass, std::vector<Float_t>& secVtx_Lxy, std::vector<Float_t>& secVtx_SigmaLxy, std::vector<Float_t>& secVtx_Chi2, std::vector<Float_t>& secVtx_Dispersion)
1430 {
1431  if(!primVtx)
1432  return;
1433 
1434  // Create ESD vertex from the existing AliVVertex
1435  Double_t vtxPos[3] = {primVtx->GetX(), primVtx->GetY(), primVtx->GetZ()};
1436  Double_t covMatrix[6] = {0};
1437  primVtx->GetCovarianceMatrix(covMatrix);
1438  AliESDVertex* esdVtx = new AliESDVertex(vtxPos, covMatrix, primVtx->GetChi2(), primVtx->GetNContributors());
1439 
1440  TClonesArray* secVertexArr = 0;
1441  vctr_pair_dbl_int arrDispersion;
1442  arrDispersion.reserve(5);
1443  //###########################################################################
1444  // ********* Calculate secondary vertices
1445  // Derived from AliAnalysisTaskEmcalJetBtagSV
1446  secVertexArr = new TClonesArray("AliAODVertex");
1447  Int_t nDauRejCount = 0;
1448  Int_t nVtx = fVtxTagger->FindVertices(jet,
1449  GetParticleContainer(0)->GetArray(),
1450  (AliAODEvent*)InputEvent(),
1451  esdVtx,
1452  InputEvent()->GetMagneticField(),
1453  secVertexArr,
1454  0,
1455  arrDispersion,
1456  nDauRejCount);
1457 
1458  if(nVtx < 0)
1459  {
1460  secVertexArr->Clear();
1461  delete secVertexArr;
1462  return;
1463  }
1464  //###########################################################################
1465 
1466  // Loop over all potential secondary vertices
1467  fSimpleSecVertices.clear();
1468  for(Int_t i=0; i<secVertexArr->GetEntriesFast(); i++)
1469  {
1470  AliAODVertex* secVtx = (AliAODVertex*)(secVertexArr->UncheckedAt(i));
1471 
1472  // Calculate vtx distance
1473  Double_t effX = secVtx->GetX() - esdVtx->GetX();
1474  Double_t effY = secVtx->GetY() - esdVtx->GetY();
1475  //Double_t effZ = secVtx->GetZ() - esdVtx->GetZ();
1476 
1477  // ##### Vertex properties
1478  // vertex dispersion
1479  Double_t dispersion = arrDispersion[i].first;
1480 
1481  // invariant mass
1482  Double_t mass = fVtxTagger->GetVertexInvariantMass(secVtx);
1483 
1484  // signed length
1485  Double_t Lxy = TMath::Sqrt(effX*effX + effY*effY);
1486  Double_t jetP[3]; jet->PxPyPz(jetP);
1487  Double_t signLxy = effX * jetP[0] + effY * jetP[1];
1488  if (signLxy < 0.) Lxy *= -1.;
1489 
1490  Double_t sigmaLxy = 0;
1491  AliAODVertex* aodVtx = (AliAODVertex*)(primVtx);
1492  if (aodVtx)
1493  sigmaLxy = aodVtx->ErrorDistanceXYToVertex(secVtx);
1494 
1495  // Add secondary vertices if they fulfill the conditions
1496  if( (dispersion > fSecondaryVertexMaxDispersion) || (TMath::Abs(secVtx->GetChi2perNDF()) > fSecondaryVertexMaxChi2) )
1497  continue;
1498 
1499  // Internally, save sec. vertices to a list
1500  // Each secondary vertex is reconstructed from 3 prongs
1502  vtx.fIndex = secVtx_X.size();
1503  vtx.fLxy = TMath::Abs(Lxy);
1504  vtx.fDaughter1 = static_cast<AliVParticle*>(secVtx->GetDaughter(0));
1505  vtx.fDaughter2 = static_cast<AliVParticle*>(secVtx->GetDaughter(1));
1506  vtx.fDaughter3 = static_cast<AliVParticle*>(secVtx->GetDaughter(2));
1507  fSimpleSecVertices.push_back(vtx);
1508 
1509  secVtx_X.push_back(secVtx->GetX()); secVtx_Y.push_back(secVtx->GetY()); secVtx_Z.push_back(secVtx->GetZ()); secVtx_Chi2.push_back(secVtx->GetChi2perNDF());
1510  secVtx_Dispersion.push_back(dispersion); secVtx_Mass.push_back(mass); secVtx_Lxy.push_back(Lxy); secVtx_SigmaLxy.push_back(sigmaLxy);
1511  }
1512 
1513  // Sort simple sec. vertices w/ descending Lxy
1514  std::sort(fSimpleSecVertices.rbegin(), fSimpleSecVertices.rend(), [](const SimpleSecondaryVertex a, const SimpleSecondaryVertex b) {return a.fLxy<b.fLxy;});
1515  if(fSimpleSecVertices.size() > 10) fSimpleSecVertices.resize(10);
1516 
1517  secVertexArr->Clear();
1518  delete secVertexArr;
1519  delete esdVtx;
1520 }
1521 
1522 //________________________________________________________________________
1524 {
1525  // Do hadron matching jet type tagging using mcparticles
1526  // ... if not explicitly deactivated
1527  if (fMCParticleArray)
1528  {
1529  for(Int_t i=0; i<fMCParticleArray->GetEntriesFast();i++)
1530  {
1531  AliAODMCParticle* part = (AliAODMCParticle*)fMCParticleArray->At(i);
1532  if(!part) continue;
1533 
1534  // Check if the particle has strangeness
1535  Int_t absPDG = TMath::Abs(part->PdgCode());
1536  if ((absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
1537  {
1538  // Check if particle is in a radius around the jet
1539  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
1541  continue;
1542  else
1543  return kTRUE;
1544  }
1545  }
1546  }
1547  // As fallback, the MC stack will be tried
1548  else if(MCEvent() && (MCEvent()->Stack()))
1549  {
1550  AliStack* stack = MCEvent()->Stack();
1551  // Go through the whole particle stack
1552  for(Int_t i=0; i<stack->GetNtrack(); i++)
1553  {
1554  TParticle *part = stack->Particle(i);
1555  if(!part) continue;
1556 
1557  // Check if the particle has strangeness
1558  Int_t absPDG = TMath::Abs(part->GetPdgCode());
1559  if ((absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
1560  {
1561  // Check if particle is in a radius around the jet
1562  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
1564  continue;
1565  else
1566  return kTRUE;
1567  }
1568  }
1569  }
1570  return kFALSE;
1571 
1572 }
1573 //________________________________________________________________________
1575 {
1576  fMultiplicity = 0;
1577  for(Int_t iCont=0; iCont<fParticleCollArray.GetEntriesFast(); iCont++)
1579 
1580  // ### Event control plots
1581  FillHistogram("hTrackCount", fMultiplicity, fCent);
1582  FillHistogram("hBackgroundPt", GetJetContainer(0)->GetRhoVal(), fCent);
1583 }
1584 
1585 //________________________________________________________________________
1587 {
1588  // ### Jet control plots
1589  AliJetContainer* jetContainer = GetJetContainer(0);
1590  FillHistogram("hJetPtRaw", jet->Pt(), fCent);
1591  FillHistogram("hJetPt", jet->Pt() - jetContainer->GetRhoVal()*jet->Area(), fCent);
1592  FillHistogram("hJetPhiEta", jet->Phi(), jet->Eta());
1593  FillHistogram("hJetArea", jet->Area(), fCent);
1594 
1595  // ### Jet constituent plots
1596  for(Int_t i = 0; i < jet->GetNumberOfParticleConstituents(); i++)
1597  {
1598  const AliVParticle* particle = jet->GetParticleConstituents()[i].GetParticle();
1599  if(!particle) continue;
1600 
1601  // Constituent eta/phi (relative to jet)
1602  Double_t deltaEta = jet->Eta() - particle->Eta();
1603  Double_t deltaPhi = TMath::Min(TMath::Abs(jet->Phi() - particle->Phi()), TMath::TwoPi() - TMath::Abs(jet->Phi() - particle->Phi()));
1604  if(!((jet->Phi() - particle->Phi() < 0) && (jet->Phi() - particle->Phi() <= TMath::Pi())))
1605  deltaPhi = -deltaPhi;
1606 
1607  FillHistogram("hConstituentPt", particle->Pt(), fCent);
1608  FillHistogram("hConstituentPhiEta", deltaPhi, deltaEta);
1609  }
1610 
1611 }
1612 
1613 //________________________________________________________________________
1615 {
1616  FillHistogram("hTrackPt", track->Pt(), fCent);
1617  FillHistogram("hTrackPhi", track->Phi(), fCent);
1618  FillHistogram("hTrackEta", track->Eta(), fCent);
1619  FillHistogram("hTrackEtaPt", track->Eta(), track->Pt());
1620  FillHistogram("hTrackPhiPt", track->Phi(), track->Pt());
1621  FillHistogram("hTrackPhiEta", track->Phi(), track->Eta());
1622 }
1623 
1624 //________________________________________________________________________
1625 void AliAnalysisTaskJetExtractor::AddPIDInformation(const AliVParticle* particle, Float_t& sigITS, Float_t& sigTPC, Float_t& sigTOF, Float_t& sigTRD, Short_t& recoPID, Int_t& truePID)
1626 {
1627  truePID = 9;
1628  if(!particle) return;
1629 
1630  // If we have AODs, retrieve particle PID signals
1631  const AliAODTrack* aodtrack = dynamic_cast<const AliAODTrack*>(particle);
1632 
1633  if(aodtrack)
1634  {
1635  // Get AOD value from reco
1636  recoPID = aodtrack->GetMostProbablePID();
1637  AliAODPid* pidObj = aodtrack->GetDetPid();
1638  if(!pidObj)
1639  return;
1640 
1641  sigITS = pidObj->GetITSsignal();
1642  sigTPC = pidObj->GetTPCsignal();
1643  sigTOF = pidObj->GetTOFsignal();
1644  sigTRD = pidObj->GetTRDsignal();
1645  }
1646 
1647  // Get truth values if we are on MC
1648  if(fMCParticleArray)
1649  {
1650  for(Int_t i=0; i<fMCParticleArray->GetEntriesFast();i++)
1651  {
1652  AliAODMCParticle* mcParticle = dynamic_cast<AliAODMCParticle*>(fMCParticleArray->At(i));
1653  if(!mcParticle) continue;
1654 
1655  if (mcParticle->GetLabel() == particle->GetLabel())
1656  {
1657  if(fSaveTrackPDGCode)
1658  {
1659  truePID = mcParticle->PdgCode();
1660  }
1661  else // Use same convention as for PID in AODs
1662  {
1663  if(TMath::Abs(mcParticle->PdgCode()) == 2212) // proton
1664  truePID = 4;
1665  else if (TMath::Abs(mcParticle->PdgCode()) == 211) // pion
1666  truePID = 2;
1667  else if (TMath::Abs(mcParticle->PdgCode()) == 321) // kaon
1668  truePID = 3;
1669  else if (TMath::Abs(mcParticle->PdgCode()) == 11) // electron
1670  truePID = 0;
1671  else if (TMath::Abs(mcParticle->PdgCode()) == 13) // muon
1672  truePID = 1;
1673  else if (TMath::Abs(mcParticle->PdgCode()) == 700201) // deuteron
1674  truePID = 5;
1675  else if (TMath::Abs(mcParticle->PdgCode()) == 700301) // triton
1676  truePID = 6;
1677  else if (TMath::Abs(mcParticle->PdgCode()) == 700302) // He3
1678  truePID = 7;
1679  else if (TMath::Abs(mcParticle->PdgCode()) == 700202) // alpha
1680  truePID = 8;
1681  else
1682  truePID = 9;
1683  }
1684 
1685  break;
1686  }
1687  }
1688  }
1689 }
1690 
1691 
1692 //________________________________________________________________________
1693 void AliAnalysisTaskJetExtractor::GetJetSplittings(AliEmcalJet* jet, std::vector<Float_t>& splittings_radiatorE, std::vector<Float_t>& splittings_kT, std::vector<Float_t>& splittings_theta, std::vector<Int_t>& splittings_secVtx_rank, std::vector<Int_t>& splittings_secVtx_index)
1694 {
1695  // ### Adapted from code in AliAnalysisTaskDmesonJetsSub ###
1696  // Define jet reclusterizer
1697  fastjet::JetAlgorithm jetAlgo(fastjet::cambridge_algorithm);
1698  Double_t jetRadius_CA = 1.0;
1699  fastjet::JetDefinition jetDefinition(jetAlgo, jetRadius_CA,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
1700 
1701  try{
1702  // Convert jet constituents to vector of fastjet::PseudoJet
1703  std::vector<fastjet::PseudoJet> particles;
1704  for(Int_t iConst=0; iConst<jet->GetNumberOfParticleConstituents(); iConst++)
1705  {
1706  const AliVParticle* constituent = jet->GetParticleConstituents()[iConst].GetParticle();
1707  Double_t p[3];
1708  constituent->PxPyPz(p);
1709  fastjet::PseudoJet pseudoJet = fastjet::PseudoJet(p[0], p[1], p[2], constituent->E());
1710 
1711  // If constituent is part of the N most significant sec. vertices, mark it
1712  for(UInt_t iVtx=0; iVtx<fSimpleSecVertices.size(); iVtx++)
1713  {
1714  if((constituent == fSimpleSecVertices[iVtx].fDaughter1) || (constituent == fSimpleSecVertices[iVtx].fDaughter2) || (constituent == fSimpleSecVertices[iVtx].fDaughter3))
1715  {
1716  pseudoJet.set_user_index(iVtx); // user_index is now index in temp sec. vtx vector
1717  break;
1718  }
1719  }
1720  particles.push_back(pseudoJet);
1721  }
1722 
1723  // Perform jet reclustering
1724  fastjet::ClusterSequence clusterSeq_CA(particles, jetDefinition);
1725  std::vector<fastjet::PseudoJet> jets_CA = clusterSeq_CA.inclusive_jets(0);
1726  jets_CA = sorted_by_pt(jets_CA);
1727  fastjet::PseudoJet radiator = jets_CA[0];
1728  fastjet::PseudoJet leadingSubJet;
1729  fastjet::PseudoJet subleadingSubJet;
1730 
1731  // Iterate through the splitting history of the CA clusterization
1732  while(radiator.has_parents(leadingSubJet,subleadingSubJet))
1733  {
1734  if(leadingSubJet.perp() < subleadingSubJet.perp())
1735  std::swap(leadingSubJet,subleadingSubJet);
1736 
1737  // Angle theta
1738  Float_t theta = leadingSubJet.delta_R(subleadingSubJet);
1739  // Radiator energy
1740  Float_t radiatorEnergy = leadingSubJet.e()+subleadingSubJet.e();
1741  // kT
1742  Float_t kT = subleadingSubJet.perp()*theta;
1743 
1744  // Go through leading subjet constituents and check if there are tracks belonging to one of the ten most significant sec. vertices
1745  Int_t secVtx_rank = -1; // rank = nth most displaced
1746  Int_t secVtx_index = -1; // index = index in sec. vertex array
1747  if(fSimpleSecVertices.size())
1748  {
1749  std::vector<fastjet::PseudoJet> leadingConsts = leadingSubJet.constituents();
1750  for(UInt_t iLeadingConst=0; iLeadingConst<leadingConsts.size(); iLeadingConst++)
1751  {
1752  if(leadingConsts[iLeadingConst].user_index()>=0)
1753  {
1754  secVtx_rank = leadingConsts[iLeadingConst].user_index();
1755  secVtx_index = fSimpleSecVertices[leadingConsts[iLeadingConst].user_index()].fIndex;
1756  }
1757  }
1758  }
1759 
1760  // Now add splitting properties to result vectors
1761  splittings_radiatorE.push_back(radiatorEnergy);
1762  splittings_theta.push_back(theta);
1763  splittings_kT.push_back(kT);
1764  splittings_secVtx_rank.push_back(secVtx_rank);
1765  splittings_secVtx_index.push_back(secVtx_index);
1766 
1767  // Continue with leadingSubJet as new radiator
1768  radiator=leadingSubJet;
1769  }
1770  } catch (fastjet::Error) { /*return -1;*/ }
1771 }
1772 
1773 
1774 //________________________________________________________________________
1776 {
1777  std::cout << "#########################################\n";
1778  std::cout << "Settings for extractor task " << GetName() << std::endl;
1779  std::cout << std::endl;
1780 
1781  std::cout << "### Event cuts ###" << std::endl;
1782  std::cout << std::endl;
1784  std::cout << Form("* Trigger track requirement: pT=%3.1f-%3.1f GeV/c, labels=%d-%d", fEventCut_TriggerTrackMinPt, fEventCut_TriggerTrackMaxPt, fEventCut_TriggerTrackMinLabel, fEventCut_TriggerTrackMaxLabel) << std::endl;
1786  std::cout << Form("* Trigger track event cut: %s", (fEventCut_TriggerTrackOrigin==1)? "From embedded event" : "Not from embedded event") << std::endl;
1787  if(fEventPercentage < 1)
1788  std::cout << Form("* Artificially lowered event efficiency: %f", fEventPercentage) << std::endl;
1790  std::cout << Form("* Random seeds explicitly set: %lu (for event/jet eff.), %lu (RCs)", fRandomSeed, fRandomSeedCones) << std::endl;
1791  std::cout << std::endl;
1792 
1793  std::cout << "### Jet properties ###" << std::endl;
1794  std::cout << "* Jet array name = " << GetJetContainer(0)->GetName() << std::endl;
1795  std::cout << "* Rho name = " << GetJetContainer(0)->GetRhoName() << std::endl;
1796  std::cout << "* Rho mass name = " << GetJetContainer(0)->GetRhoMassName() << std::endl;
1797  std::cout << std::endl;
1798 
1799  std::cout << "### Tree extraction properties ###" << std::endl;
1800  std::vector<Float_t> extractionPtBins = fJetTree->GetExtractionPercentagePtBins();
1801  std::vector<Float_t> extractionPercentages = fJetTree->GetExtractionPercentages();
1802  std::vector<Int_t> extractionHM = fJetTree->GetExtractionJetTypes_HM();
1803  std::vector<Int_t> extractionPM = fJetTree->GetExtractionJetTypes_PM();
1804  if(extractionPercentages.size())
1805  {
1806  std::cout << "* Extraction bins: (";
1807  for(Int_t i=0; i<static_cast<Int_t>(extractionPercentages.size()-1); i++)
1808  std::cout << extractionPtBins[i*2] << "-" << extractionPtBins[i*2+1] << " = " << extractionPercentages[i] << ", ";
1809  std::cout << extractionPtBins[(extractionPercentages.size()-1)*2] << "-" << extractionPtBins[(extractionPercentages.size()-1)*2+1] << " = " << extractionPercentages[(extractionPercentages.size()-1)];
1810 
1811  std::cout << ")" << std::endl;
1812  }
1813  if(extractionHM.size())
1814  {
1815  std::cout << "* Extraction of hadron-matched jets with types: (";
1816  for(Int_t i=0; i<static_cast<Int_t>(extractionHM.size()-1); i++)
1817  std::cout << extractionHM[i] << ", ";
1818  std::cout << extractionHM[extractionHM.size()-1];
1819  std::cout << ")" << std::endl;
1820  }
1821  if(extractionPM.size())
1822  {
1823  std::cout << "* Extraction of parton-matched jets with types: (";
1824  for(Int_t i=0; i<static_cast<Int_t>(extractionPM.size()-1); i++)
1825  std::cout << extractionPM[i] << ", ";
1826  std::cout << extractionPM[extractionPM.size()-1];
1827  std::cout << ")" << std::endl;
1828  }
1829  std::cout << std::endl;
1830 
1831  std::cout << "### Extracted data groups ###" << std::endl;
1832  std::cout << "* Basic event properties (ID, vertex, centrality, bgrd. densities, ...)" << std::endl;
1833  if(fSaveConstituents)
1834  std::cout << "* Jet constituents, basic properties (pt, eta, phi, charge, ...)" << std::endl;
1836  std::cout << "* Jet constituents, IPs" << std::endl;
1838  std::cout << "* Jet constituents, PID" << std::endl;
1839  if(fSaveCaloClusters)
1840  std::cout << "* Jet calorimeter clusters" << std::endl;
1841  if(fSaveMCInformation)
1842  std::cout << "* MC information (origin, matched jets, ...)" << std::endl;
1844  std::cout << "* Secondary vertices" << std::endl;
1845  if(fSaveJetShapes)
1846  std::cout << "* Jet shapes (jet mass, LeSub, pTD, ...)" << std::endl;
1847  if(fSaveJetSplittings)
1848  std::cout << "* Jet splittings (kT, theta, E from iterative CA reclustering)" << std::endl;
1849  if(fSaveTriggerTracks)
1850  std::cout << "* Trigger tracks" << std::endl;
1851  std::cout << std::endl;
1852  std::cout << "### Further settings ###" << std::endl;
1853  if(fIsEmbeddedEvent)
1854  std::cout << Form("* EMCal embedding framework will be used (at least on container has IsEmbedded() true)") << std::endl;
1856  std::cout << Form("* Detector level jet matching active with containger: %s *", GetJetContainer(1)->GetName()) << std::endl;
1858  std::cout << Form("* Particle level jet matching active with container: %s *", GetJetContainer(2)->GetName()) << std::endl;
1859  if(fMCParticleArray)
1860  std::cout << Form("* Particle level information available (for jet origin calculation, particle code): %s", fMCParticleArrayName.Data()) << std::endl;
1861  if(extractionHM.size())
1862  std::cout << Form("* Hadronic b/c-jet matching radius=%3.3f", fHadronMatchingRadius) << std::endl;
1863  if(fSaveMCInformation)
1864  std::cout << Form("* True particle label range for pt fraction=%d-%d", fTruthMinLabel, fTruthMaxLabel) << std::endl;
1866  std::cout << Form("* Particle PID codes will be PDG codes") << std::endl;
1868  std::cout << Form("* Particle PID codes will follow AliAODPid convention") << std::endl;
1870  std::cout << Form("* AliEmcalJet flavour tag is set from HF-jet tagging") << std::endl;
1871 
1872  std::cout << "#########################################\n\n";
1873 }
1874 
1875 
1876 //########################################################################
1877 // HELPERS
1878 //########################################################################
1879 
1880 //________________________________________________________________________
1881 inline Bool_t AliAnalysisTaskJetExtractor::IsTrackInCone(const AliVParticle* track, Double_t eta, Double_t phi, Double_t radius)
1882 {
1883  Double_t deltaR = GetDistance(track->Eta(), eta, track->Phi(), phi);
1884  if(deltaR <= radius)
1885  return kTRUE;
1886 
1887  return kFALSE;
1888 }
1889 
1890 //________________________________________________________________________
1891 inline Bool_t AliAnalysisTaskJetExtractor::IsClusterInCone( TLorentzVector clusterMomentum, Double_t eta, Double_t phi,Double_t radius)
1892 {
1893  Double_t deltaR = GetDistance(clusterMomentum.Eta(), eta, clusterMomentum.Phi(), phi);
1894  if(deltaR <= radius)
1895  return kTRUE;
1896 
1897  return kFALSE;
1898 }
1899 
1900 //________________________________________________________________________
1902 {
1903  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1904  if(!tmpHist)
1905  {
1906  AliError(Form("Cannot find histogram <%s> ",key)) ;
1907  return;
1908  }
1909 
1910  tmpHist->Fill(x);
1911 }
1912 
1913 //________________________________________________________________________
1915 {
1916  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1917  if(!tmpHist)
1918  {
1919  AliError(Form("Cannot find histogram <%s> ",key));
1920  return;
1921  }
1922 
1923  if (tmpHist->IsA()->GetBaseClass("TH1"))
1924  static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
1925  else if (tmpHist->IsA()->GetBaseClass("TH2"))
1926  static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
1927 }
1928 
1929 //________________________________________________________________________
1931 {
1932  TH2* tmpHist = static_cast<TH2*>(fOutput->FindObject(key));
1933  if(!tmpHist)
1934  {
1935  AliError(Form("Cannot find histogram <%s> ",key));
1936  return;
1937  }
1938 
1939  tmpHist->Fill(x,y,add);
1940 }
1941 
1942 //________________________________________________________________________
1944 {
1945  TH3* tmpHist = static_cast<TH3*>(fOutput->FindObject(key));
1946  if(!tmpHist)
1947  {
1948  AliError(Form("Cannot find histogram <%s> ",key));
1949  return;
1950  }
1951 
1952  if(add)
1953  tmpHist->Fill(x,y,z,add);
1954  else
1955  tmpHist->Fill(x,y,z);
1956 }
1957 
1958 
1959 //________________________________________________________________________
1960 template <class T> T* AliAnalysisTaskJetExtractor::AddHistogram1D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, const char* xTitle, const char* yTitle)
1961 {
1962  T* tmpHist = new T(name, title, xBins, xMin, xMax);
1963 
1964  tmpHist->GetXaxis()->SetTitle(xTitle);
1965  tmpHist->GetYaxis()->SetTitle(yTitle);
1966  tmpHist->SetOption(options);
1967  tmpHist->SetMarkerStyle(kFullCircle);
1968  tmpHist->Sumw2();
1969 
1970  fOutput->Add(tmpHist);
1971 
1972  return tmpHist;
1973 }
1974 
1975 //________________________________________________________________________
1976 template <class T> T* AliAnalysisTaskJetExtractor::AddHistogram2D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, Int_t yBins, Double_t yMin, Double_t yMax, const char* xTitle, const char* yTitle, const char* zTitle)
1977 {
1978  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax);
1979  tmpHist->GetXaxis()->SetTitle(xTitle);
1980  tmpHist->GetYaxis()->SetTitle(yTitle);
1981  tmpHist->GetZaxis()->SetTitle(zTitle);
1982  tmpHist->SetOption(options);
1983  tmpHist->SetMarkerStyle(kFullCircle);
1984  tmpHist->Sumw2();
1985 
1986  fOutput->Add(tmpHist);
1987 
1988  return tmpHist;
1989 }
1990 
1991 //________________________________________________________________________
1992 template <class T> T* AliAnalysisTaskJetExtractor::AddHistogram3D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, Int_t yBins, Double_t yMin, Double_t yMax, Int_t zBins, Double_t zMin, Double_t zMax, const char* xTitle, const char* yTitle, const char* zTitle)
1993 {
1994  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax, zBins, zMin, zMax);
1995  tmpHist->GetXaxis()->SetTitle(xTitle);
1996  tmpHist->GetYaxis()->SetTitle(yTitle);
1997  tmpHist->GetZaxis()->SetTitle(zTitle);
1998  tmpHist->SetOption(options);
1999  tmpHist->SetMarkerStyle(kFullCircle);
2000  tmpHist->Sumw2();
2001 
2002  fOutput->Add(tmpHist);
2003 
2004  return tmpHist;
2005 }
2006 
2007 //________________________________________________________________________
2009 {
2010  // Called once at the end of the analysis.
2011 }
2012 
2013 // ### ADDTASK MACRO
2014 //________________________________________________________________________
2015 AliAnalysisTaskJetExtractor* AliAnalysisTaskJetExtractor::AddTaskJetExtractor(TString trackArray, TString clusterArray, TString jetArray, TString rhoObject, Double_t jetRadius, AliRDHFJetsCutsVertex* vertexerCuts, const char* taskNameSuffix)
2016 {
2017  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
2018  Double_t minJetEta = 0.5;
2019  Double_t minJetPt = 0.15;
2020  Double_t minTrackPt = 0.15;
2021  Double_t minJetAreaPerc = 0.557;
2022  TString suffix = "";
2023  if(taskNameSuffix != 0)
2024  suffix = taskNameSuffix;
2025 
2026  // ###### Task name
2027  TString name("AliAnalysisTaskJetExtractor");
2028  if (jetArray != "") {
2029  name += "_";
2030  name += jetArray;
2031  }
2032  if (rhoObject != "") {
2033  name += "_";
2034  name += rhoObject;
2035  }
2036  if (suffix != "") {
2037  name += "_";
2038  name += suffix;
2039  }
2040 
2041  // ###### Setup task with default settings
2043  myTask->SetNeedEmcalGeom(kFALSE);
2044  myTask->SetVzRange(-10.,10.);
2045 
2046  // Particle container and track pt cut
2047  AliParticleContainer* trackCont = 0;
2048  if(trackArray == "mcparticles")
2049  trackCont = myTask->AddMCParticleContainer(trackArray);
2050  else if(trackArray =="mctracks")
2051  trackCont = myTask->AddParticleContainer(trackArray);
2052  else
2053  trackCont = myTask->AddTrackContainer(trackArray);
2054  trackCont->SetParticlePtCut(minTrackPt);
2055 
2056  // Particle container and track pt cut
2057  AliClusterContainer* clusterCont = 0;
2058  if(clusterArray != ""){
2059  clusterCont = myTask->AddClusterContainer(clusterArray);
2060  clusterCont->SetClusECut(0.);
2061  clusterCont->SetClusPtCut(0.);
2062  clusterCont->SetClusNonLinCorrEnergyCut(0.);
2063  clusterCont->SetClusHadCorrEnergyCut(0.30);
2064  clusterCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
2065  }
2066 
2067 
2068  // Secondary vertex cuts (default settings from PWGHF)
2069  // (can be overwritten by using myTask->SetVertexerCuts(cuts) from outside macro)
2070  AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts", "default");
2071  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
2072  esdTrackCuts->SetMinNClustersTPC(90);
2073  esdTrackCuts->SetMaxChi2PerClusterTPC(4);
2074  esdTrackCuts->SetRequireTPCRefit(kTRUE);
2075  esdTrackCuts->SetRequireITSRefit(kTRUE);
2076  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
2077  esdTrackCuts->SetMinDCAToVertexXY(0.);
2078  esdTrackCuts->SetEtaRange(-0.8, 0.8);
2079  esdTrackCuts->SetPtRange(1.0, 1.e10);
2080  vertexerCuts->AddTrackCuts(esdTrackCuts);
2081  vertexerCuts->SetNprongs(3);
2082  vertexerCuts->SetMinPtHardestTrack(1.0);//default 0.3
2083  vertexerCuts->SetSecVtxWithKF(kFALSE);//default with StrLinMinDist
2084  vertexerCuts->SetImpParCut(0.);//default 0
2085  vertexerCuts->SetDistPrimSec(0.);//default 0
2086  vertexerCuts->SetCospCut(-1);//default -1
2087  myTask->SetVertexerCuts(vertexerCuts);
2088 
2089  // Jet container
2090  AliJetContainer *jetCont = myTask->AddJetContainer(jetArray,6,jetRadius);
2091  if (jetCont) {
2092  jetCont->SetRhoName(rhoObject);
2093  jetCont->SetPercAreaCut(minJetAreaPerc);
2094  jetCont->SetJetPtCut(minJetPt);
2095  jetCont->SetLeadingHadronType(0);
2096  jetCont->SetPtBiasJetTrack(minTrackPt);
2097  jetCont->SetJetEtaLimits(-minJetEta, +minJetEta);
2098  jetCont->ConnectParticleContainer(trackCont);
2099  if(clusterCont)
2100  jetCont->ConnectClusterContainer(clusterCont);
2101  jetCont->SetMaxTrackPt(1000);
2102  }
2103 
2104  mgr->AddTask(myTask);
2105 
2106  // ###### Connect inputs/outputs
2107  mgr->ConnectInput (myTask, 0, mgr->GetCommonInputContainer() );
2108  mgr->ConnectOutput (myTask, 1, mgr->CreateContainer(Form("%s_histos", name.Data()), AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, Form("%s:ChargedJetsHadronCF", mgr->GetCommonFileName())) );
2109  mgr->ConnectOutput (myTask, 2, mgr->CreateContainer(Form("%s_tree", name.Data()), TTree::Class(), AliAnalysisManager::kOutputContainer, mgr->GetCommonFileName()) );
2110 
2111  return myTask;
2112 }
Float_t fBuffer_Shape_MomentumDispersion
! array buffer
void FillBuffer_Splittings(std::vector< Float_t > &splittings_radiatorE, std::vector< Float_t > &splittings_kT, std::vector< Float_t > &splittings_theta, Bool_t saveSecondaryVertices, std::vector< Int_t > &splittings_secVtx_rank, std::vector< Int_t > &splittings_secVtx_index)
void ExecOnce()
Perform steps needed to initialize the analysis.
Int_t pdg
Double_t GetFirstOrderSubtractedAngularity() const
void AddFlavourTag(Int_t tag)
Definition: AliEmcalJet.h:365
static const AliVEvent * GetEvent(const AliVEvent *inputEvent, bool isEmbedding=false)
Float_t * fBuffer_Track_ProdVtx_Z
! array buffer
Double_t Area() const
Definition: AliEmcalJet.h:130
void SetParticlePtCut(Double_t cut)
virtual AliVParticle * GetNextAcceptParticle()
Double_t GetRhoVal() const
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
Bool_t fIsEmbeddedEvent
Set to true if at least one embedding container is added to this task.
Double_t GetSecondOrderSubtractedSigma2() const
void GetTrackImpactParameters(const AliVVertex *vtx, const AliAODTrack *track, Float_t &d0, Float_t &d0cov, Float_t &z0, Float_t &z0cov)
Float_t fBuffer_Shape_LeSub_DerivCorr
! array buffer
const char * title
Definition: MakeQAPdf.C:27
Bool_t fDoPartLevelMatching
Whether or not we do particle level matching.
void GetJetType(AliEmcalJet *jet, Int_t &typeHM, Int_t &typePM, Int_t &typeIC)
Float_t GetPartonEta6() const
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:341
Bool_t fSetEmcalJetFlavour
if set, the flavour property of the AliEmcalJets will be set
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: External.C:244
Float_t GetPartonEta7() const
TTree * fJetTree
! tree structure
void FillBuffer_MonteCarlo(Int_t motherParton, Int_t motherHadron, Int_t partonInitialCollision, Float_t matchedJetDistance_Det, Float_t matchedJetPt_Det, Float_t matchedJetMass_Det, Float_t matchedJetAngularity_Det, Float_t matchedJetpTD_Det, Float_t matchedJetDistance_Part, Float_t matchedJetPt_Part, Float_t matchedJetMass_Part, Float_t matchedJetAngularity_Part, Float_t matchedJetpTD_Part, Float_t truePtFraction, Float_t truePtFraction_PartLevel, Float_t ptHard, Float_t eventWeight, Float_t impactParameter)
Double_t GetSecondOrderSubtractedConstituent() const
Double_t Eta() const
Definition: AliEmcalJet.h:121
long long Long64_t
Definition: External.C:43
void FillBuffer_TriggerTracks(std::vector< Float_t > &triggerTrackPt, std::vector< Float_t > &triggerTrackDeltaEta, std::vector< Float_t > &triggerTrackDeltaPhi)
T * AddHistogram3D(const char *name="CustomHistogram", const char *title="NO_TITLE", const char *options="", Int_t xBins=100, Double_t xMin=0.0, Double_t xMax=20.0, Int_t yBins=100, Double_t yMin=0.0, Double_t yMax=20.0, Int_t zBins=100, Double_t zMin=0.0, Double_t zMax=20.0, const char *xTitle="x axis", const char *yTitle="y axis", const char *zTitle="z axis")
Float_t fBuffer_Event_BackgroundDensity
! array buffer
void SetLeadingHadronType(Int_t t)
Float_t fBuffer_Shape_Mass_DerivCorr_1
! array buffer
Double_t Phi() const
Definition: AliEmcalJet.h:117
std::vector< Float_t > GetExtractionPercentages()
Bool_t IsClusterInCone(TLorentzVector clusterMomentum, Double_t eta, Double_t phi, Double_t radius)
Double_t mass
void SetPtBiasJetTrack(Float_t b)
Bool_t fNeedEmbedClusterContainer
If we need to get embedded cluster container (true for hybrid event)
Int_t fEventCut_TriggerTrackOrigin
Event requirement, trigger track origin (0: no cut, 1: from embedding, 2: not from embedding) ...
TSystem * gSystem
centrality
Bool_t fSaveMCInformation
save MC information
Float_t * fBuffer_Cluster_M02
! array buffer
Bool_t fSaveConstituentsIP
save arrays of constituent impact parameters
Int_t fTruthMinLabel
min track label to consider it as true particle
Double_t fSecondaryVertexMaxChi2
Max chi2 of secondary vertex (others will be discarded)
Int_t fMultiplicity
Multiplicity (number tracks, also for multiple containers)
void FillBuffer_JetShapes(AliEmcalJet *jet, Double_t leSub_noCorr, Double_t angularity, Double_t momentumDispersion, Double_t trackPtMean, Double_t trackPtMedian)
Float_t fBuffer_Jet_MC_MatchedPartLevelJet_Distance
! array buffer
Float_t GetPartonPhi7() const
Double_t fSecondaryVertexMaxDispersion
Max dispersion of secondary vertex (others will be discarded)
void AddPIDInformation(const AliVParticle *particle, Float_t &sigITS, Float_t &sigTPC, Float_t &sigTOF, Float_t &sigTRD, Short_t &recoPID, Int_t &truePID)
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
const std::vector< PWG::JETFW::AliEmcalParticleJetConstituent > & GetParticleConstituents() const
Get container with particle (track / MC particle) constituents.
Definition: AliEmcalJet.h:184
Double_t fEventCut_TriggerTrackMinPt
Event requirement, trigger track min pT.
std::vector< Float_t > GetExtractionPercentagePtBins()
std::vector< Float_t > fTriggerTracks_Pt
found trigger track pT
Int_t * fBuffer_Track_Label
! array buffer
Float_t fBuffer_Shape_Circularity_DerivCorr_1
! array buffer
Float_t fBuffer_Shape_pTD_DerivCorr_2
! array buffer
void LoadRhoMass(const AliVEvent *event)
void SetPercAreaCut(Float_t p)
Bool_t fSaveConstituentPID
save arrays of constituent PID parameters
Double_t fEventWeight
event weight for each event (implemented for JEWEL)
AliStack * stack
void SetVzRange(Double_t min, Double_t max)
Set pre-configured event cut object.
const Int_t kMaxNumConstituents
Float_t * fBuffer_Track_Eta
! array buffer
EJetType_t GetJetType() const
static AliAnalysisTaskJetExtractor * AddTaskJetExtractor(TString trackArray, TString clusterArray, TString jetArray, TString rhoObject, Double_t jetRadius, AliRDHFJetsCutsVertex *vertexerCuts, const char *taskNameSuffix)
T * AddHistogram2D(const char *name="CustomHistogram", const char *title="NO_TITLE", const char *options="", Int_t xBins=100, Double_t xMin=0.0, Double_t xMax=20.0, Int_t yBins=100, Double_t yMin=0.0, Double_t yMax=20.0, const char *xTitle="x axis", const char *yTitle="y axis", const char *zTitle="z axis")
AliClusterContainer * AddClusterContainer(const char *n)
Create new cluster container and attach it to the task.
void CalculateJetShapes(AliEmcalJet *jet, Double_t &leSub_noCorr, Double_t &angularity, Double_t &momentumDispersion, Double_t &trackPtMean, Double_t &trackPtMedian)
Int_t GetPartonFlag7() const
std::vector< Int_t > fExtractionJetTypes_HM
Jet-types that will be extracted with this tree (hadron matching)
Container for particles within the EMCAL framework.
Float_t GetPartonPhi6() const
std::vector< Float_t > fTriggerTracks_Phi
found trigger track phi
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
Int_t fEventCut_TriggerTrackMaxLabel
Event requirement, trigger track max label (can be used to selected toy particles) ...
Float_t fBuffer_Jet_MC_MatchedPartLevelJet_Pt
! array buffer
TObjArray fParticleCollArray
particle/track collection array
Float_t fBuffer_Shape_Mass_DerivCorr_2
! array buffer
Float_t fBuffer_Event_ImpactParameter
! array buffer
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
const TString & GetRhoMassName() const
TRandom3 * fRandomGenerator
! Random number generator, used for event + jet efficiency
Int_t GetPartonFlag6() const
void SetRhoName(const char *n)
void InitializeTree(Bool_t saveCaloClusters, Bool_t saveMCInformation, Bool_t saveMatchedJets_Det, Bool_t saveMatchedJets_Part, Bool_t saveConstituents, Bool_t saveConstituentsIP, Bool_t saveConstituentPID, Bool_t saveJetShapes, Bool_t saveSplittings, Bool_t saveSecondaryVertices, Bool_t saveTriggerTracks)
TRandom3 * fRandomGenerator
! random generator
std::vector< Float_t > fExtractionPercentages
Percentages which will be extracted for a given pT bin.
TClonesArray * fMCParticleArray
! Array of MC particles in event (usually mcparticles)
Float_t fBuffer_Jet_MC_MatchedPartLevelJet_Angularity
! array buffer
Float_t * fBuffer_Cluster_Time
! array buffer
Float_t fBuffer_Jet_MC_MatchedDetLevelJet_Angularity
! array buffer
Double_t fEventPercentage
percentage (0, 1] which will be extracted
Double_t fJetMatchingSharedPtFraction
Shared pT fraction required in matching.
Float_t fBuffer_Shape_Sigma2_DerivCorr_2
! array buffer
void FillHistogram3D(const char *key, Double_t x, Double_t y, Double_t z, Double_t add=0)
int Int_t
Definition: External.C:63
Int_t fBuffer_Jet_MC_MotherParton
! array buffer
Float_t fBuffer_Event_PtHard
! array buffer
void SetJetPtCut(Float_t cut)
Float_t * fBuffer_Track_Charge
! array buffer
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
void SetClusNonLinCorrEnergyCut(Double_t cut)
Float_t * fBuffer_Cluster_Phi
! array buffer
Float_t fBuffer_Event_Weight
! array buffer
T * AddHistogram1D(const char *name="CustomHistogram", const char *title="NO_TITLE", const char *options="", Int_t xBins=100, Double_t xMin=0.0, Double_t xMax=20.0, const char *xTitle="x axis", const char *yTitle="y axis")
Bool_t GetIsEmbedding() const
Get embedding status.
TRandom3 * fRandomGeneratorCones
! Random number generator, used for random cones
Float_t fBuffer_Shape_Sigma2_DerivCorr_1
! array buffer
AliParticleContainer * AddParticleContainer(const char *n)
Create new particle container and attach it to the task.
Double_t fImpactParameter
IP for each event (implemented for JEWEL)
Float_t fBuffer_Jet_MC_MatchedDetLevelJet_Mass
! array buffer
Int_t fBuffer_Event_Multiplicity
! array buffer
Float_t fBuffer_JetPhi
! array buffer
Double_t GetSecondOrderSubtractedAngularity() const
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
Float_t fBuffer_Shape_TrackPtMedian
! array buffer
AliVCluster * GetAcceptCluster(Int_t i) const
Int_t GetNClusters() const
Bool_t fSaveConstituents
save arrays of constituent basic properties
Bool_t IsTrackInCone(const AliVParticle *track, Double_t eta, Double_t phi, Double_t radius)
Int_t fTruthMaxLabel
max track label to consider it as true particle
void FillBuffer_PID(std::vector< Float_t > &trackPID_ITS, std::vector< Float_t > &trackPID_TPC, std::vector< Float_t > &trackPID_TOF, std::vector< Float_t > &trackPID_TRD, std::vector< Short_t > &trackPID_Reco, std::vector< Int_t > &trackPID_Truth)
std::vector< SimpleSecondaryVertex > fSimpleSecVertices
Vector of secondary vertices.
Float_t fBuffer_JetArea
! array buffer
std::vector< Int_t > fExtractionJetTypes_PM
Jet-types that will be extracted with this tree (parton matching)
std::vector< Float_t > fExtractionPercentagePtBins
pT-bins associated with fExtractionPercentages
Double_t fCent
!event centrality
Float_t * fBuffer_Track_Phi
! array buffer
Long64_t fBuffer_Event_ID
! array buffer
int GetNumberOfClusterConstituents() const
Get the number of cluster constituents.
Definition: AliEmcalJet.h:202
Float_t fBuffer_Event_Centrality
! array buffer
Int_t fBuffer_Jet_MC_MotherIC
! array buffer
TString fCustomStartupScript
Path to custom shell script that will be executed.
Bool_t fSaveJetSplittings
save jet splittings from iterative CA reclustering
Float_t * fBuffer_Track_ProdVtx_X
! array buffer
Float_t fBuffer_Shape_Mass_NoCorr
! array buffer
TString fMCParticleArrayName
Array name of MC particles in event (mcparticles)
Bool_t fInitialized
init state of tree
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.
Analysis task that implements AliEmcalJetTree to extract jets to a tree.
AliMCParticleContainer * AddMCParticleContainer(const char *n)
Create new container for MC particles and attach it to the task.
std::vector< Int_t > GetExtractionJetTypes_PM()
std::vector< Int_t > GetExtractionJetTypes_HM()
Float_t fBuffer_Event_BackgroundDensityMass
! array buffer
Double_t DeltaR(const AliVParticle *part) const
AliHFJetsTaggingVertex * fVtxTagger
! class for sec. vertexing
short Short_t
Definition: External.C:23
void FillBuffer_ImpactParameters(std::vector< Float_t > &trackIP_d0, std::vector< Float_t > &trackIP_z0, std::vector< Float_t > &trackIP_d0cov, std::vector< Float_t > &trackIP_z0cov)
void ConnectParticleContainer(AliParticleContainer *c)
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliRDHFJetsCutsVertex * fVertexerCuts
Cuts used for the vertexer (given in add task macro)
Float_t fBuffer_Jet_MC_TruePtFraction_PartLevel
! array buffer
Float_t fBuffer_Jet_MC_MatchedDetLevelJet_Distance
! array buffer
Double_t GetRhoVal(Int_t i=0) const
Double_t GetFirstOrderSubtractedCircularity() const
Float_t * fBuffer_Cluster_Eta
! array buffer
Float_t fBuffer_Shape_NumTracks_DerivCorr
! array buffer
Int_t fEventCut_TriggerTrackMinLabel
Event requirement, trigger track min label (can be used to selected toy particles) ...
Int_t fBuffer_NumClusters
! array buffer
const char * GetName() const
Float_t fBuffer_Jet_MC_MatchedPartLevelJet_pTD
! array buffer
Float_t fPtHard
!event -hard
void FillHistogram(const char *key, Double_t x)
void FillTrackControlHistograms(AliVTrack *track)
Float_t GetJetRadius() const
const std::vector< PWG::JETFW::AliEmcalClusterJetConstituent > & GetClusterConstituents() const
Get container with cluster constituents.
Definition: AliEmcalJet.h:190
AliEmcalList * fOutput
!output list
void SetClosestJet(AliEmcalJet *j, Double_t d)
Definition: AliEmcalJet.h:337
Float_t fBuffer_Jet_MC_MatchedDetLevelJet_Pt
! array buffer
Float_t fBuffer_JetPt
! array buffer
Definition: External.C:220
void ReconstructSecondaryVertices(const AliVVertex *primVtx, const AliEmcalJet *jet, std::vector< Float_t > &secVtx_X, std::vector< Float_t > &secVtx_Y, std::vector< Float_t > &secVtx_Z, std::vector< Float_t > &secVtx_Mass, std::vector< Float_t > &secVtx_Lxy, std::vector< Float_t > &secVtx_SigmaLxy, std::vector< Float_t > &secVtx_Chi2, std::vector< Float_t > &secVtx_Dispersion)
Float_t fBuffer_Event_Vertex_Y
! array buffer
void FillBuffer_SecVertices(std::vector< Float_t > &secVtx_X, std::vector< Float_t > &secVtx_Y, std::vector< Float_t > &secVtx_Z, std::vector< Float_t > &secVtx_Mass, std::vector< Float_t > &secVtx_Lxy, std::vector< Float_t > &secVtx_SigmaLxy, std::vector< Float_t > &secVtx_Chi2, std::vector< Float_t > &secVtx_Dispersion)
AliTrackContainer * AddTrackContainer(const char *n)
Create new track container and attach it to the task.
Int_t fBuffer_NumTracks
! array buffer
Double_t GetSecondOrderSubtractedCircularity() const
Bool_t fSaveSecondaryVertices
save reconstructed sec. vertex properties
void GetTrueJetPtFraction(AliEmcalJet *jet, Double_t &truePtFraction, Double_t &truePtFraction_mcparticles)
void SetVertexerCuts(AliRDHFJetsCutsVertex *val)
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
void SetMakeGeneralHistograms(Bool_t g)
Enable general histograms.
Float_t * fBuffer_Cluster_E
! array buffer
int GetNumberOfParticleConstituents() const
Get the number of particle constituents assigned to the given jet.
Definition: AliEmcalJet.h:196
void SetNeedEmcalGeom(Bool_t n)
Double_t fHadronMatchingRadius
Matching radius to search for beauty/charm hadrons around jet.
Float_t * fBuffer_Track_ProdVtx_Y
! array buffer
Base task in the EMCAL jet framework.
Float_t fBuffer_Shape_Angularity_DerivCorr_2
! array buffer
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
const AliEmcalPythiaInfo * GetPythiaInfo() const
Bool_t PxPyPz(Double_t p[3]) const
Definition: AliEmcalJet.h:111
Bool_t fDoDetLevelMatching
Whether or not we do det level matching.
Float_t fBuffer_Jet_MC_MatchedDetLevelJet_pTD
! array buffer
Float_t * fBuffer_Cluster_Pt
! array buffer
Float_t fBuffer_JetEta
! array buffer
void GetMatchedJetObservables(AliEmcalJet *jet, Double_t &detJetPt, Double_t &partJetPt, Double_t &detJetDistance, Double_t &partJetDistance, Double_t &detJetMass, Double_t &partJetMass, Double_t &detJetAngularity, Double_t &partJetAngularity, Double_t &detJetpTD, Double_t &partJetpTD)
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
Int_t * fBuffer_Cluster_Label
! array buffer
void GetJetSplittings(AliEmcalJet *jet, std::vector< Float_t > &splittings_radiatorE, std::vector< Float_t > &splittings_kT, std::vector< Float_t > &splittings_theta, std::vector< Int_t > &splittings_secVtx_rank, std::vector< Int_t > &splittings_secVtx_index)
Bool_t fSaveTrackPDGCode
save PDG code instead of code defined for AOD pid
Declaration of class AliEmcalPythiaInfo.
Double_t fJetMatchingRadius
Matching radius for geometrically matching jets.
const char Option_t
Definition: External.C:48
void SetClusPtCut(Double_t cut)
Bool_t fSaveTriggerTracks
save event trigger track
Float_t fBuffer_Shape_TrackPtMean
! array buffer
void UserCreateOutputObjects()
Main initialization function on the worker.
bool PerformGeometricalJetMatching(AliJetContainer &contBase, AliJetContainer &contTag, double maxDist)
Double_t GetFirstOrderSubtractedSigma2() const
Float_t fBuffer_Event_Vertex_X
! array buffer
Double_t GetDistance(Double_t eta1, Double_t eta2, Double_t phi1, Double_t phi2)
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:375
Double_t yMin
void SetClusECut(Double_t cut)
Int_t GetNAcceptedParticles() const
void SetDefaultClusterEnergy(Int_t d)
void ConnectClusterContainer(AliClusterContainer *c)
void SetRandomGenerator(TRandom3 *gen)
void SetMaxTrackPt(Float_t b)
Float_t fBuffer_Shape_Circularity_DerivCorr_2
! array buffer
void SetJetEtaLimits(Float_t min, Float_t max)
Int_t fBuffer_Jet_MC_MotherHadron
! array buffer
Double_t fEventCut_TriggerTrackMaxPt
Event requirement, trigger track max pT.
Bool_t fSaveCaloClusters
save calorimeter clusters
Container structure for EMCAL clusters.
Bool_t AddJetToTree(AliEmcalJet *jet, Bool_t saveConstituents, Bool_t saveConstituentsIP, Bool_t saveCaloClusters, Double_t *vertex, Float_t rho, Float_t rhoMass, Float_t centrality, Int_t multiplicity, Long64_t eventID, Float_t magField)
void swap(PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &first, PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &second)
std::vector< Float_t > fTriggerTracks_Eta
found trigger track eta
Double_t M() const
Definition: AliEmcalJet.h:120
void ResetCurrentID(Int_t i=-1)
Reset the iterator to a given index.
Float_t fBuffer_Shape_LeSub_NoCorr
! array buffer
Float_t fBuffer_Shape_Angularity_DerivCorr_1
! array buffer
Float_t fBuffer_Jet_MC_MatchedPartLevelJet_Mass
! array buffer
Container for jet within the EMCAL jet framework.
Float_t fBuffer_Jet_MC_TruePtFraction
! array buffer
Class managing creation of a tree containing jets.
Definition: External.C:196
Float_t fBuffer_Shape_Angularity_NoCorr
! array buffer
Float_t fBuffer_Shape_pTD_DerivCorr_1
! array buffer
Double_t yMax
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
Float_t fBuffer_Event_Vertex_Z
! array buffer
void SetClusHadCorrEnergyCut(Double_t cut)
const AliJetIterableContainer all() const
Float_t * fBuffer_Track_Pt
! array buffer
Float_t fBuffer_Event_MagneticField
! array buffer