AliPhysics  2c6b7ad (2c6b7ad)
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 <iostream>
17 #include <cstdlib>
18 #include <algorithm>
19 #include <vector>
20 #include <TClonesArray.h>
21 #include <TH1F.h>
22 #include <TH2F.h>
23 #include <TH3F.h>
24 #include <TProfile.h>
25 #include <TTree.h>
26 #include <TList.h>
27 #include <TLorentzVector.h>
28 #include <TNamed.h>
29 
30 #include "AliMCEvent.h"
31 #include "AliESDVertex.h"
32 #include "AliAODVertex.h"
33 #include "AliAODPid.h"
34 
35 #include "AliRhoParameter.h"
36 
37 #include "AliVTrack.h"
38 #include "AliVHeader.h"
39 #include "AliEmcalJet.h"
40 #include "AliLog.h"
41 #include "AliJetContainer.h"
42 #include "AliParticleContainer.h"
43 #include "AliAODTrack.h"
44 #include "AliVParticle.h"
45 #include "TRandom3.h"
46 #include "AliEmcalPythiaInfo.h"
48 #include "AliAnalysisManager.h"
49 #include "AliYAMLConfiguration.h"
50 
51 #include "AliGenHepMCEventHeader.h"
52 #include "AliGenHijingEventHeader.h"
53 #include "AliHFJetsTaggingVertex.h"
54 #include "AliRDHFJetsCutsVertex.h"
55 
57 
59 ClassImp(AliEmcalJetTree)
61 
65 
66 //________________________________________________________________________
67 AliEmcalJetTree::AliEmcalJetTree() : TNamed("CustomTree", "CustomTree"), fJetTree(0), fInitialized(0), fSaveEventProperties(0), fSaveConstituents(0), fSaveConstituentsIP(0), fSaveConstituentPID(0), fSaveJetShapes(0), fSaveMCInformation(0), fSaveSecondaryVertices(0), fSaveTriggerTracks(0), fExtractionPercentages(), fExtractionPercentagePtBins(), fExtractionJetTypes_HM(), fExtractionJetTypes_PM()
68 {
69  // For these arrays, we need to reserve memory
70  fBuffer_Const_Pt = new Float_t[kMaxNumConstituents];
71  fBuffer_Const_Eta = new Float_t[kMaxNumConstituents];
72  fBuffer_Const_Phi = new Float_t[kMaxNumConstituents];
73  fBuffer_Const_Charge = new Float_t[kMaxNumConstituents];
74  fBuffer_Const_Label = new Int_t [kMaxNumConstituents];
75  fBuffer_Const_ProdVtx_X = new Float_t[kMaxNumConstituents];
76  fBuffer_Const_ProdVtx_Y = new Float_t[kMaxNumConstituents];
77  fBuffer_Const_ProdVtx_Z = new Float_t[kMaxNumConstituents];
78 
79 }
80 
81 //________________________________________________________________________
82 AliEmcalJetTree::AliEmcalJetTree(const char* name) : TNamed(name, name), fJetTree(0), fInitialized(0), fSaveEventProperties(0), fSaveConstituents(0), fSaveConstituentsIP(0), fSaveConstituentPID(0), fSaveJetShapes(0), fSaveMCInformation(0), fSaveSecondaryVertices(0), fSaveTriggerTracks(0), fExtractionPercentages(), fExtractionPercentagePtBins(), fExtractionJetTypes_HM(), fExtractionJetTypes_PM()
83 {
84  // For these arrays, we need to reserve memory
93 }
94 
95 
96 //________________________________________________________________________
98  Int_t motherParton, Int_t motherHadron, Int_t partonInitialCollision, Float_t matchDistance, Float_t matchedPt, Float_t matchedMass, Float_t truePtFraction, Float_t ptHard, Float_t eventWeight, Float_t impactParameter,
99  Float_t* trackPID_ITS, Float_t* trackPID_TPC, Float_t* trackPID_TOF, Float_t* trackPID_TRD, Short_t* trackPID_Reco, Int_t* trackPID_Truth,
100  Int_t numTriggerTracks, Float_t* triggerTrackPt, Float_t* triggerTrackDeltaEta, Float_t* triggerTrackDeltaPhi,
101  Float_t* trackIP_d0, Float_t* trackIP_z0, Float_t* trackIP_d0cov, Float_t* trackIP_z0cov,
102  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)
103 {
104  // ############
105  // Approach: Fill buffer and then add to tree
106  //
107  if(!fInitialized)
108  AliFatal("Tree is not initialized.");
109 
110  fBuffer_JetPt = jet->Pt() - fJetContainer->GetRhoVal()*jet->Area();
111 
112  // Check if jet type is contained in extraction list
113  if( (fExtractionJetTypes_PM.size() || fExtractionJetTypes_HM.size()) &&
114  (std::find(fExtractionJetTypes_PM.begin(), fExtractionJetTypes_PM.end(), motherParton) == fExtractionJetTypes_PM.end()) &&
115  (std::find(fExtractionJetTypes_HM.begin(), fExtractionJetTypes_HM.end(), motherHadron) == fExtractionJetTypes_HM.end()) )
116  return kFALSE;
117 
118  // Check acceptance percentage for the given jet and discard statistically on demand
119  Bool_t inPtRange = kFALSE;
120  for(size_t i = 0; i<fExtractionPercentages.size(); i++)
121  {
123  {
124  inPtRange = kTRUE;
125  if(fRandomGenerator->Rndm() >= fExtractionPercentages[i])
126  return kFALSE;
127  break;
128  }
129  }
130  if(!inPtRange && fExtractionPercentagePtBins.size()) // only discard jet if fExtractionPercentagePtBins was given
131  return kFALSE;
132 
133  // Set basic properties
134  fBuffer_JetEta = jet->Eta();
135  fBuffer_JetPhi = jet->Phi();
136  fBuffer_JetArea = jet->Area();
137 
138  // Set event properties
140  {
143  fBuffer_Event_Vertex_X = vertexX;
144  fBuffer_Event_Vertex_Y = vertexY;
145  fBuffer_Event_Vertex_Z = vertexZ;
147  fBuffer_Event_ID = eventID;
148  fBuffer_Event_MagneticField = magField;
150  {
151  fBuffer_Event_PtHard = ptHard;
152  fBuffer_Event_Weight = eventWeight;
153  fBuffer_Event_ImpactParameter = impactParameter;
154  }
155  }
156 
157  // Extract basic constituent properties directly from AliEmcalJet object
160  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
161  {
162  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(i, fTrackContainer->GetArray()));
163  if(!particle) continue;
164 
166  {
167  fBuffer_Const_Pt[fBuffer_NumConstituents] = particle->Pt();
168  fBuffer_Const_Eta[fBuffer_NumConstituents] = particle->Eta();
169  fBuffer_Const_Phi[fBuffer_NumConstituents] = particle->Phi();
170  fBuffer_Const_Charge[fBuffer_NumConstituents] = particle->Charge();
171  fBuffer_Const_Label[fBuffer_NumConstituents] = particle->GetLabel();
172  }
174  {
178  }
180  }
181 
182 
183  // Set constituent arrays for impact parameters
185  {
186  fJetTree->SetBranchAddress("Jet_Const_CovIPd", trackIP_d0cov);
187  fJetTree->SetBranchAddress("Jet_Const_CovIPz", trackIP_z0cov);
188  fJetTree->SetBranchAddress("Jet_Const_IPd", trackIP_d0);
189  fJetTree->SetBranchAddress("Jet_Const_IPz", trackIP_z0);
190  }
191  // Set constituent arrays for PID parameters
193  {
194  fJetTree->SetBranchAddress("Jet_Const_PID_ITS", trackPID_ITS);
195  fJetTree->SetBranchAddress("Jet_Const_PID_TPC", trackPID_TPC);
196  fJetTree->SetBranchAddress("Jet_Const_PID_TOF", trackPID_TOF);
197  fJetTree->SetBranchAddress("Jet_Const_PID_TRD", trackPID_TRD);
198  fJetTree->SetBranchAddress("Jet_Const_PID_Reconstructed", trackPID_Reco);
200  fJetTree->SetBranchAddress("Jet_Const_PID_Truth", trackPID_Truth);
201  }
202 
203  // Set secondary vertices arrays
205  {
206  fBuffer_NumSecVertices = numSecVertices;
207  fJetTree->SetBranchAddress("Jet_SecVtx_X", secVtx_X);
208  fJetTree->SetBranchAddress("Jet_SecVtx_Y", secVtx_Y);
209  fJetTree->SetBranchAddress("Jet_SecVtx_Z", secVtx_Z);
210  fJetTree->SetBranchAddress("Jet_SecVtx_Mass", secVtx_Mass);
211  fJetTree->SetBranchAddress("Jet_SecVtx_Lxy", secVtx_Lxy);
212  fJetTree->SetBranchAddress("Jet_SecVtx_SigmaLxy", secVtx_SigmaLxy);
213  fJetTree->SetBranchAddress("Jet_SecVtx_Chi2", secVtx_Chi2);
214  fJetTree->SetBranchAddress("Jet_SecVtx_Dispersion", secVtx_Dispersion);
215  }
216 
217  // Set jet shape observables
218  fBuffer_Shape_Mass = jet->M() - fJetContainer->GetRhoMassVal()*jet->Area();
219 
220  // Set Monte Carlo information
222  {
223  fBuffer_Jet_MC_MotherParton = motherParton;
224  fBuffer_Jet_MC_MotherHadron = motherHadron;
225  fBuffer_Jet_MC_MotherIC = partonInitialCollision;
226  fBuffer_Jet_MC_MatchDistance = matchDistance;
227  fBuffer_Jet_MC_MatchedPt = matchedPt;
228  fBuffer_Jet_MC_MatchedMass = matchedMass;
229  fBuffer_Jet_MC_TruePtFraction = truePtFraction;
230  }
231 
233  {
234  fBuffer_NumTriggerTracks = numTriggerTracks;
235  fJetTree->SetBranchAddress("Jet_TriggerTrack_Pt", triggerTrackPt);
236  fJetTree->SetBranchAddress("Jet_TriggerTrack_dEta", triggerTrackDeltaEta);
237  fJetTree->SetBranchAddress("Jet_TriggerTrack_dPhi", triggerTrackDeltaPhi);
238  }
239 
240  // Add buffered jet to tree
241  fJetTree->Fill();
242 
243  return kTRUE;
244 }
245 
246 
247 //________________________________________________________________________
249 {
250  // Create the tree with active branches
251 
252  void* dummy = 0; // for some branches, do not need a buffer pointer yet
253 
254 
255  if(fInitialized)
256  AliFatal("Tree is already initialized.");
257 
258  if(!fJetContainer)
259  AliFatal("fJetContainer not set in tree");
260 
261  // ### Prepare the jet tree
262  fJetTree = new TTree(Form("JetTree_%s", GetName()), "");
263 
264  fJetTree->Branch("Jet_Pt",&fBuffer_JetPt,"Jet_Pt/F");
265  fJetTree->Branch("Jet_Phi",&fBuffer_JetPhi,"Jet_Phi/F");
266  fJetTree->Branch("Jet_Eta",&fBuffer_JetEta,"Jet_Eta/F");
267  fJetTree->Branch("Jet_Area",&fBuffer_JetArea,"Jet_Area/F");
268  fJetTree->Branch("Jet_NumConstituents",&fBuffer_NumConstituents,"Jet_NumConstituents/I");
269 
271  {
272  fJetTree->Branch("Event_BackgroundDensity",&fBuffer_Event_BackgroundDensity,"Event_BackgroundDensity/F");
273  fJetTree->Branch("Event_BackgroundDensityMass",&fBuffer_Event_BackgroundDensityMass,"Event_BackgroundDensityMass/F");
274  fJetTree->Branch("Event_Vertex_X",&fBuffer_Event_Vertex_X,"Event_Vertex_X/F");
275  fJetTree->Branch("Event_Vertex_Y",&fBuffer_Event_Vertex_Y,"Event_Vertex_Y/F");
276  fJetTree->Branch("Event_Vertex_Z",&fBuffer_Event_Vertex_Z,"Event_Vertex_Z/F");
277  fJetTree->Branch("Event_Centrality",&fBuffer_Event_Centrality,"Event_Centrality/F");
278  fJetTree->Branch("Event_ID",&fBuffer_Event_ID,"Event_ID/L");
279  fJetTree->Branch("Event_MagneticField",&fBuffer_Event_MagneticField,"Event_MagneticField/F");
280 
282  {
283  fJetTree->Branch("Event_PtHard",&fBuffer_Event_PtHard,"Event_PtHard/F");
284  fJetTree->Branch("Event_Weight",&fBuffer_Event_Weight,"Event_Weight/F");
285  fJetTree->Branch("Event_ImpactParameter",&fBuffer_Event_ImpactParameter,"Event_ImpactParameter/F");
286  }
287  }
288 
290  {
291  fJetTree->Branch("Jet_Const_Pt",fBuffer_Const_Pt,"Jet_Const_Pt[Jet_NumConstituents]/F");
292  fJetTree->Branch("Jet_Const_Phi",fBuffer_Const_Phi,"Jet_Const_Phi[Jet_NumConstituents]/F");
293  fJetTree->Branch("Jet_Const_Eta",fBuffer_Const_Eta,"Jet_Const_Eta[Jet_NumConstituents]/F");
294  fJetTree->Branch("Jet_Const_Charge",fBuffer_Const_Charge,"Jet_Const_Charge[Jet_NumConstituents]/F");
296  fJetTree->Branch("Jet_Const_Label",fBuffer_Const_Label,"Jet_Const_Label[Jet_NumConstituents]/I");
297  }
298 
300  {
301  fJetTree->Branch("Jet_Const_IPd",&dummy,"Jet_Const_IPd[Jet_NumConstituents]/F");
302  fJetTree->Branch("Jet_Const_IPz",&dummy,"Jet_Const_IPz[Jet_NumConstituents]/F");
303  fJetTree->Branch("Jet_Const_CovIPd",&dummy,"Jet_Const_CovIPd[Jet_NumConstituents]/F");
304  fJetTree->Branch("Jet_Const_CovIPz",&dummy,"Jet_Const_CovIPz[Jet_NumConstituents]/F");
305 
306  fJetTree->Branch("Jet_Const_ProdVtx_X",fBuffer_Const_ProdVtx_X,"Jet_Const_ProdVtx_X[Jet_NumConstituents]/F");
307  fJetTree->Branch("Jet_Const_ProdVtx_Y",fBuffer_Const_ProdVtx_Y,"Jet_Const_ProdVtx_Y[Jet_NumConstituents]/F");
308  fJetTree->Branch("Jet_Const_ProdVtx_Z",fBuffer_Const_ProdVtx_Z,"Jet_Const_ProdVtx_Z[Jet_NumConstituents]/F");
309  }
310 
312  {
313  fJetTree->Branch("Jet_Const_PID_ITS",&dummy,"Jet_Const_PID_ITS[Jet_NumConstituents]/F");
314  fJetTree->Branch("Jet_Const_PID_TPC",&dummy,"Jet_Const_PID_TPC[Jet_NumConstituents]/F");
315  fJetTree->Branch("Jet_Const_PID_TOF",&dummy,"Jet_Const_PID_TOF[Jet_NumConstituents]/F");
316  fJetTree->Branch("Jet_Const_PID_TRD",&dummy,"Jet_Const_PID_TRD[Jet_NumConstituents]/F");
317 
318  fJetTree->Branch("Jet_Const_PID_Reconstructed",&dummy,"Jet_Const_PID_Reconstructed[Jet_NumConstituents]/S");
320  fJetTree->Branch("Jet_Const_PID_Truth",&dummy,"Jet_Const_PID_Truth[Jet_NumConstituents]/I");
321  }
322 
323  if(fSaveJetShapes)
324  {
325  fJetTree->Branch("Jet_Shape_Mass",&fBuffer_Shape_Mass,"Jet_Shape_Mass/F");
326  }
327 
329  {
330  fJetTree->Branch("Jet_MC_MotherParton",&fBuffer_Jet_MC_MotherParton,"Jet_MC_MotherParton/I");
331  fJetTree->Branch("Jet_MC_MotherHadron",&fBuffer_Jet_MC_MotherHadron,"Jet_MC_MotherHadron/I");
332  fJetTree->Branch("Jet_MC_MotherIC",&fBuffer_Jet_MC_MotherIC,"Jet_MC_MotherIC/I");
333  fJetTree->Branch("Jet_MC_MatchDistance",&fBuffer_Jet_MC_MatchDistance,"Jet_MC_MatchDistance/F");
334  fJetTree->Branch("Jet_MC_MatchedPt",&fBuffer_Jet_MC_MatchedPt,"Jet_MC_MatchedPt/F");
335  fJetTree->Branch("Jet_MC_MatchedMass",&fBuffer_Jet_MC_MatchedMass,"Jet_MC_MatchedMass/F");
336  fJetTree->Branch("Jet_MC_TruePtFraction",&fBuffer_Jet_MC_TruePtFraction,"Jet_MC_TruePtFraction/F");
337  }
338 
340  {
341  fJetTree->Branch("Jet_NumSecVertices",&fBuffer_NumSecVertices,"Jet_NumSecVertices/I");
342 
343  fJetTree->Branch("Jet_SecVtx_X",&dummy,"Jet_SecVtx_X[Jet_NumSecVertices]/F");
344  fJetTree->Branch("Jet_SecVtx_Y",&dummy,"Jet_SecVtx_Y[Jet_NumSecVertices]/F");
345  fJetTree->Branch("Jet_SecVtx_Z",&dummy,"Jet_SecVtx_Z[Jet_NumSecVertices]/F");
346  fJetTree->Branch("Jet_SecVtx_Mass",&dummy,"Jet_SecVtx_Mass[Jet_NumSecVertices]/F");
347  fJetTree->Branch("Jet_SecVtx_Lxy",&dummy,"Jet_SecVtx_Lxy[Jet_NumSecVertices]/F");
348  fJetTree->Branch("Jet_SecVtx_SigmaLxy",&dummy,"Jet_SecVtx_SigmaLxy[Jet_NumSecVertices]/F");
349  fJetTree->Branch("Jet_SecVtx_Chi2",&dummy,"Jet_SecVtx_Chi2[Jet_NumSecVertices]/F");
350  fJetTree->Branch("Jet_SecVtx_Dispersion",&dummy,"Jet_SecVtx_Dispersion[Jet_NumSecVertices]/F");
351  }
352 
353  // Trigger track requirement active -> Save trigger track
355  {
356  fJetTree->Branch("Jet_NumTriggerTracks",&fBuffer_NumTriggerTracks,"Jet_NumTriggerTracks/I");
357  fJetTree->Branch("Jet_TriggerTrack_Pt",&dummy,"Jet_TriggerTrack_Pt[Jet_NumTriggerTracks]/F");
358  fJetTree->Branch("Jet_TriggerTrack_dEta",&dummy,"Jet_TriggerTrack_dEta[Jet_NumTriggerTracks]/F");
359  fJetTree->Branch("Jet_TriggerTrack_dPhi",&dummy,"Jet_TriggerTrack_dPhi[Jet_NumTriggerTracks]/F");
360  }
361 
362  fInitialized = kTRUE;
363 }
364 
365 //________________________________________________________________________
367  AliAnalysisTaskEmcalJet("AliAnalysisTaskJetExtractor", kTRUE),
368  fJetTree(0),
369  fVtxTagger(0),
370  fHadronMatchingRadius(0.4),
371  fTrueJetMatchingRadius(0.4),
372  fSecondaryVertexMaxChi2(1e10),
373  fSecondaryVertexMaxDispersion(0.05),
374  fCalculateSecondaryVertices(kTRUE),
375  fVertexerCuts(0),
376  fSetEmcalJetFlavour(0),
377  fEventPercentage(1.0),
378  fTruthMinLabel(0),
379  fTruthMaxLabel(100000),
380  fSaveTrackPDGCode(kTRUE),
381  fRandomSeed(0),
382  fRandomSeedCones(0),
383  fEventCut_TriggerTrackMinPt(0),
384  fEventCut_TriggerTrackMaxPt(0),
385  fEventCut_TriggerTrackMinLabel(-9999999),
386  fEventCut_TriggerTrackMaxLabel(+9999999),
387  fJetsCont(0),
388  fTracksCont(0),
389  fTruthParticleArray(0),
390  fTruthJetsArrayName(""),
391  fTruthJetsRhoName(""),
392  fTruthJetsRhoMassName(""),
393  fTruthParticleArrayName("mcparticles"),
394  fRandomGenerator(0),
395  fRandomGeneratorCones(0),
396  fEventWeight(0.),
397  fImpactParameter(0.),
398  fTriggerTracks_Pt(),
399  fTriggerTracks_Eta(),
400  fTriggerTracks_Phi()
401 {
402  fRandomGenerator = new TRandom3();
403  fRandomGeneratorCones = new TRandom3();
404 
406  fJetTree = new AliEmcalJetTree(GetName());
409  fJetTree->SetSaveJetShapes(kTRUE);
410  DefineOutput(2, TTree::Class());
411 }
412 
413 //________________________________________________________________________
415  AliAnalysisTaskEmcalJet(name, kTRUE),
416  fJetTree(0),
417  fVtxTagger(0),
423  fVertexerCuts(0),
425  fEventPercentage(1.0),
426  fTruthMinLabel(0),
427  fTruthMaxLabel(100000),
428  fSaveTrackPDGCode(kTRUE),
429  fRandomSeed(0),
430  fRandomSeedCones(0),
435  fJetsCont(0),
436  fTracksCont(0),
439  fTruthJetsRhoName(""),
441  fTruthParticleArrayName("mcparticles"),
442  fRandomGenerator(0),
444  fEventWeight(0.),
445  fImpactParameter(0.),
449 {
450  fRandomGenerator = new TRandom3();
451  fRandomGeneratorCones = new TRandom3();
452 
454  fJetTree = new AliEmcalJetTree(GetName());
457  fJetTree->SetSaveJetShapes(kTRUE);
458  DefineOutput(2, TTree::Class());
459 }
460 
461 //________________________________________________________________________
463 {
464  // Destructor.
465 }
466 
467 //________________________________________________________________________
469 {
471 
472  // ### Basic container settings
474  if(!fJetsCont)
475  AliFatal("Jet input container not found!");
476  fJetsCont->PrintCuts();
478  if(!fTracksCont)
479  AliFatal("Particle input container not found attached to jets!");
480 
481  fRandomGenerator->SetSeed(fRandomSeed);
483 
484  // Activate saving of trigger tracks if this is demanded
487 
488  // ### Initialize the jet tree (settings must all be given at this stage)
493  OpenFile(2);
494  PostData(2, fJetTree->GetTreePointer());
495 
496  // ### Add control histograms (already some created in base task)
497  AddHistogram2D<TH2D>("hTrackCount", "Number of tracks in acceptance vs. centrality", "COLZ", 500, 0., 5000., 100, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}");
498  AddHistogram2D<TH2D>("hBackgroundPt", "Background p_{T} distribution", "", 150, 0., 150., 100, 0, 100, "Background p_{T} (GeV/c)", "Centrality", "dN^{Events}/dp_{T}");
499 
500  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}");
501  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}");
502  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}");
503  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");
504  AddHistogram2D<TH2D>("hJetArea", "Jet area", "COLZ", 200, 0., 2., 100, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
505 
506  AddHistogram2D<TH2D>("hDeltaPt", "#delta p_{T} distribution", "", 400, -100., 300., 100, 0, 100, "p_{T, cone} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
507 
508  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}");
509  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");
510 
511  AddHistogram1D<TH1D>("hExtractionPercentage", "Extracted jets p_{T} distribution (background subtracted)", "e1p", 400, -100., 300., "p_{T, jet} (GeV/c)", "Extracted percentage");
512 
513  // Track QA plots
514  AddHistogram2D<TH2D>("hTrackPt", "Tracks p_{T} distribution", "", 300, 0., 300., 100, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
515  AddHistogram2D<TH2D>("hTrackPhi", "Track angular distribution in #phi", "LEGO2", 180, 0., 2*TMath::Pi(), 100, 0, 100, "#phi", "Centrality", "dN^{Tracks}/(d#phi)");
516  AddHistogram2D<TH2D>("hTrackEta", "Track angular distribution in #eta", "LEGO2", 100, -2.5, 2.5, 100, 0, 100, "#eta", "Centrality", "dN^{Tracks}/(d#eta)");
517  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");
518 
519  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})");
520  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})");
521 }
522 
523 
524 //________________________________________________________________________
526 {
528 
529  // ### Need to explicitly tell jet container to load rho mass object
530  fJetsCont->LoadRhoMass(InputEvent());
531 
532  if (fTruthParticleArrayName != "")
533  fTruthParticleArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTruthParticleArrayName.Data()));
534 
535  // ### Prepare vertexer
537  {
538  if(!fVertexerCuts)
539  AliFatal("VertexerCuts not given but secondary vertex calculation turned on.");
540  fVtxTagger = new AliHFJetsTaggingVertex();
541  fVtxTagger->SetCuts(fVertexerCuts);
542  }
543 
544  // ### Save tree extraction percentage to histogram
545  std::vector<Float_t> extractionPtBins = fJetTree->GetExtractionPercentagePtBins();
546  std::vector<Float_t> extractionPercentages = fJetTree->GetExtractionPercentages();
547 
548  for(Int_t i=0; i<static_cast<Int_t>(extractionPercentages.size()); i++)
549  {
550  Double_t percentage = extractionPercentages[i];
551  for(Int_t pt=static_cast<Int_t>(extractionPtBins[i*2]); pt<static_cast<Int_t>(extractionPtBins[i*2+1]); pt++)
552  {
553  FillHistogram("hExtractionPercentage", pt, percentage);
554  }
555  }
556 
557  // ### Add HIJING/JEWEL histograms (in case we need it)
558  if(MCEvent())
559  {
560  if(dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()))
561  AddHistogram1D<TH1D>("hHIJING_Ncoll", "N_{coll} from HIJING", "e1p", 1000, 0., 5000, "N_{coll}", "dN^{Events}/dN^{N_{coll}}");
562 
563  if(dynamic_cast<AliGenHepMCEventHeader*>(MCEvent()->GenEventHeader()))
564  {
565  AddHistogram1D<TH1D>("hJEWEL_NProduced", "NProduced from JEWEL", "e1p", 1000, 0., 1000, "Nproduced", "dN^{Events}/dN^{Nproduced}");
566  AddHistogram1D<TH1D>("hJEWEL_ImpactParameter", "Impact parameter from JEWEL", "e1p", 1000, 0., 1, "IP", "dN^{Events}/dN^{IP}");
567  TProfile* evweight = new TProfile("hJEWEL_EventWeight", "Event weight from JEWEL", 1, 0,1);
568  fOutput->Add(evweight);
569  }
570  }
571 
572  PrintConfig();
573 
574 }
575 
576 //________________________________________________________________________
578 {
579  // ################################### EVENT SELECTION
580 
581  if(!IsTriggerTrackInEvent())
582  return kFALSE;
583  if(fRandomGenerator->Rndm() >= fEventPercentage)
584  return kFALSE;
585 
586  // ################################### EVENT PROPERTIES
588 
589  // Load vertex if possible
590  Double_t vtxX = 0;
591  Double_t vtxY = 0;
592  Double_t vtxZ = 0;
593  Long64_t eventID = 0;
594  const AliVVertex* myVertex = InputEvent()->GetPrimaryVertex();
595  if(!myVertex && MCEvent())
596  myVertex = MCEvent()->GetPrimaryVertex();
597  if(myVertex)
598  {
599  vtxX = myVertex->GetX();
600  vtxY = myVertex->GetY();
601  vtxZ = myVertex->GetZ();
602  }
603  // Get event ID from header
604  AliVHeader* eventIDHeader = InputEvent()->GetHeader();
605  if(eventIDHeader)
606  eventID = eventIDHeader->GetEventIdAsLong();
607 
608  // If available, get Ncoll from HIJING
609  if(MCEvent())
610  {
611  if(AliGenHijingEventHeader* mcHeader = dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader()))
612  {
613  Double_t ncoll = mcHeader->NN() + mcHeader->NNw() + mcHeader->NwN() + mcHeader->NwNw();
614  FillHistogram("hHIJING_Ncoll", ncoll);
615  }
616  if(AliGenHepMCEventHeader* mcHeader = dynamic_cast<AliGenHepMCEventHeader*>(MCEvent()->GenEventHeader()))
617  {
618  fEventWeight = mcHeader->EventWeight();
619  fImpactParameter = mcHeader->impact_parameter();
620 
621  TProfile* tmpHist = static_cast<TProfile*>(fOutput->FindObject("hJEWEL_EventWeight"));
622  tmpHist->Fill(0., fEventWeight);
623  FillHistogram("hJEWEL_NProduced", mcHeader->NProduced());
624  FillHistogram("hJEWEL_ImpactParameter", fImpactParameter);
625  }
626  }
627 
628  // ################################### MAIN JET LOOP
629  fJetsCont->ResetCurrentID();
630  while(AliEmcalJet *jet = fJetsCont->GetNextAcceptJet())
631  {
633 
634  Double_t matchDistance = 0;
635  Double_t matchedJetPt = 0;
636  Double_t matchedJetMass = 0;
637  Double_t truePtFraction = 0;
638  Int_t currentJetType_HM = 0;
639  Int_t currentJetType_PM = 0;
640  Int_t currentJetType_IC = 0;
642  {
643  // Get jet type from MC (hadron matching, parton matching definition - for HF jets)
644  GetJetType(jet, currentJetType_HM, currentJetType_PM, currentJetType_IC);
645  // Get true estimators: for pt, jet mass, ...
646  truePtFraction = GetTrueJetPtFraction(jet);
647  matchDistance = GetMatchedTrueJetObservables(jet, matchedJetPt, matchedJetMass);
648  }
649 
650  // ### CONSTITUENT LOOP: Retrieve PID values + impact parameters
651  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<Int_t> vecTruePID;
652  std::vector<Float_t> vec_d0; std::vector<Float_t> vec_d0cov; std::vector<Float_t> vec_z0; std::vector<Float_t> vec_z0cov;
653 
655  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
656  {
657  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
658  if(!particle) continue;
660  {
661  Float_t sigITS = 0; Float_t sigTPC = 0; Float_t sigTOF = 0; Float_t sigTRD = 0; Short_t recoPID = 0; Int_t truePID = 0;
662  AddPIDInformation(particle, sigITS, sigTPC, sigTOF, sigTRD, recoPID, truePID);
663  vecSigITS.push_back(sigITS); vecSigTPC.push_back(sigTPC); vecSigTOF.push_back(sigTOF); vecSigTRD.push_back(sigTRD); vecRecoPID.push_back(recoPID); vecTruePID.push_back(truePID);
664  }
666  {
667  Float_t d0 = 0; Float_t d0cov = 0; Float_t z0 = 0; Float_t z0cov = 0;
668  GetTrackImpactParameters(myVertex, dynamic_cast<AliAODTrack*>(particle), d0, d0cov, z0, z0cov);
669  vec_d0.push_back(d0); vec_d0cov.push_back(d0cov); vec_z0.push_back(z0); vec_z0cov.push_back(z0cov);
670  }
671  }
672 
673  // Reconstruct secondary vertices
674  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;
676  ReconstructSecondaryVertices(myVertex, jet, secVtx_X, secVtx_Y, secVtx_Z, secVtx_Mass, secVtx_Lxy, secVtx_SigmaLxy, secVtx_Chi2, secVtx_Dispersion);
677 
678  // Now change trigger tracks eta/phi to dEta/dPhi relative to jet
679  std::vector<Float_t> triggerTracks_dEta(fTriggerTracks_Eta);
680  std::vector<Float_t> triggerTracks_dPhi(fTriggerTracks_Phi);
681  for(UInt_t i=0; i<triggerTracks_dEta.size(); i++)
682  {
683  triggerTracks_dEta[i] = jet->Eta() - fTriggerTracks_Eta[i];
684  triggerTracks_dPhi[i] = TMath::Min(TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]), TMath::TwoPi() - TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]));
685  if( ((TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]) <= TMath::Pi()) && (jet->Phi() - fTriggerTracks_Phi[i] > 0))
686  || ((TMath::Abs(jet->Phi() - fTriggerTracks_Phi[i]) > TMath::Pi()) && (jet->Phi() - fTriggerTracks_Phi[i] < 0)) )
687  triggerTracks_dPhi[i] = -triggerTracks_dPhi[i];
688  }
689 
690  // Fill jet to tree
691  Bool_t accepted = fJetTree->AddJetToTree(jet, vtxX, vtxY, vtxZ, fCent, eventID, InputEvent()->GetMagneticField(),
692  currentJetType_PM,currentJetType_HM,currentJetType_IC,matchDistance,matchedJetPt,matchedJetMass,truePtFraction,fPtHard,fEventWeight,fImpactParameter,
693  vecSigITS, vecSigTPC, vecSigTOF, vecSigTRD, vecRecoPID, vecTruePID,
694  fTriggerTracks_Pt, triggerTracks_dEta, triggerTracks_dPhi,
695  vec_d0, vec_z0, vec_d0cov, vec_z0cov,
696  secVtx_X, secVtx_Y, secVtx_Z, secVtx_Mass, secVtx_Lxy, secVtx_SigmaLxy, secVtx_Chi2, secVtx_Dispersion);
697 
698  if(accepted)
699  FillHistogram("hJetPtExtracted", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
700  }
701 
702  // ################################### PARTICLE PROPERTIES
703  fTracksCont->ResetCurrentID();
704  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
706 
707  return kTRUE;
708 }
709 
710 //________________________________________________________________________
712 {
713  // Cut for trigger track requirement
715  {
716  // Clear vector of trigger tracks
717  fTriggerTracks_Pt.clear();
718  fTriggerTracks_Eta.clear();
719  fTriggerTracks_Phi.clear();
720 
721  // Go through all tracks and check whether trigger tracks can be found
722  fTracksCont->ResetCurrentID();
723  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
724  {
725  if( (track->GetLabel() >= fEventCut_TriggerTrackMinLabel) && (track->GetLabel() < fEventCut_TriggerTrackMaxLabel) )
726  if( (track->Pt() >= fEventCut_TriggerTrackMinPt) && (track->Pt() < fEventCut_TriggerTrackMaxPt) )
727  {
728  fTriggerTracks_Pt.push_back(track->Pt());
729  fTriggerTracks_Eta.push_back(track->Eta());
730  fTriggerTracks_Phi.push_back(track->Phi());
731  }
732  }
733  // No particle has been found that fulfills requirement -> Do not accept event
734  if(fTriggerTracks_Pt.size() == 0)
735  return kFALSE;
736  }
737  return kTRUE;
738 }
739 
740 //________________________________________________________________________
742 {
743  // #################################################################################
744  // ##### FRACTION OF TRUE PT IN JET: Defined as "not from toy"
745  Double_t pt_nonMC = 0.;
746  Double_t pt_all = 0.;
747  Double_t truePtFraction = 0;
748 
749  for(Int_t iConst = 0; iConst < jet->GetNumberOfTracks(); iConst++)
750  {
751  // Loop over all valid jet constituents
752  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(iConst, fTracksCont->GetArray()));
753  if(!particle) continue;
754  if(particle->Pt() < 1e-6) continue;
755 
756  // Particles marked w/ labels within label range are considered from toy
757  if( (particle->GetLabel() >= fTruthMinLabel) && (particle->GetLabel() < fTruthMaxLabel))
758  pt_nonMC += particle->Pt();
759  pt_all += particle->Pt();
760  }
761  if(pt_all)
762  truePtFraction = (pt_nonMC/pt_all);
763  return truePtFraction;
764 }
765 
766 //________________________________________________________________________
768 {
769  // #################################################################################
770  // ##### OBSERVABLES FROM MATCHED JETS: Jet pt, jet mass
771  Double_t bestMatchDeltaR = 8.; // 8 is higher than maximum possible matching distance
772  if(fTruthJetsArrayName != "")
773  {
774  // "True" background for pt
775  AliRhoParameter* rho = static_cast<AliRhoParameter*>(InputEvent()->FindListObject(fTruthJetsRhoName.Data()));
776  Double_t trueRho = 0;
777  if(rho)
778  trueRho = rho->GetVal();
779 
780  // "True" background for mass
781  AliRhoParameter* rhoMass = static_cast<AliRhoParameter*>(InputEvent()->FindListObject(fTruthJetsRhoMassName.Data()));
782  Double_t trueRhoMass = 0;
783  if(rhoMass)
784  trueRhoMass = rhoMass->GetVal();
785 
786  TClonesArray* truthArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fTruthJetsArrayName.Data())));
787 
788  // Loop over all true jets to find the best match
789  matchedJetPt = 0;
790  matchedJetMass = 0;
791  if(truthArray)
792  for(Int_t i=0; i<truthArray->GetEntries(); i++)
793  {
794  AliEmcalJet* truthJet = static_cast<AliEmcalJet*>(truthArray->At(i));
795  if(truthJet->Pt() < 0.15)
796  continue;
797 
798  Double_t deltaEta = (truthJet->Eta()-jet->Eta());
799  Double_t deltaPhi = TMath::Min(TMath::Abs(truthJet->Phi()-jet->Phi()),TMath::TwoPi() - TMath::Abs(truthJet->Phi()-jet->Phi()));
800  Double_t deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
801 
802  // Cut jets too far away
803  if (deltaR > fTrueJetMatchingRadius)
804  continue;
805 
806  // Search for the best match
807  if(deltaR < bestMatchDeltaR)
808  {
809  bestMatchDeltaR = deltaR;
810  matchedJetPt = truthJet->Pt() - truthJet->Area()* trueRho;
811  matchedJetMass = truthJet->M() - truthJet->Area()* trueRhoMass;
812  }
813  }
814  }
815 
816  return bestMatchDeltaR;
817 }
818 
819 
820 //________________________________________________________________________
822 {
824 
826  return;
827 
828  typeHM = 0;
829  typePM = 0;
830 
831  AliAODMCParticle* parton[2];
832  parton[0] = (AliAODMCParticle*) fVtxTagger->IsMCJetParton(fTruthParticleArray, jet, radius); // method 1 (parton matching)
833  parton[1] = (AliAODMCParticle*) fVtxTagger->IsMCJetMeson(fTruthParticleArray, jet, radius); // method 2 (hadron matching)
834 
835  if (parton[0]) {
836  Int_t pdg = TMath::Abs(parton[0]->PdgCode());
837  typePM = pdg;
838  }
839 
840  if (!parton[1])
841  {
842  // No HF jet (according to hadron matching) -- now also separate jets in udg (1) and s-jets (3)
843  if(IsStrangeJet(jet))
844  typeHM = 3;
845  else
846  typeHM = 1;
847  }
848  else {
849  Int_t pdg = TMath::Abs(parton[1]->PdgCode());
850  if(fVtxTagger->IsDMeson(pdg)) typeHM = 4;
851  else if (fVtxTagger->IsBMeson(pdg)) typeHM = 5;
852  }
853 
854  // Set flavour of AliEmcalJet object (set ith bit while i corresponds to type)
856  jet->AddFlavourTag(static_cast<Int_t>(TMath::Power(2, typeHM)));
857 
858 
859  const AliEmcalPythiaInfo* partonsInfo = GetPythiaInfo();
860  typeIC = 0;
861  if (partonsInfo)
862  {
863  // Get primary partons directions
864  Double_t parton1phi = partonsInfo->GetPartonPhi6();
865  Double_t parton1eta = partonsInfo->GetPartonEta6();
866  Double_t parton2phi = partonsInfo->GetPartonPhi7();
867  Double_t parton2eta = partonsInfo->GetPartonEta7();
868 
869 
870  Double_t delta1Eta = (parton1eta-jet->Eta());
871  Double_t delta1Phi = TMath::Min(TMath::Abs(parton1phi-jet->Phi()),TMath::TwoPi() - TMath::Abs(parton1phi-jet->Phi()));
872  Double_t delta1R = TMath::Sqrt(delta1Eta*delta1Eta + delta1Phi*delta1Phi);
873  Double_t delta2Eta = (parton2eta-jet->Eta());
874  Double_t delta2Phi = TMath::Min(TMath::Abs(parton2phi-jet->Phi()),TMath::TwoPi() - TMath::Abs(parton2phi-jet->Phi()));
875  Double_t delta2R = TMath::Sqrt(delta2Eta*delta2Eta + delta2Phi*delta2Phi);
876 
877  // Check if one of the partons if closer than matching criterion
878  Bool_t matched = (delta1R < fJetsCont->GetJetRadius()/2.) || (delta2R < fJetsCont->GetJetRadius()/2.);
879 
880  // Matching criterion fulfilled -> Set flag to closest
881  if(matched)
882  {
883  if(delta1R < delta2R)
884  typeIC = partonsInfo->GetPartonFlag6();
885  else
886  typeIC = partonsInfo->GetPartonFlag7();
887  }
888  }
889 }
890 
891 //________________________________________________________________________
892 void AliAnalysisTaskJetExtractor::GetTrackImpactParameters(const AliVVertex* vtx, AliAODTrack* track, Float_t& d0, Float_t& d0cov, Float_t& z0, Float_t& z0cov)
893 {
894  if (track)
895  {
896  Double_t d0rphiz[2],covd0[3];
897  Bool_t isDCA=track->PropagateToDCA(vtx,InputEvent()->GetMagneticField(),3.0,d0rphiz,covd0);
898  if(isDCA)
899  {
900  d0 = d0rphiz[0];
901  z0 = d0rphiz[1];
902  d0cov = covd0[0];
903  z0cov = covd0[2];
904  }
905  }
906 }
907 
908 //________________________________________________________________________
909 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)
910 {
911  if(!primVtx)
912  return;
913 
914  // Create ESD vertex from the existing AliVVertex
915  Double_t vtxPos[3] = {primVtx->GetX(), primVtx->GetY(), primVtx->GetZ()};
916  Double_t covMatrix[6] = {0};
917  primVtx->GetCovarianceMatrix(covMatrix);
918  AliESDVertex* esdVtx = new AliESDVertex(vtxPos, covMatrix, primVtx->GetChi2(), primVtx->GetNContributors());
919 
920  TClonesArray* secVertexArr = 0;
921  vctr_pair_dbl_int arrDispersion;
922  arrDispersion.reserve(5);
924  {
925  //###########################################################################
926  // ********* Calculate secondary vertices
927  // Derived from AliAnalysisTaskEmcalJetBtagSV
928  secVertexArr = new TClonesArray("AliAODVertex");
929  Int_t nDauRejCount = 0;
930  Int_t nVtx = fVtxTagger->FindVertices(jet,
931  fTracksCont->GetArray(),
932  (AliAODEvent*)InputEvent(),
933  esdVtx,
934  InputEvent()->GetMagneticField(),
935  secVertexArr,
936  0,
937  arrDispersion,
938  nDauRejCount);
939 
940 
941  if(nVtx < 0)
942  {
943  secVertexArr->Clear();
944  delete secVertexArr;
945  return;
946  }
947  //###########################################################################
948  }
949  else // Load HF vertex branch from AOD event, if possible
950  {
951  secVertexArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("VerticesHF"));
952  if(!secVertexArr)
953  return;
954  }
955 
956  // Loop over all potential secondary vertices
957  for(Int_t i=0; i<secVertexArr->GetEntriesFast(); i++)
958  {
959  AliAODVertex* secVtx = (AliAODVertex*)(secVertexArr->UncheckedAt(i));
961  if((strcmp(secVtx->GetParent()->ClassName(), "AliAODRecoDecayHF3Prong")))
962  continue;
963 
964  // Calculate vtx distance
965  Double_t effX = secVtx->GetX() - esdVtx->GetX();
966  Double_t effY = secVtx->GetY() - esdVtx->GetY();
967  //Double_t effZ = secVtx->GetZ() - esdVtx->GetZ();
968 
969  // ##### Vertex properties
970  // vertex dispersion
971  Double_t dispersion = arrDispersion[i].first;
972 
973  // invariant mass
974  Double_t mass = fVtxTagger->GetVertexInvariantMass(secVtx);
975 
976  // signed length
977  Double_t Lxy = TMath::Sqrt(effX*effX + effY*effY);
978  Double_t jetP[3]; jet->PxPyPz(jetP);
979  Double_t signLxy = effX * jetP[0] + effY * jetP[1];
980  if (signLxy < 0.) Lxy *= -1.;
981 
982  Double_t sigmaLxy = 0;
983  AliAODVertex* aodVtx = (AliAODVertex*)(primVtx);
984  if (aodVtx)
985  sigmaLxy = aodVtx->ErrorDistanceXYToVertex(secVtx);
986 
987  // Add secondary vertices if they fulfill the conditions
988  if( (dispersion > fSecondaryVertexMaxDispersion) || (TMath::Abs(secVtx->GetChi2perNDF()) > fSecondaryVertexMaxChi2) )
989  continue;
990 
991  secVtx_X.push_back(secVtx->GetX()); secVtx_Y.push_back(secVtx->GetY()); secVtx_Z.push_back(secVtx->GetZ()); secVtx_Chi2.push_back(secVtx->GetChi2perNDF());
992  secVtx_Dispersion.push_back(dispersion); secVtx_Mass.push_back(mass); secVtx_Lxy.push_back(Lxy); secVtx_SigmaLxy.push_back(sigmaLxy);
993  }
994 
996  {
997  secVertexArr->Clear();
998  delete secVertexArr;
999  }
1000  delete esdVtx;
1001 }
1002 
1003 //________________________________________________________________________
1005 {
1006  // Do hadron matching jet type tagging using mcparticles
1007  // ... if not explicitly deactivated
1008  if (fTruthParticleArray)
1009  {
1010  for(Int_t i=0; i<fTruthParticleArray->GetEntries();i++)
1011  {
1012  AliAODMCParticle* part = (AliAODMCParticle*)fTruthParticleArray->At(i);
1013  if(!part) continue;
1014 
1015  // Check if the particle has strangeness
1016  Int_t absPDG = TMath::Abs(part->PdgCode());
1017  if ((absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
1018  {
1019  // Check if particle is in a radius around the jet
1020  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
1022  continue;
1023  else
1024  return kTRUE;
1025  }
1026  }
1027  }
1028  // As fallback, the MC stack will be tried
1029  else if(MCEvent() && (MCEvent()->Stack()))
1030  {
1031  AliStack* stack = MCEvent()->Stack();
1032  // Go through the whole particle stack
1033  for(Int_t i=0; i<stack->GetNtrack(); i++)
1034  {
1035  TParticle *part = stack->Particle(i);
1036  if(!part) continue;
1037 
1038  // Check if the particle has strangeness
1039  Int_t absPDG = TMath::Abs(part->GetPdgCode());
1040  if ((absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
1041  {
1042  // Check if particle is in a radius around the jet
1043  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
1045  continue;
1046  else
1047  return kTRUE;
1048  }
1049  }
1050  }
1051  return kFALSE;
1052 
1053 }
1054 //________________________________________________________________________
1056 {
1057  // ### Event control plots
1059  FillHistogram("hBackgroundPt", fJetsCont->GetRhoVal(), fCent);
1060 }
1061 
1062 //________________________________________________________________________
1064 {
1065  // ### Jet control plots
1066  FillHistogram("hJetPtRaw", jet->Pt(), fCent);
1067  FillHistogram("hJetPt", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
1068  FillHistogram("hJetPhiEta", jet->Phi(), jet->Eta());
1069  FillHistogram("hJetArea", jet->Area(), fCent);
1070 
1071  // ### Jet constituent plots
1072  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
1073  {
1074  AliVParticle* jetConst = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
1075  if(!jetConst) continue;
1076 
1077  // Constituent eta/phi (relative to jet)
1078  Double_t deltaEta = jet->Eta() - jetConst->Eta();
1079  Double_t deltaPhi = TMath::Min(TMath::Abs(jet->Phi() - jetConst->Phi()), TMath::TwoPi() - TMath::Abs(jet->Phi() - jetConst->Phi()));
1080  if(!((jet->Phi() - jetConst->Phi() < 0) && (jet->Phi() - jetConst->Phi() <= TMath::Pi())))
1081  deltaPhi = -deltaPhi;
1082 
1083  FillHistogram("hConstituentPt", jetConst->Pt(), fCent);
1084  FillHistogram("hConstituentPhiEta", deltaPhi, deltaEta);
1085  }
1086 
1087  // ### Random cone / delta pT plots
1088  const Int_t kNumRandomConesPerEvent = 4;
1089  for(Int_t iCone=0; iCone<kNumRandomConesPerEvent; iCone++)
1090  {
1091  // Throw random cone
1092  Double_t tmpRandConeEta = fJetsCont->GetJetEtaMin() + fRandomGeneratorCones->Rndm()*TMath::Abs(fJetsCont->GetJetEtaMax()-fJetsCont->GetJetEtaMin());
1093  Double_t tmpRandConePhi = fRandomGeneratorCones->Rndm()*TMath::TwoPi();
1094  Double_t tmpRandConePt = 0;
1095  // Fill pT that is in cone
1096  fTracksCont->ResetCurrentID();
1097  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
1098  if(IsTrackInCone(track, tmpRandConeEta, tmpRandConePhi, fJetsCont->GetJetRadius()))
1099  tmpRandConePt += track->Pt();
1100 
1101  // Fill histograms
1102  FillHistogram("hDeltaPt", tmpRandConePt - fJetsCont->GetRhoVal()*fJetsCont->GetJetRadius()*fJetsCont->GetJetRadius()*TMath::Pi(), fCent);
1103  }
1104 }
1105 
1106 //________________________________________________________________________
1108 {
1109  FillHistogram("hTrackPt", track->Pt(), fCent);
1110  FillHistogram("hTrackPhi", track->Phi(), fCent);
1111  FillHistogram("hTrackEta", track->Eta(), fCent);
1112  FillHistogram("hTrackEtaPt", track->Eta(), track->Pt());
1113  FillHistogram("hTrackPhiPt", track->Phi(), track->Pt());
1114  FillHistogram("hTrackPhiEta", track->Phi(), track->Eta());
1115 }
1116 
1117 //________________________________________________________________________
1118 void AliAnalysisTaskJetExtractor::AddPIDInformation(AliVParticle* particle, Float_t& sigITS, Float_t& sigTPC, Float_t& sigTOF, Float_t& sigTRD, Short_t& recoPID, Int_t& truePID)
1119 {
1120  truePID = 9;
1121  if(!particle) return;
1122 
1123  // If we have AODs, retrieve particle PID signals
1124  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(particle);
1125 
1126  if(aodtrack)
1127  {
1128  // Get AOD value from reco
1129  recoPID = aodtrack->GetMostProbablePID();
1130  AliAODPid* pidObj = aodtrack->GetDetPid();
1131  if(!pidObj)
1132  return;
1133 
1134  sigITS = pidObj->GetITSsignal();
1135  sigTPC = pidObj->GetTPCsignal();
1136  sigTOF = pidObj->GetTOFsignal();
1137  sigTRD = pidObj->GetTRDsignal();
1138  }
1139 
1140  // Get truth values if we are on MC
1142  {
1143  for(Int_t i=0; i<fTruthParticleArray->GetEntries();i++)
1144  {
1145  AliAODMCParticle* mcParticle = dynamic_cast<AliAODMCParticle*>(fTruthParticleArray->At(i));
1146  if(!mcParticle) continue;
1147 
1148  if (mcParticle->GetLabel() == particle->GetLabel())
1149  {
1150  if(fSaveTrackPDGCode)
1151  {
1152  truePID = mcParticle->PdgCode();
1153  }
1154  else // Use same convention as for PID in AODs
1155  {
1156  if(TMath::Abs(mcParticle->PdgCode()) == 2212) // proton
1157  truePID = 4;
1158  else if (TMath::Abs(mcParticle->PdgCode()) == 211) // pion
1159  truePID = 2;
1160  else if (TMath::Abs(mcParticle->PdgCode()) == 321) // kaon
1161  truePID = 3;
1162  else if (TMath::Abs(mcParticle->PdgCode()) == 11) // electron
1163  truePID = 0;
1164  else if (TMath::Abs(mcParticle->PdgCode()) == 13) // muon
1165  truePID = 1;
1166  else if (TMath::Abs(mcParticle->PdgCode()) == 700201) // deuteron
1167  truePID = 5;
1168  else if (TMath::Abs(mcParticle->PdgCode()) == 700301) // triton
1169  truePID = 6;
1170  else if (TMath::Abs(mcParticle->PdgCode()) == 700302) // He3
1171  truePID = 7;
1172  else if (TMath::Abs(mcParticle->PdgCode()) == 700202) // alpha
1173  truePID = 8;
1174  else
1175  truePID = 9;
1176  }
1177 
1178  break;
1179  }
1180  }
1181  }
1182 }
1183 
1184 //________________________________________________________________________
1186 {
1187  std::cout << "#########################################\n";
1188  std::cout << "Settings for extractor task " << GetName() << std::endl;
1189  std::cout << std::endl;
1190 
1191  std::cout << "### Event cuts ###" << std::endl;
1192  std::cout << std::endl;
1194  std::cout << Form("* Trigger track requirement: pT=%3.1f-%3.1f GeV/c, labels=%d-%d", fEventCut_TriggerTrackMinPt, fEventCut_TriggerTrackMaxPt, fEventCut_TriggerTrackMinLabel, fEventCut_TriggerTrackMaxLabel) << std::endl;
1195  if(fEventPercentage < 1)
1196  std::cout << Form("* Artificially lowered event efficiency: %f", fEventPercentage) << std::endl;
1198  std::cout << Form("* Random seeds explicitly set: %lu (for event/jet eff.), %lu (RCs)", fRandomSeed, fRandomSeedCones) << std::endl;
1199  std::cout << std::endl;
1200 
1201  std::cout << "### Jet properties ###" << std::endl;
1202  std::cout << "* Jet array name = " << fJetsCont->GetName() << std::endl;
1203  std::cout << "* Rho name = " << fJetsCont->GetRhoName() << std::endl;
1204  std::cout << "* Rho mass name = " << fJetsCont->GetRhoMassName() << std::endl;
1205  std::cout << std::endl;
1206 
1207  std::cout << "### Tree extraction properties ###" << std::endl;
1208  std::vector<Float_t> extractionPtBins = fJetTree->GetExtractionPercentagePtBins();
1209  std::vector<Float_t> extractionPercentages = fJetTree->GetExtractionPercentages();
1210  std::vector<Int_t> extractionHM = fJetTree->GetExtractionJetTypes_HM();
1211  std::vector<Int_t> extractionPM = fJetTree->GetExtractionJetTypes_PM();
1212  if(extractionPercentages.size())
1213  {
1214  std::cout << "* Extraction bins: (";
1215  for(Int_t i=0; i<static_cast<Int_t>(extractionPercentages.size()-1); i++)
1216  std::cout << extractionPtBins[i*2] << "-" << extractionPtBins[i*2+1] << " = " << extractionPercentages[i] << ", ";
1217  std::cout << extractionPtBins[(extractionPercentages.size()-1)*2] << "-" << extractionPtBins[(extractionPercentages.size()-1)*2+1] << " = " << extractionPercentages[(extractionPercentages.size()-1)];
1218 
1219  std::cout << ")" << std::endl;
1220  }
1221  if(extractionHM.size())
1222  {
1223  std::cout << "* Extraction of hadron-matched jets with types: (";
1224  for(Int_t i=0; i<static_cast<Int_t>(extractionHM.size()-1); i++)
1225  std::cout << extractionHM[i] << ", ";
1226  std::cout << extractionHM[extractionHM.size()-1];
1227  std::cout << ")" << std::endl;
1228  }
1229  if(extractionPM.size())
1230  {
1231  std::cout << "* Extraction of parton-matched jets with types: (";
1232  for(Int_t i=0; i<static_cast<Int_t>(extractionPM.size()-1); i++)
1233  std::cout << extractionPM[i] << ", ";
1234  std::cout << extractionPM[extractionPM.size()-1];
1235  std::cout << ")" << std::endl;
1236  }
1237  std::cout << std::endl;
1238 
1239  std::cout << "### Extracted data groups ###" << std::endl;
1241  std::cout << "* Basic event properties (ID, vertex, centrality, bgrd. densities, ...)" << std::endl;
1243  std::cout << "* Jet constituents, basic properties (pt, eta, phi, charge, ...)" << std::endl;
1245  std::cout << "* Jet constituents, IPs" << std::endl;
1247  std::cout << "* Jet constituents, PID" << std::endl;
1249  std::cout << "* MC information (origin, matched jets, ...)" << std::endl;
1251  std::cout << "* Secondary vertices" << std::endl;
1252  if(fJetTree->GetSaveJetShapes())
1253  std::cout << "* Jet shapes (jet mass)" << std::endl;
1255  std::cout << "* Trigger tracks" << std::endl;
1256  std::cout << std::endl;
1257  std::cout << "### Further settings ###" << std::endl;
1258  if(fTruthJetsArrayName != "")
1259  std::cout << Form("* True jet matching active, array=%s, rho=%s, rho_mass=%s", fTruthJetsArrayName.Data(), fTruthJetsRhoName.Data(), fTruthJetsRhoMassName.Data()) << std::endl;
1261  std::cout << Form("* Particle level information available (for jet origin calculation, particle code): %s", fTruthParticleArrayName.Data()) << std::endl;
1262  if(extractionHM.size())
1263  std::cout << Form("* Hadron--jet matching radius=%3.3f", fHadronMatchingRadius) << std::endl;
1265  std::cout << Form("* True particle label range for pt fraction=%d-%d", fTruthMinLabel, fTruthMaxLabel) << std::endl;
1267  std::cout << Form("* Particle PID codes will be PDG codes") << std::endl;
1269  std::cout << Form("* Particle PID codes will follow AliAODPid convention") << std::endl;
1271  std::cout << Form("* AliEmcalJet flavour tag is set from HF-jet tagging") << std::endl;
1272 
1273  std::cout << "#########################################\n\n";
1274 }
1275 
1276 
1277 //########################################################################
1278 // HELPERS
1279 //########################################################################
1280 
1281 //________________________________________________________________________
1282 inline Bool_t AliAnalysisTaskJetExtractor::IsTrackInCone(AliVParticle* track, Double_t eta, Double_t phi, Double_t radius)
1283 {
1284  // This is to use a full cone in phi even at the edges of phi (2pi -> 0) (0 -> 2pi)
1285  Double_t trackPhi = 0.0;
1286  if (track->Phi() > (TMath::TwoPi() - (radius-phi)))
1287  trackPhi = track->Phi() - TMath::TwoPi();
1288  else if (track->Phi() < (phi+radius - TMath::TwoPi()))
1289  trackPhi = track->Phi() + TMath::TwoPi();
1290  else
1291  trackPhi = track->Phi();
1292 
1293  if ( TMath::Abs(trackPhi-phi)*TMath::Abs(trackPhi-phi) + TMath::Abs(track->Eta()-eta)*TMath::Abs(track->Eta()-eta) <= radius*radius)
1294  return kTRUE;
1295 
1296  return kFALSE;
1297 }
1298 
1299 //________________________________________________________________________
1301 {
1302  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1303  if(!tmpHist)
1304  {
1305  AliError(Form("Cannot find histogram <%s> ",key)) ;
1306  return;
1307  }
1308 
1309  tmpHist->Fill(x);
1310 }
1311 
1312 //________________________________________________________________________
1314 {
1315  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1316  if(!tmpHist)
1317  {
1318  AliError(Form("Cannot find histogram <%s> ",key));
1319  return;
1320  }
1321 
1322  if (tmpHist->IsA()->GetBaseClass("TH1"))
1323  static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
1324  else if (tmpHist->IsA()->GetBaseClass("TH2"))
1325  static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
1326 }
1327 
1328 //________________________________________________________________________
1330 {
1331  TH2* tmpHist = static_cast<TH2*>(fOutput->FindObject(key));
1332  if(!tmpHist)
1333  {
1334  AliError(Form("Cannot find histogram <%s> ",key));
1335  return;
1336  }
1337 
1338  tmpHist->Fill(x,y,add);
1339 }
1340 
1341 //________________________________________________________________________
1343 {
1344  TH3* tmpHist = static_cast<TH3*>(fOutput->FindObject(key));
1345  if(!tmpHist)
1346  {
1347  AliError(Form("Cannot find histogram <%s> ",key));
1348  return;
1349  }
1350 
1351  if(add)
1352  tmpHist->Fill(x,y,z,add);
1353  else
1354  tmpHist->Fill(x,y,z);
1355 }
1356 
1357 
1358 //________________________________________________________________________
1359 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)
1360 {
1361  T* tmpHist = new T(name, title, xBins, xMin, xMax);
1362 
1363  tmpHist->GetXaxis()->SetTitle(xTitle);
1364  tmpHist->GetYaxis()->SetTitle(yTitle);
1365  tmpHist->SetOption(options);
1366  tmpHist->SetMarkerStyle(kFullCircle);
1367  tmpHist->Sumw2();
1368 
1369  fOutput->Add(tmpHist);
1370 
1371  return tmpHist;
1372 }
1373 
1374 //________________________________________________________________________
1375 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)
1376 {
1377  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax);
1378  tmpHist->GetXaxis()->SetTitle(xTitle);
1379  tmpHist->GetYaxis()->SetTitle(yTitle);
1380  tmpHist->GetZaxis()->SetTitle(zTitle);
1381  tmpHist->SetOption(options);
1382  tmpHist->SetMarkerStyle(kFullCircle);
1383  tmpHist->Sumw2();
1384 
1385  fOutput->Add(tmpHist);
1386 
1387  return tmpHist;
1388 }
1389 
1390 //________________________________________________________________________
1391 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)
1392 {
1393  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax, zBins, zMin, zMax);
1394  tmpHist->GetXaxis()->SetTitle(xTitle);
1395  tmpHist->GetYaxis()->SetTitle(yTitle);
1396  tmpHist->GetZaxis()->SetTitle(zTitle);
1397  tmpHist->SetOption(options);
1398  tmpHist->SetMarkerStyle(kFullCircle);
1399  tmpHist->Sumw2();
1400 
1401  fOutput->Add(tmpHist);
1402 
1403  return tmpHist;
1404 }
1405 
1406 //________________________________________________________________________
1408 {
1409  // Called once at the end of the analysis.
1410 }
1411 
1412 // ### ADDTASK MACRO
1413 //________________________________________________________________________
1414 AliAnalysisTaskJetExtractor* AliAnalysisTaskJetExtractor::AddTaskJetExtractor(TString trackArray, TString jetArray, TString rhoObject, Double_t jetRadius, TString configFile, AliRDHFJetsCutsVertex* vertexerCuts, const char* taskNameSuffix)
1415 {
1416  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1417  Double_t minJetEta = 0.5;
1418  Double_t minJetPt = 0.15;
1419  Double_t minTrackPt = 0.15;
1420  Double_t minJetAreaPerc = 0.557;
1421  TString suffix = "";
1422  if(taskNameSuffix != 0)
1423  suffix = taskNameSuffix;
1424 
1425  // ###### Load configuration from YAML files
1427  fYAMLConfig.AddConfiguration("alien:///alice/cern.ch/user/r/rhaake/ConfigFiles/JetExtractor_BaseConfig.yaml", "baseConfig");
1428  if(configFile != "")
1429  fYAMLConfig.AddConfiguration(configFile.Data(), "customConfig");
1430  fYAMLConfig.Initialize();
1431 
1432  fYAMLConfig.GetProperty("TrackArray", trackArray, kFALSE);
1433  fYAMLConfig.GetProperty("JetArray", jetArray, kFALSE);
1434  fYAMLConfig.GetProperty("RhoName", rhoObject, kFALSE);
1435  fYAMLConfig.GetProperty("JetRadius", jetRadius);
1436  fYAMLConfig.GetProperty("MinJetEta", minJetEta);
1437  fYAMLConfig.GetProperty("MinJetPt", minJetPt);
1438  fYAMLConfig.GetProperty("MinTrackPt", minTrackPt);
1439  fYAMLConfig.GetProperty("MinJetAreaPerc", minJetAreaPerc);
1440  fYAMLConfig.Print();
1441 
1442  // ###### Task name
1443  TString name("AliAnalysisTaskJetExtractor");
1444  if (jetArray != "") {
1445  name += "_";
1446  name += jetArray;
1447  }
1448  if (rhoObject != "") {
1449  name += "_";
1450  name += rhoObject;
1451  }
1452  if (suffix != "") {
1453  name += "_";
1454  name += suffix;
1455  }
1456 
1457  // ###### Setup task with default settings
1459  myTask->SetNeedEmcalGeom(kFALSE);
1460  myTask->SetVzRange(-10.,10.);
1461 
1462  // Particle container and track pt cut
1463  AliParticleContainer* trackCont = 0;
1464  if(trackArray == "mcparticles")
1465  trackCont = myTask->AddMCParticleContainer(trackArray);
1466  else if(trackArray =="mctracks")
1467  trackCont = myTask->AddParticleContainer(trackArray);
1468  else
1469  trackCont = myTask->AddTrackContainer(trackArray);
1470  trackCont->SetParticlePtCut(minTrackPt);
1471 
1472  // Secondary vertex cuts (default settings from PWGHF)
1473  // (can be overwritten by using myTask->SetVertexerCuts(cuts) from outside macro)
1474  AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts", "default");
1475  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1476  esdTrackCuts->SetMinNClustersTPC(90);
1477  esdTrackCuts->SetMaxChi2PerClusterTPC(4);
1478  esdTrackCuts->SetRequireTPCRefit(kTRUE);
1479  esdTrackCuts->SetRequireITSRefit(kTRUE);
1480  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
1481  esdTrackCuts->SetMinDCAToVertexXY(0.);
1482  esdTrackCuts->SetEtaRange(-0.8, 0.8);
1483  esdTrackCuts->SetPtRange(1.0, 1.e10);
1484  vertexerCuts->AddTrackCuts(esdTrackCuts);
1485  vertexerCuts->SetNprongs(3);
1486  vertexerCuts->SetMinPtHardestTrack(1.0);//default 0.3
1487  vertexerCuts->SetSecVtxWithKF(kFALSE);//default with StrLinMinDist
1488  vertexerCuts->SetImpParCut(0.);//default 0
1489  vertexerCuts->SetDistPrimSec(0.);//default 0
1490  vertexerCuts->SetCospCut(-1);//default -1
1491  myTask->SetVertexerCuts(vertexerCuts);
1492 
1493  // Jet container
1494  AliJetContainer *jetCont = myTask->AddJetContainer(jetArray,6,jetRadius);
1495  if (jetCont) {
1496  jetCont->SetRhoName(rhoObject);
1497  jetCont->SetPercAreaCut(minJetAreaPerc);
1498  jetCont->SetJetPtCut(minJetPt);
1499  jetCont->SetLeadingHadronType(0);
1500  jetCont->SetPtBiasJetTrack(minTrackPt);
1501  jetCont->SetJetEtaLimits(-minJetEta, +minJetEta);
1502  jetCont->ConnectParticleContainer(trackCont);
1503  jetCont->SetMaxTrackPt(1000);
1504  }
1505 
1506  mgr->AddTask(myTask);
1507 
1508  // ###### Connect inputs/outputs
1509  mgr->ConnectInput (myTask, 0, mgr->GetCommonInputContainer() );
1510  mgr->ConnectOutput (myTask, 1, mgr->CreateContainer(Form("%s_histos", name.Data()), AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, Form("%s:ChargedJetsHadronCF", mgr->GetCommonFileName())) );
1511  mgr->ConnectOutput (myTask, 2, mgr->CreateContainer(Form("%s_tree", name.Data()), TTree::Class(), AliAnalysisManager::kOutputContainer, mgr->GetCommonFileName()) );
1512 
1513  return myTask;
1514 }
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 SetParticlePtCut(Double_t cut)
virtual AliVParticle * GetNextAcceptParticle()
void SetSaveConstituents(Bool_t val)
Double_t GetRhoVal() const
const TString & GetRhoName() 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 SetTrackContainer(AliParticleContainer *cont)
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
void SetLeadingHadronType(Int_t t)
Double_t Phi() const
Definition: AliEmcalJet.h:117
std::vector< Float_t > GetExtractionPercentages()
Bool_t fSaveMCInformation
save MC information
Double_t mass
void SetPtBiasJetTrack(Float_t b)
centrality
AliParticleContainer * fTrackContainer
! track container, used to work on tracks
Int_t fTruthMinLabel
min track label to consider it as true particle
Double_t fSecondaryVertexMaxChi2
Max chi2 of secondary vertex (others will be discarded)
Double_t GetJetEtaMax() const
Float_t GetPartonPhi7() const
Double_t fSecondaryVertexMaxDispersion
Max dispersion of secondary vertex (others will be discarded)
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
Double_t fEventCut_TriggerTrackMinPt
Event requirement, trigger track min pT.
std::vector< Float_t > GetExtractionPercentagePtBins()
Bool_t fSaveEventProperties
save general event properties (bgrd. etc.)
std::vector< Float_t > fTriggerTracks_Pt
If trigger track requi. is used -> save pT.
bool GetProperty(std::vector< std::string > propertyPath, const std::string &propertyName, T &property, const bool requiredProperty) const
Bool_t fSaveJetShapes
save jet shapes
void LoadRhoMass(const AliVEvent *event)
void SetPercAreaCut(Float_t p)
Bool_t fSaveConstituentPID
save arrays of constituent PID parameters
Double_t fEventWeight
event weight for each event (only works for JEWEL)
AliStack * stack
void SetVzRange(Double_t min, Double_t max)
const Int_t kMaxNumConstituents
Bool_t IsTrackInCone(AliVParticle *track, Double_t eta, Double_t phi, Double_t radius)
Double_t GetTrueJetPtFraction(AliEmcalJet *jet)
T * AddHistogram2D(const char *name="CustomHistogram", const char *title="NO_TITLE", const char *options="", Int_t xBins=100, Double_t xMin=0.0, Double_t xMax=20.0, Int_t yBins=100, Double_t yMin=0.0, Double_t yMax=20.0, const char *xTitle="x axis", const char *yTitle="y axis", const char *zTitle="z axis")
TString fTruthJetsRhoName
Array name for particle-level rho.
Int_t GetPartonFlag7() const
std::vector< Int_t > fExtractionJetTypes_HM
Jet-types that will be extracted with this tree (hadron matching)
Container for particles within the EMCAL framework.
Float_t GetPartonPhi6() const
Double_t fTrueJetMatchingRadius
Matching radius to true jet.
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
std::vector< Float_t > fTriggerTracks_Phi
If trigger track requi. is used -> save phi.
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
Int_t fEventCut_TriggerTrackMaxLabel
Event requirement, trigger track max label (can be used to selected embedded particles) ...
Bool_t AddJetToTree(AliEmcalJet *jet, Float_t vertexX=0, Float_t vertexY=0, Float_t vertexZ=0, Float_t centrality=0, Long64_t eventID=0, Float_t magField=0, Int_t motherParton=0, Int_t motherHadron=0, Int_t partonInitialCollision=0, Float_t matchDistance=0, Float_t matchedPt=0, Float_t matchedMass=0, Float_t truePtFraction=0, Float_t ptHard=0, Float_t eventWeight=0, Float_t impactParameter=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< Int_t > &trackPID_Truth=DEFAULT_VECTOR_INT, std::vector< Float_t > &triggerTrackPt=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &triggerTrackDeltaEta=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &triggerTrackDeltaPhi=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &trackIP_d0=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &trackIP_z0=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &trackIP_d0cov=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &trackIP_z0cov=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_X=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_Y=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_Z=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_Mass=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_Lxy=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_SigmaLxy=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_Chi2=DEFAULT_VECTOR_FLOAT, std::vector< Float_t > &secVtx_Dispersion=DEFAULT_VECTOR_FLOAT)
Float_t * fBuffer_Const_Phi
! array buffer
Float_t fBuffer_Jet_MC_MatchDistance
! array buffer
Float_t fBuffer_Event_ImpactParameter
! array buffer
const TString & GetRhoMassName() const
TRandom3 * fRandomGenerator
! Random number generator, used for event + jet efficiency
Int_t GetPartonFlag6() const
void SetRhoName(const char *n)
Float_t fBuffer_Jet_MC_MatchedPt
! array buffer
Bool_t fSaveTriggerTracks
save event trigger track
AliParticleContainer * GetParticleContainer() const
TRandom3 * fRandomGenerator
! random generator
AliParticleContainer * fTracksCont
! Tracks
std::vector< Float_t > fExtractionPercentages
Percentages which will be extracted for a given pT bin.
TClonesArray * fTruthParticleArray
! Array of MC particles in event (mcparticles)
Double_t fEventPercentage
percentage (0, 1] which will be extracted
void FillHistogram3D(const char *key, Double_t x, Double_t y, Double_t z, Double_t add=0)
int Int_t
Definition: External.C:63
Int_t fBuffer_Jet_MC_MotherParton
! array buffer
Float_t fBuffer_Event_PtHard
! array buffer
void SetJetPtCut(Float_t cut)
void SetSaveJetShapes(Bool_t val)
unsigned int UInt_t
Definition: External.C:33
Float_t * fBuffer_Const_Charge
! array buffer
float Float_t
Definition: External.C:68
Float_t fBuffer_Event_Weight
! array buffer
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")
TRandom3 * fRandomGeneratorCones
! Random number generator, used for random cones
Int_t * fBuffer_Const_Label
! array buffer
AliParticleContainer * AddParticleContainer(const char *n)
Create new particle container and attach it to the task.
Double_t fImpactParameter
IP for each event (only works for JEWEL)
Float_t fBuffer_JetPhi
! array buffer
Int_t fTruthMaxLabel
max track label to consider it as true particle
Double_t GetMatchedTrueJetObservables(AliEmcalJet *jet, Double_t &matchedJetPt, Double_t &matchedJetMass)
Float_t fBuffer_JetArea
! array buffer
std::vector< Int_t > fExtractionJetTypes_PM
Jet-types that will be extracted with this tree (parton matching)
std::vector< Float_t > fExtractionPercentagePtBins
pT-bins associated with fExtractionPercentages
Double_t fCent
!event centrality
Long64_t fBuffer_Event_ID
! array buffer
Float_t fBuffer_Event_Centrality
! array buffer
Int_t fBuffer_Jet_MC_MotherIC
! array buffer
std::ostream & Print(std::ostream &in, const int index=-1) const
Double_t GetRhoMassVal() const
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.
void AddPIDInformation(AliVParticle *particle, Float_t &sigITS, Float_t &sigTPC, Float_t &sigTOF, Float_t &sigTRD, Short_t &recoPID, Int_t &truePID)
AliMCParticleContainer * AddMCParticleContainer(const char *n)
Create new container for MC particles and attach it to the task.
std::vector< Int_t > GetExtractionJetTypes_PM()
std::vector< Int_t > GetExtractionJetTypes_HM()
AliEmcalJet * GetNextAcceptJet()
Float_t fBuffer_Event_BackgroundDensityMass
! array buffer
AliHFJetsTaggingVertex * fVtxTagger
! class for sec. vertexing
short Short_t
Definition: External.C:23
TString fTruthJetsArrayName
Array name for particle-level jets.
void ConnectParticleContainer(AliParticleContainer *c)
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliRDHFJetsCutsVertex * fVertexerCuts
Cuts used for the vertexer (given in add task macro)
Int_t fEventCut_TriggerTrackMinLabel
Event requirement, trigger track min label (can be used to selected embedded particles) ...
AliJetContainer * fJetContainer
! jet container, used to work on jets
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
TString fTruthJetsRhoMassName
Array name for particle-level rho for mass.
Float_t fBuffer_Jet_MC_MatchedMass
! array buffer
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
AliTrackContainer * AddTrackContainer(const char *n)
Create new track container and attach it to the task.
Bool_t fCalculateSecondaryVertices
Calculate the secondary vertices (instead of loading)
Float_t * fBuffer_Const_ProdVtx_Y
! array buffer
Bool_t fSaveConstituents
save arrays of constituent basic properties
YAML configuration class for AliPhysics.
void SetVertexerCuts(AliRDHFJetsCutsVertex *val)
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
void SetMakeGeneralHistograms(Bool_t g)
void SetNeedEmcalGeom(Bool_t n)
Double_t fHadronMatchingRadius
Matching radius to search for beauty/charm hadrons around jet.
static AliAnalysisTaskJetExtractor * AddTaskJetExtractor(TString trackArray, TString jetArray, TString rhoObject, Double_t jetRadius, TString configFile, AliRDHFJetsCutsVertex *vertexerCuts, const char *taskNameSuffix)
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
int AddConfiguration(std::string configurationFilename, std::string configurationName="")
Bool_t fSaveTrackPDGCode
save PDG code instead of code defined for AOD pid
Declaration of class AliEmcalPythiaInfo.
const char Option_t
Definition: External.C:48
Int_t fBuffer_NumConstituents
! array buffer
void UserCreateOutputObjects()
Main initialization function on the worker.
Float_t fBuffer_Event_Vertex_X
! array buffer
bool Bool_t
Definition: External.C:53
Double_t yMin
Int_t GetNAcceptedParticles() const
void SetRandomGenerator(TRandom3 *gen)
void SetMaxTrackPt(Float_t b)
void SetJetEtaLimits(Float_t min, Float_t max)
Int_t fBuffer_Jet_MC_MotherHadron
! array buffer
Double_t fEventCut_TriggerTrackMaxPt
Event requirement, trigger track max pT.
void SetJetContainer(AliJetContainer *cont)
std::vector< Float_t > fTriggerTracks_Eta
If trigger track requi. is used -> save eta.
Float_t fBuffer_Shape_Mass
! array buffer
Double_t M() const
Definition: AliEmcalJet.h:120
Float_t * fBuffer_Const_ProdVtx_Z
! array buffer
Container for jet within the EMCAL jet framework.
Float_t fBuffer_Jet_MC_TruePtFraction
! array buffer
Class managing creation of a tree containing jets.
Definition: External.C:196
Double_t yMax
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
Float_t fBuffer_Event_Vertex_Z
! array buffer
void SetSaveTriggerTracks(Bool_t val)
Float_t fBuffer_Event_MagneticField
! array buffer