AliPhysics  608b256 (608b256)
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 
33 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
34  #include <TPython.h>
35 #endif
36 
37 #include "AliMCEvent.h"
38 #include "AliESDVertex.h"
39 #include "AliAODVertex.h"
40 #include "AliAODPid.h"
41 
42 #include "AliRhoParameter.h"
43 
44 #include "AliVTrack.h"
45 #include "AliVHeader.h"
46 #include "AliEmcalJet.h"
47 #include "AliLog.h"
48 #include "AliJetContainer.h"
49 #include "AliParticleContainer.h"
50 #include "AliClusterContainer.h"
51 #include "AliAODTrack.h"
52 #include "AliVParticle.h"
53 #include "TRandom3.h"
54 #include "AliEmcalPythiaInfo.h"
56 #include "AliAnalysisManager.h"
57 #include "AliEmcalContainerUtils.h"
58 #include "AliFJWrapper.h"
59 
60 
61 #include "AliGenHepMCEventHeader.h"
62 #include "AliGenHijingEventHeader.h"
63 #include "AliHFJetsTaggingVertex.h"
64 #include "AliRDHFJetsCutsVertex.h"
65 
67 
69 ClassImp(AliEmcalJetTree)
71 
75 
76 //________________________________________________________________________
77 AliEmcalJetTree::AliEmcalJetTree() : TNamed("CustomTree", "CustomTree"), fJetTree(0), fInitialized(0), fExtractionPercentages(), fExtractionPercentagePtBins(), fExtractionJetTypes_HM(), fExtractionJetTypes_PM()
78 {
79  // For these arrays, we need to reserve memory
80  fBuffer_Track_Pt = new Float_t[kMaxNumConstituents];
81  fBuffer_Track_Eta = new Float_t[kMaxNumConstituents];
82  fBuffer_Track_Phi = new Float_t[kMaxNumConstituents];
83  fBuffer_Track_Charge = new Float_t[kMaxNumConstituents];
84  fBuffer_Track_Label = new Int_t [kMaxNumConstituents];
85  fBuffer_Track_ProdVtx_X = new Float_t[kMaxNumConstituents];
86  fBuffer_Track_ProdVtx_Y = new Float_t[kMaxNumConstituents];
87  fBuffer_Track_ProdVtx_Z = new Float_t[kMaxNumConstituents];
88 
89  fBuffer_Cluster_Pt = new Float_t[kMaxNumConstituents];
90  fBuffer_Cluster_E = new Float_t[kMaxNumConstituents];
91  fBuffer_Cluster_Eta = new Float_t[kMaxNumConstituents];
92  fBuffer_Cluster_Phi = new Float_t[kMaxNumConstituents];
93  fBuffer_Cluster_M02 = new Float_t[kMaxNumConstituents];
94  fBuffer_Cluster_Time = new Float_t[kMaxNumConstituents];
95  fBuffer_Cluster_Label = new Int_t[kMaxNumConstituents];
96 }
97 
98 //________________________________________________________________________
99 AliEmcalJetTree::AliEmcalJetTree(const char* name) : TNamed(name, name), fJetTree(0), fInitialized(0), fExtractionPercentages(), fExtractionPercentagePtBins(), fExtractionJetTypes_HM(), fExtractionJetTypes_PM()
100 {
101  // For these arrays, we need to reserve memory
110 
118 }
119 
120 //________________________________________________________________________
121 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)
122 {
123  if(!fInitialized)
124  AliFatal("Tree is not initialized.");
125 
126  fBuffer_JetPt = jet->Pt() - rho*jet->Area();
127 
128  // Check if jet type is contained in extraction list
129  if( (fExtractionJetTypes_PM.size() || fExtractionJetTypes_HM.size()) &&
132  return kFALSE;
133 
134  // Check acceptance percentage for the given jet and discard statistically on demand
135  Bool_t inPtRange = kFALSE;
136  for(size_t i = 0; i<fExtractionPercentages.size(); i++)
137  {
139  {
140  inPtRange = kTRUE;
141  if(fRandomGenerator->Rndm() >= fExtractionPercentages[i])
142  return kFALSE;
143  break;
144  }
145  }
146  if(!inPtRange && fExtractionPercentagePtBins.size()) // only discard jet if fExtractionPercentagePtBins was given
147  return kFALSE;
148 
149  // Set basic properties
150  fBuffer_JetEta = jet->Eta();
151  fBuffer_JetPhi = jet->Phi();
152  fBuffer_JetArea = jet->Area();
153 
154  // Set event properties
157  fBuffer_Event_Vertex_X = vertex ? vertex[0] : 0;
158  fBuffer_Event_Vertex_Y = vertex ? vertex[1] : 0;
159  fBuffer_Event_Vertex_Z = vertex ? vertex[2] : 0;
161  fBuffer_Event_Multiplicity = multiplicity;
162  fBuffer_Event_ID = eventID;
163  fBuffer_Event_MagneticField = magField;
164 
165  // Extract basic constituent track properties directly from AliEmcalJet object
166  fBuffer_NumTracks = 0;
167  if(saveConstituents || saveConstituentsIP)
168  for(Int_t i = 0; i < jet->GetNumberOfParticleConstituents(); i++)
169  {
170  const AliVParticle* particle = jet->GetParticleConstituents()[i].GetParticle();
171  if(!particle) continue;
172 
173  if(saveConstituents)
174  {
175  fBuffer_Track_Pt[fBuffer_NumTracks] = particle->Pt();
176  fBuffer_Track_Eta[fBuffer_NumTracks] = particle->Eta();
177  fBuffer_Track_Phi[fBuffer_NumTracks] = particle->Phi();
178  fBuffer_Track_Charge[fBuffer_NumTracks] = particle->Charge();
179  fBuffer_Track_Label[fBuffer_NumTracks] = particle->GetLabel();
180  }
181  if(saveConstituentsIP) // track production vertices
182  {
183  fBuffer_Track_ProdVtx_X[fBuffer_NumTracks] = particle->Xv();
184  fBuffer_Track_ProdVtx_Y[fBuffer_NumTracks] = particle->Yv();
185  fBuffer_Track_ProdVtx_Z[fBuffer_NumTracks] = particle->Zv();
186  }
188  }
189 
190  // Extract basic constituent cluster properties directly from AliEmcalJet object
192  if(saveCaloClusters)
193  for(Int_t i = 0; i < jet->GetNumberOfClusterConstituents(); i++)
194  {
195  const AliVCluster* cluster = jet->GetClusterConstituents()[i].GetCluster();
196  if(!cluster) continue;
197 
198  // #### Retrieve cluster pT
199  TLorentzVector clusterMomentum;
200  cluster->GetMomentum(clusterMomentum, vertex);
201  // ####
202 
203  fBuffer_Cluster_Pt[fBuffer_NumClusters] = clusterMomentum.Perp();
204  fBuffer_Cluster_E[fBuffer_NumClusters] = cluster->E();
205  fBuffer_Cluster_Eta[fBuffer_NumClusters] = clusterMomentum.Eta();
206  fBuffer_Cluster_Phi[fBuffer_NumClusters] = clusterMomentum.Phi();
207  fBuffer_Cluster_M02[fBuffer_NumClusters] = cluster->GetM02();
208  fBuffer_Cluster_Time[fBuffer_NumClusters] = cluster->GetTOF();
209  fBuffer_Cluster_Label[fBuffer_NumClusters] = cluster->GetLabel();
211  }
212 
213 
214  // Add all buffers to tree
215  fJetTree->Fill();
216 
217  return kTRUE;
218 }
219 
220 //________________________________________________________________________
221 void AliEmcalJetTree::FillBuffer_TriggerTracks(std::vector<Float_t>& triggerTrackPt, std::vector<Float_t>& triggerTrackDeltaEta, std::vector<Float_t>& triggerTrackDeltaPhi)
222 {
223  fBuffer_NumTriggerTracks = triggerTrackPt.size();
224  fJetTree->SetBranchAddress("Jet_TriggerTrack_Pt", triggerTrackPt.data());
225  fJetTree->SetBranchAddress("Jet_TriggerTrack_dEta", triggerTrackDeltaEta.data());
226  fJetTree->SetBranchAddress("Jet_TriggerTrack_dPhi", triggerTrackDeltaPhi.data());
227 }
228 
229 //________________________________________________________________________
230 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)
231 {
232  fJetTree->SetBranchAddress("Jet_Track_CovIPd", trackIP_d0cov.data());
233  fJetTree->SetBranchAddress("Jet_Track_CovIPz", trackIP_z0cov.data());
234  fJetTree->SetBranchAddress("Jet_Track_IPd", trackIP_d0.data());
235  fJetTree->SetBranchAddress("Jet_Track_IPz", trackIP_z0.data());
236 }
237 
238 //________________________________________________________________________
239 void AliEmcalJetTree::FillBuffer_MonteCarlo(Int_t motherParton, Int_t motherHadron, Int_t partonInitialCollision, Float_t matchedJetDistance, Float_t matchedJetPt, Float_t matchedJetMass, Float_t truePtFraction, Float_t truePtFraction_mcparticles, Float_t ptHard, Float_t eventWeight, Float_t impactParameter)
240 {
241  fBuffer_Jet_MC_MotherParton = motherParton;
242  fBuffer_Jet_MC_MotherHadron = motherHadron;
243  fBuffer_Jet_MC_MotherIC = partonInitialCollision;
244  fBuffer_Jet_MC_MatchedJet_Distance = matchedJetDistance;
245  fBuffer_Jet_MC_MatchedJet_Pt = matchedJetPt;
246  fBuffer_Jet_MC_MatchedJet_Mass = matchedJetMass;
247  fBuffer_Jet_MC_TruePtFraction = truePtFraction;
248  fBuffer_Jet_MC_TruePtFraction_mcparticles = truePtFraction_mcparticles;
249 
250  fBuffer_Event_PtHard = ptHard;
251  fBuffer_Event_Weight = eventWeight;
252  fBuffer_Event_ImpactParameter = impactParameter;
253 }
254 
255 
256 //________________________________________________________________________
257 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)
258 {
259  fJetTree->SetBranchAddress("Jet_Track_PID_ITS", trackPID_ITS.data());
260  fJetTree->SetBranchAddress("Jet_Track_PID_TPC", trackPID_TPC.data());
261  fJetTree->SetBranchAddress("Jet_Track_PID_TOF", trackPID_TOF.data());
262  fJetTree->SetBranchAddress("Jet_Track_PID_TRD", trackPID_TRD.data());
263  fJetTree->SetBranchAddress("Jet_Track_PID_Reconstructed", trackPID_Reco.data());
264  if(trackPID_Truth.data())
265  fJetTree->SetBranchAddress("Jet_Track_PID_Truth", trackPID_Truth.data());
266 }
267 
268 //________________________________________________________________________
269 void AliEmcalJetTree::FillBuffer_JetShapes(AliEmcalJet* jet, Double_t leSub_noCorr, Double_t angularity, Double_t momentumDispersion, Double_t trackPtMean, Double_t trackPtMedian)
270 {
271  fBuffer_Shape_Mass_NoCorr = jet->M();
276  fBuffer_Shape_Angularity_NoCorr = angularity;
285  fBuffer_Shape_LeSub_NoCorr = leSub_noCorr;
286  fBuffer_Shape_MomentumDispersion = momentumDispersion;
287  fBuffer_Shape_TrackPtMean = trackPtMean;
288  fBuffer_Shape_TrackPtMedian = trackPtMedian;
289 }
290 
291 //________________________________________________________________________
292 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)
293 {
294  fBuffer_NumSplittings = splittings_radiatorE.size();
295  fJetTree->SetBranchAddress("Jet_Splitting_RadiatorE", splittings_radiatorE.data());
296  fJetTree->SetBranchAddress("Jet_Splitting_kT", splittings_kT.data());
297  fJetTree->SetBranchAddress("Jet_Splitting_Theta", splittings_theta.data());
298  if(saveSecondaryVertices)
299  {
300  fJetTree->SetBranchAddress("Jet_Splitting_SecVtx_Rank", splittings_secVtx_rank.data());
301  fJetTree->SetBranchAddress("Jet_Splitting_SecVtx_Index", splittings_secVtx_index.data());
302  }
303 }
304 
305 //________________________________________________________________________
306 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)
307 {
308  fBuffer_NumSecVertices = secVtx_X.size();
309  fJetTree->SetBranchAddress("Jet_SecVtx_X", secVtx_X.data());
310  fJetTree->SetBranchAddress("Jet_SecVtx_Y", secVtx_Y.data());
311  fJetTree->SetBranchAddress("Jet_SecVtx_Z", secVtx_Z.data());
312  fJetTree->SetBranchAddress("Jet_SecVtx_Mass", secVtx_Mass.data());
313  fJetTree->SetBranchAddress("Jet_SecVtx_Lxy", secVtx_Lxy.data());
314  fJetTree->SetBranchAddress("Jet_SecVtx_SigmaLxy", secVtx_SigmaLxy.data());
315  fJetTree->SetBranchAddress("Jet_SecVtx_Chi2", secVtx_Chi2.data());
316  fJetTree->SetBranchAddress("Jet_SecVtx_Dispersion", secVtx_Dispersion.data());
317 }
318 
319 //________________________________________________________________________
320 void AliEmcalJetTree::InitializeTree(Bool_t saveCaloClusters, Bool_t saveMCInformation, Bool_t saveConstituents, Bool_t saveConstituentsIP, Bool_t saveConstituentPID, Bool_t saveJetShapes, Bool_t saveSplittings, Bool_t saveSecondaryVertices, Bool_t saveTriggerTracks)
321 {
322  // Create the tree with active branches
323 
324  void* dummy = 0; // for some branches, do not need a buffer pointer yet
325 
326 
327  if(fInitialized)
328  AliFatal("Tree is already initialized.");
329 
330  // ### Prepare the jet tree
331  fJetTree = new TTree(Form("JetTree_%s", GetName()), "");
332 
333  fJetTree->Branch("Jet_Pt",&fBuffer_JetPt,"Jet_Pt/F");
334  fJetTree->Branch("Jet_Phi",&fBuffer_JetPhi,"Jet_Phi/F");
335  fJetTree->Branch("Jet_Eta",&fBuffer_JetEta,"Jet_Eta/F");
336  fJetTree->Branch("Jet_Area",&fBuffer_JetArea,"Jet_Area/F");
337  fJetTree->Branch("Jet_NumTracks",&fBuffer_NumTracks,"Jet_NumTracks/I");
338  if(saveCaloClusters)
339  fJetTree->Branch("Jet_NumClusters",&fBuffer_NumClusters,"Jet_NumClusters/I");
340 
341  fJetTree->Branch("Event_BackgroundDensity",&fBuffer_Event_BackgroundDensity,"Event_BackgroundDensity/F");
342  fJetTree->Branch("Event_BackgroundDensityMass",&fBuffer_Event_BackgroundDensityMass,"Event_BackgroundDensityMass/F");
343  fJetTree->Branch("Event_Vertex_X",&fBuffer_Event_Vertex_X,"Event_Vertex_X/F");
344  fJetTree->Branch("Event_Vertex_Y",&fBuffer_Event_Vertex_Y,"Event_Vertex_Y/F");
345  fJetTree->Branch("Event_Vertex_Z",&fBuffer_Event_Vertex_Z,"Event_Vertex_Z/F");
346  fJetTree->Branch("Event_Centrality",&fBuffer_Event_Centrality,"Event_Centrality/F");
347  fJetTree->Branch("Event_Multiplicity",&fBuffer_Event_Multiplicity,"Event_Multiplicity/I");
348  fJetTree->Branch("Event_ID",&fBuffer_Event_ID,"Event_ID/L");
349  fJetTree->Branch("Event_MagneticField",&fBuffer_Event_MagneticField,"Event_MagneticField/F");
350 
351  if(saveMCInformation)
352  {
353  fJetTree->Branch("Event_PtHard",&fBuffer_Event_PtHard,"Event_PtHard/F");
354  fJetTree->Branch("Event_Weight",&fBuffer_Event_Weight,"Event_Weight/F");
355  fJetTree->Branch("Event_ImpactParameter",&fBuffer_Event_ImpactParameter,"Event_ImpactParameter/F");
356  }
357 
358  if(saveConstituents)
359  {
360  fJetTree->Branch("Jet_Track_Pt",fBuffer_Track_Pt,"Jet_Track_Pt[Jet_NumTracks]/F");
361  fJetTree->Branch("Jet_Track_Phi",fBuffer_Track_Phi,"Jet_Track_Phi[Jet_NumTracks]/F");
362  fJetTree->Branch("Jet_Track_Eta",fBuffer_Track_Eta,"Jet_Track_Eta[Jet_NumTracks]/F");
363  fJetTree->Branch("Jet_Track_Charge",fBuffer_Track_Charge,"Jet_Track_Charge[Jet_NumTracks]/F");
364  if(saveMCInformation)
365  fJetTree->Branch("Jet_Track_Label",fBuffer_Track_Label,"Jet_Track_Label[Jet_NumTracks]/I");
366  }
367 
368  if(saveCaloClusters)
369  {
370  fJetTree->Branch("Jet_Cluster_Pt",fBuffer_Cluster_Pt,"Jet_Cluster_Pt[Jet_NumClusters]/F");
371  fJetTree->Branch("Jet_Cluster_E",fBuffer_Cluster_E,"Jet_Cluster_Pt[Jet_NumClusters]/F");
372  fJetTree->Branch("Jet_Cluster_Phi",fBuffer_Cluster_Phi,"Jet_Cluster_Phi[Jet_NumClusters]/F");
373  fJetTree->Branch("Jet_Cluster_Eta",fBuffer_Cluster_Eta,"Jet_Cluster_Eta[Jet_NumClusters]/F");
374  fJetTree->Branch("Jet_Cluster_M02",fBuffer_Cluster_M02,"Jet_Cluster_M02[Jet_NumClusters]/F");
375  fJetTree->Branch("Jet_Cluster_Time",fBuffer_Cluster_Time,"Jet_Cluster_Time[Jet_NumClusters]/F");
376  if(saveMCInformation)
377  fJetTree->Branch("Jet_Cluster_Label",fBuffer_Cluster_Label,"Jet_Cluster_Label[Jet_NumClusters]/I");
378  }
379 
380  if(saveConstituentsIP)
381  {
382  fJetTree->Branch("Jet_Track_IPd",&dummy,"Jet_Track_IPd[Jet_NumTracks]/F");
383  fJetTree->Branch("Jet_Track_IPz",&dummy,"Jet_Track_IPz[Jet_NumTracks]/F");
384  fJetTree->Branch("Jet_Track_CovIPd",&dummy,"Jet_Track_CovIPd[Jet_NumTracks]/F");
385  fJetTree->Branch("Jet_Track_CovIPz",&dummy,"Jet_Track_CovIPz[Jet_NumTracks]/F");
386 
387  fJetTree->Branch("Jet_Track_ProdVtx_X",fBuffer_Track_ProdVtx_X,"Jet_Track_ProdVtx_X[Jet_NumTracks]/F");
388  fJetTree->Branch("Jet_Track_ProdVtx_Y",fBuffer_Track_ProdVtx_Y,"Jet_Track_ProdVtx_Y[Jet_NumTracks]/F");
389  fJetTree->Branch("Jet_Track_ProdVtx_Z",fBuffer_Track_ProdVtx_Z,"Jet_Track_ProdVtx_Z[Jet_NumTracks]/F");
390  }
391 
392  if(saveConstituentPID)
393  {
394  fJetTree->Branch("Jet_Track_PID_ITS",&dummy,"Jet_Track_PID_ITS[Jet_NumTracks]/F");
395  fJetTree->Branch("Jet_Track_PID_TPC",&dummy,"Jet_Track_PID_TPC[Jet_NumTracks]/F");
396  fJetTree->Branch("Jet_Track_PID_TOF",&dummy,"Jet_Track_PID_TOF[Jet_NumTracks]/F");
397  fJetTree->Branch("Jet_Track_PID_TRD",&dummy,"Jet_Track_PID_TRD[Jet_NumTracks]/F");
398 
399  fJetTree->Branch("Jet_Track_PID_Reconstructed",&dummy,"Jet_Track_PID_Reconstructed[Jet_NumTracks]/S");
400  if(saveMCInformation)
401  fJetTree->Branch("Jet_Track_PID_Truth",&dummy,"Jet_Track_PID_Truth[Jet_NumTracks]/I");
402  }
403 
404  if(saveJetShapes)
405  {
406  fJetTree->Branch("Jet_Shape_Mass_NoCorr",&fBuffer_Shape_Mass_NoCorr,"Jet_Shape_Mass_NoCorr/F");
407  fJetTree->Branch("Jet_Shape_Mass_DerivCorr_1",&fBuffer_Shape_Mass_DerivCorr_1,"Jet_Shape_Mass_DerivCorr_1/F");
408  fJetTree->Branch("Jet_Shape_Mass_DerivCorr_2",&fBuffer_Shape_Mass_DerivCorr_2,"Jet_Shape_Mass_DerivCorr_2/F");
409  fJetTree->Branch("Jet_Shape_pTD_DerivCorr_1",&fBuffer_Shape_pTD_DerivCorr_1,"Jet_Shape_pTD_DerivCorr_1/F");
410  fJetTree->Branch("Jet_Shape_pTD_DerivCorr_2",&fBuffer_Shape_pTD_DerivCorr_2,"Jet_Shape_pTD_DerivCorr_2/F");
411  fJetTree->Branch("Jet_Shape_LeSub_NoCorr",&fBuffer_Shape_LeSub_NoCorr,"Jet_Shape_LeSub_NoCorr/F");
412  fJetTree->Branch("Jet_Shape_LeSub_DerivCorr",&fBuffer_Shape_LeSub_DerivCorr,"Jet_Shape_LeSub_DerivCorr/F");
413  fJetTree->Branch("Jet_Shape_Angularity",&fBuffer_Shape_Angularity_NoCorr,"Jet_Shape_Angularity/F");
414  fJetTree->Branch("Jet_Shape_Angularity_DerivCorr_1",&fBuffer_Shape_Angularity_DerivCorr_1,"Jet_Shape_Angularity_DerivCorr_1/F");
415  fJetTree->Branch("Jet_Shape_Angularity_DerivCorr_2",&fBuffer_Shape_Angularity_DerivCorr_2,"Jet_Shape_Angularity_DerivCorr_2/F");
416  fJetTree->Branch("Jet_Shape_Circularity_DerivCorr_1",&fBuffer_Shape_Circularity_DerivCorr_1,"Jet_Shape_Circularity_DerivCorr_1/F");
417  fJetTree->Branch("Jet_Shape_Circularity_DerivCorr_2",&fBuffer_Shape_Circularity_DerivCorr_2,"Jet_Shape_Circularity_DerivCorr_2/F");
418  fJetTree->Branch("Jet_Shape_Sigma2_DerivCorr_1",&fBuffer_Shape_Sigma2_DerivCorr_1,"Jet_Shape_Sigma2_DerivCorr_1/F");
419  fJetTree->Branch("Jet_Shape_Sigma2_DerivCorr_2",&fBuffer_Shape_Sigma2_DerivCorr_2,"Jet_Shape_Sigma2_DerivCorr_2/F");
420  fJetTree->Branch("Jet_Shape_NumTracks_DerivCorr",&fBuffer_Shape_NumTracks_DerivCorr,"Jet_Shape_NumTracks_DerivCorr/F");
421  fJetTree->Branch("Jet_Shape_MomentumDispersion",&fBuffer_Shape_MomentumDispersion,"Jet_Shape_MomentumDispersion/F");
422  fJetTree->Branch("Jet_Shape_TrackPtMean",&fBuffer_Shape_TrackPtMean,"Jet_Shape_TrackPtMean/F");
423  fJetTree->Branch("Jet_Shape_TrackPtMedian",&fBuffer_Shape_TrackPtMedian,"Jet_Shape_TrackPtMedian/F");
424  }
425 
426  if(saveSplittings)
427  {
428  fJetTree->Branch("Jet_NumSplittings",&fBuffer_NumSplittings,"Jet_NumSplittings/I");
429  fJetTree->Branch("Jet_Splitting_Theta",&dummy,"Jet_Splitting_Theta[Jet_NumSplittings]/F");
430  fJetTree->Branch("Jet_Splitting_RadiatorE",&dummy,"Jet_Splitting_RadiatorE[Jet_NumSplittings]/F");
431  fJetTree->Branch("Jet_Splitting_kT",&dummy,"Jet_Splitting_kT[Jet_NumSplittings]/F");
432  if(saveSecondaryVertices)
433  {
434  fJetTree->Branch("Jet_Splitting_SecVtx_Rank",&dummy,"Jet_Splitting_SecVtx_Rank[Jet_NumSplittings]/I");
435  fJetTree->Branch("Jet_Splitting_SecVtx_Index",&dummy,"Jet_Splitting_SecVtx_Index[Jet_NumSplittings]/I");
436  }
437  }
438 
439  if(saveMCInformation)
440  {
441  fJetTree->Branch("Jet_MC_MotherParton",&fBuffer_Jet_MC_MotherParton,"Jet_MC_MotherParton/I");
442  fJetTree->Branch("Jet_MC_MotherHadron",&fBuffer_Jet_MC_MotherHadron,"Jet_MC_MotherHadron/I");
443  fJetTree->Branch("Jet_MC_MotherIC",&fBuffer_Jet_MC_MotherIC,"Jet_MC_MotherIC/I");
444  fJetTree->Branch("Jet_MC_MatchedJet_Distance",&fBuffer_Jet_MC_MatchedJet_Distance,"Jet_MC_MatchedJet_Distance/F");
445  fJetTree->Branch("Jet_MC_MatchedJet_Pt",&fBuffer_Jet_MC_MatchedJet_Pt,"Jet_MC_MatchedJet_Pt/F");
446  fJetTree->Branch("Jet_MC_MatchedJet_Mass",&fBuffer_Jet_MC_MatchedJet_Mass,"Jet_MC_MatchedJet_Mass/F");
447  fJetTree->Branch("Jet_MC_TruePtFraction",&fBuffer_Jet_MC_TruePtFraction,"Jet_MC_TruePtFraction/F");
448  fJetTree->Branch("Jet_MC_TruePtFraction_mcparticles",&fBuffer_Jet_MC_TruePtFraction_mcparticles,"Jet_MC_TruePtFraction_mcparticles/F");
449  }
450 
451  if(saveSecondaryVertices)
452  {
453  fJetTree->Branch("Jet_NumSecVertices",&fBuffer_NumSecVertices,"Jet_NumSecVertices/I");
454 
455  fJetTree->Branch("Jet_SecVtx_X",&dummy,"Jet_SecVtx_X[Jet_NumSecVertices]/F");
456  fJetTree->Branch("Jet_SecVtx_Y",&dummy,"Jet_SecVtx_Y[Jet_NumSecVertices]/F");
457  fJetTree->Branch("Jet_SecVtx_Z",&dummy,"Jet_SecVtx_Z[Jet_NumSecVertices]/F");
458  fJetTree->Branch("Jet_SecVtx_Mass",&dummy,"Jet_SecVtx_Mass[Jet_NumSecVertices]/F");
459  fJetTree->Branch("Jet_SecVtx_Lxy",&dummy,"Jet_SecVtx_Lxy[Jet_NumSecVertices]/F");
460  fJetTree->Branch("Jet_SecVtx_SigmaLxy",&dummy,"Jet_SecVtx_SigmaLxy[Jet_NumSecVertices]/F");
461  fJetTree->Branch("Jet_SecVtx_Chi2",&dummy,"Jet_SecVtx_Chi2[Jet_NumSecVertices]/F");
462  fJetTree->Branch("Jet_SecVtx_Dispersion",&dummy,"Jet_SecVtx_Dispersion[Jet_NumSecVertices]/F");
463  }
464 
465  // Trigger track requirement active -> Save trigger track
466  if(saveTriggerTracks)
467  {
468  fJetTree->Branch("Jet_NumTriggerTracks",&fBuffer_NumTriggerTracks,"Jet_NumTriggerTracks/I");
469  fJetTree->Branch("Jet_TriggerTrack_Pt",&dummy,"Jet_TriggerTrack_Pt[Jet_NumTriggerTracks]/F");
470  fJetTree->Branch("Jet_TriggerTrack_dEta",&dummy,"Jet_TriggerTrack_dEta[Jet_NumTriggerTracks]/F");
471  fJetTree->Branch("Jet_TriggerTrack_dPhi",&dummy,"Jet_TriggerTrack_dPhi[Jet_NumTriggerTracks]/F");
472  }
473 
474  fInitialized = kTRUE;
475 }
476 
477 //________________________________________________________________________
479  AliAnalysisTaskEmcalJet("AliAnalysisTaskJetExtractor", kTRUE),
480  fSaveConstituents(0), fSaveConstituentsIP(0), fSaveConstituentPID(0), fSaveJetShapes(0), fSaveJetSplittings(0), fSaveMCInformation(0), fSaveSecondaryVertices(0), fSaveTriggerTracks(0), fSaveCaloClusters(0),
481  fEventPercentage(1.0),
482  fEventCut_TriggerTrackMinPt(0),
483  fEventCut_TriggerTrackMaxPt(0),
484  fEventCut_TriggerTrackMinLabel(-9999999),
485  fEventCut_TriggerTrackMaxLabel(+9999999),
486  fEventCut_TriggerTrackOrigin(-1),
487  fTruthMinLabel(0),
488  fTruthMaxLabel(100000),
489  fHadronMatchingRadius(0.4),
490  fMatchedJetsArrayName(""),
491  fMatchedJetsRhoName(""),
492  fMatchedJetsRhoMassName(""),
493  fMCParticleArrayName("mcparticles"),
494  fRandomSeed(0),
495  fRandomSeedCones(0),
496  fVertexerCuts(0),
497  fSecondaryVertexMaxChi2(1e10),
498  fSecondaryVertexMaxDispersion(0.05),
499  fCustomStartupScript(),
500  fSetEmcalJetFlavour(0),
501  fSaveTrackPDGCode(kTRUE),
502  fJetTree(0),
503  fEventWeight(0.),
504  fImpactParameter(0.),
505  fMultiplicity(0),
506  fTriggerTracks_Pt(),
507  fTriggerTracks_Eta(),
508  fTriggerTracks_Phi(),
509  fMCParticleArray(0),
510  fRandomGenerator(0),
511  fRandomGeneratorCones(0),
512  fVtxTagger(0),
513  fIsEmbeddedEvent(kFALSE),
514  fSimpleSecVertices()
515 {
516  fRandomGenerator = new TRandom3();
517  fRandomGeneratorCones = new TRandom3();
518 
520  fJetTree = new AliEmcalJetTree(GetName());
521  DefineOutput(2, TTree::Class());
522 }
523 
524 //________________________________________________________________________
526  AliAnalysisTaskEmcalJet(name, kTRUE),
528  fEventPercentage(1.0),
534  fTruthMinLabel(0),
535  fTruthMaxLabel(100000),
540  fMCParticleArrayName("mcparticles"),
541  fRandomSeed(0),
542  fRandomSeedCones(0),
543  fVertexerCuts(0),
548  fSaveTrackPDGCode(kTRUE),
549  fJetTree(0),
550  fEventWeight(0.),
551  fImpactParameter(0.),
552  fMultiplicity(0),
556  fMCParticleArray(0),
557  fRandomGenerator(0),
559  fVtxTagger(0),
560  fIsEmbeddedEvent(kFALSE),
562 {
563  fRandomGenerator = new TRandom3();
564  fRandomGeneratorCones = new TRandom3();
565 
567  fJetTree = new AliEmcalJetTree(GetName());
568  DefineOutput(2, TTree::Class());
569 }
570 
571 //________________________________________________________________________
573 {
574  // Destructor.
575 }
576 
577 //________________________________________________________________________
579 {
581 
582  // ### Basic container settings
583  if(!GetJetContainer(0))
584  AliFatal("Jet input container not found!");
586  if(!GetParticleContainer(0))
587  AliFatal("At least one particle input container needs to be added");
589  AliFatal("Cluster input container not found although cluster extraction demanded!");
590 
591  fRandomGenerator->SetSeed(fRandomSeed);
593 
594  // Activate saving of trigger tracks if this is demanded
596  fSaveTriggerTracks = kTRUE;
597 
598  // ### Initialize the jet tree (settings must all be given at this stage)
601  OpenFile(2);
602  PostData(2, fJetTree->GetTreePointer());
603 
604  // ### Add control histograms (already some created in base task)
605  AddHistogram2D<TH2D>("hTrackCount", "Number of tracks in acceptance vs. centrality", "COLZ", 500, 0., 5000., 100, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}");
606  AddHistogram2D<TH2D>("hBackgroundPt", "Background p_{T} distribution", "", 150, 0., 150., 100, 0, 100, "Background p_{T} (GeV/c)", "Centrality", "dN^{Events}/dp_{T}");
607 
608  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}");
609  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}");
610  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}");
611  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");
612  AddHistogram2D<TH2D>("hJetArea", "Jet area", "COLZ", 200, 0., 2., 100, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
613  AddHistogram2D<TH2D>("hDeltaPt", "#delta p_{T} distribution", "", 400, -100., 300., 100, 0, 100, "p_{T, cone} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
614  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}");
615  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");
616  AddHistogram1D<TH1D>("hExtractionPercentage", "Extracted jets p_{T} distribution (background subtracted)", "e1p", 400, -100., 300., "p_{T, jet} (GeV/c)", "Extracted percentage");
617 
618  // Track QA plots
619  AddHistogram2D<TH2D>("hTrackPt", "Tracks p_{T} distribution", "", 300, 0., 300., 100, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
620  AddHistogram2D<TH2D>("hTrackPhi", "Track angular distribution in #phi", "LEGO2", 180, 0., 2*TMath::Pi(), 100, 0, 100, "#phi", "Centrality", "dN^{Tracks}/(d#phi)");
621  AddHistogram2D<TH2D>("hTrackEta", "Track angular distribution in #eta", "LEGO2", 100, -2.5, 2.5, 100, 0, 100, "#eta", "Centrality", "dN^{Tracks}/(d#eta)");
622  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");
623 
624  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})");
625  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})");
626 }
627 
628 
629 //________________________________________________________________________
631 {
633 
634 
635  // ### If there is an embedding track container, set flag that an embedded event is used
636  for(Int_t iCont=0; iCont<fParticleCollArray.GetEntriesFast(); iCont++)
637  if (GetParticleContainer(iCont)->GetIsEmbedding())
638  fIsEmbeddedEvent = kTRUE;
639 
640  // ### Need to explicitly tell jet container to load rho mass object
641  GetJetContainer(0)->LoadRhoMass(InputEvent());
642 
643  if (fMCParticleArrayName != "")
644  {
645  if(!fIsEmbeddedEvent)
646  fMCParticleArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCParticleArrayName.Data()));
647  else
648  {
649  // In case of embedding, the MC particle array needs to be fetched differently
650  AliVEvent* event = AliEmcalContainerUtils::GetEvent(InputEvent(), kTRUE);
651  fMCParticleArray = dynamic_cast<TClonesArray*>(event->FindListObject(fMCParticleArrayName.Data()));
652  }
653  }
654 
655  // ### Prepare vertexer
657  {
658  if(!fVertexerCuts)
659  AliFatal("VertexerCuts not given but secondary vertex calculation turned on.");
660  fVtxTagger = new AliHFJetsTaggingVertex();
661  fVtxTagger->SetCuts(fVertexerCuts);
662  }
663 
664  // ### Save tree extraction percentage to histogram
665  std::vector<Float_t> extractionPtBins = fJetTree->GetExtractionPercentagePtBins();
666  std::vector<Float_t> extractionPercentages = fJetTree->GetExtractionPercentages();
667 
668  for(Int_t i=0; i<static_cast<Int_t>(extractionPercentages.size()); i++)
669  {
670  Double_t percentage = extractionPercentages[i];
671  for(Int_t pt=static_cast<Int_t>(extractionPtBins[i*2]); pt<static_cast<Int_t>(extractionPtBins[i*2+1]); pt++)
672  {
673  FillHistogram("hExtractionPercentage", pt, percentage);
674  }
675  }
676 
677  // ### Add HIJING/JEWEL histograms (in case we need it)
678  if(MCEvent())
679  {
680  if(dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()))
681  AddHistogram1D<TH1D>("hHIJING_Ncoll", "N_{coll} from HIJING", "e1p", 1000, 0., 5000, "N_{coll}", "dN^{Events}/dN^{N_{coll}}");
682 
683  if(dynamic_cast<AliGenHepMCEventHeader*>(MCEvent()->GenEventHeader()))
684  {
685  AddHistogram1D<TH1D>("hJEWEL_NProduced", "NProduced from JEWEL", "e1p", 1000, 0., 1000, "Nproduced", "dN^{Events}/dN^{Nproduced}");
686  AddHistogram1D<TH1D>("hJEWEL_ImpactParameter", "Impact parameter from JEWEL", "e1p", 1000, 0., 1, "IP", "dN^{Events}/dN^{IP}");
687  TProfile* evweight = new TProfile("hJEWEL_EventWeight", "Event weight from JEWEL", 1, 0,1);
688  fOutput->Add(evweight);
689  }
690  }
691 
692  // ### Execute shell script at startup
693  if(!fCustomStartupScript.IsNull())
694  {
695  TGrid::Connect("alien://");
696  TFile::Cp(fCustomStartupScript.Data(), "./myScript.sh");
697  gSystem->Exec("bash ./myScript.sh");
698  }
699 
700  PrintConfig();
701 
702 }
703 
704 //________________________________________________________________________
706 {
707  // ################################### EVENT SELECTION
708 
709  if(!IsTriggerTrackInEvent())
710  return kFALSE;
711  if(fRandomGenerator->Rndm() >= fEventPercentage)
712  return kFALSE;
713 
714  // ################################### EVENT PROPERTIES
716 
717  // Load vertex if possible
718  Long64_t eventID = 0;
719  const AliVVertex* myVertex = InputEvent()->GetPrimaryVertex();
720  if(!myVertex && MCEvent())
721  myVertex = MCEvent()->GetPrimaryVertex();
722  Double_t vtx[3] = {0, 0, 0};
723  if(myVertex)
724  {
725  vtx[0] = myVertex->GetX(); vtx[1] = myVertex->GetY(); vtx[2] = myVertex->GetZ();
726  }
727 
728  // Get event ID from header
729  AliVHeader* eventIDHeader = InputEvent()->GetHeader();
730  if(eventIDHeader)
731  eventID = eventIDHeader->GetEventIdAsLong();
732 
733  // If available, get Ncoll from HIJING
734  if(MCEvent())
735  {
736  if(AliGenHijingEventHeader* mcHeader = dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()))
737  {
738  Double_t ncoll = mcHeader->NN() + mcHeader->NNw() + mcHeader->NwN() + mcHeader->NwNw();
739  FillHistogram("hHIJING_Ncoll", ncoll);
740  }
741  if(AliGenHepMCEventHeader* mcHeader = dynamic_cast<AliGenHepMCEventHeader*>(MCEvent()->GenEventHeader()))
742  {
743  fEventWeight = mcHeader->EventWeight();
744  fImpactParameter = mcHeader->impact_parameter();
745 
746  TProfile* tmpHist = static_cast<TProfile*>(fOutput->FindObject("hJEWEL_EventWeight"));
747  tmpHist->Fill(0., fEventWeight);
748  FillHistogram("hJEWEL_NProduced", mcHeader->NProduced());
749  FillHistogram("hJEWEL_ImpactParameter", fImpactParameter);
750  }
751  }
752 
753  // ################################### MAIN JET LOOP
755  Int_t jetCount = 0;
756  while(AliEmcalJet *jet = GetJetContainer(0)->GetNextAcceptJet())
757  {
758  // LOCAL BUFFER (Will be set for each jet, only added to tree if jet is accepted )
759  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;
760  std::vector<Float_t> vec_d0; std::vector<Float_t> vec_d0cov; std::vector<Float_t> vec_z0; std::vector<Float_t> vec_z0cov;
761  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;
762  std::vector<Float_t> triggerTracks_dEta(fTriggerTracks_Eta);
763  std::vector<Float_t> triggerTracks_dPhi(fTriggerTracks_Phi);
764  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;
765 
768  {
769  Double_t matchedJetDistance = 0;
770  Double_t matchedJetPt = 0;
771  Double_t matchedJetMass = 0;
772  Double_t truePtFraction = 0;
773  Double_t truePtFraction_mcparticles = 0;
774  Int_t currentJetType_HM = 0;
775  Int_t currentJetType_PM = 0;
776  Int_t currentJetType_IC = 0;
777  // Get jet type from MC (hadron matching, parton matching definition - for HF jets)
778  GetJetType(jet, currentJetType_HM, currentJetType_PM, currentJetType_IC);
779  // Get true estimators: for pt, jet mass, ...
780  GetTrueJetPtFraction(jet, truePtFraction, truePtFraction_mcparticles);
781  GetMatchedJetObservables(jet, matchedJetPt, matchedJetMass, matchedJetDistance);
782  fJetTree->FillBuffer_MonteCarlo(currentJetType_PM,currentJetType_HM,currentJetType_IC,matchedJetDistance,matchedJetPt,matchedJetMass,truePtFraction,truePtFraction_mcparticles,fPtHard,fEventWeight,fImpactParameter);
783  }
784 
785  // ### CONSTITUENT LOOP: Retrieve PID values + impact parameters
787  {
788  for(Int_t i = 0; i < jet->GetNumberOfParticleConstituents(); i++)
789  {
790  const AliVParticle* particle = jet->GetParticleConstituents()[i].GetParticle();
791  if(!particle) continue;
793  {
794  Float_t sigITS = 0; Float_t sigTPC = 0; Float_t sigTOF = 0; Float_t sigTRD = 0; Short_t recoPID = 0; Int_t truePID = 0;
795  AddPIDInformation(const_cast<AliVParticle*>(particle), sigITS, sigTPC, sigTOF, sigTRD, recoPID, truePID);
796  vecSigITS.push_back(sigITS); vecSigTPC.push_back(sigTPC); vecSigTOF.push_back(sigTOF); vecSigTRD.push_back(sigTRD); vecRecoPID.push_back(recoPID); vecTruePID.push_back(truePID);
797  fJetTree->FillBuffer_PID(vecSigITS, vecSigTPC, vecSigTOF, vecSigTRD, vecRecoPID, vecTruePID);
798  }
800  {
801  Float_t d0 = 0; Float_t d0cov = 0; Float_t z0 = 0; Float_t z0cov = 0;
802  GetTrackImpactParameters(myVertex, dynamic_cast<const AliAODTrack*>(particle), d0, d0cov, z0, z0cov);
803  vec_d0.push_back(d0); vec_d0cov.push_back(d0cov); vec_z0.push_back(z0); vec_z0cov.push_back(z0cov);
804  fJetTree->FillBuffer_ImpactParameters(vec_d0, vec_z0, vec_d0cov, vec_z0cov);
805  }
806  }
807  }
808 
809  // Reconstruct secondary vertices
811  {
812  ReconstructSecondaryVertices(myVertex, jet, secVtx_X, secVtx_Y, secVtx_Z, secVtx_Mass, secVtx_Lxy, secVtx_SigmaLxy, secVtx_Chi2, secVtx_Dispersion);
813  fJetTree->FillBuffer_SecVertices(secVtx_X, secVtx_Y, secVtx_Z, secVtx_Mass, secVtx_Lxy, secVtx_SigmaLxy, secVtx_Chi2, secVtx_Dispersion);
814  }
815 
816  // Now change trigger tracks eta/phi to dEta/dPhi relative to jet and save them
818  {
819  for(UInt_t i=0; i<triggerTracks_dEta.size(); i++)
820  {
821  triggerTracks_dEta[i] = jet->Eta() - fTriggerTracks_Eta[i];
822  triggerTracks_dPhi[i] = TMath::Min(TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]), TMath::TwoPi() - TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]));
823  if( ((TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]) <= TMath::Pi()) && (jet->Phi() - fTriggerTracks_Phi[i] > 0))
824  || ((TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]) > TMath::Pi()) && (jet->Phi() - fTriggerTracks_Phi[i] < 0)) )
825  triggerTracks_dPhi[i] = -triggerTracks_dPhi[i];
826  }
827  fJetTree->FillBuffer_TriggerTracks(fTriggerTracks_Pt, triggerTracks_dEta, triggerTracks_dPhi);
828  }
829 
830  if(fSaveJetShapes)
831  {
832  // Calculate jet shapes and set them in the tree (some are retrieved in the tree itself)
833  Double_t leSub_noCorr = 0;
834  Double_t angularity = 0;
835  Double_t momentumDispersion = 0;
836  Double_t trackPtMean = 0;
837  Double_t trackPtMedian = 0;
838  CalculateJetShapes(jet, leSub_noCorr, angularity, momentumDispersion, trackPtMean, trackPtMedian);
839  fJetTree->FillBuffer_JetShapes(jet, leSub_noCorr, angularity, momentumDispersion, trackPtMean, trackPtMedian);
840  }
841 
843  {
844  GetJetSplittings(jet, splittings_radiatorE, splittings_kT, splittings_theta, splittings_secVtx_rank, splittings_secVtx_index);
845  fJetTree->FillBuffer_Splittings(splittings_radiatorE, splittings_kT, splittings_theta, fSaveSecondaryVertices, splittings_secVtx_rank, splittings_secVtx_index);
846  }
847 
848  // Fill jet to tree (here adding the minimum properties)
849  Bool_t accepted = fJetTree->AddJetToTree(jet, fSaveConstituents, fSaveConstituentsIP, fSaveCaloClusters, vtx, GetJetContainer(0)->GetRhoVal(), GetJetContainer(0)->GetRhoMassVal(), fCent, fMultiplicity, eventID, InputEvent()->GetMagneticField());
850  if(accepted)
851  FillHistogram("hJetPtExtracted", jet->Pt() - GetJetContainer(0)->GetRhoVal()*jet->Area(), fCent);
852  jetCount++;
853  }
854 
855  // ################################### PARTICLE PROPERTIES
856  for(Int_t iCont=0; iCont<fParticleCollArray.GetEntriesFast(); iCont++)
857  {
859  while(AliVTrack *track = static_cast<AliVTrack*>(GetParticleContainer(iCont)->GetNextAcceptParticle()))
861  }
862 
863  return kTRUE;
864 }
865 
866 
867 //________________________________________________________________________
869 {
870  // Cut for trigger track requirement
872  {
873  // Clear vector of trigger tracks
874  fTriggerTracks_Pt.clear();
875  fTriggerTracks_Eta.clear();
876  fTriggerTracks_Phi.clear();
877 
878  // Go through all tracks, in all attached track containers, and check whether trigger tracks can be found
879  for(Int_t iCont=0; iCont<fParticleCollArray.GetEntriesFast(); iCont++)
880  {
881  AliParticleContainer* partCont = GetParticleContainer(iCont);
882  partCont->ResetCurrentID();
883  if((fEventCut_TriggerTrackOrigin == 1) && !partCont->GetIsEmbedding())
884  continue;
885  else if((fEventCut_TriggerTrackOrigin == 2) && partCont->GetIsEmbedding())
886  continue;
887 
888  while(AliVTrack *track = static_cast<AliVTrack*>(partCont->GetNextAcceptParticle()))
889  {
890  if( (track->GetLabel() >= fEventCut_TriggerTrackMinLabel) && (track->GetLabel() < fEventCut_TriggerTrackMaxLabel) )
891  if( (track->Pt() >= fEventCut_TriggerTrackMinPt) && (track->Pt() < fEventCut_TriggerTrackMaxPt) )
892  {
893  fTriggerTracks_Pt.push_back(track->Pt());
894  fTriggerTracks_Eta.push_back(track->Eta());
895  fTriggerTracks_Phi.push_back(track->Phi());
896  }
897  }
898  }
899  // No particle has been found that fulfills requirement -> Do not accept event
900  if(fTriggerTracks_Pt.size() == 0)
901  return kFALSE;
902  }
903  return kTRUE;
904 }
905 
906 //________________________________________________________________________
907 void AliAnalysisTaskJetExtractor::CalculateJetShapes(AliEmcalJet* jet, Double_t& leSub_noCorr, Double_t& angularity, Double_t& momentumDispersion, Double_t& trackPtMean, Double_t& trackPtMedian)
908 {
909  // #### Calculate mean, median of constituents, radial moment (angularity), momentum dispersion, leSub (no correction)
910  Double_t jetLeadingHadronPt = -999.;
911  Double_t jetSubleadingHadronPt = -999.;
912  Double_t jetSummedPt = 0;
913  Double_t jetSummedPt2 = 0;
914  trackPtMean = 0;
915  trackPtMedian = 0;
916  angularity = 0;
917  momentumDispersion = 0;
918  std::vector<PWG::JETFW::AliEmcalParticleJetConstituent> tracks_sorted = jet->GetParticleConstituents();
919  std::sort(tracks_sorted.rbegin(), tracks_sorted.rend());
920  Int_t numTracks = tracks_sorted.size();
921  if(!numTracks) return;
922  Double_t* constPts = new Double_t[numTracks];
923 
924  // Loop over all constituents and do jet shape calculations
925  for (Int_t i=0;i<numTracks;i++)
926  {
927  const AliVParticle* particle = tracks_sorted[i].GetParticle();
928  trackPtMean += particle->Pt();
929  constPts[i] = particle->Pt();
930  if(particle->Pt() > jetLeadingHadronPt)
931  {
932  jetSubleadingHadronPt = jetLeadingHadronPt;
933  jetLeadingHadronPt = particle->Pt();
934  }
935  else if(particle->Pt() > jetSubleadingHadronPt)
936  jetSubleadingHadronPt = particle->Pt();
937 
938  Double_t deltaR = GetDistance(particle->Eta(), jet->Eta(), particle->Phi(), jet->Phi());
939  jetSummedPt += particle->Pt();
940  jetSummedPt2 += particle->Pt()*particle->Pt();
941  angularity += particle->Pt() * deltaR;
942  }
943 
944  if(numTracks)
945  {
946  trackPtMean /= numTracks;
947  trackPtMedian = TMath::Median(numTracks, constPts);
948  }
949 
950  if(numTracks > 1)
951  leSub_noCorr = jetLeadingHadronPt - jetSubleadingHadronPt;
952  else
953  leSub_noCorr = jetLeadingHadronPt;
954 
955  if(jetSummedPt)
956  {
957  momentumDispersion = TMath::Sqrt(jetSummedPt2)/jetSummedPt;
958  angularity /= jetSummedPt;
959  }
960 }
961 
962 //________________________________________________________________________
963 void AliAnalysisTaskJetExtractor::GetTrueJetPtFraction(AliEmcalJet* jet, Double_t& truePtFraction, Double_t& truePtFraction_mcparticles)
964 {
965  // #################################################################################
966  // ##### FRACTION OF TRUE PT IN JET: Defined as "not from toy"
967  Double_t pt_truth = 0.;
968  Double_t pt_truth_mcparticles = 0.;
969  Double_t pt_all = 0.;
970  truePtFraction = 0;
971  truePtFraction_mcparticles = 0;
972 
973  // ### Loop over all tracks constituents
974  for(Int_t iConst = 0; iConst < jet->GetNumberOfParticleConstituents(); iConst++)
975  {
976  const AliVParticle* particle = jet->GetParticleConstituents()[iConst].GetParticle();
977  if(!particle) continue;
978  if(particle->Pt() < 1e-6) continue;
979 
980  // Particles marked w/ labels within label range OR explicitly set as embedded tracks are considered to be from truth
981  if ( (fIsEmbeddedEvent && jet->GetParticleConstituents()[iConst].IsFromEmbeddedEvent()) ||
982  (!fIsEmbeddedEvent && ((particle->GetLabel() >= fTruthMinLabel) && (particle->GetLabel() < fTruthMaxLabel))) )
983  pt_truth += particle->Pt();
984  pt_all += particle->Pt();
985  }
986 
987  // ### Loop over all cluster constituents
988  for(Int_t iConst = 0; iConst < jet->GetNumberOfClusterConstituents(); iConst++)
989  {
990  const AliVCluster* cluster = jet->GetClusterConstituents()[iConst].GetCluster();
991  if(!cluster) continue;
992  if(cluster->E() < 1e-6) continue;
993 
994  // #### Retrieve cluster pT
995  Double_t vtxX = 0;
996  Double_t vtxY = 0;
997  Double_t vtxZ = 0;
998  const AliVVertex* myVertex = InputEvent()->GetPrimaryVertex();
999  if(!myVertex && MCEvent())
1000  myVertex = MCEvent()->GetPrimaryVertex();
1001  if(myVertex)
1002  {
1003  vtxX = myVertex->GetX();
1004  vtxY = myVertex->GetY();
1005  vtxZ = myVertex->GetZ();
1006  }
1007 
1008  Double_t primVtx[3] = {vtxX, vtxY, vtxZ};
1009  TLorentzVector clusterMomentum;
1010  cluster->GetMomentum(clusterMomentum, primVtx);
1011  Double_t ClusterPt = clusterMomentum.Perp();
1012  // ####
1013 
1014  // Particles marked w/ labels within label range OR explicitly set as embedded clusters are considered to be from truth
1015  if ( (fIsEmbeddedEvent && jet->GetClusterConstituents()[iConst].IsFromEmbeddedEvent()) ||
1016  (!fIsEmbeddedEvent && ((cluster->GetLabel() >= fTruthMinLabel) && (cluster->GetLabel() < fTruthMaxLabel))) )
1017  pt_truth += ClusterPt;
1018  pt_all += ClusterPt;
1019  }
1020 
1021  // ### Loop over all primary (charged) MC particles and check if they have a corresponding track/cluster
1022  // Correspondence is checked geometrically, sum of matched particles pT is truth
1024  Double_t jetRadius = GetJetContainer(0)->GetJetRadius();
1025  if(fMCParticleArray)
1026  for(Int_t iPart=0; iPart<fMCParticleArray->GetEntries();iPart++)
1027  {
1028  AliAODMCParticle* part = (AliAODMCParticle*)fMCParticleArray->At(iPart);
1029  if(!part) continue;
1030  if(!part->IsPhysicalPrimary()) continue;
1031  if(!fulljets && !part->Charge()) continue;
1032  if(part->Pt() < 1e-6) continue;
1033 
1034  if(IsTrackInCone(part, jet->Eta(), jet->Phi(), jetRadius))
1035  pt_truth_mcparticles += part->Pt();
1036  }
1037 
1038  if(pt_all)
1039  {
1040  truePtFraction = (pt_truth/pt_all);
1041  truePtFraction_mcparticles = (pt_truth_mcparticles/pt_all);
1042  }
1043 }
1044 
1045 //________________________________________________________________________
1046 void AliAnalysisTaskJetExtractor::GetMatchedJetObservables(AliEmcalJet* jet, Double_t& matchedJetPt, Double_t& matchedJetMass, Double_t& matchedJetDistance)
1047 {
1048  // #################################################################################
1049  // ##### OBSERVABLES FROM MATCHED JETS: Jet pt, jet mass
1050  matchedJetDistance = 8.; // 8 is higher than maximum possible matching distance
1051  if(fMatchedJetsArrayName != "")
1052  {
1053  // "True" background for pt
1054  AliRhoParameter* rho = static_cast<AliRhoParameter*>(InputEvent()->FindListObject(fMatchedJetsRhoName.Data()));
1055  Double_t trueRho = 0;
1056  if(rho)
1057  trueRho = rho->GetVal();
1058 
1059  // "True" background for mass
1060  AliRhoParameter* rhoMass = static_cast<AliRhoParameter*>(InputEvent()->FindListObject(fMatchedJetsRhoMassName.Data()));
1061  TClonesArray* matchedJetArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fMatchedJetsArrayName.Data())));
1062 
1063  // Loop over all true jets to find the best match
1064  matchedJetPt = 0;
1065  matchedJetMass = 0;
1066  if(matchedJetArray)
1067  {
1068  for(Int_t i=0; i<matchedJetArray->GetEntries(); i++)
1069  {
1070  AliEmcalJet* matchedJet = static_cast<AliEmcalJet*>(matchedJetArray->At(i));
1071  if(matchedJet->Pt() < 0.15)
1072  continue;
1073  Double_t deltaR = GetDistance(matchedJet->Eta(), jet->Eta(), matchedJet->Phi(), jet->Phi());
1074  // Search for the best match
1075  if(deltaR < matchedJetDistance)
1076  {
1077  matchedJetDistance = deltaR;
1078  matchedJetPt = matchedJet->Pt() - matchedJet->Area()* trueRho;
1079  if(rhoMass)
1080  matchedJetMass = matchedJet->GetShapeProperties()->GetSecondOrderSubtracted();
1081  else
1082  matchedJetMass = matchedJet->M();
1083  }
1084  }
1085  }
1086  }
1087 }
1088 
1089 
1090 //________________________________________________________________________
1092 {
1093  if(!fMCParticleArray)
1094  return;
1095 
1096  typeHM = 0;
1097  typePM = 0;
1098 
1099  AliAODMCParticle* parton[2];
1100  parton[0] = (AliAODMCParticle*) fVtxTagger->IsMCJetParton(fMCParticleArray, jet, fHadronMatchingRadius); // method 1 (parton matching)
1101  parton[1] = (AliAODMCParticle*) fVtxTagger->IsMCJetMeson(fMCParticleArray, jet, fHadronMatchingRadius); // method 2 (hadron matching)
1102 
1103  if (parton[0]) {
1104  Int_t pdg = TMath::Abs(parton[0]->PdgCode());
1105  typePM = pdg;
1106  }
1107 
1108  if (!parton[1])
1109  {
1110  // No HF jet (according to hadron matching) -- now also separate jets in udg (1) and s-jets (3)
1111  if(IsStrangeJet(jet))
1112  typeHM = 3;
1113  else
1114  typeHM = 1;
1115  }
1116  else {
1117  Int_t pdg = TMath::Abs(parton[1]->PdgCode());
1118  if(fVtxTagger->IsDMeson(pdg)) typeHM = 4;
1119  else if (fVtxTagger->IsBMeson(pdg)) typeHM = 5;
1120  }
1121 
1122  // Set flavour of AliEmcalJet object (set ith bit while i corresponds to type)
1124  jet->AddFlavourTag(static_cast<Int_t>(TMath::Power(2, typeHM)));
1125 
1126  const AliEmcalPythiaInfo* partonsInfo = GetPythiaInfo();
1127  typeIC = 0;
1128  if (partonsInfo)
1129  {
1130  // Get primary partons directions
1131  Double_t parton1phi = partonsInfo->GetPartonPhi6();
1132  Double_t parton1eta = partonsInfo->GetPartonEta6();
1133  Double_t parton2phi = partonsInfo->GetPartonPhi7();
1134  Double_t parton2eta = partonsInfo->GetPartonEta7();
1135 
1136  Double_t delta1R = GetDistance(parton1eta, jet->Eta(), parton1phi, jet->Phi());
1137  Double_t delta2R = GetDistance(parton2eta, jet->Eta(), parton2phi, jet->Phi());
1138 
1139  // Check if one of the partons if closer than matching criterion
1140  Bool_t matched = (delta1R < GetJetContainer(0)->GetJetRadius()/2.) || (delta2R < GetJetContainer(0)->GetJetRadius()/2.);
1141 
1142  // Matching criterion fulfilled -> Set flag to closest
1143  if(matched)
1144  {
1145  if(delta1R < delta2R)
1146  typeIC = partonsInfo->GetPartonFlag6();
1147  else
1148  typeIC = partonsInfo->GetPartonFlag7();
1149  }
1150  }
1151 }
1152 
1153 //________________________________________________________________________
1154 void AliAnalysisTaskJetExtractor::GetTrackImpactParameters(const AliVVertex* vtx, const AliAODTrack* track, Float_t& d0, Float_t& d0cov, Float_t& z0, Float_t& z0cov)
1155 {
1156  if (track)
1157  {
1158  AliAODTrack localTrack(*track);
1159  Double_t d0rphiz[2],covd0[3];
1160  Bool_t isDCA=localTrack.PropagateToDCA(vtx,InputEvent()->GetMagneticField(),3.0,d0rphiz,covd0);
1161  if(isDCA)
1162  {
1163  d0 = d0rphiz[0];
1164  z0 = d0rphiz[1];
1165  d0cov = covd0[0];
1166  z0cov = covd0[2];
1167  }
1168  }
1169 }
1170 
1171 //________________________________________________________________________
1172 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)
1173 {
1174  if(!primVtx)
1175  return;
1176 
1177  // Create ESD vertex from the existing AliVVertex
1178  Double_t vtxPos[3] = {primVtx->GetX(), primVtx->GetY(), primVtx->GetZ()};
1179  Double_t covMatrix[6] = {0};
1180  primVtx->GetCovarianceMatrix(covMatrix);
1181  AliESDVertex* esdVtx = new AliESDVertex(vtxPos, covMatrix, primVtx->GetChi2(), primVtx->GetNContributors());
1182 
1183  TClonesArray* secVertexArr = 0;
1184  vctr_pair_dbl_int arrDispersion;
1185  arrDispersion.reserve(5);
1186  //###########################################################################
1187  // ********* Calculate secondary vertices
1188  // Derived from AliAnalysisTaskEmcalJetBtagSV
1189  secVertexArr = new TClonesArray("AliAODVertex");
1190  Int_t nDauRejCount = 0;
1191  Int_t nVtx = fVtxTagger->FindVertices(jet,
1192  GetParticleContainer(0)->GetArray(),
1193  (AliAODEvent*)InputEvent(),
1194  esdVtx,
1195  InputEvent()->GetMagneticField(),
1196  secVertexArr,
1197  0,
1198  arrDispersion,
1199  nDauRejCount);
1200 
1201  if(nVtx < 0)
1202  {
1203  secVertexArr->Clear();
1204  delete secVertexArr;
1205  return;
1206  }
1207  //###########################################################################
1208 
1209  // Loop over all potential secondary vertices
1210  fSimpleSecVertices.clear();
1211  for(Int_t i=0; i<secVertexArr->GetEntriesFast(); i++)
1212  {
1213  AliAODVertex* secVtx = (AliAODVertex*)(secVertexArr->UncheckedAt(i));
1214 
1215  // Calculate vtx distance
1216  Double_t effX = secVtx->GetX() - esdVtx->GetX();
1217  Double_t effY = secVtx->GetY() - esdVtx->GetY();
1218  //Double_t effZ = secVtx->GetZ() - esdVtx->GetZ();
1219 
1220  // ##### Vertex properties
1221  // vertex dispersion
1222  Double_t dispersion = arrDispersion[i].first;
1223 
1224  // invariant mass
1225  Double_t mass = fVtxTagger->GetVertexInvariantMass(secVtx);
1226 
1227  // signed length
1228  Double_t Lxy = TMath::Sqrt(effX*effX + effY*effY);
1229  Double_t jetP[3]; jet->PxPyPz(jetP);
1230  Double_t signLxy = effX * jetP[0] + effY * jetP[1];
1231  if (signLxy < 0.) Lxy *= -1.;
1232 
1233  Double_t sigmaLxy = 0;
1234  AliAODVertex* aodVtx = (AliAODVertex*)(primVtx);
1235  if (aodVtx)
1236  sigmaLxy = aodVtx->ErrorDistanceXYToVertex(secVtx);
1237 
1238  // Add secondary vertices if they fulfill the conditions
1239  if( (dispersion > fSecondaryVertexMaxDispersion) || (TMath::Abs(secVtx->GetChi2perNDF()) > fSecondaryVertexMaxChi2) )
1240  continue;
1241 
1242  // Internally, save sec. vertices to a list
1243  // Each secondary vertex is reconstructed from 3 prongs
1245  vtx.fIndex = secVtx_X.size();
1246  vtx.fLxy = TMath::Abs(Lxy);
1247  vtx.fDaughter1 = static_cast<AliVParticle*>(secVtx->GetDaughter(0));
1248  vtx.fDaughter2 = static_cast<AliVParticle*>(secVtx->GetDaughter(1));
1249  vtx.fDaughter3 = static_cast<AliVParticle*>(secVtx->GetDaughter(2));
1250  fSimpleSecVertices.push_back(vtx);
1251 
1252  secVtx_X.push_back(secVtx->GetX()); secVtx_Y.push_back(secVtx->GetY()); secVtx_Z.push_back(secVtx->GetZ()); secVtx_Chi2.push_back(secVtx->GetChi2perNDF());
1253  secVtx_Dispersion.push_back(dispersion); secVtx_Mass.push_back(mass); secVtx_Lxy.push_back(Lxy); secVtx_SigmaLxy.push_back(sigmaLxy);
1254  }
1255 
1256  // Sort simple sec. vertices w/ descending Lxy
1257  std::sort(fSimpleSecVertices.rbegin(), fSimpleSecVertices.rend(), [](const SimpleSecondaryVertex a, const SimpleSecondaryVertex b) {return a.fLxy<b.fLxy;});
1258  if(fSimpleSecVertices.size() > 10) fSimpleSecVertices.resize(10);
1259 
1260  secVertexArr->Clear();
1261  delete secVertexArr;
1262  delete esdVtx;
1263 }
1264 
1265 //________________________________________________________________________
1267 {
1268  // Do hadron matching jet type tagging using mcparticles
1269  // ... if not explicitly deactivated
1270  if (fMCParticleArray)
1271  {
1272  for(Int_t i=0; i<fMCParticleArray->GetEntries();i++)
1273  {
1274  AliAODMCParticle* part = (AliAODMCParticle*)fMCParticleArray->At(i);
1275  if(!part) continue;
1276 
1277  // Check if the particle has strangeness
1278  Int_t absPDG = TMath::Abs(part->PdgCode());
1279  if ((absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
1280  {
1281  // Check if particle is in a radius around the jet
1282  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
1284  continue;
1285  else
1286  return kTRUE;
1287  }
1288  }
1289  }
1290  // As fallback, the MC stack will be tried
1291  else if(MCEvent() && (MCEvent()->Stack()))
1292  {
1293  AliStack* stack = MCEvent()->Stack();
1294  // Go through the whole particle stack
1295  for(Int_t i=0; i<stack->GetNtrack(); i++)
1296  {
1297  TParticle *part = stack->Particle(i);
1298  if(!part) continue;
1299 
1300  // Check if the particle has strangeness
1301  Int_t absPDG = TMath::Abs(part->GetPdgCode());
1302  if ((absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
1303  {
1304  // Check if particle is in a radius around the jet
1305  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
1307  continue;
1308  else
1309  return kTRUE;
1310  }
1311  }
1312  }
1313  return kFALSE;
1314 
1315 }
1316 //________________________________________________________________________
1318 {
1319  fMultiplicity = 0;
1320  for(Int_t iCont=0; iCont<fParticleCollArray.GetEntriesFast(); iCont++)
1322 
1323  // ### Event control plots
1324  FillHistogram("hTrackCount", fMultiplicity, fCent);
1325  FillHistogram("hBackgroundPt", GetJetContainer(0)->GetRhoVal(), fCent);
1326 }
1327 
1328 //________________________________________________________________________
1330 {
1331  // ### Jet control plots
1332  AliJetContainer* jetContainer = GetJetContainer(0);
1333  FillHistogram("hJetPtRaw", jet->Pt(), fCent);
1334  FillHistogram("hJetPt", jet->Pt() - jetContainer->GetRhoVal()*jet->Area(), fCent);
1335  FillHistogram("hJetPhiEta", jet->Phi(), jet->Eta());
1336  FillHistogram("hJetArea", jet->Area(), fCent);
1337 
1338  // ### Jet constituent plots
1339  for(Int_t i = 0; i < jet->GetNumberOfParticleConstituents(); i++)
1340  {
1341  const AliVParticle* particle = jet->GetParticleConstituents()[i].GetParticle();
1342  if(!particle) continue;
1343 
1344  // Constituent eta/phi (relative to jet)
1345  Double_t deltaEta = jet->Eta() - particle->Eta();
1346  Double_t deltaPhi = TMath::Min(TMath::Abs(jet->Phi() - particle->Phi()), TMath::TwoPi() - TMath::Abs(jet->Phi() - particle->Phi()));
1347  if(!((jet->Phi() - particle->Phi() < 0) && (jet->Phi() - particle->Phi() <= TMath::Pi())))
1348  deltaPhi = -deltaPhi;
1349 
1350  FillHistogram("hConstituentPt", particle->Pt(), fCent);
1351  FillHistogram("hConstituentPhiEta", deltaPhi, deltaEta);
1352  }
1353 
1354  // ### Random cone / delta pT plots
1355  const Int_t kNumRandomConesPerEvent = 4;
1356  for(Int_t iCone=0; iCone<kNumRandomConesPerEvent; iCone++)
1357  {
1358  // Throw random cone
1359  Double_t tmpRandConeEta = jetContainer->GetJetEtaMin() + fRandomGeneratorCones->Rndm()*TMath::Abs(jetContainer->GetJetEtaMax()-jetContainer->GetJetEtaMin());
1360  Double_t tmpRandConePhi = fRandomGeneratorCones->Rndm()*TMath::TwoPi();
1361  Double_t tmpRandConePt = 0;
1362  // Fill pT that is in cone
1363 
1364  for(Int_t iCont=0; iCont<fParticleCollArray.GetEntriesFast(); iCont++)
1365  {
1367  while(AliVTrack *track = static_cast<AliVTrack*>(GetParticleContainer(iCont)->GetNextAcceptParticle()))
1368  if(IsTrackInCone(track, tmpRandConeEta, tmpRandConePhi, jetContainer->GetJetRadius()))
1369  tmpRandConePt += track->Pt();
1370  }
1371 
1372  // Fill histograms
1373  FillHistogram("hDeltaPt", tmpRandConePt - jetContainer->GetRhoVal()*jetContainer->GetJetRadius()*jetContainer->GetJetRadius()*TMath::Pi(), fCent);
1374  }
1375 }
1376 
1377 //________________________________________________________________________
1379 {
1380  FillHistogram("hTrackPt", track->Pt(), fCent);
1381  FillHistogram("hTrackPhi", track->Phi(), fCent);
1382  FillHistogram("hTrackEta", track->Eta(), fCent);
1383  FillHistogram("hTrackEtaPt", track->Eta(), track->Pt());
1384  FillHistogram("hTrackPhiPt", track->Phi(), track->Pt());
1385  FillHistogram("hTrackPhiEta", track->Phi(), track->Eta());
1386 }
1387 
1388 //________________________________________________________________________
1389 void AliAnalysisTaskJetExtractor::AddPIDInformation(const AliVParticle* particle, Float_t& sigITS, Float_t& sigTPC, Float_t& sigTOF, Float_t& sigTRD, Short_t& recoPID, Int_t& truePID)
1390 {
1391  truePID = 9;
1392  if(!particle) return;
1393 
1394  // If we have AODs, retrieve particle PID signals
1395  const AliAODTrack* aodtrack = dynamic_cast<const AliAODTrack*>(particle);
1396 
1397  if(aodtrack)
1398  {
1399  // Get AOD value from reco
1400  recoPID = aodtrack->GetMostProbablePID();
1401  AliAODPid* pidObj = aodtrack->GetDetPid();
1402  if(!pidObj)
1403  return;
1404 
1405  sigITS = pidObj->GetITSsignal();
1406  sigTPC = pidObj->GetTPCsignal();
1407  sigTOF = pidObj->GetTOFsignal();
1408  sigTRD = pidObj->GetTRDsignal();
1409  }
1410 
1411  // Get truth values if we are on MC
1412  if(fMCParticleArray)
1413  {
1414  for(Int_t i=0; i<fMCParticleArray->GetEntries();i++)
1415  {
1416  AliAODMCParticle* mcParticle = dynamic_cast<AliAODMCParticle*>(fMCParticleArray->At(i));
1417  if(!mcParticle) continue;
1418 
1419  if (mcParticle->GetLabel() == particle->GetLabel())
1420  {
1421  if(fSaveTrackPDGCode)
1422  {
1423  truePID = mcParticle->PdgCode();
1424  }
1425  else // Use same convention as for PID in AODs
1426  {
1427  if(TMath::Abs(mcParticle->PdgCode()) == 2212) // proton
1428  truePID = 4;
1429  else if (TMath::Abs(mcParticle->PdgCode()) == 211) // pion
1430  truePID = 2;
1431  else if (TMath::Abs(mcParticle->PdgCode()) == 321) // kaon
1432  truePID = 3;
1433  else if (TMath::Abs(mcParticle->PdgCode()) == 11) // electron
1434  truePID = 0;
1435  else if (TMath::Abs(mcParticle->PdgCode()) == 13) // muon
1436  truePID = 1;
1437  else if (TMath::Abs(mcParticle->PdgCode()) == 700201) // deuteron
1438  truePID = 5;
1439  else if (TMath::Abs(mcParticle->PdgCode()) == 700301) // triton
1440  truePID = 6;
1441  else if (TMath::Abs(mcParticle->PdgCode()) == 700302) // He3
1442  truePID = 7;
1443  else if (TMath::Abs(mcParticle->PdgCode()) == 700202) // alpha
1444  truePID = 8;
1445  else
1446  truePID = 9;
1447  }
1448 
1449  break;
1450  }
1451  }
1452  }
1453 }
1454 
1455 
1456 //________________________________________________________________________
1457 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)
1458 {
1459  // ### Adapted from code in AliAnalysisTaskDmesonJetsSub ###
1460  // Define jet reclusterizer
1461  fastjet::JetAlgorithm jetAlgo(fastjet::cambridge_algorithm);
1462  Double_t jetRadius_CA = 1.0;
1463  fastjet::JetDefinition jetDefinition(jetAlgo, jetRadius_CA,static_cast<fastjet::RecombinationScheme>(0), fastjet::Best);
1464 
1465  try{
1466  // Convert jet constituents to vector of fastjet::PseudoJet
1467  std::vector<fastjet::PseudoJet> particles;
1468  for(Int_t iConst=0; iConst<jet->GetNumberOfParticleConstituents(); iConst++)
1469  {
1470  const AliVParticle* constituent = jet->GetParticleConstituents()[iConst].GetParticle();
1471  Double_t p[3];
1472  constituent->PxPyPz(p);
1473  fastjet::PseudoJet pseudoJet = fastjet::PseudoJet(p[0], p[1], p[2], constituent->E());
1474 
1475  // If constituent is part of the N most significant sec. vertices, mark it
1476  for(UInt_t iVtx=0; iVtx<fSimpleSecVertices.size(); iVtx++)
1477  {
1478  if((constituent == fSimpleSecVertices[iVtx].fDaughter1) || (constituent == fSimpleSecVertices[iVtx].fDaughter2) || (constituent == fSimpleSecVertices[iVtx].fDaughter3))
1479  {
1480  pseudoJet.set_user_index(iVtx); // user_index is now index in temp sec. vtx vector
1481  break;
1482  }
1483  }
1484  particles.push_back(pseudoJet);
1485  }
1486 
1487  // Perform jet reclustering
1488  fastjet::ClusterSequence clusterSeq_CA(particles, jetDefinition);
1489  std::vector<fastjet::PseudoJet> jets_CA = clusterSeq_CA.inclusive_jets(0);
1490  jets_CA = sorted_by_pt(jets_CA);
1491  fastjet::PseudoJet radiator = jets_CA[0];
1492  fastjet::PseudoJet leadingSubJet;
1493  fastjet::PseudoJet subleadingSubJet;
1494 
1495  // Iterate through the splitting history of the CA clusterization
1496  while(radiator.has_parents(leadingSubJet,subleadingSubJet))
1497  {
1498  if(leadingSubJet.perp() < subleadingSubJet.perp())
1499  std::swap(leadingSubJet,subleadingSubJet);
1500 
1501  // Angle theta
1502  Float_t theta = leadingSubJet.delta_R(subleadingSubJet);
1503  // Radiator energy
1504  Float_t radiatorEnergy = leadingSubJet.e()+subleadingSubJet.e();
1505  // kT
1506  Float_t kT = subleadingSubJet.perp()*theta;
1507 
1508  // Go through leading subjet constituents and check if there are tracks belonging to one of the ten most significant sec. vertices
1509  Int_t secVtx_rank = -1; // rank = nth most displaced
1510  Int_t secVtx_index = -1; // index = index in sec. vertex array
1511  if(fSimpleSecVertices.size())
1512  {
1513  std::vector<fastjet::PseudoJet> leadingConsts = leadingSubJet.constituents();
1514  for(UInt_t iLeadingConst=0; iLeadingConst<leadingConsts.size(); iLeadingConst++)
1515  {
1516  if(leadingConsts[iLeadingConst].user_index()>=0)
1517  {
1518  secVtx_rank = leadingConsts[iLeadingConst].user_index();
1519  secVtx_index = fSimpleSecVertices[leadingConsts[iLeadingConst].user_index()].fIndex;
1520  }
1521  }
1522  }
1523 
1524  // Now add splitting properties to result vectors
1525  splittings_radiatorE.push_back(radiatorEnergy);
1526  splittings_theta.push_back(theta);
1527  splittings_kT.push_back(kT);
1528  splittings_secVtx_rank.push_back(secVtx_rank);
1529  splittings_secVtx_index.push_back(secVtx_index);
1530 
1531  // Continue with leadingSubJet as new radiator
1532  radiator=leadingSubJet;
1533  }
1534  } catch (fastjet::Error) { /*return -1;*/ }
1535 }
1536 
1537 
1538 //________________________________________________________________________
1540 {
1541  std::cout << "#########################################\n";
1542  std::cout << "Settings for extractor task " << GetName() << std::endl;
1543  std::cout << std::endl;
1544 
1545  std::cout << "### Event cuts ###" << std::endl;
1546  std::cout << std::endl;
1548  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;
1550  std::cout << Form("* Trigger track event cut: %s", (fEventCut_TriggerTrackOrigin==1)? "From embedded event" : "Not from embedded event") << std::endl;
1551  if(fEventPercentage < 1)
1552  std::cout << Form("* Artificially lowered event efficiency: %f", fEventPercentage) << std::endl;
1554  std::cout << Form("* Random seeds explicitly set: %lu (for event/jet eff.), %lu (RCs)", fRandomSeed, fRandomSeedCones) << std::endl;
1555  std::cout << std::endl;
1556 
1557  std::cout << "### Jet properties ###" << std::endl;
1558  std::cout << "* Jet array name = " << GetJetContainer(0)->GetName() << std::endl;
1559  std::cout << "* Rho name = " << GetJetContainer(0)->GetRhoName() << std::endl;
1560  std::cout << "* Rho mass name = " << GetJetContainer(0)->GetRhoMassName() << std::endl;
1561  std::cout << std::endl;
1562 
1563  std::cout << "### Tree extraction properties ###" << std::endl;
1564  std::vector<Float_t> extractionPtBins = fJetTree->GetExtractionPercentagePtBins();
1565  std::vector<Float_t> extractionPercentages = fJetTree->GetExtractionPercentages();
1566  std::vector<Int_t> extractionHM = fJetTree->GetExtractionJetTypes_HM();
1567  std::vector<Int_t> extractionPM = fJetTree->GetExtractionJetTypes_PM();
1568  if(extractionPercentages.size())
1569  {
1570  std::cout << "* Extraction bins: (";
1571  for(Int_t i=0; i<static_cast<Int_t>(extractionPercentages.size()-1); i++)
1572  std::cout << extractionPtBins[i*2] << "-" << extractionPtBins[i*2+1] << " = " << extractionPercentages[i] << ", ";
1573  std::cout << extractionPtBins[(extractionPercentages.size()-1)*2] << "-" << extractionPtBins[(extractionPercentages.size()-1)*2+1] << " = " << extractionPercentages[(extractionPercentages.size()-1)];
1574 
1575  std::cout << ")" << std::endl;
1576  }
1577  if(extractionHM.size())
1578  {
1579  std::cout << "* Extraction of hadron-matched jets with types: (";
1580  for(Int_t i=0; i<static_cast<Int_t>(extractionHM.size()-1); i++)
1581  std::cout << extractionHM[i] << ", ";
1582  std::cout << extractionHM[extractionHM.size()-1];
1583  std::cout << ")" << std::endl;
1584  }
1585  if(extractionPM.size())
1586  {
1587  std::cout << "* Extraction of parton-matched jets with types: (";
1588  for(Int_t i=0; i<static_cast<Int_t>(extractionPM.size()-1); i++)
1589  std::cout << extractionPM[i] << ", ";
1590  std::cout << extractionPM[extractionPM.size()-1];
1591  std::cout << ")" << std::endl;
1592  }
1593  std::cout << std::endl;
1594 
1595  std::cout << "### Extracted data groups ###" << std::endl;
1596  std::cout << "* Basic event properties (ID, vertex, centrality, bgrd. densities, ...)" << std::endl;
1597  if(fSaveConstituents)
1598  std::cout << "* Jet constituents, basic properties (pt, eta, phi, charge, ...)" << std::endl;
1600  std::cout << "* Jet constituents, IPs" << std::endl;
1602  std::cout << "* Jet constituents, PID" << std::endl;
1603  if(fSaveCaloClusters)
1604  std::cout << "* Jet calorimeter clusters" << std::endl;
1605  if(fSaveMCInformation)
1606  std::cout << "* MC information (origin, matched jets, ...)" << std::endl;
1608  std::cout << "* Secondary vertices" << std::endl;
1609  if(fSaveJetShapes)
1610  std::cout << "* Jet shapes (jet mass, LeSub, pTD, ...)" << std::endl;
1611  if(fSaveJetSplittings)
1612  std::cout << "* Jet splittings (kT, theta, E from iterative CA reclustering)" << std::endl;
1613  if(fSaveTriggerTracks)
1614  std::cout << "* Trigger tracks" << std::endl;
1615  std::cout << std::endl;
1616  std::cout << "### Further settings ###" << std::endl;
1617  if(fIsEmbeddedEvent)
1618  std::cout << Form("* EMCal embedding framework will be used (at least on container has IsEmbedded() true)") << std::endl;
1619  if(fMatchedJetsArrayName != "")
1620  std::cout << Form("* Jet matching active, array=%s, rho=%s, rho_mass=%s", fMatchedJetsArrayName.Data(), fMatchedJetsRhoName.Data(), fMatchedJetsRhoMassName.Data()) << std::endl;
1621  if(fMCParticleArray)
1622  std::cout << Form("* Particle level information available (for jet origin calculation, particle code): %s", fMCParticleArrayName.Data()) << std::endl;
1623  if(extractionHM.size())
1624  std::cout << Form("* Hadronic b/c-jet matching radius=%3.3f", fHadronMatchingRadius) << std::endl;
1625  if(fSaveMCInformation)
1626  std::cout << Form("* True particle label range for pt fraction=%d-%d", fTruthMinLabel, fTruthMaxLabel) << std::endl;
1628  std::cout << Form("* Particle PID codes will be PDG codes") << std::endl;
1630  std::cout << Form("* Particle PID codes will follow AliAODPid convention") << std::endl;
1632  std::cout << Form("* AliEmcalJet flavour tag is set from HF-jet tagging") << std::endl;
1633 
1634  std::cout << "#########################################\n\n";
1635 }
1636 
1637 
1638 //########################################################################
1639 // HELPERS
1640 //########################################################################
1641 
1642 //________________________________________________________________________
1643 inline Bool_t AliAnalysisTaskJetExtractor::IsTrackInCone(const AliVParticle* track, Double_t eta, Double_t phi, Double_t radius)
1644 {
1645  Double_t deltaR = GetDistance(track->Eta(), eta, track->Phi(), phi);
1646  if(deltaR <= radius)
1647  return kTRUE;
1648 
1649  return kFALSE;
1650 }
1651 
1652 //________________________________________________________________________
1654 {
1655  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1656  if(!tmpHist)
1657  {
1658  AliError(Form("Cannot find histogram <%s> ",key)) ;
1659  return;
1660  }
1661 
1662  tmpHist->Fill(x);
1663 }
1664 
1665 //________________________________________________________________________
1667 {
1668  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1669  if(!tmpHist)
1670  {
1671  AliError(Form("Cannot find histogram <%s> ",key));
1672  return;
1673  }
1674 
1675  if (tmpHist->IsA()->GetBaseClass("TH1"))
1676  static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
1677  else if (tmpHist->IsA()->GetBaseClass("TH2"))
1678  static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
1679 }
1680 
1681 //________________________________________________________________________
1683 {
1684  TH2* tmpHist = static_cast<TH2*>(fOutput->FindObject(key));
1685  if(!tmpHist)
1686  {
1687  AliError(Form("Cannot find histogram <%s> ",key));
1688  return;
1689  }
1690 
1691  tmpHist->Fill(x,y,add);
1692 }
1693 
1694 //________________________________________________________________________
1696 {
1697  TH3* tmpHist = static_cast<TH3*>(fOutput->FindObject(key));
1698  if(!tmpHist)
1699  {
1700  AliError(Form("Cannot find histogram <%s> ",key));
1701  return;
1702  }
1703 
1704  if(add)
1705  tmpHist->Fill(x,y,z,add);
1706  else
1707  tmpHist->Fill(x,y,z);
1708 }
1709 
1710 
1711 //________________________________________________________________________
1712 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)
1713 {
1714  T* tmpHist = new T(name, title, xBins, xMin, xMax);
1715 
1716  tmpHist->GetXaxis()->SetTitle(xTitle);
1717  tmpHist->GetYaxis()->SetTitle(yTitle);
1718  tmpHist->SetOption(options);
1719  tmpHist->SetMarkerStyle(kFullCircle);
1720  tmpHist->Sumw2();
1721 
1722  fOutput->Add(tmpHist);
1723 
1724  return tmpHist;
1725 }
1726 
1727 //________________________________________________________________________
1728 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)
1729 {
1730  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax);
1731  tmpHist->GetXaxis()->SetTitle(xTitle);
1732  tmpHist->GetYaxis()->SetTitle(yTitle);
1733  tmpHist->GetZaxis()->SetTitle(zTitle);
1734  tmpHist->SetOption(options);
1735  tmpHist->SetMarkerStyle(kFullCircle);
1736  tmpHist->Sumw2();
1737 
1738  fOutput->Add(tmpHist);
1739 
1740  return tmpHist;
1741 }
1742 
1743 //________________________________________________________________________
1744 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)
1745 {
1746  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax, zBins, zMin, zMax);
1747  tmpHist->GetXaxis()->SetTitle(xTitle);
1748  tmpHist->GetYaxis()->SetTitle(yTitle);
1749  tmpHist->GetZaxis()->SetTitle(zTitle);
1750  tmpHist->SetOption(options);
1751  tmpHist->SetMarkerStyle(kFullCircle);
1752  tmpHist->Sumw2();
1753 
1754  fOutput->Add(tmpHist);
1755 
1756  return tmpHist;
1757 }
1758 
1759 //________________________________________________________________________
1761 {
1762  // Called once at the end of the analysis.
1763 }
1764 
1765 // ### ADDTASK MACRO
1766 //________________________________________________________________________
1767 AliAnalysisTaskJetExtractor* AliAnalysisTaskJetExtractor::AddTaskJetExtractor(TString trackArray, TString clusterArray, TString jetArray, TString rhoObject, Double_t jetRadius, AliRDHFJetsCutsVertex* vertexerCuts, const char* taskNameSuffix)
1768 {
1769  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1770  Double_t minJetEta = 0.5;
1771  Double_t minJetPt = 0.15;
1772  Double_t minTrackPt = 0.15;
1773  Double_t minJetAreaPerc = 0.557;
1774  TString suffix = "";
1775  if(taskNameSuffix != 0)
1776  suffix = taskNameSuffix;
1777 
1778  // ###### Task name
1779  TString name("AliAnalysisTaskJetExtractor");
1780  if (jetArray != "") {
1781  name += "_";
1782  name += jetArray;
1783  }
1784  if (rhoObject != "") {
1785  name += "_";
1786  name += rhoObject;
1787  }
1788  if (suffix != "") {
1789  name += "_";
1790  name += suffix;
1791  }
1792 
1793  // ###### Setup task with default settings
1795  myTask->SetNeedEmcalGeom(kFALSE);
1796  myTask->SetVzRange(-10.,10.);
1797 
1798  // Particle container and track pt cut
1799  AliParticleContainer* trackCont = 0;
1800  if(trackArray == "mcparticles")
1801  trackCont = myTask->AddMCParticleContainer(trackArray);
1802  else if(trackArray =="mctracks")
1803  trackCont = myTask->AddParticleContainer(trackArray);
1804  else
1805  trackCont = myTask->AddTrackContainer(trackArray);
1806  trackCont->SetParticlePtCut(minTrackPt);
1807 
1808  // Particle container and track pt cut
1809  AliClusterContainer* clusterCont = 0;
1810  if(clusterArray != "")
1811  clusterCont = myTask->AddClusterContainer(clusterArray);
1812 
1813  // Secondary vertex cuts (default settings from PWGHF)
1814  // (can be overwritten by using myTask->SetVertexerCuts(cuts) from outside macro)
1815  AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts", "default");
1816  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1817  esdTrackCuts->SetMinNClustersTPC(90);
1818  esdTrackCuts->SetMaxChi2PerClusterTPC(4);
1819  esdTrackCuts->SetRequireTPCRefit(kTRUE);
1820  esdTrackCuts->SetRequireITSRefit(kTRUE);
1821  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
1822  esdTrackCuts->SetMinDCAToVertexXY(0.);
1823  esdTrackCuts->SetEtaRange(-0.8, 0.8);
1824  esdTrackCuts->SetPtRange(1.0, 1.e10);
1825  vertexerCuts->AddTrackCuts(esdTrackCuts);
1826  vertexerCuts->SetNprongs(3);
1827  vertexerCuts->SetMinPtHardestTrack(1.0);//default 0.3
1828  vertexerCuts->SetSecVtxWithKF(kFALSE);//default with StrLinMinDist
1829  vertexerCuts->SetImpParCut(0.);//default 0
1830  vertexerCuts->SetDistPrimSec(0.);//default 0
1831  vertexerCuts->SetCospCut(-1);//default -1
1832  myTask->SetVertexerCuts(vertexerCuts);
1833 
1834  // Jet container
1835  AliJetContainer *jetCont = myTask->AddJetContainer(jetArray,6,jetRadius);
1836  if (jetCont) {
1837  jetCont->SetRhoName(rhoObject);
1838  jetCont->SetPercAreaCut(minJetAreaPerc);
1839  jetCont->SetJetPtCut(minJetPt);
1840  jetCont->SetLeadingHadronType(0);
1841  jetCont->SetPtBiasJetTrack(minTrackPt);
1842  jetCont->SetJetEtaLimits(-minJetEta, +minJetEta);
1843  jetCont->ConnectParticleContainer(trackCont);
1844  if(clusterCont)
1845  jetCont->ConnectClusterContainer(clusterCont);
1846  jetCont->SetMaxTrackPt(1000);
1847  }
1848 
1849  mgr->AddTask(myTask);
1850 
1851  // ###### Connect inputs/outputs
1852  mgr->ConnectInput (myTask, 0, mgr->GetCommonInputContainer() );
1853  mgr->ConnectOutput (myTask, 1, mgr->CreateContainer(Form("%s_histos", name.Data()), AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, Form("%s:ChargedJetsHadronCF", mgr->GetCommonFileName())) );
1854  mgr->ConnectOutput (myTask, 2, mgr->CreateContainer(Form("%s_tree", name.Data()), TTree::Class(), AliAnalysisManager::kOutputContainer, mgr->GetCommonFileName()) );
1855 
1856  return myTask;
1857 }
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_Jet_MC_TruePtFraction_mcparticles
! array buffer
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
void FillBuffer_MonteCarlo(Int_t motherParton, Int_t motherHadron, Int_t partonInitialCollision, Float_t matchJetDistance, Float_t matchedJetPt, Float_t matchedJetMass, Float_t truePtFraction, Float_t truePtFraction_mcparticles, Float_t ptHard, Float_t eventWeight, Float_t impactParameter)
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_Jet_MC_MatchedJet_Pt
! array buffer
Float_t fBuffer_Shape_LeSub_DerivCorr
! array buffer
Double_t GetJetEtaMin() const
const char * title
Definition: MakeQAPdf.C:27
void GetJetType(AliEmcalJet *jet, Int_t &typeHM, Int_t &typePM, Int_t &typeIC)
Float_t GetPartonEta6() const
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
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
Float_t fBuffer_Jet_MC_MatchedJet_Mass
! array buffer
Double_t Phi() const
Definition: AliEmcalJet.h:117
std::vector< Float_t > GetExtractionPercentages()
Double_t mass
void SetPtBiasJetTrack(Float_t b)
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)
Double_t GetJetEtaMax() const
void FillBuffer_JetShapes(AliEmcalJet *jet, Double_t leSub_noCorr, Double_t angularity, Double_t momentumDispersion, Double_t trackPtMean, Double_t trackPtMedian)
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
TString fMatchedJetsRhoMassName
Name for matched jets rho_mass object.
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) ...
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)
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_Cluster_Time
! array buffer
Double_t fEventPercentage
percentage (0, 1] which will be extracted
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
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
void InitializeTree(Bool_t saveCaloClusters, Bool_t saveMCInformation, Bool_t saveConstituents, Bool_t saveConstituentsIP, Bool_t saveConstituentPID, Bool_t saveJetShapes, Bool_t saveSplittings, Bool_t saveSecondaryVertices, Bool_t saveTriggerTracks)
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)
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
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)
void GetMatchedJetObservables(AliEmcalJet *jet, Double_t &matchedJetPt, Double_t &matchedJetMass, Double_t &matchedJetDistance)
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
TString fMatchedJetsArrayName
Array name for matched jets.
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
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)
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 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
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
Float_t * fBuffer_Cluster_Pt
! array buffer
Float_t fBuffer_JetEta
! array buffer
Int_t * fBuffer_Cluster_Label
! array buffer
TString fMatchedJetsRhoName
Name for matched jets rho object.
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.
const char Option_t
Definition: External.C:48
Bool_t fSaveTriggerTracks
save event trigger track
Float_t fBuffer_Shape_TrackPtMean
! array buffer
void UserCreateOutputObjects()
Main initialization function on the worker.
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)
bool Bool_t
Definition: External.C:53
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:375
Double_t yMin
Int_t GetNAcceptedParticles() const
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)
Float_t fBuffer_Jet_MC_MatchedJet_Distance
! array buffer
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
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
Float_t * fBuffer_Track_Pt
! array buffer
Float_t fBuffer_Event_MagneticField
! array buffer