AliPhysics  d9e9949 (d9e9949)
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 "AliParticleContainer.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  AliParticleContainer* 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  fEventPercentage(1.0),
351  fTruthMinLabel(0),
352  fTruthMaxLabel(100000),
353  fEventCut_TriggerTrackMinPt(0),
354  fEventCut_TriggerTrackMaxPt(0),
355  fEventCut_TriggerTrackMinLabel(-9999999),
356  fEventCut_TriggerTrackMaxLabel(+9999999),
357  fJetsCont(0),
358  fTracksCont(0),
359  fTruthParticleArray(0),
360  fTruthJetsArrayName(""),
361  fTruthJetsRhoName(""),
362  fTruthParticleArrayName("mcparticles"),
363  fTriggerTracks_Pt(),
364  fTriggerTracks_Eta(),
365  fTriggerTracks_Phi()
366 {
368  fJetTree = new AliEmcalJetTree(GetName());
371  fJetTree->SetSaveJetShapes(kTRUE);
372  DefineOutput(2, TTree::Class());
373 }
374 
375 //________________________________________________________________________
377  AliAnalysisTaskEmcalJet(name, kTRUE),
378  fJetTree(0),
379  fVtxTagger(0),
385  fVertexerCuts(0),
387  fEventPercentage(1.0),
388  fTruthMinLabel(0),
389  fTruthMaxLabel(100000),
394  fJetsCont(0),
395  fTracksCont(0),
398  fTruthJetsRhoName(""),
399  fTruthParticleArrayName("mcparticles"),
403 {
405  fJetTree = new AliEmcalJetTree(GetName());
408  fJetTree->SetSaveJetShapes(kTRUE);
409  DefineOutput(2, TTree::Class());
410 }
411 
412 //________________________________________________________________________
414 {
415  // Destructor.
416 }
417 
418 //________________________________________________________________________
420 {
422 
423  // ### Basic container settings
425  if(!fJetsCont)
426  AliFatal("Jet input container not found!");
427  fJetsCont->PrintCuts();
429  if(!fTracksCont)
430  AliFatal("Particle input container not found attached to jets!");
431 
432  // Activate saving of trigger tracks if this is demanded
435 
436  // ### Initialize the jet tree (settings must all be given at this stage)
438  OpenFile(2);
439  PostData(2, fJetTree->GetTreePointer());
440 
441  PrintConfig();
442 
443  // ### Add control histograms (already some created in base task)
444  AddHistogram2D<TH2D>("hTrackCount", "Number of tracks in acceptance vs. centrality", "COLZ", 500, 0., 5000., 100, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}");
445  AddHistogram2D<TH2D>("hBackgroundPt", "Background p_{T} distribution", "", 150, 0., 150., 100, 0, 100, "Background p_{T} (GeV/c)", "Centrality", "dN^{Events}/dp_{T}");
446 
447  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}");
448  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}");
449  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}");
450  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");
451  AddHistogram2D<TH2D>("hJetArea", "Jet area", "COLZ", 200, 0., 2., 100, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
452 
453  AddHistogram2D<TH2D>("hDeltaPt", "#delta p_{T} distribution", "", 400, -100., 300., 100, 0, 100, "p_{T, cone} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
454 
455  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}");
456  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");
457 
458  AddHistogram1D<TH1D>("hExtractionPercentage", "Extracted jets p_{T} distribution (background subtracted)", "e1p", 400, -100., 300., "p_{T, jet} (GeV/c)", "Extracted percentage");
459 
460  // Track QA plots
461  AddHistogram2D<TH2D>("hTrackPt", "Tracks p_{T} distribution", "", 300, 0., 300., 100, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
462  AddHistogram2D<TH2D>("hTrackPhi", "Track angular distribution in #phi", "LEGO2", 180, 0., 2*TMath::Pi(), 100, 0, 100, "#phi", "Centrality", "dN^{Tracks}/(d#phi)");
463  AddHistogram2D<TH2D>("hTrackEta", "Track angular distribution in #eta", "LEGO2", 100, -2.5, 2.5, 100, 0, 100, "#eta", "Centrality", "dN^{Tracks}/(d#eta)");
464  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");
465 
466  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})");
467  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})");
468 }
469 
470 
471 //________________________________________________________________________
473 {
475  if (fTruthParticleArrayName != "")
476  fTruthParticleArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTruthParticleArrayName.Data()));
477 
478  // ### Prepare vertexer
480  {
481  if(!fVertexerCuts)
482  AliFatal("VertexerCuts not given but secondary vertex calculation turned on.");
483  fVtxTagger = new AliHFJetsTaggingVertex();
484  fVtxTagger->SetCuts(fVertexerCuts);
485  }
486 
487  // ### Save tree extraction percentage to histogram
488  std::vector<Float_t> extractionPtBins = fJetTree->GetExtractionPercentagePtBins();
489  std::vector<Float_t> extractionPercentages = fJetTree->GetExtractionPercentages();
490 
491  for(Int_t i=0; i<static_cast<Int_t>(extractionPercentages.size()); i++)
492  {
493  Double_t percentage = extractionPercentages[i];
494  for(Int_t pt=static_cast<Int_t>(extractionPtBins[i*2]); pt<static_cast<Int_t>(extractionPtBins[i*2+1]); pt++)
495  {
496  FillHistogram("hExtractionPercentage", pt, percentage);
497  }
498  }
499 
500  // ### Add HIJING Ncoll histogram (in case we need it)
501  if(MCEvent())
502  {
503  if(dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()))
504  AddHistogram1D<TH1D>("hHijingNcoll", "N_{coll} from HIJING", "e1p", 1000, 0., 5000, "N_{coll}", "dN^{Events}/dN^{N_{coll}}");
505  }
506 }
507 
508 //________________________________________________________________________
510 {
511  // ################################### EVENT SELECTION
512 
513  if(!IsTriggerTrackInEvent())
514  return kFALSE;
515  if(gRandom->Rndm() >= fEventPercentage)
516  return kFALSE;
517 
518  // ################################### EVENT PROPERTIES
519 
521 
522  // Load vertex if possible
523  Double_t vtxX = 0;
524  Double_t vtxY = 0;
525  Double_t vtxZ = 0;
526  Long64_t eventID = 0;
527  const AliVVertex* myVertex = InputEvent()->GetPrimaryVertex();
528  if(!myVertex && MCEvent())
529  myVertex = MCEvent()->GetPrimaryVertex();
530  if(myVertex)
531  {
532  vtxX = myVertex->GetX();
533  vtxY = myVertex->GetY();
534  vtxZ = myVertex->GetZ();
535  }
536  // Get event ID from header
537  AliVHeader* eventIDHeader = InputEvent()->GetHeader();
538  if(eventIDHeader)
539  eventID = eventIDHeader->GetEventIdAsLong();
540 
541  // If available, get Ncoll from HIJING
542  if(MCEvent())
543  {
544  if(AliGenHijingEventHeader* mcHeader = dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()))
545  {
546  Double_t ncoll = mcHeader->NN() + mcHeader->NNw() + mcHeader->NwN() + mcHeader->NwNw();
547  FillHistogram("hHijingNcoll", ncoll);
548  }
549  }
550 
551  // ################################### MAIN JET LOOP
552  fJetsCont->ResetCurrentID();
553  while(AliEmcalJet *jet = fJetsCont->GetNextAcceptJet())
554  {
556 
557  Double_t matchedJetPt = 0;
558  Double_t truePtFraction = 0;
559  Int_t currentJetType_HM = 0;
560  Int_t currentJetType_PM = 0;
561  Int_t currentJetType_IC = 0;
563  {
564  // Get jet type from MC (hadron matching, parton matching definition - for HF jets)
565  GetJetType(jet, currentJetType_HM, currentJetType_PM, currentJetType_IC);
566  // Get true pT estimators
567  GetJetTruePt(jet, matchedJetPt, truePtFraction);
568  }
569 
570  // ### CONSTITUENT LOOP: Retrieve PID values + impact parameters
571  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;
572  std::vector<Float_t> vec_d0; std::vector<Float_t> vec_d0cov; std::vector<Float_t> vec_z0; std::vector<Float_t> vec_z0cov;
573 
575  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
576  {
577  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
578  if(!particle) continue;
580  {
581  Float_t sigITS = 0; Float_t sigTPC = 0; Float_t sigTOF = 0; Float_t sigTRD = 0; Short_t recoPID = 0; Short_t truePID = 0;
582  AddPIDInformation(particle, sigITS, sigTPC, sigTOF, sigTRD, recoPID, truePID);
583  vecSigITS.push_back(sigITS); vecSigTPC.push_back(sigTPC); vecSigTOF.push_back(sigTOF); vecSigTRD.push_back(sigTRD); vecRecoPID.push_back(recoPID); vecTruePID.push_back(truePID);
584  }
586  {
587  Float_t d0 = 0; Float_t d0cov = 0; Float_t z0 = 0; Float_t z0cov = 0;
588  GetTrackImpactParameters(myVertex, dynamic_cast<AliAODTrack*>(particle), d0, d0cov, z0, z0cov);
589  vec_d0.push_back(d0); vec_d0cov.push_back(d0cov); vec_z0.push_back(z0); vec_z0cov.push_back(z0cov);
590  }
591  }
592 
593  // Reconstruct secondary vertices
594  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;
596  ReconstructSecondaryVertices(myVertex, jet, secVtx_X, secVtx_Y, secVtx_Z, secVtx_Mass, secVtx_Lxy, secVtx_SigmaLxy, secVtx_Chi2, secVtx_Dispersion);
597 
598  // Now change trigger tracks eta/phi to dEta/dPhi relative to jet
599  std::vector<Float_t> triggerTracks_dEta(fTriggerTracks_Eta);
600  std::vector<Float_t> triggerTracks_dPhi(fTriggerTracks_Phi);
601  for(UInt_t i=0; i<triggerTracks_dEta.size(); i++)
602  {
603  triggerTracks_dEta[i] = jet->Eta() - fTriggerTracks_Eta[i];
604  triggerTracks_dPhi[i] = TMath::Min(TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]), TMath::TwoPi() - TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]));
605  if( ((TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]) <= TMath::Pi()) && (jet->Phi() - fTriggerTracks_Phi[i] > 0))
606  || ((TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]) > TMath::Pi()) && (jet->Phi() - fTriggerTracks_Phi[i] < 0)) )
607  triggerTracks_dPhi[i] = -triggerTracks_dPhi[i];
608  }
609 
610  // Fill jet to tree
611  Bool_t accepted = fJetTree->AddJetToTree(jet, fJetsCont->GetRhoVal(), vtxX, vtxY, vtxZ, fCent, eventID, InputEvent()->GetMagneticField(), fTracksCont,
612  currentJetType_PM,currentJetType_HM,currentJetType_IC,matchedJetPt,truePtFraction,fPtHard,
613  vecSigITS, vecSigTPC, vecSigTOF, vecSigTRD, vecRecoPID, vecTruePID,
614  fTriggerTracks_Pt, triggerTracks_dEta, triggerTracks_dPhi,
615  vec_d0, vec_z0, vec_d0cov, vec_z0cov,
616  secVtx_X, secVtx_Y, secVtx_Z, secVtx_Mass, secVtx_Lxy, secVtx_SigmaLxy, secVtx_Chi2, secVtx_Dispersion);
617 
618  if(accepted)
619  FillHistogram("hJetPtExtracted", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
620  }
621 
622  // ################################### PARTICLE PROPERTIES
623  fTracksCont->ResetCurrentID();
624  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
626 
627  return kTRUE;
628 }
629 
630 //________________________________________________________________________
632 {
633  // Cut for trigger track requirement
635  {
636  // Clear vector of trigger tracks
637  fTriggerTracks_Pt.clear();
638  fTriggerTracks_Eta.clear();
639  fTriggerTracks_Phi.clear();
640 
641  // Go through all tracks and check whether trigger tracks can be found
642  fTracksCont->ResetCurrentID();
643  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
644  {
645  if( (track->GetLabel() >= fEventCut_TriggerTrackMinLabel) && (track->GetLabel() < fEventCut_TriggerTrackMaxLabel) )
646  if( (track->Pt() >= fEventCut_TriggerTrackMinPt) && (track->Pt() < fEventCut_TriggerTrackMaxPt) )
647  {
648  fTriggerTracks_Pt.push_back(track->Pt());
649  fTriggerTracks_Eta.push_back(track->Eta());
650  fTriggerTracks_Phi.push_back(track->Phi());
651  }
652  }
653  // No particle has been found that fulfills requirement -> Do not accept event
654  if(fTriggerTracks_Pt.size() == 0)
655  return kFALSE;
656  }
657  return kTRUE;
658 }
659 
660 
661 //________________________________________________________________________
662 void AliAnalysisTaskJetExtractor::GetJetTruePt(AliEmcalJet* jet, Double_t& matchedJetPt, Double_t& truePtFraction)
663 {
664  // #################################################################################
665  // ##### METHOD 1: If fTruthJetsArrayName is set, a matching jet is searched for
666  Double_t bestMatchDeltaR = 999.;
667  if(fTruthJetsArrayName != "")
668  {
669  // "True" background
670  AliRhoParameter* rho = static_cast<AliRhoParameter*>(InputEvent()->FindListObject(fTruthJetsRhoName.Data()));
671  Double_t trueRho = 0;
672  if(rho)
673  trueRho = rho->GetVal();
674 
675  TClonesArray* truthArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fTruthJetsArrayName.Data())));
676 
677  // Loop over all true jets to find the best match
678  matchedJetPt = 0;
679  if(truthArray)
680  for(Int_t i=0; i<truthArray->GetEntries(); i++)
681  {
682  AliEmcalJet* truthJet = static_cast<AliEmcalJet*>(truthArray->At(i));
683  if(truthJet->Pt() < 0.15)
684  continue;
685 
686  Double_t deltaEta = (truthJet->Eta()-jet->Eta());
687  Double_t deltaPhi = TMath::Min(TMath::Abs(truthJet->Phi()-jet->Phi()),TMath::TwoPi() - TMath::Abs(truthJet->Phi()-jet->Phi()));
688  Double_t deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
689 
690  // Cut jets too far away
691  if (deltaR > fTrueJetMatchingRadius)
692  continue;
693 
694  // Search for the best match
695  if(deltaR < bestMatchDeltaR)
696  {
697  bestMatchDeltaR = deltaR;
698  matchedJetPt = truthJet->Pt() - truthJet->Area()* trueRho;
699  }
700  }
701  }
702 
703  // #################################################################################
704  // ##### METHOD 2: Calculate fraction of "true" pT -- pT which is not from a toy
705  Double_t pt_nonMC = 0.;
706  Double_t pt_all = 0.;
707  truePtFraction = 0;
708 
709  for(Int_t iConst = 0; iConst < jet->GetNumberOfTracks(); iConst++)
710  {
711  // Loop over all valid jet constituents
712  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(iConst, fTracksCont->GetArray()));
713  if(!particle) continue;
714  if(particle->Pt() < 1e-6) continue;
715 
716  // Particles marked w/ labels within label range are considered from toy
717  if( (particle->GetLabel() >= fTruthMinLabel) && (particle->GetLabel() < fTruthMaxLabel))
718  pt_nonMC += particle->Pt();
719  pt_all += particle->Pt();
720  }
721  if(pt_all)
722  truePtFraction = (pt_nonMC/pt_all);
723 
724 }
725 
726 
727 //________________________________________________________________________
729 {
731 
733  return;
734 
735  typeHM = 0;
736  typePM = 0;
737 
738  AliAODMCParticle* parton[2];
739  parton[0] = (AliAODMCParticle*) fVtxTagger->IsMCJetParton(fTruthParticleArray, jet, radius); // method 1 (parton matching)
740  parton[1] = (AliAODMCParticle*) fVtxTagger->IsMCJetMeson(fTruthParticleArray, jet, radius); // method 2 (hadron matching)
741 
742  if (parton[0]) {
743  Int_t pdg = TMath::Abs(parton[0]->PdgCode());
744  typePM = pdg;
745  }
746 
747  if (!parton[1])
748  {
749  // No HF jet (according to hadron matching) -- now also separate jets in udg (1) and s-jets (3)
750  if(IsStrangeJet(jet))
751  typeHM = 3;
752  else
753  typeHM = 1;
754  }
755  else {
756  Int_t pdg = TMath::Abs(parton[1]->PdgCode());
757  if(fVtxTagger->IsDMeson(pdg)) typeHM = 4;
758  else if (fVtxTagger->IsBMeson(pdg)) typeHM = 5;
759  }
760 
761  // Set flavour of AliEmcalJet object (set ith bit while i corresponds to type)
763  jet->AddFlavourTag(static_cast<Int_t>(TMath::Power(2, typeHM)));
764 
765 
766  const AliEmcalPythiaInfo* partonsInfo = GetPythiaInfo();
767  typeIC = 0;
768  if (partonsInfo)
769  {
770  // Get primary partons directions
771  Double_t parton1phi = partonsInfo->GetPartonPhi6();
772  Double_t parton1eta = partonsInfo->GetPartonEta6();
773  Double_t parton2phi = partonsInfo->GetPartonPhi7();
774  Double_t parton2eta = partonsInfo->GetPartonEta7();
775 
776 
777  Double_t delta1Eta = (parton1eta-jet->Eta());
778  Double_t delta1Phi = TMath::Min(TMath::Abs(parton1phi-jet->Phi()),TMath::TwoPi() - TMath::Abs(parton1phi-jet->Phi()));
779  Double_t delta1R = TMath::Sqrt(delta1Eta*delta1Eta + delta1Phi*delta1Phi);
780  Double_t delta2Eta = (parton2eta-jet->Eta());
781  Double_t delta2Phi = TMath::Min(TMath::Abs(parton2phi-jet->Phi()),TMath::TwoPi() - TMath::Abs(parton2phi-jet->Phi()));
782  Double_t delta2R = TMath::Sqrt(delta2Eta*delta2Eta + delta2Phi*delta2Phi);
783 
784  // Check if one of the partons if closer than matching criterion
785  Bool_t matched = (delta1R < fJetsCont->GetJetRadius()/2.) || (delta2R < fJetsCont->GetJetRadius()/2.);
786 
787  // Matching criterion fulfilled -> Set flag to closest
788  if(matched)
789  {
790  if(delta1R < delta2R)
791  typeIC = partonsInfo->GetPartonFlag6();
792  else
793  typeIC = partonsInfo->GetPartonFlag7();
794  }
795  }
796 }
797 
798 //________________________________________________________________________
799 void AliAnalysisTaskJetExtractor::GetTrackImpactParameters(const AliVVertex* vtx, AliAODTrack* track, Float_t& d0, Float_t& d0cov, Float_t& z0, Float_t& z0cov)
800 {
801  if (track)
802  {
803  Double_t d0rphiz[2],covd0[3];
804  Bool_t isDCA=track->PropagateToDCA(vtx,InputEvent()->GetMagneticField(),3.0,d0rphiz,covd0);
805  if(isDCA)
806  {
807  d0 = d0rphiz[0];
808  z0 = d0rphiz[1];
809  d0cov = covd0[0];
810  z0cov = covd0[2];
811  }
812  }
813 }
814 
815 //________________________________________________________________________
816 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)
817 {
818  if(!primVtx)
819  return;
820 
821  // Create ESD vertex from the existing AliVVertex
822  Double_t vtxPos[3] = {primVtx->GetX(), primVtx->GetY(), primVtx->GetZ()};
823  Double_t covMatrix[6] = {0};
824  primVtx->GetCovarianceMatrix(covMatrix);
825  AliESDVertex* esdVtx = new AliESDVertex(vtxPos, covMatrix, primVtx->GetChi2(), primVtx->GetNContributors());
826 
827  TClonesArray* secVertexArr = 0;
828  vctr_pair_dbl_int arrDispersion;
829  arrDispersion.reserve(5);
831  {
832  //###########################################################################
833  // ********* Calculate secondary vertices
834  // Derived from AliAnalysisTaskEmcalJetBtagSV
835  secVertexArr = new TClonesArray("AliAODVertex");
836  Int_t nDauRejCount = 0;
837  Int_t nVtx = fVtxTagger->FindVertices(jet,
838  fTracksCont->GetArray(),
839  (AliAODEvent*)InputEvent(),
840  esdVtx,
841  InputEvent()->GetMagneticField(),
842  secVertexArr,
843  0,
844  arrDispersion,
845  nDauRejCount);
846 
847 
848  if(nVtx < 0)
849  {
850  secVertexArr->Clear();
851  delete secVertexArr;
852  return;
853  }
854  //###########################################################################
855  }
856  else // Load HF vertex branch from AOD event, if possible
857  {
858  secVertexArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("VerticesHF"));
859  if(!secVertexArr)
860  return;
861  }
862 
863  // Loop over all potential secondary vertices
864  for(Int_t i=0; i<secVertexArr->GetEntriesFast(); i++)
865  {
866  AliAODVertex* secVtx = (AliAODVertex*)(secVertexArr->UncheckedAt(i));
868  if((strcmp(secVtx->GetParent()->ClassName(), "AliAODRecoDecayHF3Prong")))
869  continue;
870 
871  // Calculate vtx distance
872  Double_t effX = secVtx->GetX() - esdVtx->GetX();
873  Double_t effY = secVtx->GetY() - esdVtx->GetY();
874  //Double_t effZ = secVtx->GetZ() - esdVtx->GetZ();
875 
876  // ##### Vertex properties
877  // vertex dispersion
878  Double_t dispersion = arrDispersion[i].first;
879 
880  // invariant mass
881  Double_t mass = fVtxTagger->GetVertexInvariantMass(secVtx);
882 
883  // signed length
884  Double_t Lxy = TMath::Sqrt(effX*effX + effY*effY);
885  Double_t jetP[3]; jet->PxPyPz(jetP);
886  Double_t signLxy = effX * jetP[0] + effY * jetP[1];
887  if (signLxy < 0.) Lxy *= -1.;
888 
889  Double_t sigmaLxy = 0;
890  AliAODVertex* aodVtx = (AliAODVertex*)(primVtx);
891  if (aodVtx)
892  sigmaLxy = aodVtx->ErrorDistanceXYToVertex(secVtx);
893 
894  // Add secondary vertices if they fulfill the conditions
895  if( (dispersion > fSecondaryVertexMaxDispersion) || (TMath::Abs(secVtx->GetChi2perNDF()) > fSecondaryVertexMaxChi2) )
896  continue;
897 
898  secVtx_X.push_back(secVtx->GetX()); secVtx_Y.push_back(secVtx->GetY()); secVtx_Z.push_back(secVtx->GetZ()); secVtx_Chi2.push_back(secVtx->GetChi2perNDF());
899  secVtx_Dispersion.push_back(dispersion); secVtx_Mass.push_back(mass); secVtx_Lxy.push_back(Lxy); secVtx_SigmaLxy.push_back(sigmaLxy);
900  }
901 
903  {
904  secVertexArr->Clear();
905  delete secVertexArr;
906  }
907  delete esdVtx;
908 }
909 
910 //________________________________________________________________________
912 {
913  // Do hadron matching jet type tagging using mcparticles
914  // ... if not explicitly deactivated
916  {
917  for(Int_t i=0; i<fTruthParticleArray->GetEntries();i++)
918  {
919  AliAODMCParticle* part = (AliAODMCParticle*)fTruthParticleArray->At(i);
920  if(!part) continue;
921 
922  // Check if the particle has strangeness
923  Int_t absPDG = TMath::Abs(part->PdgCode());
924  if ((absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
925  {
926  // Check if particle is in a radius around the jet
927  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
929  continue;
930  else
931  return kTRUE;
932  }
933  }
934  }
935  // As fallback, the MC stack will be tried
936  else if(MCEvent() && (MCEvent()->Stack()))
937  {
938  AliStack* stack = MCEvent()->Stack();
939  // Go through the whole particle stack
940  for(Int_t i=0; i<stack->GetNtrack(); i++)
941  {
942  TParticle *part = stack->Particle(i);
943  if(!part) continue;
944 
945  // Check if the particle has strangeness
946  Int_t absPDG = TMath::Abs(part->GetPdgCode());
947  if ((absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
948  {
949  // Check if particle is in a radius around the jet
950  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
952  continue;
953  else
954  return kTRUE;
955  }
956  }
957  }
958  return kFALSE;
959 
960 }
961 //________________________________________________________________________
963 {
964  // ### Event control plots
966  FillHistogram("hBackgroundPt", fJetsCont->GetRhoVal(), fCent);
967 }
968 
969 //________________________________________________________________________
971 {
972  // ### Jet control plots
973  FillHistogram("hJetPtRaw", jet->Pt(), fCent);
974  FillHistogram("hJetPt", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
975  FillHistogram("hJetPhiEta", jet->Phi(), jet->Eta());
976  FillHistogram("hJetArea", jet->Area(), fCent);
977 
978  // ### Jet constituent plots
979  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
980  {
981  AliVParticle* jetConst = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
982  if(!jetConst) continue;
983 
984  // Constituent eta/phi (relative to jet)
985  Double_t deltaEta = jet->Eta() - jetConst->Eta();
986  Double_t deltaPhi = TMath::Min(TMath::Abs(jet->Phi() - jetConst->Phi()), TMath::TwoPi() - TMath::Abs(jet->Phi() - jetConst->Phi()));
987  if(!((jet->Phi() - jetConst->Phi() < 0) && (jet->Phi() - jetConst->Phi() <= TMath::Pi())))
988  deltaPhi = -deltaPhi;
989 
990  FillHistogram("hConstituentPt", jetConst->Pt(), fCent);
991  FillHistogram("hConstituentPhiEta", deltaPhi, deltaEta);
992  }
993 
994  // ### Random cone / delta pT plots
995  const Int_t kNumRandomConesPerEvent = 4;
996  for(Int_t iCone=0; iCone<kNumRandomConesPerEvent; iCone++)
997  {
998  // Throw random cone
999  Double_t tmpRandConeEta = fJetsCont->GetJetEtaMin() + gRandom->Rndm()*TMath::Abs(fJetsCont->GetJetEtaMax()-fJetsCont->GetJetEtaMin());
1000  Double_t tmpRandConePhi = gRandom->Rndm()*TMath::TwoPi();
1001  Double_t tmpRandConePt = 0;
1002  // Fill pT that is in cone
1003  fTracksCont->ResetCurrentID();
1004  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
1005  if(IsTrackInCone(track, tmpRandConeEta, tmpRandConePhi, fJetsCont->GetJetRadius()))
1006  tmpRandConePt += track->Pt();
1007 
1008  // Fill histograms
1009  FillHistogram("hDeltaPt", tmpRandConePt - fJetsCont->GetRhoVal()*fJetsCont->GetJetRadius()*fJetsCont->GetJetRadius()*TMath::Pi(), fCent);
1010  }
1011 }
1012 
1013 //________________________________________________________________________
1015 {
1016  FillHistogram("hTrackPt", track->Pt(), fCent);
1017  FillHistogram("hTrackPhi", track->Phi(), fCent);
1018  FillHistogram("hTrackEta", track->Eta(), fCent);
1019  FillHistogram("hTrackEtaPt", track->Eta(), track->Pt());
1020  FillHistogram("hTrackPhiPt", track->Phi(), track->Pt());
1021  FillHistogram("hTrackPhiEta", track->Phi(), track->Eta());
1022 }
1023 
1024 //________________________________________________________________________
1025 void AliAnalysisTaskJetExtractor::AddPIDInformation(AliVParticle* particle, Float_t& sigITS, Float_t& sigTPC, Float_t& sigTOF, Float_t& sigTRD, Short_t& recoPID, Short_t& truePID)
1026 {
1027  truePID = 9;
1028  if(!particle) return;
1029 
1030  // If we have AODs, retrieve particle PID signals
1031  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(particle);
1032 
1033  if(aodtrack)
1034  {
1035  // Get AOD value from reco
1036  recoPID = aodtrack->GetMostProbablePID();
1037  AliAODPid* pidObj = aodtrack->GetDetPid();
1038  if(!pidObj)
1039  return;
1040 
1041  sigITS = pidObj->GetITSsignal();
1042  sigTPC = pidObj->GetTPCsignal();
1043  sigTOF = pidObj->GetTOFsignal();
1044  sigTRD = pidObj->GetTRDsignal();
1045  }
1046 
1047  // Get truth values if we are on MC
1049  {
1050  for(Int_t i=0; i<fTruthParticleArray->GetEntries();i++)
1051  {
1052  AliAODMCParticle* mcParticle = dynamic_cast<AliAODMCParticle*>(fTruthParticleArray->At(i));
1053  if(!mcParticle) continue;
1054 
1055  if (mcParticle->GetLabel() == particle->GetLabel())
1056  {
1057  // Use same convention as PID in AODs
1058  if(TMath::Abs(mcParticle->PdgCode()) == 2212) // proton
1059  truePID = 4;
1060  else if (TMath::Abs(mcParticle->PdgCode()) == 211) // pion
1061  truePID = 2;
1062  else if (TMath::Abs(mcParticle->PdgCode()) == 321) // kaon
1063  truePID = 3;
1064  else if (TMath::Abs(mcParticle->PdgCode()) == 11) // electron
1065  truePID = 0;
1066  else if (TMath::Abs(mcParticle->PdgCode()) == 13) // muon
1067  truePID = 1;
1068  else if (TMath::Abs(mcParticle->PdgCode()) == 700201) // deuteron
1069  truePID = 5;
1070  else if (TMath::Abs(mcParticle->PdgCode()) == 700301) // triton
1071  truePID = 6;
1072  else if (TMath::Abs(mcParticle->PdgCode()) == 700302) // He3
1073  truePID = 7;
1074  else if (TMath::Abs(mcParticle->PdgCode()) == 700202) // alpha
1075  truePID = 8;
1076  else
1077  truePID = 9;
1078 
1079  break;
1080  }
1081  }
1082  }
1083 }
1084 
1085 //________________________________________________________________________
1087 {
1088  // Print properties for extraction
1089  // to be implemented later
1090 }
1091 
1092 
1093 //########################################################################
1094 // HELPERS
1095 //########################################################################
1096 
1097 //________________________________________________________________________
1098 inline Bool_t AliAnalysisTaskJetExtractor::IsTrackInCone(AliVParticle* track, Double_t eta, Double_t phi, Double_t radius)
1099 {
1100  // This is to use a full cone in phi even at the edges of phi (2pi -> 0) (0 -> 2pi)
1101  Double_t trackPhi = 0.0;
1102  if (track->Phi() > (TMath::TwoPi() - (radius-phi)))
1103  trackPhi = track->Phi() - TMath::TwoPi();
1104  else if (track->Phi() < (phi+radius - TMath::TwoPi()))
1105  trackPhi = track->Phi() + TMath::TwoPi();
1106  else
1107  trackPhi = track->Phi();
1108 
1109  if ( TMath::Abs(trackPhi-phi)*TMath::Abs(trackPhi-phi) + TMath::Abs(track->Eta()-eta)*TMath::Abs(track->Eta()-eta) <= radius*radius)
1110  return kTRUE;
1111 
1112  return kFALSE;
1113 }
1114 
1115 //________________________________________________________________________
1117 {
1118  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1119  if(!tmpHist)
1120  {
1121  AliError(Form("Cannot find histogram <%s> ",key)) ;
1122  return;
1123  }
1124 
1125  tmpHist->Fill(x);
1126 }
1127 
1128 //________________________________________________________________________
1130 {
1131  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1132  if(!tmpHist)
1133  {
1134  AliError(Form("Cannot find histogram <%s> ",key));
1135  return;
1136  }
1137 
1138  if (tmpHist->IsA()->GetBaseClass("TH1"))
1139  static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
1140  else if (tmpHist->IsA()->GetBaseClass("TH2"))
1141  static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
1142 }
1143 
1144 //________________________________________________________________________
1146 {
1147  TH2* tmpHist = static_cast<TH2*>(fOutput->FindObject(key));
1148  if(!tmpHist)
1149  {
1150  AliError(Form("Cannot find histogram <%s> ",key));
1151  return;
1152  }
1153 
1154  tmpHist->Fill(x,y,add);
1155 }
1156 
1157 //________________________________________________________________________
1159 {
1160  TH3* tmpHist = static_cast<TH3*>(fOutput->FindObject(key));
1161  if(!tmpHist)
1162  {
1163  AliError(Form("Cannot find histogram <%s> ",key));
1164  return;
1165  }
1166 
1167  if(add)
1168  tmpHist->Fill(x,y,z,add);
1169  else
1170  tmpHist->Fill(x,y,z);
1171 }
1172 
1173 
1174 //________________________________________________________________________
1175 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)
1176 {
1177  T* tmpHist = new T(name, title, xBins, xMin, xMax);
1178 
1179  tmpHist->GetXaxis()->SetTitle(xTitle);
1180  tmpHist->GetYaxis()->SetTitle(yTitle);
1181  tmpHist->SetOption(options);
1182  tmpHist->SetMarkerStyle(kFullCircle);
1183  tmpHist->Sumw2();
1184 
1185  fOutput->Add(tmpHist);
1186 
1187  return tmpHist;
1188 }
1189 
1190 //________________________________________________________________________
1191 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)
1192 {
1193  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax);
1194  tmpHist->GetXaxis()->SetTitle(xTitle);
1195  tmpHist->GetYaxis()->SetTitle(yTitle);
1196  tmpHist->GetZaxis()->SetTitle(zTitle);
1197  tmpHist->SetOption(options);
1198  tmpHist->SetMarkerStyle(kFullCircle);
1199  tmpHist->Sumw2();
1200 
1201  fOutput->Add(tmpHist);
1202 
1203  return tmpHist;
1204 }
1205 
1206 //________________________________________________________________________
1207 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)
1208 {
1209  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax, zBins, zMin, zMax);
1210  tmpHist->GetXaxis()->SetTitle(xTitle);
1211  tmpHist->GetYaxis()->SetTitle(yTitle);
1212  tmpHist->GetZaxis()->SetTitle(zTitle);
1213  tmpHist->SetOption(options);
1214  tmpHist->SetMarkerStyle(kFullCircle);
1215  tmpHist->Sumw2();
1216 
1217  fOutput->Add(tmpHist);
1218 
1219  return tmpHist;
1220 }
1221 
1222 //________________________________________________________________________
1224 {
1225  // Called once at the end of the analysis.
1226 }
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
virtual AliVParticle * GetNextAcceptParticle()
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
centrality
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
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
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)
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) ...
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
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)
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
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
Int_t fTruthMaxLabel
max track label to consider it as true particle
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)
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
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)
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, AliParticleContainer *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 fBuffer_Event_MagneticField
! array buffer