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