AliPhysics  64f4410 (64f4410)
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), fPhysicalPrimariesOnly(kTRUE), 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), fPhysicalPrimariesOnly(kTRUE), 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  if (fPhysicalPrimariesOnly && !aodMCParticle->IsPhysicalPrimary())
158  continue;
159  Float_t trackTheta = 2.*atan(exp(-aodMCParticle->Eta()));
160  new ((*fOutputArray)[particleCount]) AliAODTrack();
161  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetPt(aodMCParticle->Pt());
162  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetPhi(aodMCParticle->Phi());
163  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetTheta(trackTheta); // AliAODTrack cannot set eta directly
164  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetCharge(aodMCParticle->Charge());
165  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetLabel(aodMCParticle->GetLabel());
166  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetIsHybridGlobalConstrainedGlobal();
167  particleCount++;
168  }
169  }
170  }
171 
172  // ################# 2. Create a vertex if there is none (needed by some tasks)
173  if(dynamic_cast<AliESDEvent*>(InputEvent()))
174  {
175  if((dynamic_cast<AliESDEvent*>(InputEvent()))->GetPrimaryVertexTracks())
176  if(!(dynamic_cast<AliESDEvent*>(InputEvent()))->GetPrimaryVertexTracks()->GetNContributors())
177  static_cast<AliESDEvent*>(fInputEvent)->SetPrimaryVertexTracks(new AliESDVertex(0.,0., 100));
178  }
179  else if(dynamic_cast<AliAODEvent*>(InputEvent()))
180  {
181  if( (!(dynamic_cast<AliAODEvent*>(InputEvent()))->GetPrimaryVertex()) || (!(dynamic_cast<AliAODEvent*>(InputEvent()))->GetPrimaryVertex()->GetNContributors()) )
182  {
183  Double_t* p = new Double_t[3];
184  p[0] = 0.; p[1] = 0.; p[2] = 0.; // for backwards compatibility
185  AliAODVertex* vertex = new AliAODVertex(p,1.);
186  vertex->SetNContributors(100);
187  vertex->SetName("PrimaryVertex");
188  static_cast<AliAODEvent*>(fInputEvent)->AddVertex(vertex);
189  }
190  }
191 
192  // ################# 3. Create toy event
193  if(fUseMixedEvent) // get underlying event from mixed event files
194  {
195  // if input tree not loaded or index at the end, get next tree file
197  {
199  if (!fMixedEvent_Tree)
200  {
201  AliError(Form("Could not get tree %s in file!", fMixedEvent_TreeName.Data()));
202  return;
203  }
204  }
205 
206  fBuffer_NumTracks = 0; // Failsafe: Set num tracks to 0 if it is not read by GetEntry
208 
209  // Loop over tracks from event
210  for(Int_t j=0; j<fBuffer_NumTracks; j++)
211  {
212  Float_t trackTheta = 2.*atan(exp(-fBuffer_TrackEta[j]));
213 
214  // Add basic particle to event
215  new ((*fOutputArray)[particleCount]) AliAODTrack();
216  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetPt(fBuffer_TrackPt[j]);
217  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetPhi(fBuffer_TrackPhi[j]);
218  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetTheta(trackTheta); // AliAODTrack cannot set eta directly
219  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetCharge(fBuffer_TrackCharge[j]);
220  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetLabel(100000 + j);
221  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetIsHybridGlobalConstrainedGlobal();
222  particleCount++;
223  }
224 
226  }
227  else // Simple toy
228  {
229  Int_t multiplicity = (Int_t)fDistributionMultiplicity->GetRandom();
230  for(Int_t i=0;i<multiplicity; i++)
231  {
232  Double_t trackPt = fDistributionPt->GetRandom();
233  Double_t trackEta = 0;
234  Double_t trackPhi = 0;
235  static_cast<TH2*>(fDistributionEtaPhi)->GetRandom2(trackPhi, trackEta);
236  Double_t trackTheta = 2.*atan(exp(-trackEta));
237  Double_t trackCharge = fRandom->Rndm() - 0.5;
238 
239  if(trackCharge>0) trackCharge = 1; else trackCharge = -1;
240 
241  // Add flow to particle
243  trackPhi = AddFlow(trackPhi, trackPt);
244 
245 
246  // Add basic particle to event
247  new ((*fOutputArray)[particleCount]) AliAODTrack();
248  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetPt(trackPt);
249  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetPhi(trackPhi);
250  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetTheta(trackTheta); // AliAODTrack cannot set eta directly
251  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetCharge(trackCharge);
252  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetLabel(100000 + i);
253  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetIsHybridGlobalConstrainedGlobal();
254  particleCount++;
255  }
256  }
257 
258 
259 }
260 
261 //_____________________________________________________________________________________________________
263 {
264  for(Int_t iTrack=0; iTrack<fOutputArray->GetEntries(); iTrack++)
265  {
266  AliAODTrack* track = static_cast<AliAODTrack*>(fOutputArray->At(iTrack));
267  FillHistogram("hTrackPt", track->Pt());
268  FillHistogram("hTrackPhiEta", track->Phi(), track->Eta());
269  }
270  FillHistogram("hMultiplicity", fOutputArray->GetEntries());
271 }
272 
273 //_____________________________________________________________________________________________________
275 {
276  // ## Form file name
279  fileName = Form("%s%d.root", fMixedEvent_BaseFolder.Data(), fMixedEvent_CurrentFileID);
280 
281  AliInfo(Form("Opening mixed event file: %s", fileName.Data()));
282  // ## Check if file exists
283  if (fileName.BeginsWith("alien://") && !gGrid)
284  {
285  AliInfo("Trying to connect to AliEn ...");
286  TGrid::Connect("alien://");
287  }
288  if (gSystem->AccessPathName(fileName)) {
289  AliError(Form("File %s does not exist!", fileName.Data()));
290  return 0;
291  }
292 
293  // ## Open mixed event file and get tree object
295 
296  fMixedEvent_CurrentFile = TFile::Open(fileName);
298  {
299  AliError(Form("Unable to open file: %s!", fileName.Data()));
300  return 0;
301  }
302  TTree* tree = static_cast<TTree*>(fMixedEvent_CurrentFile->Get(fMixedEvent_TreeName.Data()));
303  if(tree)
304  {
305  tree->SetBranchAddress("NumTracks", &fBuffer_NumTracks);
306  tree->SetBranchAddress("Track_Pt", fBuffer_TrackPt);
307  tree->SetBranchAddress("Track_Phi", fBuffer_TrackPhi);
308  tree->SetBranchAddress("Track_Eta", fBuffer_TrackEta);
309  tree->SetBranchAddress("Track_Charge", fBuffer_TrackCharge);
310  }
311 
313 
314  return tree;
315 }
316 
317 //_____________________________________________________________________________________________________
319 {
320  // adapted from AliFlowTrackSimple
321  Double_t precisionPhi = 1e-10;
322  Int_t maxNumberOfIterations = 200;
323 
324  Double_t phi0=phi;
325  Double_t f=0.;
326  Double_t fp=0.;
327  Double_t phiprev=0.;
328  Int_t ptBin = 0;
329 
330  // Evaluate V2 for track pt/centrality
331  Double_t v2 = 0;
332  if(fDistributionV2)
333  {
334  ptBin = fDistributionV2->GetXaxis()->FindBin(pt);
335  if(ptBin>fDistributionV2->GetNbinsX())
336  v2 = fDistributionV2->GetBinContent(fDistributionV2->GetNbinsX(), fDistributionV2->GetYaxis()->FindBin(fToyCent));
337  else if(ptBin>0)
338  v2 = fDistributionV2->GetBinContent(ptBin, fDistributionV2->GetYaxis()->FindBin(fToyCent));
339  }
340 
341  // Evaluate V3 for track pt/centrality
342  Double_t v3 = 0;
343  if(fDistributionV3)
344  {
345  ptBin = fDistributionV3->GetXaxis()->FindBin(pt);
346  if(ptBin>fDistributionV3->GetNbinsX())
347  v3 = fDistributionV3->GetBinContent(fDistributionV3->GetNbinsX(), fDistributionV3->GetYaxis()->FindBin(fToyCent));
348  else if(ptBin>0)
349  v3 = fDistributionV3->GetBinContent(ptBin, fDistributionV3->GetYaxis()->FindBin(fToyCent));
350  }
351 
352  // Evaluate V4 for track pt/centrality
353  Double_t v4 = 0;
354  if(fDistributionV4)
355  {
356  ptBin = fDistributionV4->GetXaxis()->FindBin(pt);
357  if(ptBin>fDistributionV4->GetNbinsX())
358  v4 = fDistributionV4->GetBinContent(fDistributionV4->GetNbinsX(), fDistributionV4->GetYaxis()->FindBin(fToyCent));
359  else if(ptBin>0)
360  v4 = fDistributionV4->GetBinContent(ptBin, fDistributionV4->GetYaxis()->FindBin(fToyCent));
361  }
362 
363  // Evaluate V5 for track pt/centrality
364  Double_t v5 = 0;
365  if(fDistributionV5)
366  {
367  ptBin = fDistributionV5->GetXaxis()->FindBin(pt);
368  if(ptBin>fDistributionV5->GetNbinsX())
369  v5 = fDistributionV5->GetBinContent(fDistributionV5->GetNbinsX(), fDistributionV5->GetYaxis()->FindBin(fToyCent));
370  else if(ptBin>0)
371  v5 = fDistributionV5->GetBinContent(ptBin, fDistributionV5->GetYaxis()->FindBin(fToyCent));
372  }
373 
374  // Add all v's
375  for (Int_t i=0; i<maxNumberOfIterations; i++)
376  {
377  phiprev=phi; //store last value for comparison
378  f = phi-phi0
379  + v2*TMath::Sin(2.*(phi-(fEPV0+(TMath::Pi()/2.))))
380  +2./3.*v3*TMath::Sin(3.*(phi-fRandomPsi3))
381  +0.5 *v4*TMath::Sin(4.*(phi-fRandomPsi4))
382  +0.4 *v5*TMath::Sin(5.*(phi-fRandomPsi5));
383 
384  fp = 1.0+2.0*(
385  +v2*TMath::Cos(2.*(phi-(fEPV0+(TMath::Pi()/2.))))
386  +v3*TMath::Cos(3.*(phi-fRandomPsi3))
387  +v4*TMath::Cos(4.*(phi-fRandomPsi4))
388  +v5*TMath::Cos(5.*(phi-fRandomPsi5))); //first derivative
389 
390  phi -= f/fp;
391  if (TMath::AreEqualAbs(phiprev,phi,precisionPhi)) break;
392  }
393 
394  return phi;
395 }
396 
397 
398 //________________________________________________________________________
400 {
401  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
402  if(!tmpHist)
403  {
404  AliError(Form("Cannot find histogram <%s> ",key)) ;
405  return;
406  }
407 
408  tmpHist->Fill(x);
409 }
410 
411 //________________________________________________________________________
413 {
414  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
415  if(!tmpHist)
416  {
417  AliError(Form("Cannot find histogram <%s> ",key));
418  return;
419  }
420 
421  if (tmpHist->IsA()->GetBaseClass("TH1"))
422  static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
423  else if (tmpHist->IsA()->GetBaseClass("TH2"))
424  static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
425 }
426 
427 //________________________________________________________________________
429 {
430  TH2* tmpHist = static_cast<TH2*>(fOutput->FindObject(key));
431  if(!tmpHist)
432  {
433  AliError(Form("Cannot find histogram <%s> ",key));
434  return;
435  }
436 
437  tmpHist->Fill(x,y,add);
438 }
439 
440 //________________________________________________________________________
441 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)
442 {
443  T* tmpHist = new T(name, title, xBins, xMin, xMax);
444 
445  tmpHist->GetXaxis()->SetTitle(xTitle);
446  tmpHist->GetYaxis()->SetTitle(yTitle);
447  tmpHist->SetOption(options);
448  tmpHist->SetMarkerStyle(kFullCircle);
449  tmpHist->Sumw2();
450 
451  fOutput->Add(tmpHist);
452 
453  return tmpHist;
454 }
455 
456 //________________________________________________________________________
457 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)
458 {
459  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax);
460  tmpHist->GetXaxis()->SetTitle(xTitle);
461  tmpHist->GetYaxis()->SetTitle(yTitle);
462  tmpHist->GetZaxis()->SetTitle(zTitle);
463  tmpHist->SetOption(options);
464  tmpHist->SetMarkerStyle(kFullCircle);
465  tmpHist->Sumw2();
466 
467  fOutput->Add(tmpHist);
468 
469  return tmpHist;
470 }
471 
472 
473 //________________________________________________________________________
475 {
476  // Called once at the end of the analysis.
477 }
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