AliPhysics  5403132 (5403132)
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 
29 
32 
36 
37 //_____________________________________________________________________________________________________
39  AliAnalysisTaskEmcalJet("AliAnalysisTaskChargedJetsHadronToy", kTRUE), fDistributionMultiplicity(0), fDistributionPt(0), fDistributionEtaPhi(0), fMinCentrality(0), fMaxCentrality(10), fDistributionV2(0), fDistributionV3(0), fDistributionV4(0), fDistributionV5(0), fInputArrayName(""), fOutputArrayName(""), fInputArray(0), fOutputArray(0), fRandom(), fToyCent(0), fRandomPsi3(0), fRandomPsi4(0), fRandomPsi5(0)
40 {
41 // constructor
42 }
43 
44 //_____________________________________________________________________________________________________
46  AliAnalysisTaskEmcalJet(name, kTRUE), fDistributionMultiplicity(0), fDistributionPt(0), fDistributionEtaPhi(0), fMinCentrality(0), fMaxCentrality(10), fDistributionV2(0), fDistributionV3(0), fDistributionV4(0), fDistributionV5(0), fInputArrayName(""), fOutputArrayName(""), fInputArray(0), fOutputArray(0), fRandom(), fToyCent(0), fRandomPsi3(0), fRandomPsi4(0), fRandomPsi5(0)
47 {
48 // constructor
49 }
50 
51 //_____________________________________________________________________________________________________
53 {
54  // destructor
57  if(fDistributionPt)
58  delete fDistributionPt;
60  delete fDistributionEtaPhi;
61 }
62 
63 
64 //_____________________________________________________________________________________________________
66 {
68 
69  fRandom = new TRandom3(0);
70  AddHistogram1D<TH1D>("hTrackPt", "Tracks p_{T} distribution", "", 300, 0., 300., "p_{T} (GeV/c)", "dN^{Tracks}/dp_{T}");
71  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");
72  AddHistogram1D<TH1D>("hMultiplicity", "Number of tracks in acceptance vs. centrality", "", 500, 0., 5000., "N tracks","dN^{Events}/dN^{Tracks}");
73 
74  PostData(1, fOutput);
75 }
76 
77 
78 //_____________________________________________________________________________________________________
80 {
82 
83  // Check if input array can be loaded
84  if(!fInputArrayName.IsNull())
85  {
86  fInputArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fInputArrayName.Data())));
87  if(!fInputArray)
88  {
89  AliFatal(Form("Input array '%s' not found!", fInputArrayName.Data()));
90  }
91  }
92 
93  // Check if output arrays can be created
94  if((InputEvent()->FindListObject(Form("%s", fOutputArrayName.Data()))))
95  AliFatal(Form("Output array '%s' already exists in the event! Rename it.", fOutputArrayName.Data()));
96 
97  fOutputArray = new TClonesArray("AliAODTrack");
98  fOutputArray->SetName(fOutputArrayName.Data());
99  fInputEvent->AddObject(fOutputArray);
100 }
101 
102 //_____________________________________________________________________________________________________
104 {
105  AssembleEvent();
106  CreateQAPlots();
107  return kTRUE;
108 }
109 
110 //________________________________________________________________________
112 {
113  fRandomPsi3 = fRandom->Rndm()*TMath::Pi(); // once per event, create a random value dedicated for Psi3
114  fRandomPsi4 = fRandom->Rndm()*TMath::Pi(); // once per event, create a random value dedicated for Psi4
115  fRandomPsi5 = fRandom->Rndm()*TMath::Pi(); // once per event, create a random value dedicated for Psi5
116  fToyCent = fMinCentrality + (fMaxCentrality-fMinCentrality)*fRandom->Rndm(); // centrality value (flat from selected range)
117 
118  // Create the event from the several inputs and run the jet finder
119 
120  // ################# 1. Add input tracks (if available)
121  Int_t particleCount = 0;
122  if(fInputArray)
123  {
124  for(Int_t iPart=0; iPart<fInputArray->GetEntries(); iPart++)
125  {
126  AliAODTrack* inputParticle = static_cast<AliAODTrack*>(fInputArray->At(iPart));
127  new ((*fOutputArray)[particleCount]) AliAODTrack(*inputParticle);
128  particleCount++;
129  }
130  }
131 
132  // ################# 2. Create a vertex if there is none (needed by some tasks)
133  if(dynamic_cast<AliESDEvent*>(InputEvent()))
134  {
135  if(!(dynamic_cast<AliESDEvent*>(InputEvent()))->GetPrimaryVertexTracks()->GetNContributors())
136  static_cast<AliESDEvent*>(fInputEvent)->SetPrimaryVertexTracks(new AliESDVertex(0.,0., 100));
137  }
138  else if(dynamic_cast<AliAODEvent*>(InputEvent()))
139  {
140  if( (!(dynamic_cast<AliAODEvent*>(InputEvent()))->GetPrimaryVertex()) || (!(dynamic_cast<AliAODEvent*>(InputEvent()))->GetPrimaryVertex()->GetNContributors()) )
141  {
142  Double_t* p = new Double_t[3];
143  p[0] = 0.; p[1] = 0.; p[2] = 0.; // for backwards compatibility
144  AliAODVertex* vertex = new AliAODVertex(p,1.);
145  vertex->SetNContributors(100);
146  vertex->SetName("PrimaryVertex");
147  static_cast<AliAODEvent*>(fInputEvent)->AddVertex(vertex);
148  }
149  }
150 
151  // ################# 3. Create toy event
152  Int_t multiplicity = (Int_t)fDistributionMultiplicity->GetRandom();
153  for(Int_t i=0;i<multiplicity; i++)
154  {
155  Double_t trackPt = fDistributionPt->GetRandom();
156  Double_t trackEta = 0;
157  Double_t trackPhi = 0;
158  static_cast<TH2*>(fDistributionEtaPhi)->GetRandom2(trackPhi, trackEta);
159  Double_t trackTheta = 2.*atan(exp(-trackEta));
160  Double_t trackCharge = fRandom->Rndm() - 0.5;
161 
162  if(trackCharge>0) trackCharge = 1; else trackCharge = -1;
163 
164  // Add flow to particle
166  trackPhi = AddFlow(trackPhi, trackPt);
167 
168 
169  // Add basic particle to event
170  new ((*fOutputArray)[particleCount]) AliAODTrack();
171  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetPt(trackPt);
172  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetPhi(trackPhi);
173  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetTheta(trackTheta); // AliAODTrack cannot set eta directly
174  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetCharge(trackCharge);
175  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetLabel(100000 + i);
176  static_cast<AliAODTrack*>(fOutputArray->At(particleCount))->SetIsHybridGlobalConstrainedGlobal();
177  particleCount++;
178  }
179 
180 }
181 
182 //_____________________________________________________________________________________________________
184 {
185  for(Int_t iTrack=0; iTrack<fOutputArray->GetEntries(); iTrack++)
186  {
187  AliAODTrack* track = static_cast<AliAODTrack*>(fOutputArray->At(iTrack));
188  FillHistogram("hTrackPt", track->Pt());
189  FillHistogram("hTrackPhiEta", track->Phi(), track->Eta());
190  }
191  FillHistogram("hMultiplicity", fOutputArray->GetEntries());
192 }
193 
194 
195 //_____________________________________________________________________________________________________
197 {
198  // adapted from AliFlowTrackSimple
199  Double_t precisionPhi = 1e-10;
200  Int_t maxNumberOfIterations = 200;
201 
202  Double_t phi0=phi;
203  Double_t f=0.;
204  Double_t fp=0.;
205  Double_t phiprev=0.;
206  Int_t ptBin = 0;
207 
208  // Evaluate V2 for track pt/centrality
209  Double_t v2 = 0;
210  if(fDistributionV2)
211  {
212  ptBin = fDistributionV2->GetXaxis()->FindBin(pt);
213  if(ptBin>fDistributionV2->GetNbinsX())
214  v2 = fDistributionV2->GetBinContent(fDistributionV2->GetNbinsX(), fDistributionV2->GetYaxis()->FindBin(fToyCent));
215  else if(ptBin>0)
216  v2 = fDistributionV2->GetBinContent(ptBin, fDistributionV2->GetYaxis()->FindBin(fToyCent));
217  }
218 
219  // Evaluate V3 for track pt/centrality
220  Double_t v3 = 0;
221  if(fDistributionV3)
222  {
223  ptBin = fDistributionV3->GetXaxis()->FindBin(pt);
224  if(ptBin>fDistributionV3->GetNbinsX())
225  v3 = fDistributionV3->GetBinContent(fDistributionV3->GetNbinsX(), fDistributionV3->GetYaxis()->FindBin(fToyCent));
226  else if(ptBin>0)
227  v3 = fDistributionV3->GetBinContent(ptBin, fDistributionV3->GetYaxis()->FindBin(fToyCent));
228  }
229 
230  // Evaluate V4 for track pt/centrality
231  Double_t v4 = 0;
232  if(fDistributionV4)
233  {
234  ptBin = fDistributionV4->GetXaxis()->FindBin(pt);
235  if(ptBin>fDistributionV4->GetNbinsX())
236  v4 = fDistributionV4->GetBinContent(fDistributionV4->GetNbinsX(), fDistributionV4->GetYaxis()->FindBin(fToyCent));
237  else if(ptBin>0)
238  v4 = fDistributionV4->GetBinContent(ptBin, fDistributionV4->GetYaxis()->FindBin(fToyCent));
239  }
240 
241  // Evaluate V5 for track pt/centrality
242  Double_t v5 = 0;
243  if(fDistributionV5)
244  {
245  ptBin = fDistributionV5->GetXaxis()->FindBin(pt);
246  if(ptBin>fDistributionV5->GetNbinsX())
247  v5 = fDistributionV5->GetBinContent(fDistributionV5->GetNbinsX(), fDistributionV5->GetYaxis()->FindBin(fToyCent));
248  else if(ptBin>0)
249  v5 = fDistributionV5->GetBinContent(ptBin, fDistributionV5->GetYaxis()->FindBin(fToyCent));
250  }
251 
252  // Add all v's
253  for (Int_t i=0; i<maxNumberOfIterations; i++)
254  {
255  phiprev=phi; //store last value for comparison
256  f = phi-phi0
257  + v2*TMath::Sin(2.*(phi-(fEPV0+(TMath::Pi()/2.))))
258  +2./3.*v3*TMath::Sin(3.*(phi-fRandomPsi3))
259  +0.5 *v4*TMath::Sin(4.*(phi-fRandomPsi4))
260  +0.4 *v5*TMath::Sin(5.*(phi-fRandomPsi5));
261 
262  fp = 1.0+2.0*(
263  +v2*TMath::Cos(2.*(phi-(fEPV0+(TMath::Pi()/2.))))
264  +v3*TMath::Cos(3.*(phi-fRandomPsi3))
265  +v4*TMath::Cos(4.*(phi-fRandomPsi4))
266  +v5*TMath::Cos(5.*(phi-fRandomPsi5))); //first derivative
267 
268  phi -= f/fp;
269  if (TMath::AreEqualAbs(phiprev,phi,precisionPhi)) break;
270  }
271 
272  return phi;
273 }
274 
275 
276 //________________________________________________________________________
278 {
279  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
280  if(!tmpHist)
281  {
282  AliError(Form("Cannot find histogram <%s> ",key)) ;
283  return;
284  }
285 
286  tmpHist->Fill(x);
287 }
288 
289 //________________________________________________________________________
291 {
292  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
293  if(!tmpHist)
294  {
295  AliError(Form("Cannot find histogram <%s> ",key));
296  return;
297  }
298 
299  if (tmpHist->IsA()->GetBaseClass("TH1"))
300  static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
301  else if (tmpHist->IsA()->GetBaseClass("TH2"))
302  static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
303 }
304 
305 //________________________________________________________________________
307 {
308  TH2* tmpHist = static_cast<TH2*>(fOutput->FindObject(key));
309  if(!tmpHist)
310  {
311  AliError(Form("Cannot find histogram <%s> ",key));
312  return;
313  }
314 
315  tmpHist->Fill(x,y,add);
316 }
317 
318 //________________________________________________________________________
319 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)
320 {
321  T* tmpHist = new T(name, title, xBins, xMin, xMax);
322 
323  tmpHist->GetXaxis()->SetTitle(xTitle);
324  tmpHist->GetYaxis()->SetTitle(yTitle);
325  tmpHist->SetOption(options);
326  tmpHist->SetMarkerStyle(kFullCircle);
327  tmpHist->Sumw2();
328 
329  fOutput->Add(tmpHist);
330 
331  return tmpHist;
332 }
333 
334 //________________________________________________________________________
335 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)
336 {
337  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax);
338  tmpHist->GetXaxis()->SetTitle(xTitle);
339  tmpHist->GetYaxis()->SetTitle(yTitle);
340  tmpHist->GetZaxis()->SetTitle(zTitle);
341  tmpHist->SetOption(options);
342  tmpHist->SetMarkerStyle(kFullCircle);
343  tmpHist->Sumw2();
344 
345  fOutput->Add(tmpHist);
346 
347  return tmpHist;
348 }
349 
350 
351 //________________________________________________________________________
353 {
354  // Called once at the end of the analysis.
355 }
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
const char * title
Definition: MakeQAPdf.C:27
TClonesArray * fOutputArray
input array containing tracks from events
Double_t fEPV0
!event plane V0
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
Distribution for v5 in bins of pt and centrality.
TH2 * fDistributionV4
Distribution for v3 in bins of pt and centrality.
Double_t fRandomPsi5
eventwise calculated psi 4
AliEmcalList * fOutput
!output list
Definition: External.C:220
Double_t fRandomPsi3
eventwise centrality value
Base task in the EMCAL jet framework.
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.