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