AliPhysics  b0c77bb (b0c77bb)
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 "AliYAMLConfiguration.h"
58 
59 #include "AliGenHepMCEventHeader.h"
60 #include "AliGenHijingEventHeader.h"
61 #include "AliHFJetsTaggingVertex.h"
62 #include "AliRDHFJetsCutsVertex.h"
63 
65 
67 ClassImp(AliEmcalJetTree)
69 
73 
74 //________________________________________________________________________
75 AliEmcalJetTree::AliEmcalJetTree() : TNamed("CustomTree", "CustomTree"), fJetTree(0), fInitialized(0), fSaveEventProperties(0), fSaveConstituents(0), fSaveConstituentsIP(0), fSaveConstituentPID(0), fSaveJetShapes(0), fSaveMCInformation(0), fSaveSecondaryVertices(0), fSaveTriggerTracks(0), fSaveCaloClusters(0), fExtractionPercentages(), fExtractionPercentagePtBins(), fExtractionJetTypes_HM(), fExtractionJetTypes_PM()
76 {
77  // For these arrays, we need to reserve memory
78  fBuffer_Const_Pt = new Float_t[kMaxNumConstituents];
79  fBuffer_Const_Eta = new Float_t[kMaxNumConstituents];
80  fBuffer_Const_Phi = new Float_t[kMaxNumConstituents];
81  fBuffer_Const_Charge = new Float_t[kMaxNumConstituents];
82  fBuffer_Const_Label = new Int_t [kMaxNumConstituents];
83  fBuffer_Const_ProdVtx_X = new Float_t[kMaxNumConstituents];
84  fBuffer_Const_ProdVtx_Y = new Float_t[kMaxNumConstituents];
85  fBuffer_Const_ProdVtx_Z = new Float_t[kMaxNumConstituents];
86 
87  fBuffer_Cluster_Pt = new Float_t[kMaxNumConstituents];
88  fBuffer_Cluster_E = new Float_t[kMaxNumConstituents];
89  fBuffer_Cluster_Eta = new Float_t[kMaxNumConstituents];
90  fBuffer_Cluster_Phi = new Float_t[kMaxNumConstituents];
91  fBuffer_Cluster_M02 = new Float_t[kMaxNumConstituents];
92  fBuffer_Cluster_Time = new Float_t[kMaxNumConstituents];
93  fBuffer_Cluster_Label = new Int_t[kMaxNumConstituents];
94 }
95 
96 //________________________________________________________________________
97 AliEmcalJetTree::AliEmcalJetTree(const char* name) : TNamed(name, name), fJetTree(0), fInitialized(0), fSaveEventProperties(0), fSaveConstituents(0), fSaveConstituentsIP(0), fSaveConstituentPID(0), fSaveJetShapes(0), fSaveMCInformation(0), fSaveSecondaryVertices(0), fSaveTriggerTracks(0), fSaveCaloClusters(0), fExtractionPercentages(), fExtractionPercentagePtBins(), fExtractionJetTypes_HM(), fExtractionJetTypes_PM()
98 {
99  // For these arrays, we need to reserve memory
108 
116 }
117 
118 
119 //________________________________________________________________________
121  Int_t motherParton, Int_t motherHadron, Int_t partonInitialCollision, Float_t matchDistance, Float_t matchedPt, Float_t matchedMass, Float_t truePtFraction, Float_t ptHard, Float_t eventWeight, Float_t impactParameter,
122  Float_t* trackPID_ITS, Float_t* trackPID_TPC, Float_t* trackPID_TOF, Float_t* trackPID_TRD, Short_t* trackPID_Reco, Int_t* trackPID_Truth,
123  Int_t numTriggerTracks, Float_t* triggerTrackPt, Float_t* triggerTrackDeltaEta, Float_t* triggerTrackDeltaPhi,
124  Float_t* trackIP_d0, Float_t* trackIP_z0, Float_t* trackIP_d0cov, Float_t* trackIP_z0cov,
125  Int_t numSecVertices, Float_t* secVtx_X, Float_t* secVtx_Y, Float_t* secVtx_Z, Float_t* secVtx_Mass, Float_t* secVtx_Lxy, Float_t* secVtx_SigmaLxy, Float_t* secVtx_Chi2, Float_t* secVtx_Dispersion)
126 {
127  // ############
128  // Approach: Fill buffer and then add to tree
129  //
130  if(!fInitialized)
131  AliFatal("Tree is not initialized.");
132 
133  fBuffer_JetPt = jet->Pt() - fJetContainer->GetRhoVal()*jet->Area();
134 
135  // Check if jet type is contained in extraction list
136  if( (fExtractionJetTypes_PM.size() || fExtractionJetTypes_HM.size()) &&
137  (std::find(fExtractionJetTypes_PM.begin(), fExtractionJetTypes_PM.end(), motherParton) == fExtractionJetTypes_PM.end()) &&
138  (std::find(fExtractionJetTypes_HM.begin(), fExtractionJetTypes_HM.end(), motherHadron) == fExtractionJetTypes_HM.end()) )
139  return kFALSE;
140 
141  // Check acceptance percentage for the given jet and discard statistically on demand
142  Bool_t inPtRange = kFALSE;
143  for(size_t i = 0; i<fExtractionPercentages.size(); i++)
144  {
146  {
147  inPtRange = kTRUE;
148  if(fRandomGenerator->Rndm() >= fExtractionPercentages[i])
149  return kFALSE;
150  break;
151  }
152  }
153  if(!inPtRange && fExtractionPercentagePtBins.size()) // only discard jet if fExtractionPercentagePtBins was given
154  return kFALSE;
155 
156  // Set basic properties
157  fBuffer_JetEta = jet->Eta();
158  fBuffer_JetPhi = jet->Phi();
159  fBuffer_JetArea = jet->Area();
160 
161  // Set event properties
163  {
166  fBuffer_Event_Vertex_X = vertexX;
167  fBuffer_Event_Vertex_Y = vertexY;
168  fBuffer_Event_Vertex_Z = vertexZ;
170  fBuffer_Event_ID = eventID;
171  fBuffer_Event_MagneticField = magField;
173  {
174  fBuffer_Event_PtHard = ptHard;
175  fBuffer_Event_Weight = eventWeight;
176  fBuffer_Event_ImpactParameter = impactParameter;
177  }
178  }
179 
180  // Extract basic constituent properties directly from AliEmcalJet object
183  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
184  {
185  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(i, fTrackContainer->GetArray()));
186  if(!particle) continue;
187 
189  {
190  fBuffer_Const_Pt[fBuffer_NumConstituents] = particle->Pt();
191  fBuffer_Const_Eta[fBuffer_NumConstituents] = particle->Eta();
192  fBuffer_Const_Phi[fBuffer_NumConstituents] = particle->Phi();
193  fBuffer_Const_Charge[fBuffer_NumConstituents] = particle->Charge();
194  fBuffer_Const_Label[fBuffer_NumConstituents] = particle->GetLabel();
195  }
197  {
201  }
203  }
204 
205  // Extract basic CALOCLUSTER properties directly from AliEmcalJet object
208  for(Int_t i = 0; i < jet->GetNumberOfClusters(); i++)
209  {
210  AliVCluster* cluster = static_cast<AliVCluster*>(jet->ClusterAt(i, fClusterContainer->GetArray()));
211  if(!cluster) continue;
212 
213  // #### Retrieve cluster pT
214  Double_t primVtx[3] = {vertexX, vertexY, vertexZ};
215  TLorentzVector clusterMomentum;
216  cluster->GetMomentum(clusterMomentum, primVtx);
217  // ####
218  fBuffer_Cluster_Pt[fBuffer_NumClusters] = clusterMomentum.Perp();
219  fBuffer_Cluster_E[fBuffer_NumClusters] = cluster->E();
220  fBuffer_Cluster_Eta[fBuffer_NumClusters] = clusterMomentum.Eta();
221  fBuffer_Cluster_Phi[fBuffer_NumClusters] = clusterMomentum.Phi();
222  fBuffer_Cluster_M02[fBuffer_NumClusters] = cluster->GetM02();
223  fBuffer_Cluster_Time[fBuffer_NumClusters] = cluster->GetTOF();
224  fBuffer_Cluster_Label[fBuffer_NumClusters] = cluster->GetLabel();
225 
227  }
228 
229 
230  // Set constituent arrays for impact parameters
232  {
233  fJetTree->SetBranchAddress("Jet_Const_CovIPd", trackIP_d0cov);
234  fJetTree->SetBranchAddress("Jet_Const_CovIPz", trackIP_z0cov);
235  fJetTree->SetBranchAddress("Jet_Const_IPd", trackIP_d0);
236  fJetTree->SetBranchAddress("Jet_Const_IPz", trackIP_z0);
237  }
238  // Set constituent arrays for PID parameters
240  {
241  fJetTree->SetBranchAddress("Jet_Const_PID_ITS", trackPID_ITS);
242  fJetTree->SetBranchAddress("Jet_Const_PID_TPC", trackPID_TPC);
243  fJetTree->SetBranchAddress("Jet_Const_PID_TOF", trackPID_TOF);
244  fJetTree->SetBranchAddress("Jet_Const_PID_TRD", trackPID_TRD);
245  fJetTree->SetBranchAddress("Jet_Const_PID_Reconstructed", trackPID_Reco);
247  fJetTree->SetBranchAddress("Jet_Const_PID_Truth", trackPID_Truth);
248  }
249 
250  // Set secondary vertices arrays
252  {
253  fBuffer_NumSecVertices = numSecVertices;
254  fJetTree->SetBranchAddress("Jet_SecVtx_X", secVtx_X);
255  fJetTree->SetBranchAddress("Jet_SecVtx_Y", secVtx_Y);
256  fJetTree->SetBranchAddress("Jet_SecVtx_Z", secVtx_Z);
257  fJetTree->SetBranchAddress("Jet_SecVtx_Mass", secVtx_Mass);
258  fJetTree->SetBranchAddress("Jet_SecVtx_Lxy", secVtx_Lxy);
259  fJetTree->SetBranchAddress("Jet_SecVtx_SigmaLxy", secVtx_SigmaLxy);
260  fJetTree->SetBranchAddress("Jet_SecVtx_Chi2", secVtx_Chi2);
261  fJetTree->SetBranchAddress("Jet_SecVtx_Dispersion", secVtx_Dispersion);
262  }
263 
264  // Set jet shape observables (some others are set in SetJetShapesInBuffer)
265  if(fSaveJetShapes)
266  {
267  fBuffer_Shape_Mass_NoCorr = jet->M();
280  }
281 
282  // Set Monte Carlo information
284  {
285  fBuffer_Jet_MC_MotherParton = motherParton;
286  fBuffer_Jet_MC_MotherHadron = motherHadron;
287  fBuffer_Jet_MC_MotherIC = partonInitialCollision;
288  fBuffer_Jet_MC_MatchDistance = matchDistance;
289  fBuffer_Jet_MC_MatchedPt = matchedPt;
290  fBuffer_Jet_MC_MatchedMass = matchedMass;
291  fBuffer_Jet_MC_TruePtFraction = truePtFraction;
292  }
293 
295  {
296  fBuffer_NumTriggerTracks = numTriggerTracks;
297  fJetTree->SetBranchAddress("Jet_TriggerTrack_Pt", triggerTrackPt);
298  fJetTree->SetBranchAddress("Jet_TriggerTrack_dEta", triggerTrackDeltaEta);
299  fJetTree->SetBranchAddress("Jet_TriggerTrack_dPhi", triggerTrackDeltaPhi);
300  }
301 
302  // Add buffered jet to tree
303  fJetTree->Fill();
304 
305  return kTRUE;
306 }
307 
308 //________________________________________________________________________
309 void AliEmcalJetTree::SetJetShapesInBuffer(Double_t leSub_noCorr, Double_t radialMoment, Double_t momentumDispersion, Double_t constPtMean, Double_t constPtMedian)
310 {
311  fBuffer_Shape_LeSub_NoCorr = leSub_noCorr;
312  fBuffer_Shape_RadialMoment = radialMoment;
313  fBuffer_Shape_MomentumDispersion = momentumDispersion;
314  fBuffer_Shape_ConstPtMean = constPtMean;
315  fBuffer_Shape_ConstPtMedian = constPtMedian;
316 }
317 
318 //________________________________________________________________________
320 {
321  // Create the tree with active branches
322 
323  void* dummy = 0; // for some branches, do not need a buffer pointer yet
324 
325 
326  if(fInitialized)
327  AliFatal("Tree is already initialized.");
328 
329  if(!fJetContainer)
330  AliFatal("fJetContainer not set in tree");
331 
332  // ### Prepare the jet tree
333  fJetTree = new TTree(Form("JetTree_%s", GetName()), "");
334 
335  fJetTree->Branch("Jet_Pt",&fBuffer_JetPt,"Jet_Pt/F");
336  fJetTree->Branch("Jet_Phi",&fBuffer_JetPhi,"Jet_Phi/F");
337  fJetTree->Branch("Jet_Eta",&fBuffer_JetEta,"Jet_Eta/F");
338  fJetTree->Branch("Jet_Area",&fBuffer_JetArea,"Jet_Area/F");
339  fJetTree->Branch("Jet_NumConstituents",&fBuffer_NumConstituents,"Jet_NumConstituents/I");
341  fJetTree->Branch("Jet_NumClusters",&fBuffer_NumClusters,"Jet_NumClusters/I");
342 
344  {
345  fJetTree->Branch("Event_BackgroundDensity",&fBuffer_Event_BackgroundDensity,"Event_BackgroundDensity/F");
346  fJetTree->Branch("Event_BackgroundDensityMass",&fBuffer_Event_BackgroundDensityMass,"Event_BackgroundDensityMass/F");
347  fJetTree->Branch("Event_Vertex_X",&fBuffer_Event_Vertex_X,"Event_Vertex_X/F");
348  fJetTree->Branch("Event_Vertex_Y",&fBuffer_Event_Vertex_Y,"Event_Vertex_Y/F");
349  fJetTree->Branch("Event_Vertex_Z",&fBuffer_Event_Vertex_Z,"Event_Vertex_Z/F");
350  fJetTree->Branch("Event_Centrality",&fBuffer_Event_Centrality,"Event_Centrality/F");
351  fJetTree->Branch("Event_ID",&fBuffer_Event_ID,"Event_ID/L");
352  fJetTree->Branch("Event_MagneticField",&fBuffer_Event_MagneticField,"Event_MagneticField/F");
353 
355  {
356  fJetTree->Branch("Event_PtHard",&fBuffer_Event_PtHard,"Event_PtHard/F");
357  fJetTree->Branch("Event_Weight",&fBuffer_Event_Weight,"Event_Weight/F");
358  fJetTree->Branch("Event_ImpactParameter",&fBuffer_Event_ImpactParameter,"Event_ImpactParameter/F");
359  }
360  }
361 
363  {
364  fJetTree->Branch("Jet_Const_Pt",fBuffer_Const_Pt,"Jet_Const_Pt[Jet_NumConstituents]/F");
365  fJetTree->Branch("Jet_Const_Phi",fBuffer_Const_Phi,"Jet_Const_Phi[Jet_NumConstituents]/F");
366  fJetTree->Branch("Jet_Const_Eta",fBuffer_Const_Eta,"Jet_Const_Eta[Jet_NumConstituents]/F");
367  fJetTree->Branch("Jet_Const_Charge",fBuffer_Const_Charge,"Jet_Const_Charge[Jet_NumConstituents]/F");
369  fJetTree->Branch("Jet_Const_Label",fBuffer_Const_Label,"Jet_Const_Label[Jet_NumConstituents]/I");
370  }
371 
373  {
374  fJetTree->Branch("Jet_Cluster_Pt",fBuffer_Cluster_Pt,"Jet_Cluster_Pt[Jet_NumClusters]/F");
375  fJetTree->Branch("Jet_Cluster_E",fBuffer_Cluster_E,"Jet_Cluster_Pt[Jet_NumClusters]/F");
376  fJetTree->Branch("Jet_Cluster_Phi",fBuffer_Cluster_Phi,"Jet_Cluster_Phi[Jet_NumClusters]/F");
377  fJetTree->Branch("Jet_Cluster_Eta",fBuffer_Cluster_Eta,"Jet_Cluster_Eta[Jet_NumClusters]/F");
378  fJetTree->Branch("Jet_Cluster_M02",fBuffer_Cluster_M02,"Jet_Cluster_M02[Jet_NumClusters]/F");
379  fJetTree->Branch("Jet_Cluster_Time",fBuffer_Cluster_Time,"Jet_Cluster_Time[Jet_NumClusters]/F");
381  fJetTree->Branch("Jet_Cluster_Label",fBuffer_Cluster_Label,"Jet_Cluster_Label[Jet_NumClusters]/I");
382  }
383 
385  {
386  fJetTree->Branch("Jet_Const_IPd",&dummy,"Jet_Const_IPd[Jet_NumConstituents]/F");
387  fJetTree->Branch("Jet_Const_IPz",&dummy,"Jet_Const_IPz[Jet_NumConstituents]/F");
388  fJetTree->Branch("Jet_Const_CovIPd",&dummy,"Jet_Const_CovIPd[Jet_NumConstituents]/F");
389  fJetTree->Branch("Jet_Const_CovIPz",&dummy,"Jet_Const_CovIPz[Jet_NumConstituents]/F");
390 
391  fJetTree->Branch("Jet_Const_ProdVtx_X",fBuffer_Const_ProdVtx_X,"Jet_Const_ProdVtx_X[Jet_NumConstituents]/F");
392  fJetTree->Branch("Jet_Const_ProdVtx_Y",fBuffer_Const_ProdVtx_Y,"Jet_Const_ProdVtx_Y[Jet_NumConstituents]/F");
393  fJetTree->Branch("Jet_Const_ProdVtx_Z",fBuffer_Const_ProdVtx_Z,"Jet_Const_ProdVtx_Z[Jet_NumConstituents]/F");
394  }
395 
397  {
398  fJetTree->Branch("Jet_Const_PID_ITS",&dummy,"Jet_Const_PID_ITS[Jet_NumConstituents]/F");
399  fJetTree->Branch("Jet_Const_PID_TPC",&dummy,"Jet_Const_PID_TPC[Jet_NumConstituents]/F");
400  fJetTree->Branch("Jet_Const_PID_TOF",&dummy,"Jet_Const_PID_TOF[Jet_NumConstituents]/F");
401  fJetTree->Branch("Jet_Const_PID_TRD",&dummy,"Jet_Const_PID_TRD[Jet_NumConstituents]/F");
402 
403  fJetTree->Branch("Jet_Const_PID_Reconstructed",&dummy,"Jet_Const_PID_Reconstructed[Jet_NumConstituents]/S");
405  fJetTree->Branch("Jet_Const_PID_Truth",&dummy,"Jet_Const_PID_Truth[Jet_NumConstituents]/I");
406  }
407 
408  if(fSaveJetShapes)
409  {
410  fJetTree->Branch("Jet_Shape_Mass_NoCorr",&fBuffer_Shape_Mass_NoCorr,"Jet_Shape_Mass_NoCorr/F");
411  fJetTree->Branch("Jet_Shape_Mass_DerivCorr_1",&fBuffer_Shape_Mass_DerivCorr_1,"Jet_Shape_Mass_DerivCorr_1/F");
412  fJetTree->Branch("Jet_Shape_Mass_DerivCorr_2",&fBuffer_Shape_Mass_DerivCorr_2,"Jet_Shape_Mass_DerivCorr_2/F");
413  fJetTree->Branch("Jet_Shape_pTD_DerivCorr_1",&fBuffer_Shape_pTD_DerivCorr_1,"Jet_Shape_pTD_DerivCorr_1/F");
414  fJetTree->Branch("Jet_Shape_pTD_DerivCorr_2",&fBuffer_Shape_pTD_DerivCorr_2,"Jet_Shape_pTD_DerivCorr_2/F");
415  fJetTree->Branch("Jet_Shape_LeSub_NoCorr",&fBuffer_Shape_LeSub_NoCorr,"Jet_Shape_LeSub_NoCorr/F");
416  fJetTree->Branch("Jet_Shape_LeSub_DerivCorr",&fBuffer_Shape_LeSub_DerivCorr,"Jet_Shape_LeSub_DerivCorr/F");
417  fJetTree->Branch("Jet_Shape_Angularity_DerivCorr_1",&fBuffer_Shape_Angularity_DerivCorr_1,"Jet_Shape_Angularity_DerivCorr_1/F");
418  fJetTree->Branch("Jet_Shape_Angularity_DerivCorr_2",&fBuffer_Shape_Angularity_DerivCorr_2,"Jet_Shape_Angularity_DerivCorr_2/F");
419  fJetTree->Branch("Jet_Shape_Circularity_DerivCorr_1",&fBuffer_Shape_Circularity_DerivCorr_1,"Jet_Shape_Circularity_DerivCorr_1/F");
420  fJetTree->Branch("Jet_Shape_Circularity_DerivCorr_2",&fBuffer_Shape_Circularity_DerivCorr_2,"Jet_Shape_Circularity_DerivCorr_2/F");
421  fJetTree->Branch("Jet_Shape_Sigma2_DerivCorr_1",&fBuffer_Shape_Sigma2_DerivCorr_1,"Jet_Shape_Sigma2_DerivCorr_1/F");
422  fJetTree->Branch("Jet_Shape_Sigma2_DerivCorr_2",&fBuffer_Shape_Sigma2_DerivCorr_2,"Jet_Shape_Sigma2_DerivCorr_2/F");
423  fJetTree->Branch("Jet_Shape_NumConst_DerivCorr",&fBuffer_Shape_NumConst_DerivCorr,"Jet_Shape_NumConst_DerivCorr/F");
424  fJetTree->Branch("Jet_Shape_RadialMoment",&fBuffer_Shape_RadialMoment,"Jet_Shape_RadialMoment/F");
425  fJetTree->Branch("Jet_Shape_MomentumDispersion",&fBuffer_Shape_MomentumDispersion,"Jet_Shape_MomentumDispersion/F");
426  fJetTree->Branch("Jet_Shape_ConstPtMean",&fBuffer_Shape_ConstPtMean,"Jet_Shape_ConstPtMean/F");
427  fJetTree->Branch("Jet_Shape_ConstPtMedian",&fBuffer_Shape_ConstPtMedian,"Jet_Shape_ConstPtMedian/F");
428 
429  }
430 
432  {
433  fJetTree->Branch("Jet_MC_MotherParton",&fBuffer_Jet_MC_MotherParton,"Jet_MC_MotherParton/I");
434  fJetTree->Branch("Jet_MC_MotherHadron",&fBuffer_Jet_MC_MotherHadron,"Jet_MC_MotherHadron/I");
435  fJetTree->Branch("Jet_MC_MotherIC",&fBuffer_Jet_MC_MotherIC,"Jet_MC_MotherIC/I");
436  fJetTree->Branch("Jet_MC_MatchDistance",&fBuffer_Jet_MC_MatchDistance,"Jet_MC_MatchDistance/F");
437  fJetTree->Branch("Jet_MC_MatchedPt",&fBuffer_Jet_MC_MatchedPt,"Jet_MC_MatchedPt/F");
438  fJetTree->Branch("Jet_MC_MatchedMass",&fBuffer_Jet_MC_MatchedMass,"Jet_MC_MatchedMass/F");
439  fJetTree->Branch("Jet_MC_TruePtFraction",&fBuffer_Jet_MC_TruePtFraction,"Jet_MC_TruePtFraction/F");
440  }
441 
443  {
444  fJetTree->Branch("Jet_NumSecVertices",&fBuffer_NumSecVertices,"Jet_NumSecVertices/I");
445 
446  fJetTree->Branch("Jet_SecVtx_X",&dummy,"Jet_SecVtx_X[Jet_NumSecVertices]/F");
447  fJetTree->Branch("Jet_SecVtx_Y",&dummy,"Jet_SecVtx_Y[Jet_NumSecVertices]/F");
448  fJetTree->Branch("Jet_SecVtx_Z",&dummy,"Jet_SecVtx_Z[Jet_NumSecVertices]/F");
449  fJetTree->Branch("Jet_SecVtx_Mass",&dummy,"Jet_SecVtx_Mass[Jet_NumSecVertices]/F");
450  fJetTree->Branch("Jet_SecVtx_Lxy",&dummy,"Jet_SecVtx_Lxy[Jet_NumSecVertices]/F");
451  fJetTree->Branch("Jet_SecVtx_SigmaLxy",&dummy,"Jet_SecVtx_SigmaLxy[Jet_NumSecVertices]/F");
452  fJetTree->Branch("Jet_SecVtx_Chi2",&dummy,"Jet_SecVtx_Chi2[Jet_NumSecVertices]/F");
453  fJetTree->Branch("Jet_SecVtx_Dispersion",&dummy,"Jet_SecVtx_Dispersion[Jet_NumSecVertices]/F");
454  }
455 
456  // Trigger track requirement active -> Save trigger track
458  {
459  fJetTree->Branch("Jet_NumTriggerTracks",&fBuffer_NumTriggerTracks,"Jet_NumTriggerTracks/I");
460  fJetTree->Branch("Jet_TriggerTrack_Pt",&dummy,"Jet_TriggerTrack_Pt[Jet_NumTriggerTracks]/F");
461  fJetTree->Branch("Jet_TriggerTrack_dEta",&dummy,"Jet_TriggerTrack_dEta[Jet_NumTriggerTracks]/F");
462  fJetTree->Branch("Jet_TriggerTrack_dPhi",&dummy,"Jet_TriggerTrack_dPhi[Jet_NumTriggerTracks]/F");
463  }
464 
465  fInitialized = kTRUE;
466 }
467 
468 //________________________________________________________________________
470  AliAnalysisTaskEmcalJet("AliAnalysisTaskJetExtractor", kTRUE),
471  #if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
472  fPythonCLI(),
473  #endif
474  fJetTree(0),
475  fVtxTagger(0),
476  fHadronMatchingRadius(0.4),
477  fTrueJetMatchingRadius(0.4),
478  fSecondaryVertexMaxChi2(1e10),
479  fSecondaryVertexMaxDispersion(0.05),
480  fCalculateSecondaryVertices(kTRUE),
481  fCalculateModelBackground(kFALSE),
482  fBackgroundModelFileName(),
483  fBackgroundModelInputParameters(),
484  fBackgroundModelPtOnly(kTRUE),
485  fCustomStartupScript(),
486  fVertexerCuts(0),
487  fSetEmcalJetFlavour(0),
488  fEventPercentage(1.0),
489  fTruthMinLabel(0),
490  fTruthMaxLabel(100000),
491  fSaveTrackPDGCode(kTRUE),
492  fRandomSeed(0),
493  fRandomSeedCones(0),
494  fEventCut_TriggerTrackMinPt(0),
495  fEventCut_TriggerTrackMaxPt(0),
496  fEventCut_TriggerTrackMinLabel(-9999999),
497  fEventCut_TriggerTrackMaxLabel(+9999999),
498  fJetsCont(0),
499  fTracksCont(0),
500  fClustersCont(0),
501  fJetOutputArray(0),
502  fTruthParticleArray(0),
503  fTruthJetsArrayName(""),
504  fTruthJetsRhoName(""),
505  fTruthJetsRhoMassName(""),
506  fTruthParticleArrayName("mcparticles"),
507  fRandomGenerator(0),
508  fRandomGeneratorCones(0),
509  fEventWeight(0.),
510  fImpactParameter(0.),
511  fTriggerTracks_Pt(),
512  fTriggerTracks_Eta(),
513  fTriggerTracks_Phi()
514 {
515  fRandomGenerator = new TRandom3();
516  fRandomGeneratorCones = new TRandom3();
517 
519  fJetTree = new AliEmcalJetTree(GetName());
522  fJetTree->SetSaveJetShapes(kTRUE);
523  DefineOutput(2, TTree::Class());
524 }
525 
526 //________________________________________________________________________
528  AliAnalysisTaskEmcalJet(name, kTRUE),
529  #if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
530  fPythonCLI(),
531  #endif
532  fJetTree(0),
533  fVtxTagger(0),
542  fBackgroundModelPtOnly(kTRUE),
544  fVertexerCuts(0),
546  fEventPercentage(1.0),
547  fTruthMinLabel(0),
548  fTruthMaxLabel(100000),
549  fSaveTrackPDGCode(kTRUE),
550  fRandomSeed(0),
551  fRandomSeedCones(0),
556  fJetsCont(0),
557  fTracksCont(0),
558  fClustersCont(0),
559  fJetOutputArray(0),
562  fTruthJetsRhoName(""),
564  fTruthParticleArrayName("mcparticles"),
565  fRandomGenerator(0),
567  fEventWeight(0.),
568  fImpactParameter(0.),
572 {
573  fRandomGenerator = new TRandom3();
574  fRandomGeneratorCones = new TRandom3();
575 
577  fJetTree = new AliEmcalJetTree(GetName());
580  fJetTree->SetSaveJetShapes(kTRUE);
581  DefineOutput(2, TTree::Class());
582 }
583 
584 //________________________________________________________________________
586 {
587  // Destructor.
588 }
589 
590 //________________________________________________________________________
592 {
594 
595  // ### Basic container settings
597  if(!fJetsCont)
598  AliFatal("Jet input container not found!");
599  fJetsCont->PrintCuts();
601  if(!fTracksCont)
602  AliFatal("Particle input container not found attached to jets!");
605  AliFatal("Cluster input container not found attached to jets!");
606 
607  fRandomGenerator->SetSeed(fRandomSeed);
609 
610  // Activate saving of trigger tracks if this is demanded
613 
614  // ### Initialize the jet tree (settings must all be given at this stage)
620  OpenFile(2);
621  PostData(2, fJetTree->GetTreePointer());
622 
623  // ### Add control histograms (already some created in base task)
624  AddHistogram2D<TH2D>("hTrackCount", "Number of tracks in acceptance vs. centrality", "COLZ", 500, 0., 5000., 100, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}");
625  AddHistogram2D<TH2D>("hBackgroundPt", "Background p_{T} distribution", "", 150, 0., 150., 100, 0, 100, "Background p_{T} (GeV/c)", "Centrality", "dN^{Events}/dp_{T}");
626 
627  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}");
628  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}");
629  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}");
630  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");
631  AddHistogram2D<TH2D>("hJetArea", "Jet area", "COLZ", 200, 0., 2., 100, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
632 
633  AddHistogram2D<TH2D>("hJetPtModel", "Jet p_{T} distribution (model background)", "COLZ", 400, -100., 300., 100, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
634  AddHistogram2D<TH2D>("hJetMassModel", "Jet mass distribution (model background)", "COLZ", 400, -100., 300., 100, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
635  AddHistogram2D<TH2D>("hDeltaPt", "#delta p_{T} distribution", "", 400, -100., 300., 100, 0, 100, "p_{T, cone} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
636 
637  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}");
638  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");
639 
640  AddHistogram1D<TH1D>("hExtractionPercentage", "Extracted jets p_{T} distribution (background subtracted)", "e1p", 400, -100., 300., "p_{T, jet} (GeV/c)", "Extracted percentage");
641 
642  // Track QA plots
643  AddHistogram2D<TH2D>("hTrackPt", "Tracks p_{T} distribution", "", 300, 0., 300., 100, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
644  AddHistogram2D<TH2D>("hTrackPhi", "Track angular distribution in #phi", "LEGO2", 180, 0., 2*TMath::Pi(), 100, 0, 100, "#phi", "Centrality", "dN^{Tracks}/(d#phi)");
645  AddHistogram2D<TH2D>("hTrackEta", "Track angular distribution in #eta", "LEGO2", 100, -2.5, 2.5, 100, 0, 100, "#eta", "Centrality", "dN^{Tracks}/(d#eta)");
646  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");
647 
648  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})");
649  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})");
650 }
651 
652 
653 //________________________________________________________________________
655 {
657 
658  // ### Need to explicitly tell jet container to load rho mass object
659  fJetsCont->LoadRhoMass(InputEvent());
660 
661  if (fTruthParticleArrayName != "")
662  fTruthParticleArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTruthParticleArrayName.Data()));
663 
664  // ### Prepare vertexer
666  {
667  if(!fVertexerCuts)
668  AliFatal("VertexerCuts not given but secondary vertex calculation turned on.");
669  fVtxTagger = new AliHFJetsTaggingVertex();
670  fVtxTagger->SetCuts(fVertexerCuts);
671  }
672 
673  // ### Save tree extraction percentage to histogram
674  std::vector<Float_t> extractionPtBins = fJetTree->GetExtractionPercentagePtBins();
675  std::vector<Float_t> extractionPercentages = fJetTree->GetExtractionPercentages();
676 
677  for(Int_t i=0; i<static_cast<Int_t>(extractionPercentages.size()); i++)
678  {
679  Double_t percentage = extractionPercentages[i];
680  for(Int_t pt=static_cast<Int_t>(extractionPtBins[i*2]); pt<static_cast<Int_t>(extractionPtBins[i*2+1]); pt++)
681  {
682  FillHistogram("hExtractionPercentage", pt, percentage);
683  }
684  }
685 
686  // ### Add HIJING/JEWEL histograms (in case we need it)
687  if(MCEvent())
688  {
689  if(dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()))
690  AddHistogram1D<TH1D>("hHIJING_Ncoll", "N_{coll} from HIJING", "e1p", 1000, 0., 5000, "N_{coll}", "dN^{Events}/dN^{N_{coll}}");
691 
692  if(dynamic_cast<AliGenHepMCEventHeader*>(MCEvent()->GenEventHeader()))
693  {
694  AddHistogram1D<TH1D>("hJEWEL_NProduced", "NProduced from JEWEL", "e1p", 1000, 0., 1000, "Nproduced", "dN^{Events}/dN^{Nproduced}");
695  AddHistogram1D<TH1D>("hJEWEL_ImpactParameter", "Impact parameter from JEWEL", "e1p", 1000, 0., 1, "IP", "dN^{Events}/dN^{IP}");
696  TProfile* evweight = new TProfile("hJEWEL_EventWeight", "Event weight from JEWEL", 1, 0,1);
697  fOutput->Add(evweight);
698  }
699  }
700 
701  // ### Execute shell script at startup
702  if(!fCustomStartupScript.IsNull())
703  {
704  TGrid::Connect("alien://");
705  TFile::Cp(fCustomStartupScript.Data(), "./myScript.sh");
706  gSystem->Exec("bash ./myScript.sh");
707  }
708 
709  // ### Model background
711  {
712  // # Prepare PYTHON cmd, load estimator
713  #if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
714  TGrid::Connect("alien://");
715 
716  if (gSystem->AccessPathName(fBackgroundModelFileName.Data()))
717  AliFatal(Form("Background model %s does not exist!", fBackgroundModelFileName.Data()));
718  TFile::Cp(fBackgroundModelFileName.Data(), "./Model.pkl");
719 
720  fPythonCLI = new TPython();
721  fPythonCLI->Exec("import sys");
722  fPythonCLI->Exec("sys.path.insert(0, './my-local-python/lib/python2.7/site-packages/')");
723  fPythonCLI->Exec("import sklearn, numpy");
724  fPythonCLI->Exec("estimator = sklearn.externals.joblib.load(\"./Model.pkl\")");
725  #endif
726 
727  if (!(InputEvent()->FindListObject(Form("%s_RhoMVA", fJetsCont->GetArrayName().Data())))) {
728  fJetOutputArray = new TClonesArray("AliEmcalJet");
729  fJetOutputArray->SetName(Form("%s_RhoMVA", fJetsCont->GetArrayName().Data()));
730  InputEvent()->AddObject(fJetOutputArray);
731  }
732  else {
733  AliFatal(Form("Jet output array with name %s_RhoMVA already in event! Returning", fJetsCont->GetArrayName().Data()));
734  return;
735  }
736  }
737 
738  PrintConfig();
739 
740 }
741 
742 //________________________________________________________________________
744 {
745  // ################################### EVENT SELECTION
746 
747  if(!IsTriggerTrackInEvent())
748  return kFALSE;
749  if(fRandomGenerator->Rndm() >= fEventPercentage)
750  return kFALSE;
751 
752  // ################################### EVENT PROPERTIES
754 
755  // Load vertex if possible
756  Double_t vtxX = 0;
757  Double_t vtxY = 0;
758  Double_t vtxZ = 0;
759  Long64_t eventID = 0;
760  const AliVVertex* myVertex = InputEvent()->GetPrimaryVertex();
761  if(!myVertex && MCEvent())
762  myVertex = MCEvent()->GetPrimaryVertex();
763  if(myVertex)
764  {
765  vtxX = myVertex->GetX();
766  vtxY = myVertex->GetY();
767  vtxZ = myVertex->GetZ();
768  }
769  // Get event ID from header
770  AliVHeader* eventIDHeader = InputEvent()->GetHeader();
771  if(eventIDHeader)
772  eventID = eventIDHeader->GetEventIdAsLong();
773 
774  // If available, get Ncoll from HIJING
775  if(MCEvent())
776  {
777  if(AliGenHijingEventHeader* mcHeader = dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()))
778  {
779  Double_t ncoll = mcHeader->NN() + mcHeader->NNw() + mcHeader->NwN() + mcHeader->NwNw();
780  FillHistogram("hHIJING_Ncoll", ncoll);
781  }
782  if(AliGenHepMCEventHeader* mcHeader = dynamic_cast<AliGenHepMCEventHeader*>(MCEvent()->GenEventHeader()))
783  {
784  fEventWeight = mcHeader->EventWeight();
785  fImpactParameter = mcHeader->impact_parameter();
786 
787  TProfile* tmpHist = static_cast<TProfile*>(fOutput->FindObject("hJEWEL_EventWeight"));
788  tmpHist->Fill(0., fEventWeight);
789  FillHistogram("hJEWEL_NProduced", mcHeader->NProduced());
790  FillHistogram("hJEWEL_ImpactParameter", fImpactParameter);
791  }
792  }
793 
794  // ################################### MAIN JET LOOP
795  fJetsCont->ResetCurrentID();
796  Int_t jetCount = 0;
797  while(AliEmcalJet *jet = fJetsCont->GetNextAcceptJet())
798  {
800 
801  Double_t matchDistance = 0;
802  Double_t matchedJetPt = 0;
803  Double_t matchedJetMass = 0;
804  Double_t truePtFraction = 0;
805  Int_t currentJetType_HM = 0;
806  Int_t currentJetType_PM = 0;
807  Int_t currentJetType_IC = 0;
809  {
810  // Get jet type from MC (hadron matching, parton matching definition - for HF jets)
811  GetJetType(jet, currentJetType_HM, currentJetType_PM, currentJetType_IC);
812  // Get true estimators: for pt, jet mass, ...
813  truePtFraction = GetTrueJetPtFraction(jet);
814  matchDistance = GetMatchedTrueJetObservables(jet, matchedJetPt, matchedJetMass);
815  }
816 
817  // ### CONSTITUENT LOOP: Retrieve PID values + impact parameters
818  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;
819  std::vector<Float_t> vec_d0; std::vector<Float_t> vec_d0cov; std::vector<Float_t> vec_z0; std::vector<Float_t> vec_z0cov;
820 
822  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
823  {
824  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
825  if(!particle) continue;
827  {
828  Float_t sigITS = 0; Float_t sigTPC = 0; Float_t sigTOF = 0; Float_t sigTRD = 0; Short_t recoPID = 0; Int_t truePID = 0;
829  AddPIDInformation(particle, sigITS, sigTPC, sigTOF, sigTRD, recoPID, truePID);
830  vecSigITS.push_back(sigITS); vecSigTPC.push_back(sigTPC); vecSigTOF.push_back(sigTOF); vecSigTRD.push_back(sigTRD); vecRecoPID.push_back(recoPID); vecTruePID.push_back(truePID);
831  }
833  {
834  Float_t d0 = 0; Float_t d0cov = 0; Float_t z0 = 0; Float_t z0cov = 0;
835  GetTrackImpactParameters(myVertex, dynamic_cast<AliAODTrack*>(particle), d0, d0cov, z0, z0cov);
836  vec_d0.push_back(d0); vec_d0cov.push_back(d0cov); vec_z0.push_back(z0); vec_z0cov.push_back(z0cov);
837  }
838  }
839 
840  // Reconstruct secondary vertices
841  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;
843  ReconstructSecondaryVertices(myVertex, jet, secVtx_X, secVtx_Y, secVtx_Z, secVtx_Mass, secVtx_Lxy, secVtx_SigmaLxy, secVtx_Chi2, secVtx_Dispersion);
844 
845  // Now change trigger tracks eta/phi to dEta/dPhi relative to jet
846  std::vector<Float_t> triggerTracks_dEta(fTriggerTracks_Eta);
847  std::vector<Float_t> triggerTracks_dPhi(fTriggerTracks_Phi);
848  for(UInt_t i=0; i<triggerTracks_dEta.size(); i++)
849  {
850  triggerTracks_dEta[i] = jet->Eta() - fTriggerTracks_Eta[i];
851  triggerTracks_dPhi[i] = TMath::Min(TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]), TMath::TwoPi() - TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]));
852  if( ((TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]) <= TMath::Pi()) && (jet->Phi() - fTriggerTracks_Phi[i] > 0))
853  || ((TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]) > TMath::Pi()) && (jet->Phi() - fTriggerTracks_Phi[i] < 0)) )
854  triggerTracks_dPhi[i] = -triggerTracks_dPhi[i];
855  }
856 
858  {
859  Float_t pt_ML = 0;
860  Float_t mass_ML = 0;
861  GetPtAndMassFromModel(jet, pt_ML, mass_ML);
862  AliEmcalJet *outJet = new ((*fJetOutputArray)[jetCount]) AliEmcalJet(pt_ML, jet->Eta(), jet->Phi(), mass_ML);
863  outJet->SetLabel(jet->Label());
864  outJet->SetArea(jet->Area());
865  outJet->SetJetAcceptanceType(jet->GetJetAcceptanceType());
866  FillHistogram("hJetPtModel", pt_ML, fCent);
867  FillHistogram("hJetMassModel", mass_ML, fCent);
868 
869  }
870 
872  {
873  // Calculate jet shapes and set them in the tree (some are retrieved in the tree itself)
874  Double_t leSub_noCorr = 0;
875  Double_t radialMoment = 0;
876  Double_t momentumDispersion = 0;
877  Double_t constPtMean = 0;
878  Double_t constPtMedian = 0;
879  CalculateJetShapes(jet, leSub_noCorr, radialMoment, momentumDispersion, constPtMean, constPtMedian);
880  fJetTree->SetJetShapesInBuffer(leSub_noCorr, radialMoment, momentumDispersion, constPtMean, constPtMedian);
881  }
882  // Fill jet to tree
883  Bool_t accepted = fJetTree->AddJetToTree(jet, vtxX, vtxY, vtxZ, fCent, eventID, InputEvent()->GetMagneticField(),
884  currentJetType_PM,currentJetType_HM,currentJetType_IC,matchDistance,matchedJetPt,matchedJetMass,truePtFraction,fPtHard,fEventWeight,fImpactParameter,
885  vecSigITS, vecSigTPC, vecSigTOF, vecSigTRD, vecRecoPID, vecTruePID,
886  fTriggerTracks_Pt, triggerTracks_dEta, triggerTracks_dPhi,
887  vec_d0, vec_z0, vec_d0cov, vec_z0cov,
888  secVtx_X, secVtx_Y, secVtx_Z, secVtx_Mass, secVtx_Lxy, secVtx_SigmaLxy, secVtx_Chi2, secVtx_Dispersion);
889 
890  if(accepted)
891  FillHistogram("hJetPtExtracted", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
892  jetCount++;
893  }
894 
895  // ################################### PARTICLE PROPERTIES
896  fTracksCont->ResetCurrentID();
897  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
899 
900  return kTRUE;
901 }
902 
903 //________________________________________________________________________
905 {
906  #if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
907  TString arrayStr = GetBackgroundModelArrayString(jet);
908  // ####### Run Python script that does inference on jet
909  fPythonCLI->Exec(Form("data_inference = numpy.array([[%s]])", arrayStr.Data()));
910  fPythonCLI->Exec("result = estimator.predict(data_inference)");
912  {
913  pt_ML = fPythonCLI->Eval("result[0][0]");
914  mass_ML = fPythonCLI->Eval("result[0][1]");
915  }
916  else
917  {
918  pt_ML = fPythonCLI->Eval("result[0]");
919  mass_ML = 0.0;
920  }
921  #endif
922 }
923 
924 //________________________________________________________________________
926 {
927  // ####### Calculate inference input parameters
928  std::vector<Int_t> index_sorted_list = jet->GetPtSortedTrackConstituentIndexes(fJetsCont->GetParticleContainer()->GetArray());
929  Int_t numConst = index_sorted_list.size();
930  Double_t* constPts = new Double_t[TMath::Max(Int_t(index_sorted_list.size()), 10)];
931  for(Int_t i = 0; i < TMath::Max(Int_t(index_sorted_list.size()), 10); i++)
932  constPts[i] = 0;
933  for(Int_t i = 0; i < numConst; i++)
934  {
935  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(index_sorted_list.at(i), fJetsCont->GetParticleContainer()->GetArray()));
936  constPts[i] = particle->Pt();
937  }
938 
939  // Calculate jet shapes that could be demanded
940  Double_t leSub_noCorr = 0;
941  Double_t radialMoment = 0;
942  Double_t momentumDispersion = 0;
943  Double_t constPtMean = 0;
944  Double_t constPtMedian = 0;
945  CalculateJetShapes(jet, leSub_noCorr, radialMoment, momentumDispersion, constPtMean, constPtMedian);
946 
947  TString resultStr = "";
948  TObjArray* data_tokens = fBackgroundModelInputParameters.Tokenize(",");
949  for(Int_t iToken = 0; iToken < data_tokens->GetEntries(); iToken++)
950  {
951  TString token = ((TObjString *)(data_tokens->At(iToken)))->String();
952  if(token == "Jet_Pt")
953  resultStr += Form("%E", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
954  else if(token == "Jet_NumConstituents")
955  resultStr += Form("%E", (Double_t)numConst);
956  else if(token == "Jet_Shape_Mass_NoCorr")
957  resultStr += Form("%E", jet->M());
958  else if(token == "Jet_Shape_Mass_DerivCorr_1")
959  resultStr += Form("%E", jet->GetShapeProperties()->GetFirstOrderSubtracted());
960  else if(token == "Jet_Shape_Mass_DerivCorr_2")
961  resultStr += Form("%E", jet->GetShapeProperties()->GetSecondOrderSubtracted());
962  else if(token == "Jet_Shape_pTD_DerivCorr_1")
963  resultStr += Form("%E", jet->GetShapeProperties()->GetFirstOrderSubtractedpTD());
964  else if(token == "Jet_Shape_pTD_DerivCorr_2")
965  resultStr += Form("%E", jet->GetShapeProperties()->GetSecondOrderSubtractedpTD());
966  else if(token == "Jet_Shape_LeSub_NoCorr")
967  resultStr += Form("%E", leSub_noCorr);
968  else if(token == "Jet_Shape_LeSub_DerivCorr")
969  resultStr += Form("%E", jet->GetShapeProperties()->GetSecondOrderSubtractedLeSub());
970  else if(token == "Jet_Shape_Sigma2_DerivCorr_1")
971  resultStr += Form("%E", jet->GetShapeProperties()->GetFirstOrderSubtractedSigma2());
972  else if(token == "Jet_Shape_Sigma2_DerivCorr_2")
973  resultStr += Form("%E", jet->GetShapeProperties()->GetSecondOrderSubtractedSigma2());
974  else if(token == "Jet_Shape_Angularity_DerivCorr_1")
975  resultStr += Form("%E", jet->GetShapeProperties()->GetFirstOrderSubtractedAngularity());
976  else if(token == "Jet_Shape_Angularity_DerivCorr_2")
977  resultStr += Form("%E", jet->GetShapeProperties()->GetSecondOrderSubtractedAngularity());
978  else if(token == "Jet_Shape_Circularity_DerivCorr_1")
979  resultStr += Form("%E", jet->GetShapeProperties()->GetFirstOrderSubtractedCircularity());
980  else if(token == "Jet_Shape_Circularity_DerivCorr_2")
981  resultStr += Form("%E", jet->GetShapeProperties()->GetSecondOrderSubtractedCircularity());
982  else if(token == "Jet_Shape_NumConst_DerivCorr")
983  resultStr += Form("%E", jet->GetShapeProperties()->GetSecondOrderSubtractedConstituent());
984  else if(token == "Jet_Shape_RadialMoment")
985  resultStr += Form("%E", radialMoment);
986  else if(token == "Jet_Shape_RadialMoment_Shrinked")
987  {
988  if (radialMoment < 0)
989  resultStr += Form("%E", 0.0);
990  else if (radialMoment > 1)
991  resultStr += Form("%E", 1.0);
992  else
993  resultStr += Form("%E", radialMoment);
994  }
995  else if(token == "Jet_Shape_MomentumDispersion")
996  resultStr += Form("%E", momentumDispersion);
997  else if(token == "Jet_Shape_ConstPtMean")
998  resultStr += Form("%E", constPtMean);
999  else if(token == "Jet_Shape_ConstPtMedian")
1000  resultStr += Form("%E", constPtMedian);
1001  else if(token == "Event_BackgroundDensity")
1002  resultStr += Form("%E", fJetsCont->GetRhoVal());
1003  else if(token == "Event_BackgroundDensityMass")
1004  resultStr += Form("%E", fJetsCont->GetRhoMassVal());
1005  else if(token == "Jet_Area")
1006  resultStr += Form("%E", jet->Area());
1007  else if(token.BeginsWith("Jet_Shape_ConstPt"))
1008  {
1009  TString num = token(17,(token.Length()-17));
1010  resultStr += Form("%E", constPts[num.Atoi()]);
1011  }
1012 
1013  // Add comma after numbers
1014  if((iToken < data_tokens->GetEntries()-1))
1015  resultStr += ", ";
1016  }
1017 
1018  data_tokens->SetOwner();
1019  delete data_tokens;
1020  delete[] constPts;
1021  return resultStr;
1022 }
1023 
1024 
1025 //________________________________________________________________________
1027 {
1028  // Cut for trigger track requirement
1030  {
1031  // Clear vector of trigger tracks
1032  fTriggerTracks_Pt.clear();
1033  fTriggerTracks_Eta.clear();
1034  fTriggerTracks_Phi.clear();
1035 
1036  // Go through all tracks and check whether trigger tracks can be found
1037  fTracksCont->ResetCurrentID();
1038  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
1039  {
1040  if( (track->GetLabel() >= fEventCut_TriggerTrackMinLabel) && (track->GetLabel() < fEventCut_TriggerTrackMaxLabel) )
1041  if( (track->Pt() >= fEventCut_TriggerTrackMinPt) && (track->Pt() < fEventCut_TriggerTrackMaxPt) )
1042  {
1043  fTriggerTracks_Pt.push_back(track->Pt());
1044  fTriggerTracks_Eta.push_back(track->Eta());
1045  fTriggerTracks_Phi.push_back(track->Phi());
1046  }
1047  }
1048  // No particle has been found that fulfills requirement -> Do not accept event
1049  if(fTriggerTracks_Pt.size() == 0)
1050  return kFALSE;
1051  }
1052  return kTRUE;
1053 }
1054 
1055 //________________________________________________________________________
1056 void AliAnalysisTaskJetExtractor::CalculateJetShapes(AliEmcalJet* jet, Double_t& leSub_noCorr, Double_t& radialMoment, Double_t& momentumDispersion, Double_t& constPtMean, Double_t& constPtMedian)
1057 {
1058  // #### Calculate mean, median of constituents, radial moment, momentum dispersion, leSub (no correction)
1059  Double_t jetCorrectedPt = jet->Pt() - jet->Area() * fJetsCont->GetRhoVal();
1060  Double_t jetLeadingHadronPt = -999.;
1061  Double_t jetSubleadingHadronPt = -999.;
1062  Double_t jetSum = 0;
1063  Double_t jetSum2 = 0;
1064  constPtMean = 0;
1065  constPtMedian = 0;
1066  radialMoment = 0;
1067  momentumDispersion = 0;
1068  std::vector<Int_t> index_sorted_list = jet->GetPtSortedTrackConstituentIndexes(fJetsCont->GetParticleContainer()->GetArray());
1069  Int_t numConst = index_sorted_list.size();
1070  Double_t* constPts = new Double_t[TMath::Max(Int_t(index_sorted_list.size()), 10)];
1071  for(Int_t i = 0; i < TMath::Max(Int_t(index_sorted_list.size()), 10); i++)
1072  constPts[i] = 0;
1073 
1074  // Loop over all constituents and do jet shape calculations
1075  for (Int_t i=0;i<numConst;i++)
1076  {
1077  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(index_sorted_list.at(i), fJetsCont->GetParticleContainer()->GetArray()));
1078  constPtMean += particle->Pt();
1079  constPts[i] = particle->Pt();
1080  if(particle->Pt() > jetLeadingHadronPt)
1081  {
1082  jetSubleadingHadronPt = jetLeadingHadronPt;
1083  jetLeadingHadronPt = particle->Pt();
1084  }
1085  else if(particle->Pt() > jetSubleadingHadronPt)
1086  jetSubleadingHadronPt = particle->Pt();
1087 
1088  Double_t deltaPhi = TMath::Min(TMath::Abs(jet->Phi()-particle->Phi()),TMath::TwoPi() - TMath::Abs(jet->Phi()-particle->Phi()));
1089  Double_t deltaR = TMath::Sqrt( (jet->Eta() - particle->Eta())*(jet->Eta() - particle->Eta()) + deltaPhi*deltaPhi );
1090 
1091  jetSum += particle->Pt();
1092  jetSum2 += particle->Pt()*particle->Pt();
1093  radialMoment += particle->Pt() * deltaR;
1094  }
1095 
1096  if(jetCorrectedPt)
1097  radialMoment /= jetCorrectedPt;
1098  if(numConst)
1099  {
1100  constPtMean /= numConst;
1101  constPtMedian = TMath::Median(numConst, constPts);
1102  }
1103 
1104  if(numConst > 1)
1105  leSub_noCorr = jetLeadingHadronPt - jetSubleadingHadronPt;
1106  else
1107  leSub_noCorr = jetLeadingHadronPt;
1108 
1109  if(jetSum)
1110  momentumDispersion = TMath::Sqrt(jetSum2)/jetSum;
1111 
1112 }
1113 
1114 //________________________________________________________________________
1116 {
1117  // #################################################################################
1118  // ##### FRACTION OF TRUE PT IN JET: Defined as "not from toy"
1119  Double_t pt_nonMC = 0.;
1120  Double_t pt_all = 0.;
1121  Double_t truePtFraction = 0;
1122 
1123  // Loop over all jet tracks+clusters
1124  for(Int_t iConst = 0; iConst < jet->GetNumberOfTracks(); iConst++)
1125  {
1126  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(iConst, fTracksCont->GetArray()));
1127  if(!particle) continue;
1128  if(particle->Pt() < 1e-6) continue;
1129 
1130  // Particles marked w/ labels within label range are considered from toy
1131  if( (particle->GetLabel() >= fTruthMinLabel) && (particle->GetLabel() < fTruthMaxLabel))
1132  pt_nonMC += particle->Pt();
1133  pt_all += particle->Pt();
1134  }
1135  for(Int_t iConst = 0; iConst < jet->GetNumberOfClusters(); iConst++)
1136  {
1137  AliVCluster* cluster = static_cast<AliVCluster*>(jet->ClusterAt(iConst, fClustersCont->GetArray()));
1138  if(!cluster) continue;
1139  if(cluster->E() < 1e-6) continue;
1140 
1141  // #### Retrieve cluster pT
1142  Double_t vtxX = 0;
1143  Double_t vtxY = 0;
1144  Double_t vtxZ = 0;
1145  const AliVVertex* myVertex = InputEvent()->GetPrimaryVertex();
1146  if(!myVertex && MCEvent())
1147  myVertex = MCEvent()->GetPrimaryVertex();
1148  if(myVertex)
1149  {
1150  vtxX = myVertex->GetX();
1151  vtxY = myVertex->GetY();
1152  vtxZ = myVertex->GetZ();
1153  }
1154 
1155  Double_t primVtx[3] = {vtxX, vtxY, vtxZ};
1156  TLorentzVector clusterMomentum;
1157  cluster->GetMomentum(clusterMomentum, primVtx);
1158  Double_t ClusterPt = clusterMomentum.Perp();
1159  // ####
1160 
1161  // Particles marked w/ labels within label range are considered from toy
1162  if( (cluster->GetLabel() >= fTruthMinLabel) && (cluster->GetLabel() < fTruthMaxLabel))
1163  pt_nonMC += ClusterPt;
1164  pt_all += ClusterPt;
1165  }
1166 
1167  if(pt_all)
1168  truePtFraction = (pt_nonMC/pt_all);
1169  return truePtFraction;
1170 }
1171 
1172 //________________________________________________________________________
1174 {
1175  // #################################################################################
1176  // ##### OBSERVABLES FROM MATCHED JETS: Jet pt, jet mass
1177  Double_t bestMatchDeltaR = 8.; // 8 is higher than maximum possible matching distance
1178  if(fTruthJetsArrayName != "")
1179  {
1180  // "True" background for pt
1181  AliRhoParameter* rho = static_cast<AliRhoParameter*>(InputEvent()->FindListObject(fTruthJetsRhoName.Data()));
1182  Double_t trueRho = 0;
1183  if(rho)
1184  trueRho = rho->GetVal();
1185 
1186  // "True" background for mass
1187  AliRhoParameter* rhoMass = static_cast<AliRhoParameter*>(InputEvent()->FindListObject(fTruthJetsRhoMassName.Data()));
1188  TClonesArray* truthArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fTruthJetsArrayName.Data())));
1189 
1190  // Loop over all true jets to find the best match
1191  matchedJetPt = 0;
1192  matchedJetMass = 0;
1193  if(truthArray)
1194  for(Int_t i=0; i<truthArray->GetEntries(); i++)
1195  {
1196  AliEmcalJet* truthJet = static_cast<AliEmcalJet*>(truthArray->At(i));
1197  if(truthJet->Pt() < 0.15)
1198  continue;
1199 
1200  Double_t deltaEta = (truthJet->Eta()-jet->Eta());
1201  Double_t deltaPhi = TMath::Min(TMath::Abs(truthJet->Phi()-jet->Phi()),TMath::TwoPi() - TMath::Abs(truthJet->Phi()-jet->Phi()));
1202  Double_t deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
1203 
1204  // Cut jets too far away
1205  if (deltaR > fTrueJetMatchingRadius)
1206  continue;
1207 
1208  // Search for the best match
1209  if(deltaR < bestMatchDeltaR)
1210  {
1211  bestMatchDeltaR = deltaR;
1212  matchedJetPt = truthJet->Pt() - truthJet->Area()* trueRho;
1213  if(rhoMass)
1214  matchedJetMass = truthJet->GetShapeProperties()->GetSecondOrderSubtracted();
1215  else
1216  matchedJetMass = truthJet->M();
1217  }
1218  }
1219  }
1220 
1221  return bestMatchDeltaR;
1222 }
1223 
1224 
1225 //________________________________________________________________________
1227 {
1229 
1230  if(!fTruthParticleArray)
1231  return;
1232 
1233  typeHM = 0;
1234  typePM = 0;
1235 
1236  AliAODMCParticle* parton[2];
1237  parton[0] = (AliAODMCParticle*) fVtxTagger->IsMCJetParton(fTruthParticleArray, jet, radius); // method 1 (parton matching)
1238  parton[1] = (AliAODMCParticle*) fVtxTagger->IsMCJetMeson(fTruthParticleArray, jet, radius); // method 2 (hadron matching)
1239 
1240  if (parton[0]) {
1241  Int_t pdg = TMath::Abs(parton[0]->PdgCode());
1242  typePM = pdg;
1243  }
1244 
1245  if (!parton[1])
1246  {
1247  // No HF jet (according to hadron matching) -- now also separate jets in udg (1) and s-jets (3)
1248  if(IsStrangeJet(jet))
1249  typeHM = 3;
1250  else
1251  typeHM = 1;
1252  }
1253  else {
1254  Int_t pdg = TMath::Abs(parton[1]->PdgCode());
1255  if(fVtxTagger->IsDMeson(pdg)) typeHM = 4;
1256  else if (fVtxTagger->IsBMeson(pdg)) typeHM = 5;
1257  }
1258 
1259  // Set flavour of AliEmcalJet object (set ith bit while i corresponds to type)
1261  jet->AddFlavourTag(static_cast<Int_t>(TMath::Power(2, typeHM)));
1262 
1263 
1264  const AliEmcalPythiaInfo* partonsInfo = GetPythiaInfo();
1265  typeIC = 0;
1266  if (partonsInfo)
1267  {
1268  // Get primary partons directions
1269  Double_t parton1phi = partonsInfo->GetPartonPhi6();
1270  Double_t parton1eta = partonsInfo->GetPartonEta6();
1271  Double_t parton2phi = partonsInfo->GetPartonPhi7();
1272  Double_t parton2eta = partonsInfo->GetPartonEta7();
1273 
1274 
1275  Double_t delta1Eta = (parton1eta-jet->Eta());
1276  Double_t delta1Phi = TMath::Min(TMath::Abs(parton1phi-jet->Phi()),TMath::TwoPi() - TMath::Abs(parton1phi-jet->Phi()));
1277  Double_t delta1R = TMath::Sqrt(delta1Eta*delta1Eta + delta1Phi*delta1Phi);
1278  Double_t delta2Eta = (parton2eta-jet->Eta());
1279  Double_t delta2Phi = TMath::Min(TMath::Abs(parton2phi-jet->Phi()),TMath::TwoPi() - TMath::Abs(parton2phi-jet->Phi()));
1280  Double_t delta2R = TMath::Sqrt(delta2Eta*delta2Eta + delta2Phi*delta2Phi);
1281 
1282  // Check if one of the partons if closer than matching criterion
1283  Bool_t matched = (delta1R < fJetsCont->GetJetRadius()/2.) || (delta2R < fJetsCont->GetJetRadius()/2.);
1284 
1285  // Matching criterion fulfilled -> Set flag to closest
1286  if(matched)
1287  {
1288  if(delta1R < delta2R)
1289  typeIC = partonsInfo->GetPartonFlag6();
1290  else
1291  typeIC = partonsInfo->GetPartonFlag7();
1292  }
1293  }
1294 }
1295 
1296 //________________________________________________________________________
1297 void AliAnalysisTaskJetExtractor::GetTrackImpactParameters(const AliVVertex* vtx, AliAODTrack* track, Float_t& d0, Float_t& d0cov, Float_t& z0, Float_t& z0cov)
1298 {
1299  if (track)
1300  {
1301  Double_t d0rphiz[2],covd0[3];
1302  Bool_t isDCA=track->PropagateToDCA(vtx,InputEvent()->GetMagneticField(),3.0,d0rphiz,covd0);
1303  if(isDCA)
1304  {
1305  d0 = d0rphiz[0];
1306  z0 = d0rphiz[1];
1307  d0cov = covd0[0];
1308  z0cov = covd0[2];
1309  }
1310  }
1311 }
1312 
1313 //________________________________________________________________________
1314 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)
1315 {
1316  if(!primVtx)
1317  return;
1318 
1319  // Create ESD vertex from the existing AliVVertex
1320  Double_t vtxPos[3] = {primVtx->GetX(), primVtx->GetY(), primVtx->GetZ()};
1321  Double_t covMatrix[6] = {0};
1322  primVtx->GetCovarianceMatrix(covMatrix);
1323  AliESDVertex* esdVtx = new AliESDVertex(vtxPos, covMatrix, primVtx->GetChi2(), primVtx->GetNContributors());
1324 
1325  TClonesArray* secVertexArr = 0;
1326  vctr_pair_dbl_int arrDispersion;
1327  arrDispersion.reserve(5);
1329  {
1330  //###########################################################################
1331  // ********* Calculate secondary vertices
1332  // Derived from AliAnalysisTaskEmcalJetBtagSV
1333  secVertexArr = new TClonesArray("AliAODVertex");
1334  Int_t nDauRejCount = 0;
1335  Int_t nVtx = fVtxTagger->FindVertices(jet,
1336  fTracksCont->GetArray(),
1337  (AliAODEvent*)InputEvent(),
1338  esdVtx,
1339  InputEvent()->GetMagneticField(),
1340  secVertexArr,
1341  0,
1342  arrDispersion,
1343  nDauRejCount);
1344 
1345 
1346  if(nVtx < 0)
1347  {
1348  secVertexArr->Clear();
1349  delete secVertexArr;
1350  return;
1351  }
1352  //###########################################################################
1353  }
1354  else // Load HF vertex branch from AOD event, if possible
1355  {
1356  secVertexArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("VerticesHF"));
1357  if(!secVertexArr)
1358  return;
1359  }
1360 
1361  // Loop over all potential secondary vertices
1362  for(Int_t i=0; i<secVertexArr->GetEntriesFast(); i++)
1363  {
1364  AliAODVertex* secVtx = (AliAODVertex*)(secVertexArr->UncheckedAt(i));
1366  if((strcmp(secVtx->GetParent()->ClassName(), "AliAODRecoDecayHF3Prong")))
1367  continue;
1368 
1369  // Calculate vtx distance
1370  Double_t effX = secVtx->GetX() - esdVtx->GetX();
1371  Double_t effY = secVtx->GetY() - esdVtx->GetY();
1372  //Double_t effZ = secVtx->GetZ() - esdVtx->GetZ();
1373 
1374  // ##### Vertex properties
1375  // vertex dispersion
1376  Double_t dispersion = arrDispersion[i].first;
1377 
1378  // invariant mass
1379  Double_t mass = fVtxTagger->GetVertexInvariantMass(secVtx);
1380 
1381  // signed length
1382  Double_t Lxy = TMath::Sqrt(effX*effX + effY*effY);
1383  Double_t jetP[3]; jet->PxPyPz(jetP);
1384  Double_t signLxy = effX * jetP[0] + effY * jetP[1];
1385  if (signLxy < 0.) Lxy *= -1.;
1386 
1387  Double_t sigmaLxy = 0;
1388  AliAODVertex* aodVtx = (AliAODVertex*)(primVtx);
1389  if (aodVtx)
1390  sigmaLxy = aodVtx->ErrorDistanceXYToVertex(secVtx);
1391 
1392  // Add secondary vertices if they fulfill the conditions
1393  if( (dispersion > fSecondaryVertexMaxDispersion) || (TMath::Abs(secVtx->GetChi2perNDF()) > fSecondaryVertexMaxChi2) )
1394  continue;
1395 
1396  secVtx_X.push_back(secVtx->GetX()); secVtx_Y.push_back(secVtx->GetY()); secVtx_Z.push_back(secVtx->GetZ()); secVtx_Chi2.push_back(secVtx->GetChi2perNDF());
1397  secVtx_Dispersion.push_back(dispersion); secVtx_Mass.push_back(mass); secVtx_Lxy.push_back(Lxy); secVtx_SigmaLxy.push_back(sigmaLxy);
1398  }
1399 
1401  {
1402  secVertexArr->Clear();
1403  delete secVertexArr;
1404  }
1405  delete esdVtx;
1406 }
1407 
1408 //________________________________________________________________________
1410 {
1411  // Do hadron matching jet type tagging using mcparticles
1412  // ... if not explicitly deactivated
1413  if (fTruthParticleArray)
1414  {
1415  for(Int_t i=0; i<fTruthParticleArray->GetEntries();i++)
1416  {
1417  AliAODMCParticle* part = (AliAODMCParticle*)fTruthParticleArray->At(i);
1418  if(!part) continue;
1419 
1420  // Check if the particle has strangeness
1421  Int_t absPDG = TMath::Abs(part->PdgCode());
1422  if ((absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
1423  {
1424  // Check if particle is in a radius around the jet
1425  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
1427  continue;
1428  else
1429  return kTRUE;
1430  }
1431  }
1432  }
1433  // As fallback, the MC stack will be tried
1434  else if(MCEvent() && (MCEvent()->Stack()))
1435  {
1436  AliStack* stack = MCEvent()->Stack();
1437  // Go through the whole particle stack
1438  for(Int_t i=0; i<stack->GetNtrack(); i++)
1439  {
1440  TParticle *part = stack->Particle(i);
1441  if(!part) continue;
1442 
1443  // Check if the particle has strangeness
1444  Int_t absPDG = TMath::Abs(part->GetPdgCode());
1445  if ((absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
1446  {
1447  // Check if particle is in a radius around the jet
1448  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
1450  continue;
1451  else
1452  return kTRUE;
1453  }
1454  }
1455  }
1456  return kFALSE;
1457 
1458 }
1459 //________________________________________________________________________
1461 {
1462  // ### Event control plots
1464  FillHistogram("hBackgroundPt", fJetsCont->GetRhoVal(), fCent);
1465 }
1466 
1467 //________________________________________________________________________
1469 {
1470  // ### Jet control plots
1471  FillHistogram("hJetPtRaw", jet->Pt(), fCent);
1472  FillHistogram("hJetPt", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
1473  FillHistogram("hJetPhiEta", jet->Phi(), jet->Eta());
1474  FillHistogram("hJetArea", jet->Area(), fCent);
1475 
1476  // ### Jet constituent plots
1477  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
1478  {
1479  AliVParticle* jetConst = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
1480  if(!jetConst) continue;
1481 
1482  // Constituent eta/phi (relative to jet)
1483  Double_t deltaEta = jet->Eta() - jetConst->Eta();
1484  Double_t deltaPhi = TMath::Min(TMath::Abs(jet->Phi() - jetConst->Phi()), TMath::TwoPi() - TMath::Abs(jet->Phi() - jetConst->Phi()));
1485  if(!((jet->Phi() - jetConst->Phi() < 0) && (jet->Phi() - jetConst->Phi() <= TMath::Pi())))
1486  deltaPhi = -deltaPhi;
1487 
1488  FillHistogram("hConstituentPt", jetConst->Pt(), fCent);
1489  FillHistogram("hConstituentPhiEta", deltaPhi, deltaEta);
1490  }
1491 
1492  // ### Random cone / delta pT plots
1493  const Int_t kNumRandomConesPerEvent = 4;
1494  for(Int_t iCone=0; iCone<kNumRandomConesPerEvent; iCone++)
1495  {
1496  // Throw random cone
1497  Double_t tmpRandConeEta = fJetsCont->GetJetEtaMin() + fRandomGeneratorCones->Rndm()*TMath::Abs(fJetsCont->GetJetEtaMax()-fJetsCont->GetJetEtaMin());
1498  Double_t tmpRandConePhi = fRandomGeneratorCones->Rndm()*TMath::TwoPi();
1499  Double_t tmpRandConePt = 0;
1500  // Fill pT that is in cone
1501  fTracksCont->ResetCurrentID();
1502  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
1503  if(IsTrackInCone(track, tmpRandConeEta, tmpRandConePhi, fJetsCont->GetJetRadius()))
1504  tmpRandConePt += track->Pt();
1505 
1506  // Fill histograms
1507  FillHistogram("hDeltaPt", tmpRandConePt - fJetsCont->GetRhoVal()*fJetsCont->GetJetRadius()*fJetsCont->GetJetRadius()*TMath::Pi(), fCent);
1508  }
1509 }
1510 
1511 //________________________________________________________________________
1513 {
1514  FillHistogram("hTrackPt", track->Pt(), fCent);
1515  FillHistogram("hTrackPhi", track->Phi(), fCent);
1516  FillHistogram("hTrackEta", track->Eta(), fCent);
1517  FillHistogram("hTrackEtaPt", track->Eta(), track->Pt());
1518  FillHistogram("hTrackPhiPt", track->Phi(), track->Pt());
1519  FillHistogram("hTrackPhiEta", track->Phi(), track->Eta());
1520 }
1521 
1522 //________________________________________________________________________
1523 void AliAnalysisTaskJetExtractor::AddPIDInformation(AliVParticle* particle, Float_t& sigITS, Float_t& sigTPC, Float_t& sigTOF, Float_t& sigTRD, Short_t& recoPID, Int_t& truePID)
1524 {
1525  truePID = 9;
1526  if(!particle) return;
1527 
1528  // If we have AODs, retrieve particle PID signals
1529  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(particle);
1530 
1531  if(aodtrack)
1532  {
1533  // Get AOD value from reco
1534  recoPID = aodtrack->GetMostProbablePID();
1535  AliAODPid* pidObj = aodtrack->GetDetPid();
1536  if(!pidObj)
1537  return;
1538 
1539  sigITS = pidObj->GetITSsignal();
1540  sigTPC = pidObj->GetTPCsignal();
1541  sigTOF = pidObj->GetTOFsignal();
1542  sigTRD = pidObj->GetTRDsignal();
1543  }
1544 
1545  // Get truth values if we are on MC
1547  {
1548  for(Int_t i=0; i<fTruthParticleArray->GetEntries();i++)
1549  {
1550  AliAODMCParticle* mcParticle = dynamic_cast<AliAODMCParticle*>(fTruthParticleArray->At(i));
1551  if(!mcParticle) continue;
1552 
1553  if (mcParticle->GetLabel() == particle->GetLabel())
1554  {
1555  if(fSaveTrackPDGCode)
1556  {
1557  truePID = mcParticle->PdgCode();
1558  }
1559  else // Use same convention as for PID in AODs
1560  {
1561  if(TMath::Abs(mcParticle->PdgCode()) == 2212) // proton
1562  truePID = 4;
1563  else if (TMath::Abs(mcParticle->PdgCode()) == 211) // pion
1564  truePID = 2;
1565  else if (TMath::Abs(mcParticle->PdgCode()) == 321) // kaon
1566  truePID = 3;
1567  else if (TMath::Abs(mcParticle->PdgCode()) == 11) // electron
1568  truePID = 0;
1569  else if (TMath::Abs(mcParticle->PdgCode()) == 13) // muon
1570  truePID = 1;
1571  else if (TMath::Abs(mcParticle->PdgCode()) == 700201) // deuteron
1572  truePID = 5;
1573  else if (TMath::Abs(mcParticle->PdgCode()) == 700301) // triton
1574  truePID = 6;
1575  else if (TMath::Abs(mcParticle->PdgCode()) == 700302) // He3
1576  truePID = 7;
1577  else if (TMath::Abs(mcParticle->PdgCode()) == 700202) // alpha
1578  truePID = 8;
1579  else
1580  truePID = 9;
1581  }
1582 
1583  break;
1584  }
1585  }
1586  }
1587 }
1588 
1589 //________________________________________________________________________
1591 {
1592  std::cout << "#########################################\n";
1593  std::cout << "Settings for extractor task " << GetName() << std::endl;
1594  std::cout << std::endl;
1595 
1596  std::cout << "### Event cuts ###" << std::endl;
1597  std::cout << std::endl;
1599  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;
1600  if(fEventPercentage < 1)
1601  std::cout << Form("* Artificially lowered event efficiency: %f", fEventPercentage) << std::endl;
1603  std::cout << Form("* Random seeds explicitly set: %lu (for event/jet eff.), %lu (RCs)", fRandomSeed, fRandomSeedCones) << std::endl;
1604  std::cout << std::endl;
1605 
1606  std::cout << "### Jet properties ###" << std::endl;
1607  std::cout << "* Jet array name = " << fJetsCont->GetName() << std::endl;
1608  std::cout << "* Rho name = " << fJetsCont->GetRhoName() << std::endl;
1609  std::cout << "* Rho mass name = " << fJetsCont->GetRhoMassName() << std::endl;
1610  std::cout << std::endl;
1611 
1612  std::cout << "### Tree extraction properties ###" << std::endl;
1613  std::vector<Float_t> extractionPtBins = fJetTree->GetExtractionPercentagePtBins();
1614  std::vector<Float_t> extractionPercentages = fJetTree->GetExtractionPercentages();
1615  std::vector<Int_t> extractionHM = fJetTree->GetExtractionJetTypes_HM();
1616  std::vector<Int_t> extractionPM = fJetTree->GetExtractionJetTypes_PM();
1617  if(extractionPercentages.size())
1618  {
1619  std::cout << "* Extraction bins: (";
1620  for(Int_t i=0; i<static_cast<Int_t>(extractionPercentages.size()-1); i++)
1621  std::cout << extractionPtBins[i*2] << "-" << extractionPtBins[i*2+1] << " = " << extractionPercentages[i] << ", ";
1622  std::cout << extractionPtBins[(extractionPercentages.size()-1)*2] << "-" << extractionPtBins[(extractionPercentages.size()-1)*2+1] << " = " << extractionPercentages[(extractionPercentages.size()-1)];
1623 
1624  std::cout << ")" << std::endl;
1625  }
1626  if(extractionHM.size())
1627  {
1628  std::cout << "* Extraction of hadron-matched jets with types: (";
1629  for(Int_t i=0; i<static_cast<Int_t>(extractionHM.size()-1); i++)
1630  std::cout << extractionHM[i] << ", ";
1631  std::cout << extractionHM[extractionHM.size()-1];
1632  std::cout << ")" << std::endl;
1633  }
1634  if(extractionPM.size())
1635  {
1636  std::cout << "* Extraction of parton-matched jets with types: (";
1637  for(Int_t i=0; i<static_cast<Int_t>(extractionPM.size()-1); i++)
1638  std::cout << extractionPM[i] << ", ";
1639  std::cout << extractionPM[extractionPM.size()-1];
1640  std::cout << ")" << std::endl;
1641  }
1642  std::cout << std::endl;
1643 
1644  std::cout << "### Extracted data groups ###" << std::endl;
1646  std::cout << "* Basic event properties (ID, vertex, centrality, bgrd. densities, ...)" << std::endl;
1648  std::cout << "* Jet constituents, basic properties (pt, eta, phi, charge, ...)" << std::endl;
1650  std::cout << "* Jet constituents, IPs" << std::endl;
1652  std::cout << "* Jet constituents, PID" << std::endl;
1654  std::cout << "* Jet calorimeter clusters" << std::endl;
1656  std::cout << "* MC information (origin, matched jets, ...)" << std::endl;
1658  std::cout << "* Secondary vertices" << std::endl;
1659  if(fJetTree->GetSaveJetShapes())
1660  std::cout << "* Jet shapes (jet mass, LeSub, pTD, ...)" << std::endl;
1662  std::cout << "* Trigger tracks" << std::endl;
1663  std::cout << std::endl;
1664  std::cout << "### Further settings ###" << std::endl;
1665  if(fTruthJetsArrayName != "")
1666  std::cout << Form("* True jet matching active, array=%s, rho=%s, rho_mass=%s", fTruthJetsArrayName.Data(), fTruthJetsRhoName.Data(), fTruthJetsRhoMassName.Data()) << std::endl;
1668  std::cout << Form("* Particle level information available (for jet origin calculation, particle code): %s", fTruthParticleArrayName.Data()) << std::endl;
1669  if(extractionHM.size())
1670  std::cout << Form("* Hadron--jet matching radius=%3.3f", fHadronMatchingRadius) << std::endl;
1672  std::cout << Form("* True particle label range for pt fraction=%d-%d", fTruthMinLabel, fTruthMaxLabel) << std::endl;
1674  std::cout << Form("* Particle PID codes will be PDG codes") << std::endl;
1676  std::cout << Form("* Particle PID codes will follow AliAODPid convention") << std::endl;
1678  std::cout << Form("* AliEmcalJet flavour tag is set from HF-jet tagging") << std::endl;
1679 
1680  std::cout << "#########################################\n\n";
1681 }
1682 
1683 
1684 //########################################################################
1685 // HELPERS
1686 //########################################################################
1687 
1688 //________________________________________________________________________
1689 inline Bool_t AliAnalysisTaskJetExtractor::IsTrackInCone(AliVParticle* track, Double_t eta, Double_t phi, Double_t radius)
1690 {
1691  // This is to use a full cone in phi even at the edges of phi (2pi -> 0) (0 -> 2pi)
1692  Double_t trackPhi = 0.0;
1693  if (track->Phi() > (TMath::TwoPi() - (radius-phi)))
1694  trackPhi = track->Phi() - TMath::TwoPi();
1695  else if (track->Phi() < (phi+radius - TMath::TwoPi()))
1696  trackPhi = track->Phi() + TMath::TwoPi();
1697  else
1698  trackPhi = track->Phi();
1699 
1700  if ( TMath::Abs(trackPhi-phi)*TMath::Abs(trackPhi-phi) + TMath::Abs(track->Eta()-eta)*TMath::Abs(track->Eta()-eta) <= radius*radius)
1701  return kTRUE;
1702 
1703  return kFALSE;
1704 }
1705 
1706 //________________________________________________________________________
1708 {
1709  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1710  if(!tmpHist)
1711  {
1712  AliError(Form("Cannot find histogram <%s> ",key)) ;
1713  return;
1714  }
1715 
1716  tmpHist->Fill(x);
1717 }
1718 
1719 //________________________________________________________________________
1721 {
1722  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1723  if(!tmpHist)
1724  {
1725  AliError(Form("Cannot find histogram <%s> ",key));
1726  return;
1727  }
1728 
1729  if (tmpHist->IsA()->GetBaseClass("TH1"))
1730  static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
1731  else if (tmpHist->IsA()->GetBaseClass("TH2"))
1732  static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
1733 }
1734 
1735 //________________________________________________________________________
1737 {
1738  TH2* tmpHist = static_cast<TH2*>(fOutput->FindObject(key));
1739  if(!tmpHist)
1740  {
1741  AliError(Form("Cannot find histogram <%s> ",key));
1742  return;
1743  }
1744 
1745  tmpHist->Fill(x,y,add);
1746 }
1747 
1748 //________________________________________________________________________
1750 {
1751  TH3* tmpHist = static_cast<TH3*>(fOutput->FindObject(key));
1752  if(!tmpHist)
1753  {
1754  AliError(Form("Cannot find histogram <%s> ",key));
1755  return;
1756  }
1757 
1758  if(add)
1759  tmpHist->Fill(x,y,z,add);
1760  else
1761  tmpHist->Fill(x,y,z);
1762 }
1763 
1764 
1765 //________________________________________________________________________
1766 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)
1767 {
1768  T* tmpHist = new T(name, title, xBins, xMin, xMax);
1769 
1770  tmpHist->GetXaxis()->SetTitle(xTitle);
1771  tmpHist->GetYaxis()->SetTitle(yTitle);
1772  tmpHist->SetOption(options);
1773  tmpHist->SetMarkerStyle(kFullCircle);
1774  tmpHist->Sumw2();
1775 
1776  fOutput->Add(tmpHist);
1777 
1778  return tmpHist;
1779 }
1780 
1781 //________________________________________________________________________
1782 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)
1783 {
1784  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax);
1785  tmpHist->GetXaxis()->SetTitle(xTitle);
1786  tmpHist->GetYaxis()->SetTitle(yTitle);
1787  tmpHist->GetZaxis()->SetTitle(zTitle);
1788  tmpHist->SetOption(options);
1789  tmpHist->SetMarkerStyle(kFullCircle);
1790  tmpHist->Sumw2();
1791 
1792  fOutput->Add(tmpHist);
1793 
1794  return tmpHist;
1795 }
1796 
1797 //________________________________________________________________________
1798 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)
1799 {
1800  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax, zBins, zMin, zMax);
1801  tmpHist->GetXaxis()->SetTitle(xTitle);
1802  tmpHist->GetYaxis()->SetTitle(yTitle);
1803  tmpHist->GetZaxis()->SetTitle(zTitle);
1804  tmpHist->SetOption(options);
1805  tmpHist->SetMarkerStyle(kFullCircle);
1806  tmpHist->Sumw2();
1807 
1808  fOutput->Add(tmpHist);
1809 
1810  return tmpHist;
1811 }
1812 
1813 //________________________________________________________________________
1815 {
1816  // Called once at the end of the analysis.
1817 }
1818 
1819 // ### ADDTASK MACRO
1820 //________________________________________________________________________
1821 AliAnalysisTaskJetExtractor* AliAnalysisTaskJetExtractor::AddTaskJetExtractor(TString trackArray, TString clusterArray, TString jetArray, TString rhoObject, Double_t jetRadius, TString configFile, AliRDHFJetsCutsVertex* vertexerCuts, const char* taskNameSuffix)
1822 {
1823  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1824  Double_t minJetEta = 0.5;
1825  Double_t minJetPt = 0.15;
1826  Double_t minTrackPt = 0.15;
1827  Double_t minJetAreaPerc = 0.557;
1828  TString suffix = "";
1829  if(taskNameSuffix != 0)
1830  suffix = taskNameSuffix;
1831 
1832  // ###### Load configuration from YAML files
1833  if(gGrid)
1834  {
1836  fYAMLConfig.AddConfiguration("alien:///alice/cern.ch/user/r/rhaake/ConfigFiles/JetExtractor_BaseConfig.yaml", "baseConfig");
1837  if(configFile != "")
1838  fYAMLConfig.AddConfiguration(configFile.Data(), "customConfig");
1839  fYAMLConfig.Initialize();
1840 
1841  fYAMLConfig.GetProperty("TrackArray", trackArray, kFALSE);
1842  fYAMLConfig.GetProperty("ClusterArray", clusterArray, kFALSE);
1843  fYAMLConfig.GetProperty("JetArray", jetArray, kFALSE);
1844  fYAMLConfig.GetProperty("RhoName", rhoObject, kFALSE);
1845  fYAMLConfig.GetProperty("JetRadius", jetRadius);
1846  fYAMLConfig.GetProperty("MinJetEta", minJetEta);
1847  fYAMLConfig.GetProperty("MinJetPt", minJetPt);
1848  fYAMLConfig.GetProperty("MinTrackPt", minTrackPt);
1849  fYAMLConfig.GetProperty("MinJetAreaPerc", minJetAreaPerc);
1850  fYAMLConfig.Print();
1851  }
1852 
1853  // ###### Task name
1854  TString name("AliAnalysisTaskJetExtractor");
1855  if (jetArray != "") {
1856  name += "_";
1857  name += jetArray;
1858  }
1859  if (rhoObject != "") {
1860  name += "_";
1861  name += rhoObject;
1862  }
1863  if (suffix != "") {
1864  name += "_";
1865  name += suffix;
1866  }
1867 
1868  // ###### Setup task with default settings
1870  myTask->SetNeedEmcalGeom(kFALSE);
1871  myTask->SetVzRange(-10.,10.);
1872 
1873  // Particle container and track pt cut
1874  AliParticleContainer* trackCont = 0;
1875  if(trackArray == "mcparticles")
1876  trackCont = myTask->AddMCParticleContainer(trackArray);
1877  else if(trackArray =="mctracks")
1878  trackCont = myTask->AddParticleContainer(trackArray);
1879  else
1880  trackCont = myTask->AddTrackContainer(trackArray);
1881  trackCont->SetParticlePtCut(minTrackPt);
1882 
1883  // Particle container and track pt cut
1884  AliClusterContainer* clusterCont = 0;
1885  if(clusterArray != "")
1886  clusterCont = myTask->AddClusterContainer(clusterArray);
1887 
1888  // Secondary vertex cuts (default settings from PWGHF)
1889  // (can be overwritten by using myTask->SetVertexerCuts(cuts) from outside macro)
1890  AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts", "default");
1891  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1892  esdTrackCuts->SetMinNClustersTPC(90);
1893  esdTrackCuts->SetMaxChi2PerClusterTPC(4);
1894  esdTrackCuts->SetRequireTPCRefit(kTRUE);
1895  esdTrackCuts->SetRequireITSRefit(kTRUE);
1896  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
1897  esdTrackCuts->SetMinDCAToVertexXY(0.);
1898  esdTrackCuts->SetEtaRange(-0.8, 0.8);
1899  esdTrackCuts->SetPtRange(1.0, 1.e10);
1900  vertexerCuts->AddTrackCuts(esdTrackCuts);
1901  vertexerCuts->SetNprongs(3);
1902  vertexerCuts->SetMinPtHardestTrack(1.0);//default 0.3
1903  vertexerCuts->SetSecVtxWithKF(kFALSE);//default with StrLinMinDist
1904  vertexerCuts->SetImpParCut(0.);//default 0
1905  vertexerCuts->SetDistPrimSec(0.);//default 0
1906  vertexerCuts->SetCospCut(-1);//default -1
1907  myTask->SetVertexerCuts(vertexerCuts);
1908 
1909  // Jet container
1910  AliJetContainer *jetCont = myTask->AddJetContainer(jetArray,6,jetRadius);
1911  if (jetCont) {
1912  jetCont->SetRhoName(rhoObject);
1913  jetCont->SetPercAreaCut(minJetAreaPerc);
1914  jetCont->SetJetPtCut(minJetPt);
1915  jetCont->SetLeadingHadronType(0);
1916  jetCont->SetPtBiasJetTrack(minTrackPt);
1917  jetCont->SetJetEtaLimits(-minJetEta, +minJetEta);
1918  jetCont->ConnectParticleContainer(trackCont);
1919  if(clusterCont)
1920  jetCont->ConnectClusterContainer(clusterCont);
1921  jetCont->SetMaxTrackPt(1000);
1922  }
1923 
1924  mgr->AddTask(myTask);
1925 
1926  // ###### Connect inputs/outputs
1927  mgr->ConnectInput (myTask, 0, mgr->GetCommonInputContainer() );
1928  mgr->ConnectOutput (myTask, 1, mgr->CreateContainer(Form("%s_histos", name.Data()), AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, Form("%s:ChargedJetsHadronCF", mgr->GetCommonFileName())) );
1929  mgr->ConnectOutput (myTask, 2, mgr->CreateContainer(Form("%s_tree", name.Data()), TTree::Class(), AliAnalysisManager::kOutputContainer, mgr->GetCommonFileName()) );
1930 
1931  return myTask;
1932 }
Float_t fBuffer_Shape_MomentumDispersion
! array buffer
Bool_t fSaveConstituentsIP
save arrays of constituent impact parameters
void ExecOnce()
Perform steps needed to initialize the analysis.
Int_t pdg
Bool_t fSaveSecondaryVertices
save reconstructed sec. vertex properties
Double_t GetFirstOrderSubtractedAngularity() const
void AddFlavourTag(Int_t tag)
Definition: AliEmcalJet.h:365
Double_t Area() const
Definition: AliEmcalJet.h:130
void SetParticlePtCut(Double_t cut)
virtual AliVParticle * GetNextAcceptParticle()
void SetSaveConstituents(Bool_t val)
Double_t GetRhoVal() const
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
TString fTruthParticleArrayName
Array name of MC particles in event (mcparticles)
Double_t GetSecondOrderSubtractedSigma2() const
Float_t fBuffer_Shape_LeSub_DerivCorr
! array buffer
Double_t GetJetEtaMin() const
const char * title
Definition: MakeQAPdf.C:27
void SetTrackContainer(AliParticleContainer *cont)
void GetJetType(AliEmcalJet *jet, Int_t &typeHM, Int_t &typePM, Int_t &typeIC)
Bool_t fSaveCaloClusters
save calorimeter clusters
Float_t GetPartonEta6() const
void GetTrackImpactParameters(const AliVVertex *vtx, AliAODTrack *track, Float_t &d0, Float_t &d0cov, Float_t &z0, Float_t &z0cov)
Bool_t fSetEmcalJetFlavour
if set, the flavour property of the AliEmcalJets will be set
AliJetContainer * GetJetContainer(Int_t i=0) const
TString fBackgroundModelFileName
MVA model file name.
Definition: External.C:244
Float_t GetPartonEta7() const
Bool_t fBackgroundModelPtOnly
MVA model, use approximation for pT only (not for mass)
TTree * fJetTree
! tree structure
Double_t GetSecondOrderSubtractedConstituent() const
Double_t Eta() const
Definition: AliEmcalJet.h:121
long long Long64_t
Definition: External.C:43
T * AddHistogram3D(const char *name="CustomHistogram", const char *title="NO_TITLE", const char *options="", Int_t xBins=100, Double_t xMin=0.0, Double_t xMax=20.0, Int_t yBins=100, Double_t yMin=0.0, Double_t yMax=20.0, Int_t zBins=100, Double_t zMin=0.0, Double_t zMax=20.0, const char *xTitle="x axis", const char *yTitle="y axis", const char *zTitle="z axis")
Float_t fBuffer_Event_BackgroundDensity
! array buffer
void SetLeadingHadronType(Int_t t)
Float_t fBuffer_Shape_Mass_DerivCorr_1
! array buffer
Double_t Phi() const
Definition: AliEmcalJet.h:117
std::vector< Float_t > GetExtractionPercentages()
Bool_t fSaveMCInformation
save MC information
Double_t mass
void SetJetShapesInBuffer(Double_t leSub_noCorr, Double_t radialMoment, Double_t momentumDispersion, Double_t constPtMean, Double_t constPtMedian)
void SetPtBiasJetTrack(Float_t b)
TSystem * gSystem
centrality
Float_t * fBuffer_Cluster_M02
! array buffer
AliParticleContainer * fTrackContainer
! track container, used to work on tracks
Int_t ClusterAt(Int_t idx) const
Definition: AliEmcalJet.h:137
Int_t fTruthMinLabel
min track label to consider it as true particle
Double_t fSecondaryVertexMaxChi2
Max chi2 of secondary vertex (others will be discarded)
Double_t GetJetEtaMax() const
AliClusterContainer * GetClusterContainer() const
void SetArea(Double_t a)
Definition: AliEmcalJet.h:256
Float_t GetPartonPhi7() const
Double_t fSecondaryVertexMaxDispersion
Max dispersion of secondary vertex (others will be discarded)
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
Double_t fEventCut_TriggerTrackMinPt
Event requirement, trigger track min pT.
std::vector< Float_t > GetExtractionPercentagePtBins()
Bool_t fSaveEventProperties
save general event properties (bgrd. etc.)
std::vector< Float_t > fTriggerTracks_Pt
If trigger track requi. is used -> save pT.
Float_t fBuffer_Shape_Circularity_DerivCorr_1
! array buffer
bool GetProperty(std::vector< std::string > propertyPath, const std::string &propertyName, T &property, const bool requiredProperty) const
Bool_t fSaveJetShapes
save jet shapes
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 (only works for JEWEL)
AliStack * stack
void SetVzRange(Double_t min, Double_t max)
Set pre-configured event cut object.
const Int_t kMaxNumConstituents
Bool_t IsTrackInCone(AliVParticle *track, Double_t eta, Double_t phi, Double_t radius)
Double_t GetTrueJetPtFraction(AliEmcalJet *jet)
void SetLabel(Int_t l)
Definition: AliEmcalJet.h:255
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.
TString fTruthJetsRhoName
Array name for particle-level rho.
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
Double_t fTrueJetMatchingRadius
Matching radius to true jet.
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
std::vector< Float_t > fTriggerTracks_Phi
If trigger track requi. is used -> save phi.
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
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 embedded particles) ...
Bool_t AddJetToTree(AliEmcalJet *jet, Float_t vertexX=0, Float_t vertexY=0, Float_t vertexZ=0, Float_t centrality=0, Long64_t eventID=0, Float_t magField=0, Int_t motherParton=0, Int_t motherHadron=0, Int_t partonInitialCollision=0, Float_t matchDistance=0, Float_t matchedPt=0, Float_t matchedMass=0, Float_t truePtFraction=0, Float_t ptHard=0, Float_t eventWeight=0, Float_t impactParameter=0, std::vector< Float_t > &trackPID_ITS=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &trackPID_TPC=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &trackPID_TOF=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &trackPID_TRD=DEFAULT_VECTOR_FLOAT, std::vector< Short_t > &trackPID_Reco=DEFAULT_VECTOR_SHORT, std::vector< Int_t > &trackPID_Truth=DEFAULT_VECTOR_INT, std::vector< Float_t > &triggerTrackPt=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &triggerTrackDeltaEta=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &triggerTrackDeltaPhi=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &trackIP_d0=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &trackIP_z0=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &trackIP_d0cov=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &trackIP_z0cov=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_X=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_Y=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_Z=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_Mass=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_Lxy=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_SigmaLxy=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_Chi2=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_Dispersion=DEFAULT_VECTOR_FLOAT)
Float_t * fBuffer_Const_Phi
! array buffer
Float_t fBuffer_Jet_MC_MatchDistance
! array buffer
Float_t fBuffer_Shape_Mass_DerivCorr_2
! array buffer
Float_t fBuffer_Event_ImpactParameter
! array buffer
const TString & GetRhoMassName() const
TRandom3 * fRandomGenerator
! Random number generator, used for event + jet efficiency
Int_t GetPartonFlag6() const
void SetRhoName(const char *n)
Float_t fBuffer_Jet_MC_MatchedPt
! array buffer
Bool_t fSaveTriggerTracks
save event trigger track
AliParticleContainer * GetParticleContainer() const
TRandom3 * fRandomGenerator
! random generator
void SetJetAcceptanceType(UInt_t type)
Definition: AliEmcalJet.h:380
AliParticleContainer * fTracksCont
! Tracks
std::vector< Float_t > fExtractionPercentages
Percentages which will be extracted for a given pT bin.
TClonesArray * fTruthParticleArray
! Array of MC particles in event (mcparticles)
Float_t * fBuffer_Cluster_Time
! array buffer
Double_t fEventPercentage
percentage (0, 1] which will be extracted
TString GetBackgroundModelArrayString(AliEmcalJet *jet)
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
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:138
void SetJetPtCut(Float_t cut)
void SetSaveJetShapes(Bool_t val)
TClonesArray * fJetOutputArray
! Array of corr. jets, attached to event
unsigned int UInt_t
Definition: External.C:33
TString fBackgroundModelInputParameters
MVA model input parameters (comma-separated)
Float_t * fBuffer_Const_Charge
! array buffer
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")
void GetPtAndMassFromModel(AliEmcalJet *jet, Float_t &pt_ML, Float_t &mass_ML)
TRandom3 * fRandomGeneratorCones
! Random number generator, used for random cones
Int_t * fBuffer_Const_Label
! array buffer
Float_t fBuffer_Shape_Sigma2_DerivCorr_1
! array buffer
AliParticleContainer * AddParticleContainer(const char *n)
Create new particle container and attach it to the task.
Double_t fImpactParameter
IP for each event (only works for JEWEL)
Float_t fBuffer_JetPhi
! array buffer
Double_t GetSecondOrderSubtractedAngularity() const
Int_t fTruthMaxLabel
max track label to consider it as true particle
Double_t GetMatchedTrueJetObservables(AliEmcalJet *jet, Double_t &matchedJetPt, Double_t &matchedJetMass)
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
Long64_t fBuffer_Event_ID
! array buffer
Float_t fBuffer_Event_Centrality
! array buffer
Int_t fBuffer_Jet_MC_MotherIC
! array buffer
std::ostream & Print(std::ostream &in, const int index=-1) const
TString fCustomStartupScript
Path to custom shell script that will be executed.
Double_t GetRhoMassVal() const
Float_t fBuffer_Shape_Mass_NoCorr
! array buffer
Float_t * fBuffer_Const_ProdVtx_X
! array buffer
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.
void AddPIDInformation(AliVParticle *particle, Float_t &sigITS, Float_t &sigTPC, Float_t &sigTOF, Float_t &sigTRD, Short_t &recoPID, Int_t &truePID)
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()
AliEmcalJet * GetNextAcceptJet()
Float_t fBuffer_Event_BackgroundDensityMass
! array buffer
AliHFJetsTaggingVertex * fVtxTagger
! class for sec. vertexing
short Short_t
Definition: External.C:23
TString fTruthJetsArrayName
Array name for particle-level jets.
Bool_t fCalculateModelBackground
Calculate MVA model background and attach to event.
void ConnectParticleContainer(AliParticleContainer *c)
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliRDHFJetsCutsVertex * fVertexerCuts
Cuts used for the vertexer (given in add task macro)
static AliAnalysisTaskJetExtractor * AddTaskJetExtractor(TString trackArray, TString clusterArray, TString jetArray, TString rhoObject, Double_t jetRadius, TString configFile, AliRDHFJetsCutsVertex *vertexerCuts, const char *taskNameSuffix)
Double_t GetFirstOrderSubtractedCircularity() const
Float_t * fBuffer_Cluster_Eta
! array buffer
Int_t fEventCut_TriggerTrackMinLabel
Event requirement, trigger track min label (can be used to selected embedded particles) ...
AliJetContainer * fJetContainer
! jet container, used to work on jets
Int_t fBuffer_NumClusters
! array buffer
void SetSaveEventProperties(Bool_t val)
Float_t fPtHard
!event -hard
void FillHistogram(const char *key, Double_t x)
void FillTrackControlHistograms(AliVTrack *track)
Float_t GetJetRadius() const
AliEmcalList * fOutput
!output list
TString fTruthJetsRhoMassName
Array name for particle-level rho for mass.
Float_t fBuffer_Jet_MC_MatchedMass
! array buffer
Float_t fBuffer_JetPt
! array buffer
Definition: External.C:220
Float_t fBuffer_Shape_ConstPtMedian
! array buffer
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
AliTrackContainer * AddTrackContainer(const char *n)
Create new track container and attach it to the task.
void SetClusterContainer(AliClusterContainer *cont)
Bool_t fCalculateSecondaryVertices
Calculate the secondary vertices (instead of loading)
Double_t GetSecondOrderSubtractedCircularity() const
void CalculateJetShapes(AliEmcalJet *jet, Double_t &leSub_noCorr, Double_t &radialMoment, Double_t &momentumDispersion, Double_t &constPtMean, Double_t &constPtMedian)
Float_t fBuffer_Shape_NumConst_DerivCorr
! array buffer
Float_t * fBuffer_Const_ProdVtx_Y
! array buffer
Bool_t fSaveConstituents
save arrays of constituent basic properties
YAML configuration class for AliPhysics.
Float_t fBuffer_Shape_RadialMoment
! array buffer
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)
Float_t * fBuffer_Cluster_E
! array buffer
void SetNeedEmcalGeom(Bool_t n)
Double_t fHadronMatchingRadius
Matching radius to search for beauty/charm hadrons around jet.
Base task in the EMCAL jet framework.
Float_t fBuffer_Shape_Angularity_DerivCorr_2
! array buffer
Float_t * fBuffer_Const_Pt
! 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
AliClusterContainer * fClustersCont
! Clusters
Float_t * fBuffer_Cluster_Pt
! array buffer
Float_t * fBuffer_Const_Eta
! array buffer
Float_t fBuffer_JetEta
! array buffer
int AddConfiguration(std::string configurationFilename, std::string configurationName="")
Int_t * fBuffer_Cluster_Label
! array buffer
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
Int_t fBuffer_NumConstituents
! array buffer
void UserCreateOutputObjects()
Main initialization function on the worker.
Double_t GetFirstOrderSubtractedSigma2() const
Float_t fBuffer_Event_Vertex_X
! array buffer
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.
void SetJetContainer(AliJetContainer *cont)
Container structure for EMCAL clusters.
std::vector< Float_t > fTriggerTracks_Eta
If trigger track requi. is used -> save eta.
Double_t M() const
Definition: AliEmcalJet.h:120
Float_t fBuffer_Shape_ConstPtMean
! array buffer
Float_t fBuffer_Shape_LeSub_NoCorr
! array buffer
Float_t * fBuffer_Const_ProdVtx_Z
! 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_pTD_DerivCorr_1
! array buffer
std::vector< int > GetPtSortedTrackConstituentIndexes(TClonesArray *tracks) const
Double_t yMax
AliClusterContainer * fClusterContainer
! cluster container, used to work on clusters
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
Float_t fBuffer_Event_Vertex_Z
! array buffer
void SetSaveTriggerTracks(Bool_t val)
Float_t fBuffer_Event_MagneticField
! array buffer