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