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