AliPhysics  41af4b0 (41af4b0)
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 <cstdlib>
17 #include <algorithm>
18 #include <vector>
19 #include <TClonesArray.h>
20 #include <TH1F.h>
21 #include <TH2F.h>
22 #include <TH3F.h>
23 #include <TTree.h>
24 #include <TList.h>
25 #include <TLorentzVector.h>
26 #include <TNamed.h>
27 
28 #include "AliMCEvent.h"
29 #include "AliESDVertex.h"
30 #include "AliAODVertex.h"
31 #include "AliAODPid.h"
32 
33 #include "AliRhoParameter.h"
34 
35 #include "AliVTrack.h"
36 #include "AliVHeader.h"
37 #include "AliEmcalJet.h"
38 #include "AliLog.h"
39 #include "AliJetContainer.h"
40 #include "AliTrackContainer.h"
41 #include "AliAODTrack.h"
42 #include "AliVParticle.h"
43 #include "TRandom3.h"
44 #include "AliEmcalPythiaInfo.h"
46 
47 #include "AliGenHijingEventHeader.h"
48 #include "AliHFJetsTaggingVertex.h"
49 
50 
52 
54 ClassImp(AliEmcalJetTree)
56 
60 
61 //________________________________________________________________________
62 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), fExtractionPercentages(), fExtractionPercentagePtBins(), fExtractionJetTypes_HM(), fExtractionJetTypes_PM()
63 {
64  // For these arrays, we need to reserve memory
65  fBuffer_Const_Pt = new Float_t[kMaxNumConstituents];
66  fBuffer_Const_Eta = new Float_t[kMaxNumConstituents];
67  fBuffer_Const_Phi = new Float_t[kMaxNumConstituents];
68  fBuffer_Const_Charge = new Float_t[kMaxNumConstituents];
69  fBuffer_Const_ProdVtx_X = new Float_t[kMaxNumConstituents];
70  fBuffer_Const_ProdVtx_Y = new Float_t[kMaxNumConstituents];
71  fBuffer_Const_ProdVtx_Z = new Float_t[kMaxNumConstituents];
72 
73 }
74 
75 //________________________________________________________________________
76 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), fExtractionPercentages(), fExtractionPercentagePtBins(), fExtractionJetTypes_HM(), fExtractionJetTypes_PM()
77 {
78  // For these arrays, we need to reserve memory
86 }
87 
88 
89 //________________________________________________________________________
90 Bool_t AliEmcalJetTree::AddJetToTree(AliEmcalJet* jet, Float_t bgrdDensity, Float_t vertexX, Float_t vertexY, Float_t vertexZ, Float_t centrality, Long64_t eventID, Float_t magField,
91  AliTrackContainer* trackCont, Int_t motherParton, Int_t motherHadron, Int_t partonInitialCollision, Float_t matchedPt, Float_t truePtFraction, Float_t ptHard,
92  Float_t* trackPID_ITS, Float_t* trackPID_TPC, Float_t* trackPID_TOF, Float_t* trackPID_TRD, Short_t* trackPID_Reco, Short_t* trackPID_Truth,
93  Int_t numTriggerTracks, Float_t* triggerTrackPt, Float_t* triggerTrackDeltaEta, Float_t* triggerTrackDeltaPhi,
94  Float_t* trackIP_d0, Float_t* trackIP_z0, Float_t* trackIP_d0cov, Float_t* trackIP_z0cov,
95  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)
96 {
97  // ############
98  // Approach: Fill buffer and then add to tree
99  //
100  if(!fInitialized)
101  AliFatal("Tree is not initialized.");
102 
103  fBuffer_JetPt = jet->Pt() - bgrdDensity*jet->Area();
104 
105  // Check if jet type is contained in extraction list
106  if( (fExtractionJetTypes_PM.size() || fExtractionJetTypes_HM.size()) &&
107  (std::find(fExtractionJetTypes_PM.begin(), fExtractionJetTypes_PM.end(), motherParton) == fExtractionJetTypes_PM.end()) &&
108  (std::find(fExtractionJetTypes_HM.begin(), fExtractionJetTypes_HM.end(), motherHadron) == fExtractionJetTypes_HM.end()) )
109  return kFALSE;
110 
111  // Check acceptance percentage for the given jet and discard statistically on demand
112  Bool_t inPtRange = kFALSE;
113  for(size_t i = 0; i<fExtractionPercentages.size(); i++)
114  {
116  {
117  inPtRange = kTRUE;
118  if(gRandom->Rndm() >= fExtractionPercentages[i])
119  return kFALSE;
120  break;
121  }
122  }
123  if(!inPtRange && fExtractionPercentagePtBins.size()) // only discard jet if fExtractionPercentagePtBins was given
124  return kFALSE;
125 
126  // Set basic properties
127  fBuffer_JetEta = jet->Eta();
128  fBuffer_JetPhi = jet->Phi();
129  fBuffer_JetArea = jet->Area();
130 
131  // Set event properties
133  {
134  fBuffer_Event_BackgroundDensity = bgrdDensity;
135  fBuffer_Event_Vertex_X = vertexX;
136  fBuffer_Event_Vertex_Y = vertexY;
137  fBuffer_Event_Vertex_Z = vertexZ;
139  fBuffer_Event_ID = eventID;
140  fBuffer_Event_MagneticField = magField;
142  fBuffer_Event_PtHard = ptHard;
143  }
144 
145  // Extract basic constituent properties directly from AliEmcalJet object
147  if(trackCont && (fSaveConstituents || fSaveConstituentsIP))
148  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
149  {
150  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(i, trackCont->GetArray()));
151  if(!particle) continue;
152 
154  {
155  fBuffer_Const_Pt[fBuffer_NumConstituents] = particle->Pt();
156  fBuffer_Const_Eta[fBuffer_NumConstituents] = particle->Eta();
157  fBuffer_Const_Phi[fBuffer_NumConstituents] = particle->Phi();
158  fBuffer_Const_Charge[fBuffer_NumConstituents] = particle->Charge();
159  }
161  {
165  }
167  }
168 
169 
170  // Set constituent arrays for impact parameters
172  {
173  fJetTree->SetBranchAddress("Jet_Const_CovIPd", trackIP_d0cov);
174  fJetTree->SetBranchAddress("Jet_Const_CovIPz", trackIP_z0cov);
175  fJetTree->SetBranchAddress("Jet_Const_IPd", trackIP_d0);
176  fJetTree->SetBranchAddress("Jet_Const_IPz", trackIP_z0);
177  }
178  // Set constituent arrays for PID parameters
180  {
181  fJetTree->SetBranchAddress("Jet_Const_PID_ITS", trackPID_ITS);
182  fJetTree->SetBranchAddress("Jet_Const_PID_TPC", trackPID_TPC);
183  fJetTree->SetBranchAddress("Jet_Const_PID_TOF", trackPID_TOF);
184  fJetTree->SetBranchAddress("Jet_Const_PID_TRD", trackPID_TRD);
185  fJetTree->SetBranchAddress("Jet_Const_PID_Reconstructed", trackPID_Reco);
187  fJetTree->SetBranchAddress("Jet_Const_PID_Truth", trackPID_Truth);
188  }
189 
190  // Set secondary vertices arrays
192  {
193  fBuffer_NumSecVertices = numSecVertices;
194  fJetTree->SetBranchAddress("Jet_SecVtx_X", secVtx_X);
195  fJetTree->SetBranchAddress("Jet_SecVtx_Y", secVtx_Y);
196  fJetTree->SetBranchAddress("Jet_SecVtx_Z", secVtx_Z);
197  fJetTree->SetBranchAddress("Jet_SecVtx_Mass", secVtx_Mass);
198  fJetTree->SetBranchAddress("Jet_SecVtx_Lxy", secVtx_Lxy);
199  fJetTree->SetBranchAddress("Jet_SecVtx_SigmaLxy", secVtx_SigmaLxy);
200  fJetTree->SetBranchAddress("Jet_SecVtx_Chi2", secVtx_Chi2);
201  fJetTree->SetBranchAddress("Jet_SecVtx_Dispersion", secVtx_Dispersion);
202  }
203 
204  // Set jet shape observables
205  fBuffer_Shape_Mass = jet->M();
206 
207  // Set Monte Carlo information
209  {
210  fBuffer_Jet_MC_MotherParton = motherParton;
211  fBuffer_Jet_MC_MotherHadron = motherHadron;
212  fBuffer_Jet_MC_MotherIC = partonInitialCollision;
213  fBuffer_Jet_MC_MatchedPt = matchedPt;
214  fBuffer_Jet_MC_TruePtFraction = truePtFraction;
215  }
216 
218  {
219  fBuffer_NumTriggerTracks = numTriggerTracks;
220  fJetTree->SetBranchAddress("Jet_TriggerTrack_Pt", triggerTrackPt);
221  fJetTree->SetBranchAddress("Jet_TriggerTrack_dEta", triggerTrackDeltaEta);
222  fJetTree->SetBranchAddress("Jet_TriggerTrack_dPhi", triggerTrackDeltaPhi);
223  }
224 
225  // Add buffered jet to tree
226  fJetTree->Fill();
227 
228  return kTRUE;
229 }
230 
231 
232 //________________________________________________________________________
234 {
235  // Create the tree with active branches
236 
237  void* dummy = 0; // for some branches, do not need a buffer pointer yet
238 
239 
240  if(fInitialized)
241  AliFatal("Tree is already initialized.");
242 
243  // ### Prepare the jet tree
244  fJetTree = new TTree(Form("JetTree_%s", GetName()), "");
245 
246  fJetTree->Branch("Jet_Pt",&fBuffer_JetPt,"Jet_Pt/F");
247  fJetTree->Branch("Jet_Phi",&fBuffer_JetPhi,"Jet_Phi/F");
248  fJetTree->Branch("Jet_Eta",&fBuffer_JetEta,"Jet_Eta/F");
249  fJetTree->Branch("Jet_Area",&fBuffer_JetArea,"Jet_Area/F");
250  fJetTree->Branch("Jet_NumConstituents",&fBuffer_NumConstituents,"Jet_NumConstituents/I");
251 
253  {
254  fJetTree->Branch("Event_BackgroundDensity",&fBuffer_Event_BackgroundDensity,"Event_BackgroundDensity/F");
255  fJetTree->Branch("Event_Vertex_X",&fBuffer_Event_Vertex_X,"Event_Vertex_X/F");
256  fJetTree->Branch("Event_Vertex_Y",&fBuffer_Event_Vertex_Y,"Event_Vertex_Y/F");
257  fJetTree->Branch("Event_Vertex_Z",&fBuffer_Event_Vertex_Z,"Event_Vertex_Z/F");
258  fJetTree->Branch("Event_Centrality",&fBuffer_Event_Centrality,"Event_Centrality/F");
259  fJetTree->Branch("Event_ID",&fBuffer_Event_ID,"Event_ID/L");
260  fJetTree->Branch("Event_MagneticField",&fBuffer_Event_MagneticField,"Event_MagneticField/F");
261 
263  fJetTree->Branch("Event_PtHard",&fBuffer_Event_PtHard,"Event_PtHard/F");
264  }
265 
267  {
268  fJetTree->Branch("Jet_Const_Pt",fBuffer_Const_Pt,"Jet_Const_Pt[Jet_NumConstituents]/F");
269  fJetTree->Branch("Jet_Const_Phi",fBuffer_Const_Phi,"Jet_Const_Phi[Jet_NumConstituents]/F");
270  fJetTree->Branch("Jet_Const_Eta",fBuffer_Const_Eta,"Jet_Const_Eta[Jet_NumConstituents]/F");
271  fJetTree->Branch("Jet_Const_Charge",fBuffer_Const_Charge,"Jet_Const_Charge[Jet_NumConstituents]/F");
272  }
273 
275  {
276  fJetTree->Branch("Jet_Const_IPd",&dummy,"Jet_Const_IPd[Jet_NumConstituents]/F");
277  fJetTree->Branch("Jet_Const_IPz",&dummy,"Jet_Const_IPz[Jet_NumConstituents]/F");
278  fJetTree->Branch("Jet_Const_CovIPd",&dummy,"Jet_Const_CovIPd[Jet_NumConstituents]/F");
279  fJetTree->Branch("Jet_Const_CovIPz",&dummy,"Jet_Const_CovIPz[Jet_NumConstituents]/F");
280 
281  fJetTree->Branch("Jet_Const_ProdVtx_X",fBuffer_Const_ProdVtx_X,"Jet_Const_ProdVtx_X[Jet_NumConstituents]/F");
282  fJetTree->Branch("Jet_Const_ProdVtx_Y",fBuffer_Const_ProdVtx_Y,"Jet_Const_ProdVtx_Y[Jet_NumConstituents]/F");
283  fJetTree->Branch("Jet_Const_ProdVtx_Z",fBuffer_Const_ProdVtx_Z,"Jet_Const_ProdVtx_Z[Jet_NumConstituents]/F");
284  }
285 
287  {
288  fJetTree->Branch("Jet_Const_PID_ITS",&dummy,"Jet_Const_PID_ITS[Jet_NumConstituents]/F");
289  fJetTree->Branch("Jet_Const_PID_TPC",&dummy,"Jet_Const_PID_TPC[Jet_NumConstituents]/F");
290  fJetTree->Branch("Jet_Const_PID_TOF",&dummy,"Jet_Const_PID_TOF[Jet_NumConstituents]/F");
291  fJetTree->Branch("Jet_Const_PID_TRD",&dummy,"Jet_Const_PID_TRD[Jet_NumConstituents]/F");
292 
293  fJetTree->Branch("Jet_Const_PID_Reconstructed",&dummy,"Jet_Const_PID_Reconstructed[Jet_NumConstituents]/S");
295  fJetTree->Branch("Jet_Const_PID_Truth",&dummy,"Jet_Const_PID_Truth[Jet_NumConstituents]/S");
296  }
297 
298  if(fSaveJetShapes)
299  {
300  fJetTree->Branch("Jet_Shape_Mass",&fBuffer_Shape_Mass,"Jet_Shape_Mass/F");
301  }
302 
304  {
305  fJetTree->Branch("Jet_MC_MotherParton",&fBuffer_Jet_MC_MotherParton,"Jet_MC_MotherParton/I");
306  fJetTree->Branch("Jet_MC_MotherHadron",&fBuffer_Jet_MC_MotherHadron,"Jet_MC_MotherHadron/I");
307  fJetTree->Branch("Jet_MC_MotherIC",&fBuffer_Jet_MC_MotherIC,"Jet_MC_MotherIC/I");
308  fJetTree->Branch("Jet_MC_MatchedPt",&fBuffer_Jet_MC_MatchedPt,"Jet_MC_MatchedPt/F");
309  fJetTree->Branch("Jet_MC_TruePtFraction",&fBuffer_Jet_MC_TruePtFraction,"Jet_MC_TruePtFraction/F");
310  }
311 
313  {
314  fJetTree->Branch("Jet_NumSecVertices",&fBuffer_NumSecVertices,"Jet_NumSecVertices/I");
315 
316  fJetTree->Branch("Jet_SecVtx_X",&dummy,"Jet_SecVtx_X[Jet_NumSecVertices]/F");
317  fJetTree->Branch("Jet_SecVtx_Y",&dummy,"Jet_SecVtx_Y[Jet_NumSecVertices]/F");
318  fJetTree->Branch("Jet_SecVtx_Z",&dummy,"Jet_SecVtx_Z[Jet_NumSecVertices]/F");
319  fJetTree->Branch("Jet_SecVtx_Mass",&dummy,"Jet_SecVtx_Mass[Jet_NumSecVertices]/F");
320  fJetTree->Branch("Jet_SecVtx_Lxy",&dummy,"Jet_SecVtx_Lxy[Jet_NumSecVertices]/F");
321  fJetTree->Branch("Jet_SecVtx_SigmaLxy",&dummy,"Jet_SecVtx_SigmaLxy[Jet_NumSecVertices]/F");
322  fJetTree->Branch("Jet_SecVtx_Chi2",&dummy,"Jet_SecVtx_Chi2[Jet_NumSecVertices]/F");
323  fJetTree->Branch("Jet_SecVtx_Dispersion",&dummy,"Jet_SecVtx_Dispersion[Jet_NumSecVertices]/F");
324  }
325 
326  // Trigger track requirement active -> Save trigger track
328  {
329  fJetTree->Branch("Jet_NumTriggerTracks",&fBuffer_NumTriggerTracks,"Jet_NumTriggerTracks/I");
330  fJetTree->Branch("Jet_TriggerTrack_Pt",&dummy,"Jet_TriggerTrack_Pt[Jet_NumTriggerTracks]/F");
331  fJetTree->Branch("Jet_TriggerTrack_dEta",&dummy,"Jet_TriggerTrack_dEta[Jet_NumTriggerTracks]/F");
332  fJetTree->Branch("Jet_TriggerTrack_dPhi",&dummy,"Jet_TriggerTrack_dPhi[Jet_NumTriggerTracks]/F");
333  }
334 
335  fInitialized = kTRUE;
336 }
337 
338 //________________________________________________________________________
340  AliAnalysisTaskEmcalJet("AliAnalysisTaskJetExtractor", kTRUE),
341  fJetTree(0),
342  fVtxTagger(0),
343  fHadronMatchingRadius(0.4),
344  fTrueJetMatchingRadius(0.4),
345  fSecondaryVertexMaxChi2(1e10),
346  fSecondaryVertexMaxDispersion(0.05),
347  fCalculateSecondaryVertices(kTRUE),
348  fVertexerCuts(0),
349  fSetEmcalJetFlavour(0),
350  fEventCut_TriggerTrackMinPt(0),
351  fEventCut_TriggerTrackMaxPt(0),
352  fEventCut_TriggerTrackMinLabel(0),
353  fJetsCont(0),
354  fTracksCont(0),
355  fTruthParticleArray(0),
356  fTruthJetsArrayName(""),
357  fTruthJetsRhoName(""),
358  fTruthParticleArrayName("mcparticles"),
359  fTriggerTracks_Pt(),
360  fTriggerTracks_Eta(),
361  fTriggerTracks_Phi()
362 {
364  fJetTree = new AliEmcalJetTree(GetName());
367  fJetTree->SetSaveJetShapes(kTRUE);
368  DefineOutput(2, TTree::Class());
369 }
370 
371 //________________________________________________________________________
373  AliAnalysisTaskEmcalJet(name, kTRUE),
374  fJetTree(0),
375  fVtxTagger(0),
381  fVertexerCuts(0),
386  fJetsCont(0),
387  fTracksCont(0),
390  fTruthJetsRhoName(""),
391  fTruthParticleArrayName("mcparticles"),
395 {
397  fJetTree = new AliEmcalJetTree(GetName());
400  fJetTree->SetSaveJetShapes(kTRUE);
401  DefineOutput(2, TTree::Class());
402 }
403 
404 //________________________________________________________________________
406 {
407  // Destructor.
408 }
409 
410 //________________________________________________________________________
412 {
414 
415  // ### Basic container settings
417  if(!fJetsCont)
418  AliFatal("Jet input container not found!");
419  fJetsCont->PrintCuts();
421  if(!fTracksCont)
422  AliFatal("Particle input container not found attached to jets!");
423 
424  // Activate saving of trigger tracks if this is demanded
427 
428  // ### Initialize the jet tree (settings must all be given at this stage)
430  OpenFile(2);
431  PostData(2, fJetTree->GetTreePointer());
432 
433  // ### Add control histograms (already some created in base task)
434  AddHistogram2D<TH2D>("hTrackCount", "Number of tracks in acceptance vs. centrality", "COLZ", 500, 0., 5000., 100, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}");
435  AddHistogram2D<TH2D>("hBackgroundPt", "Background p_{T} distribution", "", 150, 0., 150., 100, 0, 100, "Background p_{T} (GeV/c)", "Centrality", "dN^{Events}/dp_{T}");
436 
437  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}");
438  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}");
439  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}");
440  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");
441  AddHistogram2D<TH2D>("hJetArea", "Jet area", "COLZ", 200, 0., 2., 100, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
442 
443  AddHistogram2D<TH2D>("hDeltaPt", "#delta p_{T} distribution", "", 400, -100., 300., 100, 0, 100, "p_{T, cone} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
444 
445  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}");
446  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");
447 
448  AddHistogram1D<TH1D>("hExtractionPercentage", "Extracted jets p_{T} distribution (background subtracted)", "e1p", 400, -100., 300., "p_{T, jet} (GeV/c)", "Extracted percentage");
449 
450  // Track QA plots
451  AddHistogram2D<TH2D>("hTrackPt", "Tracks p_{T} distribution", "", 300, 0., 300., 100, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
452  AddHistogram2D<TH2D>("hTrackPhi", "Track angular distribution in #phi", "LEGO2", 180, 0., 2*TMath::Pi(), 100, 0, 100, "#phi", "Centrality", "dN^{Tracks}/(d#phi)");
453  AddHistogram2D<TH2D>("hTrackEta", "Track angular distribution in #eta", "LEGO2", 100, -2.5, 2.5, 100, 0, 100, "#eta", "Centrality", "dN^{Tracks}/(d#eta)");
454  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");
455 
456  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})");
457  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})");
458 }
459 
460 
461 //________________________________________________________________________
463 {
465  if (fTruthParticleArrayName != "")
466  fTruthParticleArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTruthParticleArrayName.Data()));
467 
468  // ### Prepare vertexer
470  {
471  if(!fVertexerCuts)
472  AliFatal("VertexerCuts not given but secondary vertex calculation turned on.");
473  fVtxTagger = new AliHFJetsTaggingVertex();
474  fVtxTagger->SetCuts(fVertexerCuts);
475  }
476 
477  // ### Save tree extraction percentage to histogram
478  std::vector<Float_t> extractionPtBins = fJetTree->GetExtractionPercentagePtBins();
479  std::vector<Float_t> extractionPercentages = fJetTree->GetExtractionPercentages();
480 
481  for(Int_t i=0; i<static_cast<Int_t>(extractionPercentages.size()); i++)
482  {
483  Double_t percentage = extractionPercentages[i];
484  for(Int_t pt=static_cast<Int_t>(extractionPtBins[i*2]); pt<static_cast<Int_t>(extractionPtBins[i*2+1]); pt++)
485  {
486  FillHistogram("hExtractionPercentage", pt, percentage);
487  }
488  }
489 
490  // ### Add HIJING Ncoll histogram (in case we need it)
491  if(MCEvent())
492  {
493  if(dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()))
494  AddHistogram1D<TH1D>("hHijingNcoll", "N_{coll} from HIJING", "e1p", 1000, 0., 5000, "N_{coll}", "dN^{Events}/dN^{N_{coll}}");
495  }
496 }
497 
498 //________________________________________________________________________
500 {
501  // ################################### EVENT PROPERTIES
502  if(!IsEventSelected())
503  return kFALSE;
504 
506 
507  // Load vertex if possible
508  Double_t vtxX = 0;
509  Double_t vtxY = 0;
510  Double_t vtxZ = 0;
511  Long64_t eventID = 0;
512  const AliVVertex* myVertex = InputEvent()->GetPrimaryVertex();
513  if(!myVertex && MCEvent())
514  myVertex = MCEvent()->GetPrimaryVertex();
515  if(myVertex)
516  {
517  vtxX = myVertex->GetX();
518  vtxY = myVertex->GetY();
519  vtxZ = myVertex->GetZ();
520  }
521  // Get event ID from header
522  AliVHeader* eventIDHeader = InputEvent()->GetHeader();
523  if(eventIDHeader)
524  eventID = eventIDHeader->GetEventIdAsLong();
525 
526  // If available, get Ncoll from HIJING
527  if(MCEvent())
528  {
529  if(AliGenHijingEventHeader* mcHeader = dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()))
530  {
531  Double_t ncoll = mcHeader->NN() + mcHeader->NNw() + mcHeader->NwN() + mcHeader->NwNw();
532  FillHistogram("hHijingNcoll", ncoll);
533  }
534  }
535 
536  // ################################### MAIN JET LOOP
537  fJetsCont->ResetCurrentID();
538  while(AliEmcalJet *jet = fJetsCont->GetNextAcceptJet())
539  {
541 
542  Double_t matchedJetPt = 0;
543  Double_t truePtFraction = 0;
544  Int_t currentJetType_HM = 0;
545  Int_t currentJetType_PM = 0;
546  Int_t currentJetType_IC = 0;
548  {
549  // Get jet type from MC (hadron matching, parton matching definition - for HF jets)
550  GetJetType(jet, currentJetType_HM, currentJetType_PM, currentJetType_IC);
551  // Get true pT estimators
552  GetJetTruePt(jet, matchedJetPt, truePtFraction);
553  }
554 
555  // ### CONSTITUENT LOOP: Retrieve PID values + impact parameters
556  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<Short_t> vecTruePID;
557  std::vector<Float_t> vec_d0; std::vector<Float_t> vec_d0cov; std::vector<Float_t> vec_z0; std::vector<Float_t> vec_z0cov;
558 
560  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
561  {
562  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
563  if(!particle) continue;
565  {
566  Float_t sigITS = 0; Float_t sigTPC = 0; Float_t sigTOF = 0; Float_t sigTRD = 0; Short_t recoPID = 0; Short_t truePID = 0;
567  AddPIDInformation(particle, sigITS, sigTPC, sigTOF, sigTRD, recoPID, truePID);
568  vecSigITS.push_back(sigITS); vecSigTPC.push_back(sigTPC); vecSigTOF.push_back(sigTOF); vecSigTRD.push_back(sigTRD); vecRecoPID.push_back(recoPID); vecTruePID.push_back(truePID);
569  }
571  {
572  Float_t d0 = 0; Float_t d0cov = 0; Float_t z0 = 0; Float_t z0cov = 0;
573  GetTrackImpactParameters(myVertex, dynamic_cast<AliAODTrack*>(particle), d0, d0cov, z0, z0cov);
574  vec_d0.push_back(d0); vec_d0cov.push_back(d0cov); vec_z0.push_back(z0); vec_z0cov.push_back(z0cov);
575  }
576  }
577 
578  // Reconstruct secondary vertices
579  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;
581  ReconstructSecondaryVertices(myVertex, jet, secVtx_X, secVtx_Y, secVtx_Z, secVtx_Mass, secVtx_Lxy, secVtx_SigmaLxy, secVtx_Chi2, secVtx_Dispersion);
582 
583  // Now change trigger tracks eta/phi to dEta/dPhi relative to jet
584  std::vector<Float_t> triggerTracks_dEta(fTriggerTracks_Eta);
585  std::vector<Float_t> triggerTracks_dPhi(fTriggerTracks_Phi);
586  for(UInt_t i=0; i<triggerTracks_dEta.size(); i++)
587  {
588  triggerTracks_dEta[i] = jet->Eta() - triggerTracks_dEta[i];
589  triggerTracks_dPhi[i] = TMath::Min(TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]), TMath::TwoPi() - TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]));
590  if(!((jet->Phi() - fTriggerTracks_Phi[i] < 0) && (jet->Phi() - fTriggerTracks_Phi[i] <= TMath::Pi())))
591  triggerTracks_dPhi[i] = -triggerTracks_dPhi[i];
592  }
593 
594  // Fill jet to tree
595  Bool_t accepted = fJetTree->AddJetToTree(jet, fJetsCont->GetRhoVal(), vtxX, vtxY, vtxZ, fCent, eventID, InputEvent()->GetMagneticField(), fTracksCont,
596  currentJetType_PM,currentJetType_HM,currentJetType_IC,matchedJetPt,truePtFraction,fPtHard,
597  vecSigITS, vecSigTPC, vecSigTOF, vecSigTRD, vecRecoPID, vecTruePID,
598  fTriggerTracks_Pt, triggerTracks_dEta, triggerTracks_dPhi,
599  vec_d0, vec_z0, vec_d0cov, vec_z0cov,
600  secVtx_X, secVtx_Y, secVtx_Z, secVtx_Mass, secVtx_Lxy, secVtx_SigmaLxy, secVtx_Chi2, secVtx_Dispersion);
601 
602  if(accepted)
603  FillHistogram("hJetPtExtracted", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
604  }
605 
606  // ################################### PARTICLE PROPERTIES
607  fTracksCont->ResetCurrentID();
608  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
610 
611  return kTRUE;
612 }
613 
614 //________________________________________________________________________
616 {
617  // Cut for trigger track requirement
619  {
620  // Clear vector of trigger tracks
621  fTriggerTracks_Pt.clear();
622  fTriggerTracks_Eta.clear();
623  fTriggerTracks_Phi.clear();
624 
625  // Go through all tracks and check whether trigger tracks can be found
626  fTracksCont->ResetCurrentID();
627  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
628  if( (track->GetLabel() >= fEventCut_TriggerTrackMinLabel) && (track->Pt() >= fEventCut_TriggerTrackMinPt) && (track->Pt() < fEventCut_TriggerTrackMaxPt) )
629  {
630  fTriggerTracks_Pt.push_back(track->Pt());
631  fTriggerTracks_Eta.push_back(track->Eta());
632  fTriggerTracks_Phi.push_back(track->Phi());
633  }
634 
635  // No particle has been found that fulfills requirement -> Do not accept event
636  if(fTriggerTracks_Pt.size() == 0)
637  return kFALSE;
638  }
639  return kTRUE;
640 }
641 
642 
643 //________________________________________________________________________
644 void AliAnalysisTaskJetExtractor::GetJetTruePt(AliEmcalJet* jet, Double_t& matchedJetPt, Double_t& truePtFraction)
645 {
646  // #################################################################################
647  // ##### METHOD 1: If fTruthJetsArrayName is set, a matching jet is searched for
648  Double_t bestMatchDeltaR = 999.;
649  if(fTruthJetsArrayName != "")
650  {
651  // "True" background
652  AliRhoParameter* rho = static_cast<AliRhoParameter*>(InputEvent()->FindListObject(fTruthJetsRhoName.Data()));
653  Double_t trueRho = 0;
654  if(rho)
655  trueRho = rho->GetVal();
656 
657  TClonesArray* truthArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fTruthJetsArrayName.Data())));
658 
659  // Loop over all true jets to find the best match
660  matchedJetPt = 0;
661  if(truthArray)
662  for(Int_t i=0; i<truthArray->GetEntries(); i++)
663  {
664  AliEmcalJet* truthJet = static_cast<AliEmcalJet*>(truthArray->At(i));
665  if(truthJet->Pt() < 0.15)
666  continue;
667 
668  Double_t deltaEta = (truthJet->Eta()-jet->Eta());
669  Double_t deltaPhi = TMath::Min(TMath::Abs(truthJet->Phi()-jet->Phi()),TMath::TwoPi() - TMath::Abs(truthJet->Phi()-jet->Phi()));
670  Double_t deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
671 
672  // Cut jets too far away
673  if (deltaR > fTrueJetMatchingRadius)
674  continue;
675 
676  // Search for the best match
677  if(deltaR < bestMatchDeltaR)
678  {
679  bestMatchDeltaR = deltaR;
680  matchedJetPt = truthJet->Pt() - truthJet->Area()* trueRho;
681  }
682  }
683  }
684 
685  // #################################################################################
686  // ##### METHOD 2: Calculate fraction of "true" pT -- pT which is not from a toy
687  Double_t pt_nonMC = 0.;
688  Double_t pt_all = 0.;
689  truePtFraction = 0;
690 
691  for(Int_t iConst = 0; iConst < jet->GetNumberOfTracks(); iConst++)
692  {
693  // Loop over all valid jet constituents
694  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(iConst, fTracksCont->GetArray()));
695  if(!particle) continue;
696  if(particle->Pt() < 1e-6) continue;
697 
698  // Particles marked w/ labels above 100000 are considered from toy
699  if(particle->GetLabel() < 100000)
700  pt_nonMC += particle->Pt();
701  pt_all += particle->Pt();
702  }
703  if(pt_all)
704  truePtFraction = (pt_nonMC/pt_all);
705 
706 }
707 
708 
709 //________________________________________________________________________
711 {
713 
715  return;
716 
717  typeHM = 0;
718  typePM = 0;
719 
720  AliAODMCParticle* parton[2];
721  parton[0] = (AliAODMCParticle*) fVtxTagger->IsMCJetParton(fTruthParticleArray, jet, radius); // method 1 (parton matching)
722  parton[1] = (AliAODMCParticle*) fVtxTagger->IsMCJetMeson(fTruthParticleArray, jet, radius); // method 2 (hadron matching)
723 
724  if (parton[0]) {
725  Int_t pdg = TMath::Abs(parton[0]->PdgCode());
726  typePM = pdg;
727  }
728 
729  if (!parton[1])
730  {
731  // No HF jet (according to hadron matching) -- now also separate jets in udg (1) and s-jets (3)
732  if(IsStrangeJet(jet))
733  typeHM = 3;
734  else
735  typeHM = 1;
736  }
737  else {
738  Int_t pdg = TMath::Abs(parton[1]->PdgCode());
739  if(fVtxTagger->IsDMeson(pdg)) typeHM = 4;
740  else if (fVtxTagger->IsBMeson(pdg)) typeHM = 5;
741  }
742 
743  // Set flavour of AliEmcalJet object (set ith bit while i corresponds to type)
745  jet->AddFlavourTag(static_cast<Int_t>(TMath::Power(2, typeHM)));
746 
747 
748  const AliEmcalPythiaInfo* partonsInfo = GetPythiaInfo();
749  typeIC = 0;
750  if (partonsInfo)
751  {
752  // Get primary partons directions
753  Double_t parton1phi = partonsInfo->GetPartonPhi6();
754  Double_t parton1eta = partonsInfo->GetPartonEta6();
755  Double_t parton2phi = partonsInfo->GetPartonPhi7();
756  Double_t parton2eta = partonsInfo->GetPartonEta7();
757 
758 
759  Double_t delta1Eta = (parton1eta-jet->Eta());
760  Double_t delta1Phi = TMath::Min(TMath::Abs(parton1phi-jet->Phi()),TMath::TwoPi() - TMath::Abs(parton1phi-jet->Phi()));
761  Double_t delta1R = TMath::Sqrt(delta1Eta*delta1Eta + delta1Phi*delta1Phi);
762  Double_t delta2Eta = (parton2eta-jet->Eta());
763  Double_t delta2Phi = TMath::Min(TMath::Abs(parton2phi-jet->Phi()),TMath::TwoPi() - TMath::Abs(parton2phi-jet->Phi()));
764  Double_t delta2R = TMath::Sqrt(delta2Eta*delta2Eta + delta2Phi*delta2Phi);
765 
766  // Check if one of the partons if closer than matching criterion
767  Bool_t matched = (delta1R < fJetsCont->GetJetRadius()/2.) || (delta2R < fJetsCont->GetJetRadius()/2.);
768 
769  // Matching criterion fulfilled -> Set flag to closest
770  if(matched)
771  {
772  if(delta1R < delta2R)
773  typeIC = partonsInfo->GetPartonFlag6();
774  else
775  typeIC = partonsInfo->GetPartonFlag7();
776  }
777  }
778 }
779 
780 //________________________________________________________________________
781 void AliAnalysisTaskJetExtractor::GetTrackImpactParameters(const AliVVertex* vtx, AliAODTrack* track, Float_t& d0, Float_t& d0cov, Float_t& z0, Float_t& z0cov)
782 {
783  if (track)
784  {
785  Double_t d0rphiz[2],covd0[3];
786  Bool_t isDCA=track->PropagateToDCA(vtx,InputEvent()->GetMagneticField(),3.0,d0rphiz,covd0);
787  if(isDCA)
788  {
789  d0 = d0rphiz[0];
790  z0 = d0rphiz[1];
791  d0cov = covd0[0];
792  z0cov = covd0[2];
793  }
794  }
795 }
796 
797 //________________________________________________________________________
798 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)
799 {
800  if(!primVtx)
801  return;
802 
803  // Create ESD vertex from the existing AliVVertex
804  Double_t vtxPos[3] = {primVtx->GetX(), primVtx->GetY(), primVtx->GetZ()};
805  Double_t covMatrix[6] = {0};
806  primVtx->GetCovarianceMatrix(covMatrix);
807  AliESDVertex* esdVtx = new AliESDVertex(vtxPos, covMatrix, primVtx->GetChi2(), primVtx->GetNContributors());
808 
809  TClonesArray* secVertexArr = 0;
810  vctr_pair_dbl_int arrDispersion;
811  arrDispersion.reserve(5);
813  {
814  //###########################################################################
815  // ********* Calculate secondary vertices
816  // Derived from AliAnalysisTaskEmcalJetBtagSV
817  secVertexArr = new TClonesArray("AliAODVertex");
818  Int_t nDauRejCount = 0;
819  Int_t nVtx = fVtxTagger->FindVertices(jet,
820  fTracksCont->GetArray(),
821  (AliAODEvent*)InputEvent(),
822  esdVtx,
823  InputEvent()->GetMagneticField(),
824  secVertexArr,
825  0,
826  arrDispersion,
827  nDauRejCount);
828 
829 
830  if(nVtx < 0)
831  {
832  secVertexArr->Clear();
833  delete secVertexArr;
834  return;
835  }
836  //###########################################################################
837  }
838  else // Load HF vertex branch from AOD event, if possible
839  {
840  secVertexArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("VerticesHF"));
841  if(!secVertexArr)
842  return;
843  }
844 
845  // Loop over all potential secondary vertices
846  for(Int_t i=0; i<secVertexArr->GetEntriesFast(); i++)
847  {
848  AliAODVertex* secVtx = (AliAODVertex*)(secVertexArr->UncheckedAt(i));
850  if((strcmp(secVtx->GetParent()->ClassName(), "AliAODRecoDecayHF3Prong")))
851  continue;
852 
853  // Calculate vtx distance
854  Double_t effX = secVtx->GetX() - esdVtx->GetX();
855  Double_t effY = secVtx->GetY() - esdVtx->GetY();
856  //Double_t effZ = secVtx->GetZ() - esdVtx->GetZ();
857 
858  // ##### Vertex properties
859  // vertex dispersion
860  Double_t dispersion = arrDispersion[i].first;
861 
862  // invariant mass
863  Double_t mass = fVtxTagger->GetVertexInvariantMass(secVtx);
864 
865  // signed length
866  Double_t Lxy = TMath::Sqrt(effX*effX + effY*effY);
867  Double_t jetP[3]; jet->PxPyPz(jetP);
868  Double_t signLxy = effX * jetP[0] + effY * jetP[1];
869  if (signLxy < 0.) Lxy *= -1.;
870 
871  Double_t sigmaLxy = 0;
872  AliAODVertex* aodVtx = (AliAODVertex*)(primVtx);
873  if (aodVtx)
874  sigmaLxy = aodVtx->ErrorDistanceXYToVertex(secVtx);
875 
876  // Add secondary vertices if they fulfill the conditions
877  if( (dispersion > fSecondaryVertexMaxDispersion) || (TMath::Abs(secVtx->GetChi2perNDF()) > fSecondaryVertexMaxChi2) )
878  continue;
879 
880  secVtx_X.push_back(secVtx->GetX()); secVtx_Y.push_back(secVtx->GetY()); secVtx_Z.push_back(secVtx->GetZ()); secVtx_Chi2.push_back(secVtx->GetChi2perNDF());
881  secVtx_Dispersion.push_back(dispersion); secVtx_Mass.push_back(mass); secVtx_Lxy.push_back(Lxy); secVtx_SigmaLxy.push_back(sigmaLxy);
882  }
883 
885  {
886  secVertexArr->Clear();
887  delete secVertexArr;
888  }
889  delete esdVtx;
890 }
891 
892 //________________________________________________________________________
894 {
895  // Do hadron matching jet type tagging using mcparticles
896  // ... if not explicitly deactivated
898  {
899  for(Int_t i=0; i<fTruthParticleArray->GetEntries();i++)
900  {
901  AliAODMCParticle* part = (AliAODMCParticle*)fTruthParticleArray->At(i);
902  if(!part) continue;
903 
904  // Check if the particle has strangeness
905  Int_t absPDG = TMath::Abs(part->PdgCode());
906  if ((absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
907  {
908  // Check if particle is in a radius around the jet
909  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
911  continue;
912  else
913  return kTRUE;
914  }
915  }
916  }
917  // As fallback, the MC stack will be tried
918  else if(MCEvent() && (MCEvent()->Stack()))
919  {
920  AliStack* stack = MCEvent()->Stack();
921  // Go through the whole particle stack
922  for(Int_t i=0; i<stack->GetNtrack(); i++)
923  {
924  TParticle *part = stack->Particle(i);
925  if(!part) continue;
926 
927  // Check if the particle has strangeness
928  Int_t absPDG = TMath::Abs(part->GetPdgCode());
929  if ((absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
930  {
931  // Check if particle is in a radius around the jet
932  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
934  continue;
935  else
936  return kTRUE;
937  }
938  }
939  }
940  return kFALSE;
941 
942 }
943 //________________________________________________________________________
945 {
946  // ### Event control plots
948  FillHistogram("hBackgroundPt", fJetsCont->GetRhoVal(), fCent);
949 }
950 
951 //________________________________________________________________________
953 {
954  // ### Jet control plots
955  FillHistogram("hJetPtRaw", jet->Pt(), fCent);
956  FillHistogram("hJetPt", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
957  FillHistogram("hJetPhiEta", jet->Phi(), jet->Eta());
958  FillHistogram("hJetArea", jet->Area(), fCent);
959 
960  // ### Jet constituent plots
961  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
962  {
963  AliVParticle* jetConst = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
964  if(!jetConst) continue;
965 
966  // Constituent eta/phi (relative to jet)
967  Double_t deltaEta = jet->Eta() - jetConst->Eta();
968  Double_t deltaPhi = TMath::Min(TMath::Abs(jet->Phi() - jetConst->Phi()), TMath::TwoPi() - TMath::Abs(jet->Phi() - jetConst->Phi()));
969  if(!((jet->Phi() - jetConst->Phi() < 0) && (jet->Phi() - jetConst->Phi() <= TMath::Pi())))
970  deltaPhi = -deltaPhi;
971 
972  FillHistogram("hConstituentPt", jetConst->Pt(), fCent);
973  FillHistogram("hConstituentPhiEta", deltaPhi, deltaEta);
974  }
975 
976  // ### Random cone / delta pT plots
977  const Int_t kNumRandomConesPerEvent = 4;
978  for(Int_t iCone=0; iCone<kNumRandomConesPerEvent; iCone++)
979  {
980  // Throw random cone
981  Double_t tmpRandConeEta = fJetsCont->GetJetEtaMin() + gRandom->Rndm()*TMath::Abs(fJetsCont->GetJetEtaMax()-fJetsCont->GetJetEtaMin());
982  Double_t tmpRandConePhi = gRandom->Rndm()*TMath::TwoPi();
983  Double_t tmpRandConePt = 0;
984  // Fill pT that is in cone
985  fTracksCont->ResetCurrentID();
986  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
987  if(IsTrackInCone(track, tmpRandConeEta, tmpRandConePhi, fJetsCont->GetJetRadius()))
988  tmpRandConePt += track->Pt();
989 
990  // Fill histograms
991  FillHistogram("hDeltaPt", tmpRandConePt - fJetsCont->GetRhoVal()*fJetsCont->GetJetRadius()*fJetsCont->GetJetRadius()*TMath::Pi(), fCent);
992  }
993 }
994 
995 //________________________________________________________________________
997 {
998  FillHistogram("hTrackPt", track->Pt(), fCent);
999  FillHistogram("hTrackPhi", track->Phi(), fCent);
1000  FillHistogram("hTrackEta", track->Eta(), fCent);
1001  FillHistogram("hTrackEtaPt", track->Eta(), track->Pt());
1002  FillHistogram("hTrackPhiPt", track->Phi(), track->Pt());
1003  FillHistogram("hTrackPhiEta", track->Phi(), track->Eta());
1004 }
1005 
1006 //________________________________________________________________________
1007 void AliAnalysisTaskJetExtractor::AddPIDInformation(AliVParticle* particle, Float_t& sigITS, Float_t& sigTPC, Float_t& sigTOF, Float_t& sigTRD, Short_t& recoPID, Short_t& truePID)
1008 {
1009  truePID = 9;
1010  if(!particle) return;
1011 
1012  // If we have AODs, retrieve particle PID signals
1013  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(particle);
1014 
1015  if(aodtrack)
1016  {
1017  // Get AOD value from reco
1018  recoPID = aodtrack->GetMostProbablePID();
1019  AliAODPid* pidObj = aodtrack->GetDetPid();
1020  if(!pidObj)
1021  return;
1022 
1023  sigITS = pidObj->GetITSsignal();
1024  sigTPC = pidObj->GetTPCsignal();
1025  sigTOF = pidObj->GetTOFsignal();
1026  sigTRD = pidObj->GetTRDsignal();
1027  }
1028 
1029  // Get truth values if we are on MC
1031  {
1032  for(Int_t i=0; i<fTruthParticleArray->GetEntries();i++)
1033  {
1034  AliAODMCParticle* mcParticle = dynamic_cast<AliAODMCParticle*>(fTruthParticleArray->At(i));
1035  if(!mcParticle) continue;
1036 
1037  if (mcParticle->GetLabel() == particle->GetLabel())
1038  {
1039  // Use same convention as PID in AODs
1040  if(TMath::Abs(mcParticle->PdgCode()) == 2212) // proton
1041  truePID = 4;
1042  else if (TMath::Abs(mcParticle->PdgCode()) == 211) // pion
1043  truePID = 2;
1044  else if (TMath::Abs(mcParticle->PdgCode()) == 321) // kaon
1045  truePID = 3;
1046  else if (TMath::Abs(mcParticle->PdgCode()) == 11) // electron
1047  truePID = 0;
1048  else if (TMath::Abs(mcParticle->PdgCode()) == 13) // muon
1049  truePID = 1;
1050  else if (TMath::Abs(mcParticle->PdgCode()) == 700201) // deuteron
1051  truePID = 5;
1052  else if (TMath::Abs(mcParticle->PdgCode()) == 700301) // triton
1053  truePID = 6;
1054  else if (TMath::Abs(mcParticle->PdgCode()) == 700302) // He3
1055  truePID = 7;
1056  else if (TMath::Abs(mcParticle->PdgCode()) == 700202) // alpha
1057  truePID = 8;
1058  else
1059  truePID = 9;
1060 
1061  break;
1062  }
1063  }
1064  }
1065 }
1066 
1067 
1068 //########################################################################
1069 // HELPERS
1070 //########################################################################
1071 
1072 //________________________________________________________________________
1073 inline Bool_t AliAnalysisTaskJetExtractor::IsTrackInCone(AliVParticle* track, Double_t eta, Double_t phi, Double_t radius)
1074 {
1075  // This is to use a full cone in phi even at the edges of phi (2pi -> 0) (0 -> 2pi)
1076  Double_t trackPhi = 0.0;
1077  if (track->Phi() > (TMath::TwoPi() - (radius-phi)))
1078  trackPhi = track->Phi() - TMath::TwoPi();
1079  else if (track->Phi() < (phi+radius - TMath::TwoPi()))
1080  trackPhi = track->Phi() + TMath::TwoPi();
1081  else
1082  trackPhi = track->Phi();
1083 
1084  if ( TMath::Abs(trackPhi-phi)*TMath::Abs(trackPhi-phi) + TMath::Abs(track->Eta()-eta)*TMath::Abs(track->Eta()-eta) <= radius*radius)
1085  return kTRUE;
1086 
1087  return kFALSE;
1088 }
1089 
1090 //________________________________________________________________________
1092 {
1093  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1094  if(!tmpHist)
1095  {
1096  AliError(Form("Cannot find histogram <%s> ",key)) ;
1097  return;
1098  }
1099 
1100  tmpHist->Fill(x);
1101 }
1102 
1103 //________________________________________________________________________
1105 {
1106  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1107  if(!tmpHist)
1108  {
1109  AliError(Form("Cannot find histogram <%s> ",key));
1110  return;
1111  }
1112 
1113  if (tmpHist->IsA()->GetBaseClass("TH1"))
1114  static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
1115  else if (tmpHist->IsA()->GetBaseClass("TH2"))
1116  static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
1117 }
1118 
1119 //________________________________________________________________________
1121 {
1122  TH2* tmpHist = static_cast<TH2*>(fOutput->FindObject(key));
1123  if(!tmpHist)
1124  {
1125  AliError(Form("Cannot find histogram <%s> ",key));
1126  return;
1127  }
1128 
1129  tmpHist->Fill(x,y,add);
1130 }
1131 
1132 //________________________________________________________________________
1134 {
1135  TH3* tmpHist = static_cast<TH3*>(fOutput->FindObject(key));
1136  if(!tmpHist)
1137  {
1138  AliError(Form("Cannot find histogram <%s> ",key));
1139  return;
1140  }
1141 
1142  if(add)
1143  tmpHist->Fill(x,y,z,add);
1144  else
1145  tmpHist->Fill(x,y,z);
1146 }
1147 
1148 
1149 //________________________________________________________________________
1150 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)
1151 {
1152  T* tmpHist = new T(name, title, xBins, xMin, xMax);
1153 
1154  tmpHist->GetXaxis()->SetTitle(xTitle);
1155  tmpHist->GetYaxis()->SetTitle(yTitle);
1156  tmpHist->SetOption(options);
1157  tmpHist->SetMarkerStyle(kFullCircle);
1158  tmpHist->Sumw2();
1159 
1160  fOutput->Add(tmpHist);
1161 
1162  return tmpHist;
1163 }
1164 
1165 //________________________________________________________________________
1166 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)
1167 {
1168  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax);
1169  tmpHist->GetXaxis()->SetTitle(xTitle);
1170  tmpHist->GetYaxis()->SetTitle(yTitle);
1171  tmpHist->GetZaxis()->SetTitle(zTitle);
1172  tmpHist->SetOption(options);
1173  tmpHist->SetMarkerStyle(kFullCircle);
1174  tmpHist->Sumw2();
1175 
1176  fOutput->Add(tmpHist);
1177 
1178  return tmpHist;
1179 }
1180 
1181 //________________________________________________________________________
1182 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)
1183 {
1184  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax, zBins, zMin, zMax);
1185  tmpHist->GetXaxis()->SetTitle(xTitle);
1186  tmpHist->GetYaxis()->SetTitle(yTitle);
1187  tmpHist->GetZaxis()->SetTitle(zTitle);
1188  tmpHist->SetOption(options);
1189  tmpHist->SetMarkerStyle(kFullCircle);
1190  tmpHist->Sumw2();
1191 
1192  fOutput->Add(tmpHist);
1193 
1194  return tmpHist;
1195 }
1196 
1197 //________________________________________________________________________
1199 {
1200  // Called once at the end of the analysis.
1201 }
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 SetSaveConstituents(Bool_t val)
Double_t GetRhoVal() 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 GetJetType(AliEmcalJet *jet, Int_t &typeHM, Int_t &typePM, Int_t &typeIC)
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
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
Double_t Phi() const
Definition: AliEmcalJet.h:117
std::vector< Float_t > GetExtractionPercentages()
Bool_t fSaveMCInformation
save MC information
Double_t mass
Container with name, TClonesArray and cuts for particles.
centrality
Bool_t IsEventSelected()
Performing event selection.
Double_t fSecondaryVertexMaxChi2
Max chi2 of secondary vertex (others will be discarded)
Double_t GetJetEtaMax() const
Float_t GetPartonPhi7() const
Double_t fSecondaryVertexMaxDispersion
Max dispersion of secondary vertex (others will be discarded)
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_t fSaveJetShapes
save jet shapes
Bool_t fSaveConstituentPID
save arrays of constituent PID parameters
AliTrackContainer * fTracksCont
! Tracks
AliStack * stack
const Int_t kMaxNumConstituents
Bool_t IsTrackInCone(AliVParticle *track, Double_t eta, Double_t phi, Double_t radius)
TRandom * gRandom
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")
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)
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
Float_t * fBuffer_Const_Phi
! array buffer
Int_t GetPartonFlag6() const
Float_t fBuffer_Jet_MC_MatchedPt
! array buffer
Bool_t fSaveTriggerTracks
save event trigger track
AliParticleContainer * GetParticleContainer() const
std::vector< Float_t > fExtractionPercentages
Percentages which will be extracted for a given pT bin.
TClonesArray * fTruthParticleArray
! Array of MC particles in event (mcparticles)
void FillHistogram3D(const char *key, Double_t x, Double_t y, Double_t z, Double_t add=0)
int Int_t
Definition: External.C:63
Int_t fBuffer_Jet_MC_MotherParton
! array buffer
Float_t fBuffer_Event_PtHard
! array buffer
void SetSaveJetShapes(Bool_t val)
unsigned int UInt_t
Definition: External.C:33
Float_t * fBuffer_Const_Charge
! array buffer
float Float_t
Definition: External.C:68
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")
Float_t fBuffer_JetPhi
! array buffer
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
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.
AliEmcalJet * GetNextAcceptJet()
AliHFJetsTaggingVertex * fVtxTagger
! class for sec. vertexing
short Short_t
Definition: External.C:23
TString fTruthJetsArrayName
Array name for particle-level jets.
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliRDHFJetsCutsVertex * fVertexerCuts
Cuts used for the vertexer (given in add task macro)
Int_t fEventCut_TriggerTrackMinLabel
Event requirement, trigger track min label (can be used to selected embedded particles) ...
void SetSaveEventProperties(Bool_t val)
Float_t fPtHard
!event -hard
void FillHistogram(const char *key, Double_t x)
void FillTrackControlHistograms(AliVTrack *track)
Bool_t AddJetToTree(AliEmcalJet *jet, Float_t bgrdDensity=0, Float_t vertexX=0, Float_t vertexY=0, Float_t vertexZ=0, Float_t centrality=0, Long64_t eventID=0, Float_t magField=0, AliTrackContainer *trackCont=0, Int_t motherParton=0, Int_t motherHadron=0, Int_t partonInitialCollision=0, Float_t matchedPt=0, Float_t truePtFraction=0, Float_t ptHard=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< Short_t > &trackPID_Truth=DEFAULT_VECTOR_SHORT, 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 GetJetRadius() const
AliEmcalList * fOutput
!output list
Float_t fBuffer_JetPt
! array buffer
Definition: External.C:220
void ReconstructSecondaryVertices(const AliVVertex *primVtx, const AliEmcalJet *jet, std::vector< Float_t > &secVtx_X, std::vector< Float_t > &secVtx_Y, std::vector< Float_t > &secVtx_Z, std::vector< Float_t > &secVtx_Mass, std::vector< Float_t > &secVtx_Lxy, std::vector< Float_t > &secVtx_SigmaLxy, std::vector< Float_t > &secVtx_Chi2, std::vector< Float_t > &secVtx_Dispersion)
Float_t fBuffer_Event_Vertex_Y
! array buffer
Bool_t fCalculateSecondaryVertices
Calculate the secondary vertices (instead of loading)
void AddPIDInformation(AliVParticle *particle, Float_t &sigITS, Float_t &sigTPC, Float_t &sigTOF, Float_t &sigTRD, Short_t &recoPID, Short_t &truePID)
Float_t * fBuffer_Const_ProdVtx_Y
! array buffer
void GetJetTruePt(AliEmcalJet *jet, Double_t &matchedJetPt, Double_t &truePtFraction)
Bool_t fSaveConstituents
save arrays of constituent basic properties
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
void SetMakeGeneralHistograms(Bool_t g)
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
Float_t * fBuffer_Const_Eta
! array buffer
Float_t fBuffer_JetEta
! array buffer
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
virtual AliVParticle * GetNextAcceptParticle()
bool Bool_t
Definition: External.C:53
Double_t yMin
Int_t GetNAcceptedParticles() const
Int_t fBuffer_Jet_MC_MotherHadron
! array buffer
Double_t fEventCut_TriggerTrackMaxPt
Event requirement, trigger track max pT.
std::vector< Float_t > fTriggerTracks_Eta
If trigger track requi. is used -> save eta.
Float_t fBuffer_Shape_Mass
! array buffer
Double_t M() const
Definition: AliEmcalJet.h:120
void FillJetControlHistograms(AliEmcalJet *jet)
Float_t * fBuffer_Const_ProdVtx_Z
! array buffer
Float_t fBuffer_Jet_MC_TruePtFraction
! array buffer
Class managing creation of a tree containing jets.
Definition: External.C:196
Double_t yMax
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