AliPhysics  a3be53f (a3be53f)
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
265  // Jet mass
266  fBuffer_Shape_Mass_NoCorr = jet->M();
269 
270  // Set Monte Carlo information
272  {
273  fBuffer_Jet_MC_MotherParton = motherParton;
274  fBuffer_Jet_MC_MotherHadron = motherHadron;
275  fBuffer_Jet_MC_MotherIC = partonInitialCollision;
276  fBuffer_Jet_MC_MatchDistance = matchDistance;
277  fBuffer_Jet_MC_MatchedPt = matchedPt;
278  fBuffer_Jet_MC_MatchedMass = matchedMass;
279  fBuffer_Jet_MC_TruePtFraction = truePtFraction;
280  }
281 
283  {
284  fBuffer_NumTriggerTracks = numTriggerTracks;
285  fJetTree->SetBranchAddress("Jet_TriggerTrack_Pt", triggerTrackPt);
286  fJetTree->SetBranchAddress("Jet_TriggerTrack_dEta", triggerTrackDeltaEta);
287  fJetTree->SetBranchAddress("Jet_TriggerTrack_dPhi", triggerTrackDeltaPhi);
288  }
289 
290  // Add buffered jet to tree
291  fJetTree->Fill();
292 
293  return kTRUE;
294 }
295 
296 
297 //________________________________________________________________________
299 {
300  // Create the tree with active branches
301 
302  void* dummy = 0; // for some branches, do not need a buffer pointer yet
303 
304 
305  if(fInitialized)
306  AliFatal("Tree is already initialized.");
307 
308  if(!fJetContainer)
309  AliFatal("fJetContainer not set in tree");
310 
311  // ### Prepare the jet tree
312  fJetTree = new TTree(Form("JetTree_%s", GetName()), "");
313 
314  fJetTree->Branch("Jet_Pt",&fBuffer_JetPt,"Jet_Pt/F");
315  fJetTree->Branch("Jet_Phi",&fBuffer_JetPhi,"Jet_Phi/F");
316  fJetTree->Branch("Jet_Eta",&fBuffer_JetEta,"Jet_Eta/F");
317  fJetTree->Branch("Jet_Area",&fBuffer_JetArea,"Jet_Area/F");
318  fJetTree->Branch("Jet_NumConstituents",&fBuffer_NumConstituents,"Jet_NumConstituents/I");
320  fJetTree->Branch("Jet_NumClusters",&fBuffer_NumClusters,"Jet_NumClusters/I");
321 
323  {
324  fJetTree->Branch("Event_BackgroundDensity",&fBuffer_Event_BackgroundDensity,"Event_BackgroundDensity/F");
325  fJetTree->Branch("Event_BackgroundDensityMass",&fBuffer_Event_BackgroundDensityMass,"Event_BackgroundDensityMass/F");
326  fJetTree->Branch("Event_Vertex_X",&fBuffer_Event_Vertex_X,"Event_Vertex_X/F");
327  fJetTree->Branch("Event_Vertex_Y",&fBuffer_Event_Vertex_Y,"Event_Vertex_Y/F");
328  fJetTree->Branch("Event_Vertex_Z",&fBuffer_Event_Vertex_Z,"Event_Vertex_Z/F");
329  fJetTree->Branch("Event_Centrality",&fBuffer_Event_Centrality,"Event_Centrality/F");
330  fJetTree->Branch("Event_ID",&fBuffer_Event_ID,"Event_ID/L");
331  fJetTree->Branch("Event_MagneticField",&fBuffer_Event_MagneticField,"Event_MagneticField/F");
332 
334  {
335  fJetTree->Branch("Event_PtHard",&fBuffer_Event_PtHard,"Event_PtHard/F");
336  fJetTree->Branch("Event_Weight",&fBuffer_Event_Weight,"Event_Weight/F");
337  fJetTree->Branch("Event_ImpactParameter",&fBuffer_Event_ImpactParameter,"Event_ImpactParameter/F");
338  }
339  }
340 
342  {
343  fJetTree->Branch("Jet_Const_Pt",fBuffer_Const_Pt,"Jet_Const_Pt[Jet_NumConstituents]/F");
344  fJetTree->Branch("Jet_Const_Phi",fBuffer_Const_Phi,"Jet_Const_Phi[Jet_NumConstituents]/F");
345  fJetTree->Branch("Jet_Const_Eta",fBuffer_Const_Eta,"Jet_Const_Eta[Jet_NumConstituents]/F");
346  fJetTree->Branch("Jet_Const_Charge",fBuffer_Const_Charge,"Jet_Const_Charge[Jet_NumConstituents]/F");
348  fJetTree->Branch("Jet_Const_Label",fBuffer_Const_Label,"Jet_Const_Label[Jet_NumConstituents]/I");
349  }
350 
352  {
353  fJetTree->Branch("Jet_Cluster_Pt",fBuffer_Cluster_Pt,"Jet_Cluster_Pt[Jet_NumClusters]/F");
354  fJetTree->Branch("Jet_Cluster_E",fBuffer_Cluster_E,"Jet_Cluster_Pt[Jet_NumClusters]/F");
355  fJetTree->Branch("Jet_Cluster_Phi",fBuffer_Cluster_Phi,"Jet_Cluster_Phi[Jet_NumClusters]/F");
356  fJetTree->Branch("Jet_Cluster_Eta",fBuffer_Cluster_Eta,"Jet_Cluster_Eta[Jet_NumClusters]/F");
357  fJetTree->Branch("Jet_Cluster_M02",fBuffer_Cluster_M02,"Jet_Cluster_M02[Jet_NumClusters]/F");
358  fJetTree->Branch("Jet_Cluster_Time",fBuffer_Cluster_Time,"Jet_Cluster_Time[Jet_NumClusters]/F");
360  fJetTree->Branch("Jet_Cluster_Label",fBuffer_Cluster_Label,"Jet_Cluster_Label[Jet_NumClusters]/I");
361  }
362 
364  {
365  fJetTree->Branch("Jet_Const_IPd",&dummy,"Jet_Const_IPd[Jet_NumConstituents]/F");
366  fJetTree->Branch("Jet_Const_IPz",&dummy,"Jet_Const_IPz[Jet_NumConstituents]/F");
367  fJetTree->Branch("Jet_Const_CovIPd",&dummy,"Jet_Const_CovIPd[Jet_NumConstituents]/F");
368  fJetTree->Branch("Jet_Const_CovIPz",&dummy,"Jet_Const_CovIPz[Jet_NumConstituents]/F");
369 
370  fJetTree->Branch("Jet_Const_ProdVtx_X",fBuffer_Const_ProdVtx_X,"Jet_Const_ProdVtx_X[Jet_NumConstituents]/F");
371  fJetTree->Branch("Jet_Const_ProdVtx_Y",fBuffer_Const_ProdVtx_Y,"Jet_Const_ProdVtx_Y[Jet_NumConstituents]/F");
372  fJetTree->Branch("Jet_Const_ProdVtx_Z",fBuffer_Const_ProdVtx_Z,"Jet_Const_ProdVtx_Z[Jet_NumConstituents]/F");
373  }
374 
376  {
377  fJetTree->Branch("Jet_Const_PID_ITS",&dummy,"Jet_Const_PID_ITS[Jet_NumConstituents]/F");
378  fJetTree->Branch("Jet_Const_PID_TPC",&dummy,"Jet_Const_PID_TPC[Jet_NumConstituents]/F");
379  fJetTree->Branch("Jet_Const_PID_TOF",&dummy,"Jet_Const_PID_TOF[Jet_NumConstituents]/F");
380  fJetTree->Branch("Jet_Const_PID_TRD",&dummy,"Jet_Const_PID_TRD[Jet_NumConstituents]/F");
381 
382  fJetTree->Branch("Jet_Const_PID_Reconstructed",&dummy,"Jet_Const_PID_Reconstructed[Jet_NumConstituents]/S");
384  fJetTree->Branch("Jet_Const_PID_Truth",&dummy,"Jet_Const_PID_Truth[Jet_NumConstituents]/I");
385  }
386 
387  if(fSaveJetShapes)
388  {
389  fJetTree->Branch("Jet_Shape_Mass_NoCorr",&fBuffer_Shape_Mass_NoCorr,"Jet_Shape_Mass_NoCorr/F");
390  fJetTree->Branch("Jet_Shape_Mass_DerivCorr_1",&fBuffer_Shape_Mass_DerivCorr_1,"Jet_Shape_Mass_DerivCorr_1/F");
391  fJetTree->Branch("Jet_Shape_Mass_DerivCorr_2",&fBuffer_Shape_Mass_DerivCorr_2,"Jet_Shape_Mass_DerivCorr_2/F");
392  }
393 
395  {
396  fJetTree->Branch("Jet_MC_MotherParton",&fBuffer_Jet_MC_MotherParton,"Jet_MC_MotherParton/I");
397  fJetTree->Branch("Jet_MC_MotherHadron",&fBuffer_Jet_MC_MotherHadron,"Jet_MC_MotherHadron/I");
398  fJetTree->Branch("Jet_MC_MotherIC",&fBuffer_Jet_MC_MotherIC,"Jet_MC_MotherIC/I");
399  fJetTree->Branch("Jet_MC_MatchDistance",&fBuffer_Jet_MC_MatchDistance,"Jet_MC_MatchDistance/F");
400  fJetTree->Branch("Jet_MC_MatchedPt",&fBuffer_Jet_MC_MatchedPt,"Jet_MC_MatchedPt/F");
401  fJetTree->Branch("Jet_MC_MatchedMass",&fBuffer_Jet_MC_MatchedMass,"Jet_MC_MatchedMass/F");
402  fJetTree->Branch("Jet_MC_TruePtFraction",&fBuffer_Jet_MC_TruePtFraction,"Jet_MC_TruePtFraction/F");
403  }
404 
406  {
407  fJetTree->Branch("Jet_NumSecVertices",&fBuffer_NumSecVertices,"Jet_NumSecVertices/I");
408 
409  fJetTree->Branch("Jet_SecVtx_X",&dummy,"Jet_SecVtx_X[Jet_NumSecVertices]/F");
410  fJetTree->Branch("Jet_SecVtx_Y",&dummy,"Jet_SecVtx_Y[Jet_NumSecVertices]/F");
411  fJetTree->Branch("Jet_SecVtx_Z",&dummy,"Jet_SecVtx_Z[Jet_NumSecVertices]/F");
412  fJetTree->Branch("Jet_SecVtx_Mass",&dummy,"Jet_SecVtx_Mass[Jet_NumSecVertices]/F");
413  fJetTree->Branch("Jet_SecVtx_Lxy",&dummy,"Jet_SecVtx_Lxy[Jet_NumSecVertices]/F");
414  fJetTree->Branch("Jet_SecVtx_SigmaLxy",&dummy,"Jet_SecVtx_SigmaLxy[Jet_NumSecVertices]/F");
415  fJetTree->Branch("Jet_SecVtx_Chi2",&dummy,"Jet_SecVtx_Chi2[Jet_NumSecVertices]/F");
416  fJetTree->Branch("Jet_SecVtx_Dispersion",&dummy,"Jet_SecVtx_Dispersion[Jet_NumSecVertices]/F");
417  }
418 
419  // Trigger track requirement active -> Save trigger track
421  {
422  fJetTree->Branch("Jet_NumTriggerTracks",&fBuffer_NumTriggerTracks,"Jet_NumTriggerTracks/I");
423  fJetTree->Branch("Jet_TriggerTrack_Pt",&dummy,"Jet_TriggerTrack_Pt[Jet_NumTriggerTracks]/F");
424  fJetTree->Branch("Jet_TriggerTrack_dEta",&dummy,"Jet_TriggerTrack_dEta[Jet_NumTriggerTracks]/F");
425  fJetTree->Branch("Jet_TriggerTrack_dPhi",&dummy,"Jet_TriggerTrack_dPhi[Jet_NumTriggerTracks]/F");
426  }
427 
428  fInitialized = kTRUE;
429 }
430 
431 //________________________________________________________________________
433  AliAnalysisTaskEmcalJet("AliAnalysisTaskJetExtractor", kTRUE),
434  #if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
435  fPythonCLI(),
436  #endif
437  fJetTree(0),
438  fVtxTagger(0),
439  fHadronMatchingRadius(0.4),
440  fTrueJetMatchingRadius(0.4),
441  fSecondaryVertexMaxChi2(1e10),
442  fSecondaryVertexMaxDispersion(0.05),
443  fCalculateSecondaryVertices(kTRUE),
444  fCalculateModelBackground(kFALSE),
445  fBackgroundModelFileName(),
446  fVertexerCuts(0),
447  fSetEmcalJetFlavour(0),
448  fEventPercentage(1.0),
449  fTruthMinLabel(0),
450  fTruthMaxLabel(100000),
451  fSaveTrackPDGCode(kTRUE),
452  fRandomSeed(0),
453  fRandomSeedCones(0),
454  fEventCut_TriggerTrackMinPt(0),
455  fEventCut_TriggerTrackMaxPt(0),
456  fEventCut_TriggerTrackMinLabel(-9999999),
457  fEventCut_TriggerTrackMaxLabel(+9999999),
458  fJetsCont(0),
459  fTracksCont(0),
460  fClustersCont(0),
461  fJetOutputArray(0),
462  fTruthParticleArray(0),
463  fTruthJetsArrayName(""),
464  fTruthJetsRhoName(""),
465  fTruthJetsRhoMassName(""),
466  fTruthParticleArrayName("mcparticles"),
467  fRandomGenerator(0),
468  fRandomGeneratorCones(0),
469  fEventWeight(0.),
470  fImpactParameter(0.),
471  fTriggerTracks_Pt(),
472  fTriggerTracks_Eta(),
473  fTriggerTracks_Phi()
474 {
475  fRandomGenerator = new TRandom3();
476  fRandomGeneratorCones = new TRandom3();
477 
479  fJetTree = new AliEmcalJetTree(GetName());
482  fJetTree->SetSaveJetShapes(kTRUE);
483  DefineOutput(2, TTree::Class());
484 }
485 
486 //________________________________________________________________________
488  AliAnalysisTaskEmcalJet(name, kTRUE),
489  #if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
490  fPythonCLI(),
491  #endif
492  fJetTree(0),
493  fVtxTagger(0),
501  fVertexerCuts(0),
503  fEventPercentage(1.0),
504  fTruthMinLabel(0),
505  fTruthMaxLabel(100000),
506  fSaveTrackPDGCode(kTRUE),
507  fRandomSeed(0),
508  fRandomSeedCones(0),
513  fJetsCont(0),
514  fTracksCont(0),
515  fClustersCont(0),
516  fJetOutputArray(0),
519  fTruthJetsRhoName(""),
521  fTruthParticleArrayName("mcparticles"),
522  fRandomGenerator(0),
524  fEventWeight(0.),
525  fImpactParameter(0.),
529 {
530  fRandomGenerator = new TRandom3();
531  fRandomGeneratorCones = new TRandom3();
532 
534  fJetTree = new AliEmcalJetTree(GetName());
537  fJetTree->SetSaveJetShapes(kTRUE);
538  DefineOutput(2, TTree::Class());
539 }
540 
541 //________________________________________________________________________
543 {
544  // Destructor.
545 }
546 
547 //________________________________________________________________________
549 {
551 
552  // ### Basic container settings
554  if(!fJetsCont)
555  AliFatal("Jet input container not found!");
556  fJetsCont->PrintCuts();
558  if(!fTracksCont)
559  AliFatal("Particle input container not found attached to jets!");
562  AliFatal("Cluster input container not found attached to jets!");
563 
564  fRandomGenerator->SetSeed(fRandomSeed);
566 
567  // Activate saving of trigger tracks if this is demanded
570 
571  // ### Initialize the jet tree (settings must all be given at this stage)
577  OpenFile(2);
578  PostData(2, fJetTree->GetTreePointer());
579 
580  // ### Add control histograms (already some created in base task)
581  AddHistogram2D<TH2D>("hTrackCount", "Number of tracks in acceptance vs. centrality", "COLZ", 500, 0., 5000., 100, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}");
582  AddHistogram2D<TH2D>("hBackgroundPt", "Background p_{T} distribution", "", 150, 0., 150., 100, 0, 100, "Background p_{T} (GeV/c)", "Centrality", "dN^{Events}/dp_{T}");
583 
584  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}");
585  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}");
586  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}");
587  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");
588  AddHistogram2D<TH2D>("hJetArea", "Jet area", "COLZ", 200, 0., 2., 100, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
589 
590  AddHistogram2D<TH2D>("hDeltaPt", "#delta p_{T} distribution", "", 400, -100., 300., 100, 0, 100, "p_{T, cone} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
591 
592  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}");
593  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");
594 
595  AddHistogram1D<TH1D>("hExtractionPercentage", "Extracted jets p_{T} distribution (background subtracted)", "e1p", 400, -100., 300., "p_{T, jet} (GeV/c)", "Extracted percentage");
596 
597  // Track QA plots
598  AddHistogram2D<TH2D>("hTrackPt", "Tracks p_{T} distribution", "", 300, 0., 300., 100, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
599  AddHistogram2D<TH2D>("hTrackPhi", "Track angular distribution in #phi", "LEGO2", 180, 0., 2*TMath::Pi(), 100, 0, 100, "#phi", "Centrality", "dN^{Tracks}/(d#phi)");
600  AddHistogram2D<TH2D>("hTrackEta", "Track angular distribution in #eta", "LEGO2", 100, -2.5, 2.5, 100, 0, 100, "#eta", "Centrality", "dN^{Tracks}/(d#eta)");
601  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");
602 
603  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})");
604  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})");
605 }
606 
607 
608 //________________________________________________________________________
610 {
612 
613  // ### Need to explicitly tell jet container to load rho mass object
614  fJetsCont->LoadRhoMass(InputEvent());
615 
616  if (fTruthParticleArrayName != "")
617  fTruthParticleArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTruthParticleArrayName.Data()));
618 
619  // ### Prepare vertexer
621  {
622  if(!fVertexerCuts)
623  AliFatal("VertexerCuts not given but secondary vertex calculation turned on.");
624  fVtxTagger = new AliHFJetsTaggingVertex();
625  fVtxTagger->SetCuts(fVertexerCuts);
626  }
627 
628  // ### Save tree extraction percentage to histogram
629  std::vector<Float_t> extractionPtBins = fJetTree->GetExtractionPercentagePtBins();
630  std::vector<Float_t> extractionPercentages = fJetTree->GetExtractionPercentages();
631 
632  for(Int_t i=0; i<static_cast<Int_t>(extractionPercentages.size()); i++)
633  {
634  Double_t percentage = extractionPercentages[i];
635  for(Int_t pt=static_cast<Int_t>(extractionPtBins[i*2]); pt<static_cast<Int_t>(extractionPtBins[i*2+1]); pt++)
636  {
637  FillHistogram("hExtractionPercentage", pt, percentage);
638  }
639  }
640 
641  // ### Add HIJING/JEWEL histograms (in case we need it)
642  if(MCEvent())
643  {
644  if(dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()))
645  AddHistogram1D<TH1D>("hHIJING_Ncoll", "N_{coll} from HIJING", "e1p", 1000, 0., 5000, "N_{coll}", "dN^{Events}/dN^{N_{coll}}");
646 
647  if(dynamic_cast<AliGenHepMCEventHeader*>(MCEvent()->GenEventHeader()))
648  {
649  AddHistogram1D<TH1D>("hJEWEL_NProduced", "NProduced from JEWEL", "e1p", 1000, 0., 1000, "Nproduced", "dN^{Events}/dN^{Nproduced}");
650  AddHistogram1D<TH1D>("hJEWEL_ImpactParameter", "Impact parameter from JEWEL", "e1p", 1000, 0., 1, "IP", "dN^{Events}/dN^{IP}");
651  TProfile* evweight = new TProfile("hJEWEL_EventWeight", "Event weight from JEWEL", 1, 0,1);
652  fOutput->Add(evweight);
653  }
654  }
655 
656  // ### Model background
658  {
659  // # Prepare PYTHON cmd, load estimator
660  #if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
661  TGrid::Connect("alien://");
662 
663  if (gSystem->AccessPathName(fBackgroundModelFileName.Data()))
664  AliFatal(Form("Background model %s does not exist!", fBackgroundModelFileName.Data()));
665  TFile::Cp(fBackgroundModelFileName.Data(), "./Model.pkl");
666 
667  fPythonCLI = new TPython();
668  fPythonCLI->Exec("import sklearn, numpy");
669  fPythonCLI->Exec("estimator = sklearn.externals.joblib.load(\"./Model.pkl\")");
670  #endif
671 
672  if (!(InputEvent()->FindListObject(Form("%s_RhoMVA", fJetsCont->GetArrayName().Data())))) {
673  fJetOutputArray = new TClonesArray("AliEmcalJet");
674  fJetOutputArray->SetName(Form("%s_RhoMVA", fJetsCont->GetArrayName().Data()));
675  InputEvent()->AddObject(fJetOutputArray);
676  }
677  else {
678  AliFatal(Form("Jet output array with name %s_RhoMVA already in event! Returning", fJetsCont->GetArrayName().Data()));
679  return;
680  }
681  }
682 
683  PrintConfig();
684 
685 }
686 
687 //________________________________________________________________________
689 {
690  // ################################### EVENT SELECTION
691 
692  if(!IsTriggerTrackInEvent())
693  return kFALSE;
694  if(fRandomGenerator->Rndm() >= fEventPercentage)
695  return kFALSE;
696 
697  // ################################### EVENT PROPERTIES
699 
700  // Load vertex if possible
701  Double_t vtxX = 0;
702  Double_t vtxY = 0;
703  Double_t vtxZ = 0;
704  Long64_t eventID = 0;
705  const AliVVertex* myVertex = InputEvent()->GetPrimaryVertex();
706  if(!myVertex && MCEvent())
707  myVertex = MCEvent()->GetPrimaryVertex();
708  if(myVertex)
709  {
710  vtxX = myVertex->GetX();
711  vtxY = myVertex->GetY();
712  vtxZ = myVertex->GetZ();
713  }
714  // Get event ID from header
715  AliVHeader* eventIDHeader = InputEvent()->GetHeader();
716  if(eventIDHeader)
717  eventID = eventIDHeader->GetEventIdAsLong();
718 
719  // If available, get Ncoll from HIJING
720  if(MCEvent())
721  {
722  if(AliGenHijingEventHeader* mcHeader = dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()))
723  {
724  Double_t ncoll = mcHeader->NN() + mcHeader->NNw() + mcHeader->NwN() + mcHeader->NwNw();
725  FillHistogram("hHIJING_Ncoll", ncoll);
726  }
727  if(AliGenHepMCEventHeader* mcHeader = dynamic_cast<AliGenHepMCEventHeader*>(MCEvent()->GenEventHeader()))
728  {
729  fEventWeight = mcHeader->EventWeight();
730  fImpactParameter = mcHeader->impact_parameter();
731 
732  TProfile* tmpHist = static_cast<TProfile*>(fOutput->FindObject("hJEWEL_EventWeight"));
733  tmpHist->Fill(0., fEventWeight);
734  FillHistogram("hJEWEL_NProduced", mcHeader->NProduced());
735  FillHistogram("hJEWEL_ImpactParameter", fImpactParameter);
736  }
737  }
738 
739  // ################################### MAIN JET LOOP
740  fJetsCont->ResetCurrentID();
741  Int_t jetCount = 0;
742  while(AliEmcalJet *jet = fJetsCont->GetNextAcceptJet())
743  {
745 
746  Double_t matchDistance = 0;
747  Double_t matchedJetPt = 0;
748  Double_t matchedJetMass = 0;
749  Double_t truePtFraction = 0;
750  Int_t currentJetType_HM = 0;
751  Int_t currentJetType_PM = 0;
752  Int_t currentJetType_IC = 0;
754  {
755  // Get jet type from MC (hadron matching, parton matching definition - for HF jets)
756  GetJetType(jet, currentJetType_HM, currentJetType_PM, currentJetType_IC);
757  // Get true estimators: for pt, jet mass, ...
758  truePtFraction = GetTrueJetPtFraction(jet);
759  matchDistance = GetMatchedTrueJetObservables(jet, matchedJetPt, matchedJetMass);
760  }
761 
762  // ### CONSTITUENT LOOP: Retrieve PID values + impact parameters
763  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;
764  std::vector<Float_t> vec_d0; std::vector<Float_t> vec_d0cov; std::vector<Float_t> vec_z0; std::vector<Float_t> vec_z0cov;
765 
767  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
768  {
769  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
770  if(!particle) continue;
772  {
773  Float_t sigITS = 0; Float_t sigTPC = 0; Float_t sigTOF = 0; Float_t sigTRD = 0; Short_t recoPID = 0; Int_t truePID = 0;
774  AddPIDInformation(particle, sigITS, sigTPC, sigTOF, sigTRD, recoPID, truePID);
775  vecSigITS.push_back(sigITS); vecSigTPC.push_back(sigTPC); vecSigTOF.push_back(sigTOF); vecSigTRD.push_back(sigTRD); vecRecoPID.push_back(recoPID); vecTruePID.push_back(truePID);
776  }
778  {
779  Float_t d0 = 0; Float_t d0cov = 0; Float_t z0 = 0; Float_t z0cov = 0;
780  GetTrackImpactParameters(myVertex, dynamic_cast<AliAODTrack*>(particle), d0, d0cov, z0, z0cov);
781  vec_d0.push_back(d0); vec_d0cov.push_back(d0cov); vec_z0.push_back(z0); vec_z0cov.push_back(z0cov);
782  }
783  }
784 
785  // Reconstruct secondary vertices
786  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;
788  ReconstructSecondaryVertices(myVertex, jet, secVtx_X, secVtx_Y, secVtx_Z, secVtx_Mass, secVtx_Lxy, secVtx_SigmaLxy, secVtx_Chi2, secVtx_Dispersion);
789 
790  // Now change trigger tracks eta/phi to dEta/dPhi relative to jet
791  std::vector<Float_t> triggerTracks_dEta(fTriggerTracks_Eta);
792  std::vector<Float_t> triggerTracks_dPhi(fTriggerTracks_Phi);
793  for(UInt_t i=0; i<triggerTracks_dEta.size(); i++)
794  {
795  triggerTracks_dEta[i] = jet->Eta() - fTriggerTracks_Eta[i];
796  triggerTracks_dPhi[i] = TMath::Min(TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]), TMath::TwoPi() - TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]));
797  if( ((TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]) <= TMath::Pi()) && (jet->Phi() - fTriggerTracks_Phi[i] > 0))
798  || ((TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]) > TMath::Pi()) && (jet->Phi() - fTriggerTracks_Phi[i] < 0)) )
799  triggerTracks_dPhi[i] = -triggerTracks_dPhi[i];
800  }
801 
803  {
804  Float_t pt_ML = 0;
805  Float_t mass_ML = 0;
806  GetPtAndMassFromModel(jet, pt_ML, mass_ML);
807  AliEmcalJet *outJet = new ((*fJetOutputArray)[jetCount]) AliEmcalJet(pt_ML, jet->Eta(), jet->Phi(), mass_ML);
808  outJet->SetLabel(jet->Label());
809  outJet->SetArea(jet->Area());
810  }
811 
812  // Fill jet to tree
813  Bool_t accepted = fJetTree->AddJetToTree(jet, vtxX, vtxY, vtxZ, fCent, eventID, InputEvent()->GetMagneticField(),
814  currentJetType_PM,currentJetType_HM,currentJetType_IC,matchDistance,matchedJetPt,matchedJetMass,truePtFraction,fPtHard,fEventWeight,fImpactParameter,
815  vecSigITS, vecSigTPC, vecSigTOF, vecSigTRD, vecRecoPID, vecTruePID,
816  fTriggerTracks_Pt, triggerTracks_dEta, triggerTracks_dPhi,
817  vec_d0, vec_z0, vec_d0cov, vec_z0cov,
818  secVtx_X, secVtx_Y, secVtx_Z, secVtx_Mass, secVtx_Lxy, secVtx_SigmaLxy, secVtx_Chi2, secVtx_Dispersion);
819 
820  if(accepted)
821  FillHistogram("hJetPtExtracted", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
822  jetCount++;
823  }
824 
825  // ################################### PARTICLE PROPERTIES
826  fTracksCont->ResetCurrentID();
827  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
829 
830  return kTRUE;
831 }
832 
833 //________________________________________________________________________
835 {
836  #if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
837  // ####### Calculate inference input parameters
838  Double_t constPtMean = 0;
839  Double_t constPtMedian = 0;
840  std::vector<Int_t> index_sorted_list = jet->GetPtSortedTrackConstituentIndexes(fJetsCont->GetParticleContainer()->GetArray());
841  Int_t numConst = index_sorted_list.size();
842  Double_t* constPts = new Double_t[TMath::Max(Int_t(index_sorted_list.size()), 10)];
843 
844  // Calculate mean, median of constituents
845  for(Int_t i = 0; i < numConst; i++)
846  {
847  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(index_sorted_list.at(i), fJetsCont->GetParticleContainer()->GetArray()));
848  constPtMean += particle->Pt();
849  constPts[i] = particle->Pt();
850  }
851  if(numConst)
852  {
853  constPtMean /= numConst;
854  constPtMedian = TMath::Median(numConst, constPts);
855  }
856 
857  // ####### Run Python script that does inference on jet
858  fPythonCLI->Exec(Form("data_inference = numpy.array([[%E, %E, %E, %E, %E, %E, %E, %E, %E, %E, %E, %E, %E, %E, %E, %E, %E, %E]])", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), (Double_t)numConst, jet->GetShapeProperties()->GetSecondOrderSubtracted(), fJetsCont->GetRhoVal(), jet->M()-jet->GetShapeProperties()->GetSecondOrderSubtracted(), jet->Area(), constPtMean, constPtMedian, constPts[0], constPts[1], constPts[2], constPts[3], constPts[4], constPts[5], constPts[6], constPts[7], constPts[8], constPts[9]));
859  fPythonCLI->Exec("result = estimator.predict(data_inference)");
860  pt_ML = fPythonCLI->Eval("result[0][0]");
861  mass_ML = fPythonCLI->Eval("result[0][1]");
862 
863  delete[] constPts;
864  #endif
865 }
866 
867 //________________________________________________________________________
869 {
870  // Cut for trigger track requirement
872  {
873  // Clear vector of trigger tracks
874  fTriggerTracks_Pt.clear();
875  fTriggerTracks_Eta.clear();
876  fTriggerTracks_Phi.clear();
877 
878  // Go through all tracks and check whether trigger tracks can be found
879  fTracksCont->ResetCurrentID();
880  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
881  {
882  if( (track->GetLabel() >= fEventCut_TriggerTrackMinLabel) && (track->GetLabel() < fEventCut_TriggerTrackMaxLabel) )
883  if( (track->Pt() >= fEventCut_TriggerTrackMinPt) && (track->Pt() < fEventCut_TriggerTrackMaxPt) )
884  {
885  fTriggerTracks_Pt.push_back(track->Pt());
886  fTriggerTracks_Eta.push_back(track->Eta());
887  fTriggerTracks_Phi.push_back(track->Phi());
888  }
889  }
890  // No particle has been found that fulfills requirement -> Do not accept event
891  if(fTriggerTracks_Pt.size() == 0)
892  return kFALSE;
893  }
894  return kTRUE;
895 }
896 
897 //________________________________________________________________________
899 {
900  // #################################################################################
901  // ##### FRACTION OF TRUE PT IN JET: Defined as "not from toy"
902  Double_t pt_nonMC = 0.;
903  Double_t pt_all = 0.;
904  Double_t truePtFraction = 0;
905 
906  // Loop over all jet tracks+clusters
907  for(Int_t iConst = 0; iConst < jet->GetNumberOfTracks(); iConst++)
908  {
909  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(iConst, fTracksCont->GetArray()));
910  if(!particle) continue;
911  if(particle->Pt() < 1e-6) continue;
912 
913  // Particles marked w/ labels within label range are considered from toy
914  if( (particle->GetLabel() >= fTruthMinLabel) && (particle->GetLabel() < fTruthMaxLabel))
915  pt_nonMC += particle->Pt();
916  pt_all += particle->Pt();
917  }
918  for(Int_t iConst = 0; iConst < jet->GetNumberOfClusters(); iConst++)
919  {
920  AliVCluster* cluster = static_cast<AliVCluster*>(jet->ClusterAt(iConst, fClustersCont->GetArray()));
921  if(!cluster) continue;
922  if(cluster->E() < 1e-6) continue;
923 
924  // #### Retrieve cluster pT
925  Double_t vtxX = 0;
926  Double_t vtxY = 0;
927  Double_t vtxZ = 0;
928  const AliVVertex* myVertex = InputEvent()->GetPrimaryVertex();
929  if(!myVertex && MCEvent())
930  myVertex = MCEvent()->GetPrimaryVertex();
931  if(myVertex)
932  {
933  vtxX = myVertex->GetX();
934  vtxY = myVertex->GetY();
935  vtxZ = myVertex->GetZ();
936  }
937 
938  Double_t primVtx[3] = {vtxX, vtxY, vtxZ};
939  TLorentzVector clusterMomentum;
940  cluster->GetMomentum(clusterMomentum, primVtx);
941  Double_t ClusterPt = clusterMomentum.Perp();
942  // ####
943 
944  // Particles marked w/ labels within label range are considered from toy
945  if( (cluster->GetLabel() >= fTruthMinLabel) && (cluster->GetLabel() < fTruthMaxLabel))
946  pt_nonMC += ClusterPt;
947  pt_all += ClusterPt;
948  }
949 
950  if(pt_all)
951  truePtFraction = (pt_nonMC/pt_all);
952  return truePtFraction;
953 }
954 
955 //________________________________________________________________________
957 {
958  // #################################################################################
959  // ##### OBSERVABLES FROM MATCHED JETS: Jet pt, jet mass
960  Double_t bestMatchDeltaR = 8.; // 8 is higher than maximum possible matching distance
961  if(fTruthJetsArrayName != "")
962  {
963  // "True" background for pt
964  AliRhoParameter* rho = static_cast<AliRhoParameter*>(InputEvent()->FindListObject(fTruthJetsRhoName.Data()));
965  Double_t trueRho = 0;
966  if(rho)
967  trueRho = rho->GetVal();
968 
969  // "True" background for mass
970  AliRhoParameter* rhoMass = static_cast<AliRhoParameter*>(InputEvent()->FindListObject(fTruthJetsRhoMassName.Data()));
971  TClonesArray* truthArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fTruthJetsArrayName.Data())));
972 
973  // Loop over all true jets to find the best match
974  matchedJetPt = 0;
975  matchedJetMass = 0;
976  if(truthArray)
977  for(Int_t i=0; i<truthArray->GetEntries(); i++)
978  {
979  AliEmcalJet* truthJet = static_cast<AliEmcalJet*>(truthArray->At(i));
980  if(truthJet->Pt() < 0.15)
981  continue;
982 
983  Double_t deltaEta = (truthJet->Eta()-jet->Eta());
984  Double_t deltaPhi = TMath::Min(TMath::Abs(truthJet->Phi()-jet->Phi()),TMath::TwoPi() - TMath::Abs(truthJet->Phi()-jet->Phi()));
985  Double_t deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
986 
987  // Cut jets too far away
988  if (deltaR > fTrueJetMatchingRadius)
989  continue;
990 
991  // Search for the best match
992  if(deltaR < bestMatchDeltaR)
993  {
994  bestMatchDeltaR = deltaR;
995  matchedJetPt = truthJet->Pt() - truthJet->Area()* trueRho;
996  if(rhoMass)
997  matchedJetMass = truthJet->GetShapeProperties()->GetSecondOrderSubtracted();
998  else
999  matchedJetMass = truthJet->M();
1000  }
1001  }
1002  }
1003 
1004  return bestMatchDeltaR;
1005 }
1006 
1007 
1008 //________________________________________________________________________
1010 {
1012 
1013  if(!fTruthParticleArray)
1014  return;
1015 
1016  typeHM = 0;
1017  typePM = 0;
1018 
1019  AliAODMCParticle* parton[2];
1020  parton[0] = (AliAODMCParticle*) fVtxTagger->IsMCJetParton(fTruthParticleArray, jet, radius); // method 1 (parton matching)
1021  parton[1] = (AliAODMCParticle*) fVtxTagger->IsMCJetMeson(fTruthParticleArray, jet, radius); // method 2 (hadron matching)
1022 
1023  if (parton[0]) {
1024  Int_t pdg = TMath::Abs(parton[0]->PdgCode());
1025  typePM = pdg;
1026  }
1027 
1028  if (!parton[1])
1029  {
1030  // No HF jet (according to hadron matching) -- now also separate jets in udg (1) and s-jets (3)
1031  if(IsStrangeJet(jet))
1032  typeHM = 3;
1033  else
1034  typeHM = 1;
1035  }
1036  else {
1037  Int_t pdg = TMath::Abs(parton[1]->PdgCode());
1038  if(fVtxTagger->IsDMeson(pdg)) typeHM = 4;
1039  else if (fVtxTagger->IsBMeson(pdg)) typeHM = 5;
1040  }
1041 
1042  // Set flavour of AliEmcalJet object (set ith bit while i corresponds to type)
1044  jet->AddFlavourTag(static_cast<Int_t>(TMath::Power(2, typeHM)));
1045 
1046 
1047  const AliEmcalPythiaInfo* partonsInfo = GetPythiaInfo();
1048  typeIC = 0;
1049  if (partonsInfo)
1050  {
1051  // Get primary partons directions
1052  Double_t parton1phi = partonsInfo->GetPartonPhi6();
1053  Double_t parton1eta = partonsInfo->GetPartonEta6();
1054  Double_t parton2phi = partonsInfo->GetPartonPhi7();
1055  Double_t parton2eta = partonsInfo->GetPartonEta7();
1056 
1057 
1058  Double_t delta1Eta = (parton1eta-jet->Eta());
1059  Double_t delta1Phi = TMath::Min(TMath::Abs(parton1phi-jet->Phi()),TMath::TwoPi() - TMath::Abs(parton1phi-jet->Phi()));
1060  Double_t delta1R = TMath::Sqrt(delta1Eta*delta1Eta + delta1Phi*delta1Phi);
1061  Double_t delta2Eta = (parton2eta-jet->Eta());
1062  Double_t delta2Phi = TMath::Min(TMath::Abs(parton2phi-jet->Phi()),TMath::TwoPi() - TMath::Abs(parton2phi-jet->Phi()));
1063  Double_t delta2R = TMath::Sqrt(delta2Eta*delta2Eta + delta2Phi*delta2Phi);
1064 
1065  // Check if one of the partons if closer than matching criterion
1066  Bool_t matched = (delta1R < fJetsCont->GetJetRadius()/2.) || (delta2R < fJetsCont->GetJetRadius()/2.);
1067 
1068  // Matching criterion fulfilled -> Set flag to closest
1069  if(matched)
1070  {
1071  if(delta1R < delta2R)
1072  typeIC = partonsInfo->GetPartonFlag6();
1073  else
1074  typeIC = partonsInfo->GetPartonFlag7();
1075  }
1076  }
1077 }
1078 
1079 //________________________________________________________________________
1080 void AliAnalysisTaskJetExtractor::GetTrackImpactParameters(const AliVVertex* vtx, AliAODTrack* track, Float_t& d0, Float_t& d0cov, Float_t& z0, Float_t& z0cov)
1081 {
1082  if (track)
1083  {
1084  Double_t d0rphiz[2],covd0[3];
1085  Bool_t isDCA=track->PropagateToDCA(vtx,InputEvent()->GetMagneticField(),3.0,d0rphiz,covd0);
1086  if(isDCA)
1087  {
1088  d0 = d0rphiz[0];
1089  z0 = d0rphiz[1];
1090  d0cov = covd0[0];
1091  z0cov = covd0[2];
1092  }
1093  }
1094 }
1095 
1096 //________________________________________________________________________
1097 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)
1098 {
1099  if(!primVtx)
1100  return;
1101 
1102  // Create ESD vertex from the existing AliVVertex
1103  Double_t vtxPos[3] = {primVtx->GetX(), primVtx->GetY(), primVtx->GetZ()};
1104  Double_t covMatrix[6] = {0};
1105  primVtx->GetCovarianceMatrix(covMatrix);
1106  AliESDVertex* esdVtx = new AliESDVertex(vtxPos, covMatrix, primVtx->GetChi2(), primVtx->GetNContributors());
1107 
1108  TClonesArray* secVertexArr = 0;
1109  vctr_pair_dbl_int arrDispersion;
1110  arrDispersion.reserve(5);
1112  {
1113  //###########################################################################
1114  // ********* Calculate secondary vertices
1115  // Derived from AliAnalysisTaskEmcalJetBtagSV
1116  secVertexArr = new TClonesArray("AliAODVertex");
1117  Int_t nDauRejCount = 0;
1118  Int_t nVtx = fVtxTagger->FindVertices(jet,
1119  fTracksCont->GetArray(),
1120  (AliAODEvent*)InputEvent(),
1121  esdVtx,
1122  InputEvent()->GetMagneticField(),
1123  secVertexArr,
1124  0,
1125  arrDispersion,
1126  nDauRejCount);
1127 
1128 
1129  if(nVtx < 0)
1130  {
1131  secVertexArr->Clear();
1132  delete secVertexArr;
1133  return;
1134  }
1135  //###########################################################################
1136  }
1137  else // Load HF vertex branch from AOD event, if possible
1138  {
1139  secVertexArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("VerticesHF"));
1140  if(!secVertexArr)
1141  return;
1142  }
1143 
1144  // Loop over all potential secondary vertices
1145  for(Int_t i=0; i<secVertexArr->GetEntriesFast(); i++)
1146  {
1147  AliAODVertex* secVtx = (AliAODVertex*)(secVertexArr->UncheckedAt(i));
1149  if((strcmp(secVtx->GetParent()->ClassName(), "AliAODRecoDecayHF3Prong")))
1150  continue;
1151 
1152  // Calculate vtx distance
1153  Double_t effX = secVtx->GetX() - esdVtx->GetX();
1154  Double_t effY = secVtx->GetY() - esdVtx->GetY();
1155  //Double_t effZ = secVtx->GetZ() - esdVtx->GetZ();
1156 
1157  // ##### Vertex properties
1158  // vertex dispersion
1159  Double_t dispersion = arrDispersion[i].first;
1160 
1161  // invariant mass
1162  Double_t mass = fVtxTagger->GetVertexInvariantMass(secVtx);
1163 
1164  // signed length
1165  Double_t Lxy = TMath::Sqrt(effX*effX + effY*effY);
1166  Double_t jetP[3]; jet->PxPyPz(jetP);
1167  Double_t signLxy = effX * jetP[0] + effY * jetP[1];
1168  if (signLxy < 0.) Lxy *= -1.;
1169 
1170  Double_t sigmaLxy = 0;
1171  AliAODVertex* aodVtx = (AliAODVertex*)(primVtx);
1172  if (aodVtx)
1173  sigmaLxy = aodVtx->ErrorDistanceXYToVertex(secVtx);
1174 
1175  // Add secondary vertices if they fulfill the conditions
1176  if( (dispersion > fSecondaryVertexMaxDispersion) || (TMath::Abs(secVtx->GetChi2perNDF()) > fSecondaryVertexMaxChi2) )
1177  continue;
1178 
1179  secVtx_X.push_back(secVtx->GetX()); secVtx_Y.push_back(secVtx->GetY()); secVtx_Z.push_back(secVtx->GetZ()); secVtx_Chi2.push_back(secVtx->GetChi2perNDF());
1180  secVtx_Dispersion.push_back(dispersion); secVtx_Mass.push_back(mass); secVtx_Lxy.push_back(Lxy); secVtx_SigmaLxy.push_back(sigmaLxy);
1181  }
1182 
1184  {
1185  secVertexArr->Clear();
1186  delete secVertexArr;
1187  }
1188  delete esdVtx;
1189 }
1190 
1191 //________________________________________________________________________
1193 {
1194  // Do hadron matching jet type tagging using mcparticles
1195  // ... if not explicitly deactivated
1196  if (fTruthParticleArray)
1197  {
1198  for(Int_t i=0; i<fTruthParticleArray->GetEntries();i++)
1199  {
1200  AliAODMCParticle* part = (AliAODMCParticle*)fTruthParticleArray->At(i);
1201  if(!part) continue;
1202 
1203  // Check if the particle has strangeness
1204  Int_t absPDG = TMath::Abs(part->PdgCode());
1205  if ((absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
1206  {
1207  // Check if particle is in a radius around the jet
1208  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
1210  continue;
1211  else
1212  return kTRUE;
1213  }
1214  }
1215  }
1216  // As fallback, the MC stack will be tried
1217  else if(MCEvent() && (MCEvent()->Stack()))
1218  {
1219  AliStack* stack = MCEvent()->Stack();
1220  // Go through the whole particle stack
1221  for(Int_t i=0; i<stack->GetNtrack(); i++)
1222  {
1223  TParticle *part = stack->Particle(i);
1224  if(!part) continue;
1225 
1226  // Check if the particle has strangeness
1227  Int_t absPDG = TMath::Abs(part->GetPdgCode());
1228  if ((absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
1229  {
1230  // Check if particle is in a radius around the jet
1231  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
1233  continue;
1234  else
1235  return kTRUE;
1236  }
1237  }
1238  }
1239  return kFALSE;
1240 
1241 }
1242 //________________________________________________________________________
1244 {
1245  // ### Event control plots
1247  FillHistogram("hBackgroundPt", fJetsCont->GetRhoVal(), fCent);
1248 }
1249 
1250 //________________________________________________________________________
1252 {
1253  // ### Jet control plots
1254  FillHistogram("hJetPtRaw", jet->Pt(), fCent);
1255  FillHistogram("hJetPt", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
1256  FillHistogram("hJetPhiEta", jet->Phi(), jet->Eta());
1257  FillHistogram("hJetArea", jet->Area(), fCent);
1258 
1259  // ### Jet constituent plots
1260  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
1261  {
1262  AliVParticle* jetConst = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
1263  if(!jetConst) continue;
1264 
1265  // Constituent eta/phi (relative to jet)
1266  Double_t deltaEta = jet->Eta() - jetConst->Eta();
1267  Double_t deltaPhi = TMath::Min(TMath::Abs(jet->Phi() - jetConst->Phi()), TMath::TwoPi() - TMath::Abs(jet->Phi() - jetConst->Phi()));
1268  if(!((jet->Phi() - jetConst->Phi() < 0) && (jet->Phi() - jetConst->Phi() <= TMath::Pi())))
1269  deltaPhi = -deltaPhi;
1270 
1271  FillHistogram("hConstituentPt", jetConst->Pt(), fCent);
1272  FillHistogram("hConstituentPhiEta", deltaPhi, deltaEta);
1273  }
1274 
1275  // ### Random cone / delta pT plots
1276  const Int_t kNumRandomConesPerEvent = 4;
1277  for(Int_t iCone=0; iCone<kNumRandomConesPerEvent; iCone++)
1278  {
1279  // Throw random cone
1280  Double_t tmpRandConeEta = fJetsCont->GetJetEtaMin() + fRandomGeneratorCones->Rndm()*TMath::Abs(fJetsCont->GetJetEtaMax()-fJetsCont->GetJetEtaMin());
1281  Double_t tmpRandConePhi = fRandomGeneratorCones->Rndm()*TMath::TwoPi();
1282  Double_t tmpRandConePt = 0;
1283  // Fill pT that is in cone
1284  fTracksCont->ResetCurrentID();
1285  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
1286  if(IsTrackInCone(track, tmpRandConeEta, tmpRandConePhi, fJetsCont->GetJetRadius()))
1287  tmpRandConePt += track->Pt();
1288 
1289  // Fill histograms
1290  FillHistogram("hDeltaPt", tmpRandConePt - fJetsCont->GetRhoVal()*fJetsCont->GetJetRadius()*fJetsCont->GetJetRadius()*TMath::Pi(), fCent);
1291  }
1292 }
1293 
1294 //________________________________________________________________________
1296 {
1297  FillHistogram("hTrackPt", track->Pt(), fCent);
1298  FillHistogram("hTrackPhi", track->Phi(), fCent);
1299  FillHistogram("hTrackEta", track->Eta(), fCent);
1300  FillHistogram("hTrackEtaPt", track->Eta(), track->Pt());
1301  FillHistogram("hTrackPhiPt", track->Phi(), track->Pt());
1302  FillHistogram("hTrackPhiEta", track->Phi(), track->Eta());
1303 }
1304 
1305 //________________________________________________________________________
1306 void AliAnalysisTaskJetExtractor::AddPIDInformation(AliVParticle* particle, Float_t& sigITS, Float_t& sigTPC, Float_t& sigTOF, Float_t& sigTRD, Short_t& recoPID, Int_t& truePID)
1307 {
1308  truePID = 9;
1309  if(!particle) return;
1310 
1311  // If we have AODs, retrieve particle PID signals
1312  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(particle);
1313 
1314  if(aodtrack)
1315  {
1316  // Get AOD value from reco
1317  recoPID = aodtrack->GetMostProbablePID();
1318  AliAODPid* pidObj = aodtrack->GetDetPid();
1319  if(!pidObj)
1320  return;
1321 
1322  sigITS = pidObj->GetITSsignal();
1323  sigTPC = pidObj->GetTPCsignal();
1324  sigTOF = pidObj->GetTOFsignal();
1325  sigTRD = pidObj->GetTRDsignal();
1326  }
1327 
1328  // Get truth values if we are on MC
1330  {
1331  for(Int_t i=0; i<fTruthParticleArray->GetEntries();i++)
1332  {
1333  AliAODMCParticle* mcParticle = dynamic_cast<AliAODMCParticle*>(fTruthParticleArray->At(i));
1334  if(!mcParticle) continue;
1335 
1336  if (mcParticle->GetLabel() == particle->GetLabel())
1337  {
1338  if(fSaveTrackPDGCode)
1339  {
1340  truePID = mcParticle->PdgCode();
1341  }
1342  else // Use same convention as for PID in AODs
1343  {
1344  if(TMath::Abs(mcParticle->PdgCode()) == 2212) // proton
1345  truePID = 4;
1346  else if (TMath::Abs(mcParticle->PdgCode()) == 211) // pion
1347  truePID = 2;
1348  else if (TMath::Abs(mcParticle->PdgCode()) == 321) // kaon
1349  truePID = 3;
1350  else if (TMath::Abs(mcParticle->PdgCode()) == 11) // electron
1351  truePID = 0;
1352  else if (TMath::Abs(mcParticle->PdgCode()) == 13) // muon
1353  truePID = 1;
1354  else if (TMath::Abs(mcParticle->PdgCode()) == 700201) // deuteron
1355  truePID = 5;
1356  else if (TMath::Abs(mcParticle->PdgCode()) == 700301) // triton
1357  truePID = 6;
1358  else if (TMath::Abs(mcParticle->PdgCode()) == 700302) // He3
1359  truePID = 7;
1360  else if (TMath::Abs(mcParticle->PdgCode()) == 700202) // alpha
1361  truePID = 8;
1362  else
1363  truePID = 9;
1364  }
1365 
1366  break;
1367  }
1368  }
1369  }
1370 }
1371 
1372 //________________________________________________________________________
1374 {
1375  std::cout << "#########################################\n";
1376  std::cout << "Settings for extractor task " << GetName() << std::endl;
1377  std::cout << std::endl;
1378 
1379  std::cout << "### Event cuts ###" << std::endl;
1380  std::cout << std::endl;
1382  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;
1383  if(fEventPercentage < 1)
1384  std::cout << Form("* Artificially lowered event efficiency: %f", fEventPercentage) << std::endl;
1386  std::cout << Form("* Random seeds explicitly set: %lu (for event/jet eff.), %lu (RCs)", fRandomSeed, fRandomSeedCones) << std::endl;
1387  std::cout << std::endl;
1388 
1389  std::cout << "### Jet properties ###" << std::endl;
1390  std::cout << "* Jet array name = " << fJetsCont->GetName() << std::endl;
1391  std::cout << "* Rho name = " << fJetsCont->GetRhoName() << std::endl;
1392  std::cout << "* Rho mass name = " << fJetsCont->GetRhoMassName() << std::endl;
1393  std::cout << std::endl;
1394 
1395  std::cout << "### Tree extraction properties ###" << std::endl;
1396  std::vector<Float_t> extractionPtBins = fJetTree->GetExtractionPercentagePtBins();
1397  std::vector<Float_t> extractionPercentages = fJetTree->GetExtractionPercentages();
1398  std::vector<Int_t> extractionHM = fJetTree->GetExtractionJetTypes_HM();
1399  std::vector<Int_t> extractionPM = fJetTree->GetExtractionJetTypes_PM();
1400  if(extractionPercentages.size())
1401  {
1402  std::cout << "* Extraction bins: (";
1403  for(Int_t i=0; i<static_cast<Int_t>(extractionPercentages.size()-1); i++)
1404  std::cout << extractionPtBins[i*2] << "-" << extractionPtBins[i*2+1] << " = " << extractionPercentages[i] << ", ";
1405  std::cout << extractionPtBins[(extractionPercentages.size()-1)*2] << "-" << extractionPtBins[(extractionPercentages.size()-1)*2+1] << " = " << extractionPercentages[(extractionPercentages.size()-1)];
1406 
1407  std::cout << ")" << std::endl;
1408  }
1409  if(extractionHM.size())
1410  {
1411  std::cout << "* Extraction of hadron-matched jets with types: (";
1412  for(Int_t i=0; i<static_cast<Int_t>(extractionHM.size()-1); i++)
1413  std::cout << extractionHM[i] << ", ";
1414  std::cout << extractionHM[extractionHM.size()-1];
1415  std::cout << ")" << std::endl;
1416  }
1417  if(extractionPM.size())
1418  {
1419  std::cout << "* Extraction of parton-matched jets with types: (";
1420  for(Int_t i=0; i<static_cast<Int_t>(extractionPM.size()-1); i++)
1421  std::cout << extractionPM[i] << ", ";
1422  std::cout << extractionPM[extractionPM.size()-1];
1423  std::cout << ")" << std::endl;
1424  }
1425  std::cout << std::endl;
1426 
1427  std::cout << "### Extracted data groups ###" << std::endl;
1429  std::cout << "* Basic event properties (ID, vertex, centrality, bgrd. densities, ...)" << std::endl;
1431  std::cout << "* Jet constituents, basic properties (pt, eta, phi, charge, ...)" << std::endl;
1433  std::cout << "* Jet constituents, IPs" << std::endl;
1435  std::cout << "* Jet constituents, PID" << std::endl;
1437  std::cout << "* Jet calorimeter clusters" << std::endl;
1439  std::cout << "* MC information (origin, matched jets, ...)" << std::endl;
1441  std::cout << "* Secondary vertices" << std::endl;
1442  if(fJetTree->GetSaveJetShapes())
1443  std::cout << "* Jet shapes (jet mass)" << std::endl;
1445  std::cout << "* Trigger tracks" << std::endl;
1446  std::cout << std::endl;
1447  std::cout << "### Further settings ###" << std::endl;
1448  if(fTruthJetsArrayName != "")
1449  std::cout << Form("* True jet matching active, array=%s, rho=%s, rho_mass=%s", fTruthJetsArrayName.Data(), fTruthJetsRhoName.Data(), fTruthJetsRhoMassName.Data()) << std::endl;
1451  std::cout << Form("* Particle level information available (for jet origin calculation, particle code): %s", fTruthParticleArrayName.Data()) << std::endl;
1452  if(extractionHM.size())
1453  std::cout << Form("* Hadron--jet matching radius=%3.3f", fHadronMatchingRadius) << std::endl;
1455  std::cout << Form("* True particle label range for pt fraction=%d-%d", fTruthMinLabel, fTruthMaxLabel) << std::endl;
1457  std::cout << Form("* Particle PID codes will be PDG codes") << std::endl;
1459  std::cout << Form("* Particle PID codes will follow AliAODPid convention") << std::endl;
1461  std::cout << Form("* AliEmcalJet flavour tag is set from HF-jet tagging") << std::endl;
1462 
1463  std::cout << "#########################################\n\n";
1464 }
1465 
1466 
1467 //########################################################################
1468 // HELPERS
1469 //########################################################################
1470 
1471 //________________________________________________________________________
1472 inline Bool_t AliAnalysisTaskJetExtractor::IsTrackInCone(AliVParticle* track, Double_t eta, Double_t phi, Double_t radius)
1473 {
1474  // This is to use a full cone in phi even at the edges of phi (2pi -> 0) (0 -> 2pi)
1475  Double_t trackPhi = 0.0;
1476  if (track->Phi() > (TMath::TwoPi() - (radius-phi)))
1477  trackPhi = track->Phi() - TMath::TwoPi();
1478  else if (track->Phi() < (phi+radius - TMath::TwoPi()))
1479  trackPhi = track->Phi() + TMath::TwoPi();
1480  else
1481  trackPhi = track->Phi();
1482 
1483  if ( TMath::Abs(trackPhi-phi)*TMath::Abs(trackPhi-phi) + TMath::Abs(track->Eta()-eta)*TMath::Abs(track->Eta()-eta) <= radius*radius)
1484  return kTRUE;
1485 
1486  return kFALSE;
1487 }
1488 
1489 //________________________________________________________________________
1491 {
1492  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1493  if(!tmpHist)
1494  {
1495  AliError(Form("Cannot find histogram <%s> ",key)) ;
1496  return;
1497  }
1498 
1499  tmpHist->Fill(x);
1500 }
1501 
1502 //________________________________________________________________________
1504 {
1505  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1506  if(!tmpHist)
1507  {
1508  AliError(Form("Cannot find histogram <%s> ",key));
1509  return;
1510  }
1511 
1512  if (tmpHist->IsA()->GetBaseClass("TH1"))
1513  static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
1514  else if (tmpHist->IsA()->GetBaseClass("TH2"))
1515  static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
1516 }
1517 
1518 //________________________________________________________________________
1520 {
1521  TH2* tmpHist = static_cast<TH2*>(fOutput->FindObject(key));
1522  if(!tmpHist)
1523  {
1524  AliError(Form("Cannot find histogram <%s> ",key));
1525  return;
1526  }
1527 
1528  tmpHist->Fill(x,y,add);
1529 }
1530 
1531 //________________________________________________________________________
1533 {
1534  TH3* tmpHist = static_cast<TH3*>(fOutput->FindObject(key));
1535  if(!tmpHist)
1536  {
1537  AliError(Form("Cannot find histogram <%s> ",key));
1538  return;
1539  }
1540 
1541  if(add)
1542  tmpHist->Fill(x,y,z,add);
1543  else
1544  tmpHist->Fill(x,y,z);
1545 }
1546 
1547 
1548 //________________________________________________________________________
1549 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)
1550 {
1551  T* tmpHist = new T(name, title, xBins, xMin, xMax);
1552 
1553  tmpHist->GetXaxis()->SetTitle(xTitle);
1554  tmpHist->GetYaxis()->SetTitle(yTitle);
1555  tmpHist->SetOption(options);
1556  tmpHist->SetMarkerStyle(kFullCircle);
1557  tmpHist->Sumw2();
1558 
1559  fOutput->Add(tmpHist);
1560 
1561  return tmpHist;
1562 }
1563 
1564 //________________________________________________________________________
1565 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)
1566 {
1567  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax);
1568  tmpHist->GetXaxis()->SetTitle(xTitle);
1569  tmpHist->GetYaxis()->SetTitle(yTitle);
1570  tmpHist->GetZaxis()->SetTitle(zTitle);
1571  tmpHist->SetOption(options);
1572  tmpHist->SetMarkerStyle(kFullCircle);
1573  tmpHist->Sumw2();
1574 
1575  fOutput->Add(tmpHist);
1576 
1577  return tmpHist;
1578 }
1579 
1580 //________________________________________________________________________
1581 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)
1582 {
1583  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax, zBins, zMin, zMax);
1584  tmpHist->GetXaxis()->SetTitle(xTitle);
1585  tmpHist->GetYaxis()->SetTitle(yTitle);
1586  tmpHist->GetZaxis()->SetTitle(zTitle);
1587  tmpHist->SetOption(options);
1588  tmpHist->SetMarkerStyle(kFullCircle);
1589  tmpHist->Sumw2();
1590 
1591  fOutput->Add(tmpHist);
1592 
1593  return tmpHist;
1594 }
1595 
1596 //________________________________________________________________________
1598 {
1599  // Called once at the end of the analysis.
1600 }
1601 
1602 // ### ADDTASK MACRO
1603 //________________________________________________________________________
1604 AliAnalysisTaskJetExtractor* AliAnalysisTaskJetExtractor::AddTaskJetExtractor(TString trackArray, TString clusterArray, TString jetArray, TString rhoObject, Double_t jetRadius, TString configFile, AliRDHFJetsCutsVertex* vertexerCuts, const char* taskNameSuffix)
1605 {
1606  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1607  Double_t minJetEta = 0.5;
1608  Double_t minJetPt = 0.15;
1609  Double_t minTrackPt = 0.15;
1610  Double_t minJetAreaPerc = 0.557;
1611  TString suffix = "";
1612  if(taskNameSuffix != 0)
1613  suffix = taskNameSuffix;
1614 
1615  // ###### Load configuration from YAML files
1617  fYAMLConfig.AddConfiguration("alien:///alice/cern.ch/user/r/rhaake/ConfigFiles/JetExtractor_BaseConfig.yaml", "baseConfig");
1618  if(configFile != "")
1619  fYAMLConfig.AddConfiguration(configFile.Data(), "customConfig");
1620  fYAMLConfig.Initialize();
1621 
1622  fYAMLConfig.GetProperty("TrackArray", trackArray, kFALSE);
1623  fYAMLConfig.GetProperty("ClusterArray", clusterArray, kFALSE);
1624  fYAMLConfig.GetProperty("JetArray", jetArray, kFALSE);
1625  fYAMLConfig.GetProperty("RhoName", rhoObject, kFALSE);
1626  fYAMLConfig.GetProperty("JetRadius", jetRadius);
1627  fYAMLConfig.GetProperty("MinJetEta", minJetEta);
1628  fYAMLConfig.GetProperty("MinJetPt", minJetPt);
1629  fYAMLConfig.GetProperty("MinTrackPt", minTrackPt);
1630  fYAMLConfig.GetProperty("MinJetAreaPerc", minJetAreaPerc);
1631  fYAMLConfig.Print();
1632 
1633  // ###### Task name
1634  TString name("AliAnalysisTaskJetExtractor");
1635  if (jetArray != "") {
1636  name += "_";
1637  name += jetArray;
1638  }
1639  if (rhoObject != "") {
1640  name += "_";
1641  name += rhoObject;
1642  }
1643  if (suffix != "") {
1644  name += "_";
1645  name += suffix;
1646  }
1647 
1648  // ###### Setup task with default settings
1650  myTask->SetNeedEmcalGeom(kFALSE);
1651  myTask->SetVzRange(-10.,10.);
1652 
1653  // Particle container and track pt cut
1654  AliParticleContainer* trackCont = 0;
1655  if(trackArray == "mcparticles")
1656  trackCont = myTask->AddMCParticleContainer(trackArray);
1657  else if(trackArray =="mctracks")
1658  trackCont = myTask->AddParticleContainer(trackArray);
1659  else
1660  trackCont = myTask->AddTrackContainer(trackArray);
1661  trackCont->SetParticlePtCut(minTrackPt);
1662 
1663  // Particle container and track pt cut
1664  AliClusterContainer* clusterCont = 0;
1665  if(clusterArray != "")
1666  clusterCont = myTask->AddClusterContainer(clusterArray);
1667 
1668  // Secondary vertex cuts (default settings from PWGHF)
1669  // (can be overwritten by using myTask->SetVertexerCuts(cuts) from outside macro)
1670  AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts", "default");
1671  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1672  esdTrackCuts->SetMinNClustersTPC(90);
1673  esdTrackCuts->SetMaxChi2PerClusterTPC(4);
1674  esdTrackCuts->SetRequireTPCRefit(kTRUE);
1675  esdTrackCuts->SetRequireITSRefit(kTRUE);
1676  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
1677  esdTrackCuts->SetMinDCAToVertexXY(0.);
1678  esdTrackCuts->SetEtaRange(-0.8, 0.8);
1679  esdTrackCuts->SetPtRange(1.0, 1.e10);
1680  vertexerCuts->AddTrackCuts(esdTrackCuts);
1681  vertexerCuts->SetNprongs(3);
1682  vertexerCuts->SetMinPtHardestTrack(1.0);//default 0.3
1683  vertexerCuts->SetSecVtxWithKF(kFALSE);//default with StrLinMinDist
1684  vertexerCuts->SetImpParCut(0.);//default 0
1685  vertexerCuts->SetDistPrimSec(0.);//default 0
1686  vertexerCuts->SetCospCut(-1);//default -1
1687  myTask->SetVertexerCuts(vertexerCuts);
1688 
1689  // Jet container
1690  AliJetContainer *jetCont = myTask->AddJetContainer(jetArray,6,jetRadius);
1691  if (jetCont) {
1692  jetCont->SetRhoName(rhoObject);
1693  jetCont->SetPercAreaCut(minJetAreaPerc);
1694  jetCont->SetJetPtCut(minJetPt);
1695  jetCont->SetLeadingHadronType(0);
1696  jetCont->SetPtBiasJetTrack(minTrackPt);
1697  jetCont->SetJetEtaLimits(-minJetEta, +minJetEta);
1698  jetCont->ConnectParticleContainer(trackCont);
1699  if(clusterCont)
1700  jetCont->ConnectClusterContainer(clusterCont);
1701  jetCont->SetMaxTrackPt(1000);
1702  }
1703 
1704  mgr->AddTask(myTask);
1705 
1706  // ###### Connect inputs/outputs
1707  mgr->ConnectInput (myTask, 0, mgr->GetCommonInputContainer() );
1708  mgr->ConnectOutput (myTask, 1, mgr->CreateContainer(Form("%s_histos", name.Data()), AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, Form("%s:ChargedJetsHadronCF", mgr->GetCommonFileName())) );
1709  mgr->ConnectOutput (myTask, 2, mgr->CreateContainer(Form("%s_tree", name.Data()), TTree::Class(), AliAnalysisManager::kOutputContainer, mgr->GetCommonFileName()) );
1710 
1711  return myTask;
1712 }
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
void AddFlavourTag(Int_t tag)
Definition: AliEmcalJet.h:351
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 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
TTree * fJetTree
! tree structure
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 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.
bool GetProperty(std::vector< std::string > propertyPath, const std::string &propertyName, T &property, const bool requiredProperty) const
Bool_t fSaveJetShapes
save jet shapes
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
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
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
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
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
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
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)
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
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)
Float_t * fBuffer_Const_ProdVtx_Y
! array buffer
Bool_t fSaveConstituents
save arrays of constituent basic properties
YAML configuration class for AliPhysics.
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_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.
Float_t fBuffer_Event_Vertex_X
! array buffer
bool Bool_t
Definition: External.C:53
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:361
Double_t yMin
Int_t GetNAcceptedParticles() const
void ConnectClusterContainer(AliClusterContainer *c)
void SetRandomGenerator(TRandom3 *gen)
void SetMaxTrackPt(Float_t b)
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_Const_ProdVtx_Z
! 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
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