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