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