AliPhysics  59e0e03 (59e0e03)
AliAnalysisTaskChargedJetsHadronToy.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2016, 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 #include <iostream>
16 #include <TRandom3.h>
17 #include <AliLog.h>
18 #include <TString.h>
19 #include <TMath.h>
20 #include <TF1.h>
21 #include <AliEmcalJet.h>
22 #include <AliESDEvent.h>
23 #include <AliAODEvent.h>
24 #include "AliESDVertex.h"
25 #include "AliAODVertex.h"
26 #include <AliAODTrack.h>
27 #include <TClonesArray.h>
28 #include <AliAODMCParticle.h>
29 
30 #include <TFile.h>
31 #include <TGrid.h>
32 #include <TSystem.h>
33 
34 
37 
41 //_____________________________________________________________________________________________________
43  AliAnalysisTaskEmcalJet("AliAnalysisTaskChargedJetsHadronToy", kTRUE), fDistributionMultiplicity(0), fDistributionPt(0), fDistributionEtaPhi(0), fMinCentrality(0), fMaxCentrality(10), fDistributionV2(0), fDistributionV3(0), fDistributionV4(0), fDistributionV5(0), fUseMixedEvent(0), fMixedEvent_Tree(0), fMixedEvent_CurrentFile(0), fMixedEvent_CurrentFileID(-1), fMixedEvent_BaseFolder(""), fMixedEvent_TreeName("ME_tree"), fMixedEvent_CurrentEventID(0), fMixedEvent_NumTotalFiles(30), fBuffer_NumTracks(0), fBuffer_TrackPt(0), fBuffer_TrackPhi(0), fBuffer_TrackEta(0), fBuffer_TrackCharge(0), fInputArrayName(""), fOutputArrayName(""), fInputArray(0), fOutputArray(0), fRandom(), fToyCent(0), fRandomPsi3(0), fRandomPsi4(0), fRandomPsi5(0)
44 {
45  // constructor
46  fBuffer_TrackPt = new Float_t[10000];
47  fBuffer_TrackPhi = new Float_t[10000];
48  fBuffer_TrackEta = new Float_t[10000];
49  fBuffer_TrackCharge = new Short_t[10000];
50 
51 }
52 
53 
54 //_____________________________________________________________________________________________________
56  AliAnalysisTaskEmcalJet(name, kTRUE), fDistributionMultiplicity(0), fDistributionPt(0), fDistributionEtaPhi(0), fMinCentrality(0), fMaxCentrality(10), fDistributionV2(0), fDistributionV3(0), fDistributionV4(0), fDistributionV5(0), fUseMixedEvent(0), fMixedEvent_Tree(0), fMixedEvent_CurrentFile(0), fMixedEvent_CurrentFileID(-1), fMixedEvent_BaseFolder(""), fMixedEvent_TreeName("ME_tree"), fMixedEvent_CurrentEventID(0), fMixedEvent_NumTotalFiles(30), fBuffer_NumTracks(0), fBuffer_TrackPt(0), fBuffer_TrackPhi(0), fBuffer_TrackEta(0), fBuffer_TrackCharge(0), fInputArrayName(""), fOutputArrayName(""), fInputArray(0), fOutputArray(0), fRandom(), fToyCent(0), fRandomPsi3(0), fRandomPsi4(0), fRandomPsi5(0)
57 {
58  // constructor
59  fBuffer_TrackPt = new Float_t[10000];
60  fBuffer_TrackPhi = new Float_t[10000];
61  fBuffer_TrackEta = new Float_t[10000];
62  fBuffer_TrackCharge = new Short_t[10000];
63 }
64 
65 //_____________________________________________________________________________________________________
67 {
68  // destructor
71  if(fDistributionPt)
72  delete fDistributionPt;
74  delete fDistributionEtaPhi;
75 
76  delete fBuffer_TrackPt;
77  delete fBuffer_TrackPhi;
78  delete fBuffer_TrackEta;
79  delete fBuffer_TrackCharge;
80 
81 }
82 
83 
84 //_____________________________________________________________________________________________________
86 {
88 
89  fRandom = new TRandom3(0);
90  AddHistogram1D<TH1D>("hTrackPt", "Tracks p_{T} distribution", "", 300, 0., 300., "p_{T} (GeV/c)", "dN^{Tracks}/dp_{T}");
91  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");
92  AddHistogram1D<TH1D>("hMultiplicity", "Number of tracks in acceptance vs. centrality", "", 500, 0., 5000., "N tracks","dN^{Events}/dN^{Tracks}");
93 
94  PostData(1, fOutput);
95 }
96 
97 
98 //_____________________________________________________________________________________________________
100 {
102 
103  // Check if input array can be loaded
104  if(!fInputArrayName.IsNull())
105  {
106  fInputArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fInputArrayName.Data())));
107  if(!fInputArray)
108  {
109  AliFatal(Form("Input array '%s' not found!", fInputArrayName.Data()));
110  }
111  }
112 
113  // Check if output arrays can be created
114  if((InputEvent()->FindListObject(Form("%s", fOutputArrayName.Data()))))
115  AliFatal(Form("Output array '%s' already exists in the event! Rename it.", fOutputArrayName.Data()));
116 
117  fOutputArray = new TClonesArray("AliAODTrack");
118  fOutputArray->SetName(fOutputArrayName.Data());
119  fInputEvent->AddObject(fOutputArray);
120 }
121 
122 //_____________________________________________________________________________________________________
124 {
125  AssembleEvent();
126  CreateQAPlots();
127  return kTRUE;
128 }
129 
130 //________________________________________________________________________
132 {
133  fRandomPsi3 = fRandom->Rndm()*TMath::Pi(); // once per event, create a random value dedicated for Psi3
134  fRandomPsi4 = fRandom->Rndm()*TMath::Pi(); // once per event, create a random value dedicated for Psi4
135  fRandomPsi5 = fRandom->Rndm()*TMath::Pi(); // once per event, create a random value dedicated for Psi5
136  fToyCent = fMinCentrality + (fMaxCentrality-fMinCentrality)*fRandom->Rndm(); // centrality value (flat from selected range)
137 
138  // Create the event from the several inputs and run the jet finder
139 
140  // ################# 1. Add input tracks (if available)
141  Int_t particleCount = 0;
142  if(fInputArray)
143  {
144  for(Int_t iPart=0; iPart<fInputArray->GetEntries(); iPart++)
145  {
146  // Load AOD tracks or AOD MC particles (for MC gen train)
147  AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(fInputArray->At(iPart));
148  AliAODMCParticle* aodMCParticle = dynamic_cast<AliAODMCParticle*>(fInputArray->At(iPart));
149 
150  if(aodTrack)
151  {
152  new ((*fOutputArray)[particleCount]) AliAODTrack(*aodTrack);
153  particleCount++;
154  }
155  else if(aodMCParticle)
156  {
157  Float_t trackTheta = 2.*atan(exp(-aodMCParticle->Eta()));
158  new ((*fOutputArray)[particleCount]) AliAODTrack();
159  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetPt(aodMCParticle->Pt());
160  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetPhi(aodMCParticle->Phi());
161  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetTheta(trackTheta); // AliAODTrack cannot set eta directly
162  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetCharge(aodMCParticle->Charge());
163  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetLabel(aodMCParticle->GetLabel());
164  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetIsHybridGlobalConstrainedGlobal();
165  particleCount++;
166  }
167  }
168  }
169 
170  // ################# 2. Create a vertex if there is none (needed by some tasks)
171  if(dynamic_cast<AliESDEvent*>(InputEvent()))
172  {
173  if((dynamic_cast<AliESDEvent*>(InputEvent()))->GetPrimaryVertexTracks())
174  if(!(dynamic_cast<AliESDEvent*>(InputEvent()))->GetPrimaryVertexTracks()->GetNContributors())
175  static_cast<AliESDEvent*>(fInputEvent)->SetPrimaryVertexTracks(new AliESDVertex(0.,0., 100));
176  }
177  else if(dynamic_cast<AliAODEvent*>(InputEvent()))
178  {
179  if( (!(dynamic_cast<AliAODEvent*>(InputEvent()))->GetPrimaryVertex()) || (!(dynamic_cast<AliAODEvent*>(InputEvent()))->GetPrimaryVertex()->GetNContributors()) )
180  {
181  Double_t* p = new Double_t[3];
182  p[0] = 0.; p[1] = 0.; p[2] = 0.; // for backwards compatibility
183  AliAODVertex* vertex = new AliAODVertex(p,1.);
184  vertex->SetNContributors(100);
185  vertex->SetName("PrimaryVertex");
186  static_cast<AliAODEvent*>(fInputEvent)->AddVertex(vertex);
187  }
188  }
189 
190  // ################# 3. Create toy event
191  if(fUseMixedEvent) // get underlying event from mixed event files
192  {
193  // if input tree not loaded or index at the end, get next tree file
195  {
197  if (!fMixedEvent_Tree)
198  {
199  AliError(Form("Could not get tree %s in file!", fMixedEvent_TreeName.Data()));
200  return;
201  }
202  }
203 
204  fBuffer_NumTracks = 0; // Failsafe: Set num tracks to 0 if it is not read by GetEntry
206 
207  // Loop over tracks from event
208  for(Int_t j=0; j<fBuffer_NumTracks; j++)
209  {
210  Float_t trackTheta = 2.*atan(exp(-fBuffer_TrackEta[j]));
211 
212  // Add basic particle to event
213  new ((*fOutputArray)[particleCount]) AliAODTrack();
214  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetPt(fBuffer_TrackPt[j]);
215  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetPhi(fBuffer_TrackPhi[j]);
216  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetTheta(trackTheta); // AliAODTrack cannot set eta directly
217  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetCharge(fBuffer_TrackCharge[j]);
218  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetLabel(100000 + j);
219  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetIsHybridGlobalConstrainedGlobal();
220  particleCount++;
221  }
222 
224  }
225  else // Simple toy
226  {
227  Int_t multiplicity = (Int_t)fDistributionMultiplicity->GetRandom();
228  for(Int_t i=0;i<multiplicity; i++)
229  {
230  Double_t trackPt = fDistributionPt->GetRandom();
231  Double_t trackEta = 0;
232  Double_t trackPhi = 0;
233  static_cast<TH2*>(fDistributionEtaPhi)->GetRandom2(trackPhi, trackEta);
234  Double_t trackTheta = 2.*atan(exp(-trackEta));
235  Double_t trackCharge = fRandom->Rndm() - 0.5;
236 
237  if(trackCharge>0) trackCharge = 1; else trackCharge = -1;
238 
239  // Add flow to particle
241  trackPhi = AddFlow(trackPhi, trackPt);
242 
243 
244  // Add basic particle to event
245  new ((*fOutputArray)[particleCount]) AliAODTrack();
246  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetPt(trackPt);
247  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetPhi(trackPhi);
248  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetTheta(trackTheta); // AliAODTrack cannot set eta directly
249  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetCharge(trackCharge);
250  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetLabel(100000 + i);
251  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetIsHybridGlobalConstrainedGlobal();
252  particleCount++;
253  }
254  }
255 
256 
257 }
258 
259 //_____________________________________________________________________________________________________
261 {
262  for(Int_t iTrack=0; iTrack<fOutputArray->GetEntries(); iTrack++)
263  {
264  AliAODTrack* track = static_cast<AliAODTrack*>(fOutputArray->At(iTrack));
265  FillHistogram("hTrackPt", track->Pt());
266  FillHistogram("hTrackPhiEta", track->Phi(), track->Eta());
267  }
268  FillHistogram("hMultiplicity", fOutputArray->GetEntries());
269 }
270 
271 //_____________________________________________________________________________________________________
273 {
274  // ## Form file name
277  fileName = Form("%s%d.root", fMixedEvent_BaseFolder.Data(), fMixedEvent_CurrentFileID);
278 
279  AliInfo(Form("Opening mixed event file: %s", fileName.Data()));
280  // ## Check if file exists
281  if (fileName.BeginsWith("alien://") && !gGrid)
282  {
283  AliInfo("Trying to connect to AliEn ...");
284  TGrid::Connect("alien://");
285  }
286  if (gSystem->AccessPathName(fileName)) {
287  AliError(Form("File %s does not exist!", fileName.Data()));
288  return 0;
289  }
290 
291  // ## Open mixed event file and get tree object
293 
294  fMixedEvent_CurrentFile = TFile::Open(fileName);
296  {
297  AliError(Form("Unable to open file: %s!", fileName.Data()));
298  return 0;
299  }
300  TTree* tree = static_cast<TTree*>(fMixedEvent_CurrentFile->Get(fMixedEvent_TreeName.Data()));
301  if(tree)
302  {
303  tree->SetBranchAddress("NumTracks", &fBuffer_NumTracks);
304  tree->SetBranchAddress("Track_Pt", fBuffer_TrackPt);
305  tree->SetBranchAddress("Track_Phi", fBuffer_TrackPhi);
306  tree->SetBranchAddress("Track_Eta", fBuffer_TrackEta);
307  tree->SetBranchAddress("Track_Charge", fBuffer_TrackCharge);
308  }
309 
311 
312  return tree;
313 }
314 
315 //_____________________________________________________________________________________________________
317 {
318  // adapted from AliFlowTrackSimple
319  Double_t precisionPhi = 1e-10;
320  Int_t maxNumberOfIterations = 200;
321 
322  Double_t phi0=phi;
323  Double_t f=0.;
324  Double_t fp=0.;
325  Double_t phiprev=0.;
326  Int_t ptBin = 0;
327 
328  // Evaluate V2 for track pt/centrality
329  Double_t v2 = 0;
330  if(fDistributionV2)
331  {
332  ptBin = fDistributionV2->GetXaxis()->FindBin(pt);
333  if(ptBin>fDistributionV2->GetNbinsX())
334  v2 = fDistributionV2->GetBinContent(fDistributionV2->GetNbinsX(), fDistributionV2->GetYaxis()->FindBin(fToyCent));
335  else if(ptBin>0)
336  v2 = fDistributionV2->GetBinContent(ptBin, fDistributionV2->GetYaxis()->FindBin(fToyCent));
337  }
338 
339  // Evaluate V3 for track pt/centrality
340  Double_t v3 = 0;
341  if(fDistributionV3)
342  {
343  ptBin = fDistributionV3->GetXaxis()->FindBin(pt);
344  if(ptBin>fDistributionV3->GetNbinsX())
345  v3 = fDistributionV3->GetBinContent(fDistributionV3->GetNbinsX(), fDistributionV3->GetYaxis()->FindBin(fToyCent));
346  else if(ptBin>0)
347  v3 = fDistributionV3->GetBinContent(ptBin, fDistributionV3->GetYaxis()->FindBin(fToyCent));
348  }
349 
350  // Evaluate V4 for track pt/centrality
351  Double_t v4 = 0;
352  if(fDistributionV4)
353  {
354  ptBin = fDistributionV4->GetXaxis()->FindBin(pt);
355  if(ptBin>fDistributionV4->GetNbinsX())
356  v4 = fDistributionV4->GetBinContent(fDistributionV4->GetNbinsX(), fDistributionV4->GetYaxis()->FindBin(fToyCent));
357  else if(ptBin>0)
358  v4 = fDistributionV4->GetBinContent(ptBin, fDistributionV4->GetYaxis()->FindBin(fToyCent));
359  }
360 
361  // Evaluate V5 for track pt/centrality
362  Double_t v5 = 0;
363  if(fDistributionV5)
364  {
365  ptBin = fDistributionV5->GetXaxis()->FindBin(pt);
366  if(ptBin>fDistributionV5->GetNbinsX())
367  v5 = fDistributionV5->GetBinContent(fDistributionV5->GetNbinsX(), fDistributionV5->GetYaxis()->FindBin(fToyCent));
368  else if(ptBin>0)
369  v5 = fDistributionV5->GetBinContent(ptBin, fDistributionV5->GetYaxis()->FindBin(fToyCent));
370  }
371 
372  // Add all v's
373  for (Int_t i=0; i<maxNumberOfIterations; i++)
374  {
375  phiprev=phi; //store last value for comparison
376  f = phi-phi0
377  + v2*TMath::Sin(2.*(phi-(fEPV0+(TMath::Pi()/2.))))
378  +2./3.*v3*TMath::Sin(3.*(phi-fRandomPsi3))
379  +0.5 *v4*TMath::Sin(4.*(phi-fRandomPsi4))
380  +0.4 *v5*TMath::Sin(5.*(phi-fRandomPsi5));
381 
382  fp = 1.0+2.0*(
383  +v2*TMath::Cos(2.*(phi-(fEPV0+(TMath::Pi()/2.))))
384  +v3*TMath::Cos(3.*(phi-fRandomPsi3))
385  +v4*TMath::Cos(4.*(phi-fRandomPsi4))
386  +v5*TMath::Cos(5.*(phi-fRandomPsi5))); //first derivative
387 
388  phi -= f/fp;
389  if (TMath::AreEqualAbs(phiprev,phi,precisionPhi)) break;
390  }
391 
392  return phi;
393 }
394 
395 
396 //________________________________________________________________________
398 {
399  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
400  if(!tmpHist)
401  {
402  AliError(Form("Cannot find histogram <%s> ",key)) ;
403  return;
404  }
405 
406  tmpHist->Fill(x);
407 }
408 
409 //________________________________________________________________________
411 {
412  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
413  if(!tmpHist)
414  {
415  AliError(Form("Cannot find histogram <%s> ",key));
416  return;
417  }
418 
419  if (tmpHist->IsA()->GetBaseClass("TH1"))
420  static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
421  else if (tmpHist->IsA()->GetBaseClass("TH2"))
422  static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
423 }
424 
425 //________________________________________________________________________
427 {
428  TH2* tmpHist = static_cast<TH2*>(fOutput->FindObject(key));
429  if(!tmpHist)
430  {
431  AliError(Form("Cannot find histogram <%s> ",key));
432  return;
433  }
434 
435  tmpHist->Fill(x,y,add);
436 }
437 
438 //________________________________________________________________________
439 template <class T> T* AliAnalysisTaskChargedJetsHadronToy::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)
440 {
441  T* tmpHist = new T(name, title, xBins, xMin, xMax);
442 
443  tmpHist->GetXaxis()->SetTitle(xTitle);
444  tmpHist->GetYaxis()->SetTitle(yTitle);
445  tmpHist->SetOption(options);
446  tmpHist->SetMarkerStyle(kFullCircle);
447  tmpHist->Sumw2();
448 
449  fOutput->Add(tmpHist);
450 
451  return tmpHist;
452 }
453 
454 //________________________________________________________________________
455 template <class T> T* AliAnalysisTaskChargedJetsHadronToy::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)
456 {
457  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax);
458  tmpHist->GetXaxis()->SetTitle(xTitle);
459  tmpHist->GetYaxis()->SetTitle(yTitle);
460  tmpHist->GetZaxis()->SetTitle(zTitle);
461  tmpHist->SetOption(options);
462  tmpHist->SetMarkerStyle(kFullCircle);
463  tmpHist->Sumw2();
464 
465  fOutput->Add(tmpHist);
466 
467  return tmpHist;
468 }
469 
470 
471 //________________________________________________________________________
473 {
474  // Called once at the end of the analysis.
475 }
TH2 * fDistributionV5
Distribution for v4 in bins of pt and centrality.
TH2 * fDistributionV3
Distribution for v2 in bins of pt and centrality.
double Double_t
Definition: External.C:58
Double_t fRandomPsi4
eventwise calculated psi 3
Float_t * fBuffer_TrackPhi
buffer for track pt array
const char * title
Definition: MakeQAPdf.C:27
Short_t * fBuffer_TrackCharge
buffer for track eta array
TString fileName
TSystem * gSystem
TClonesArray * fOutputArray
input array containing tracks from events
Double_t fEPV0
!event plane V0
Bool_t fUseMixedEvent
Distribution for v5 in bins of pt and centrality.
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")
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
TRandom3 * fRandom
input array containing tracks from events
int Int_t
Definition: External.C:63
TString fInputArrayName
buffer for track charge array
float Float_t
Definition: External.C:68
TH2 * fDistributionV4
Distribution for v3 in bins of pt and centrality.
Double_t fRandomPsi5
eventwise calculated psi 4
short Short_t
Definition: External.C:23
AliEmcalList * fOutput
!output list
Definition: External.C:220
Double_t fRandomPsi3
eventwise centrality value
Base task in the EMCAL jet framework.
TFile * fMixedEvent_CurrentFile
ME: The tree of the mixed event.
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")
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
Double_t yMin
Definition: External.C:196
Double_t yMax
void ExecOnce()
Perform steps needed to initialize the analysis.
Float_t * fBuffer_TrackEta
buffer for track phi array