AliPhysics  781d0c7 (781d0c7)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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>
29 
31 
32 //_____________________________________________________________________________________________________
34  AliAnalysisTaskSE("AliAnalysisTaskChargedJetsHadronToy"), fCreateUE(1), fCreateJets(1), fUEMultDistribution(0), fUEDistribution(0), fUEMultiplicity(1000), fGeneratedJetParticleDistribution(0), fGeneratedJetPtDistribution(0), fGeneratedJetCount(1), fGeneratedJetPtMin(20.), fGeneratedJetPtMax(30.), fGeneratedJetWidthPhi(0.2), fGeneratedJetWidthEta(0.2), fGeneratedJetMinEta(-0.9), fGeneratedJetMaxEta(0.9), fInputArrTracks(0), fInputArrTracksName(""), fOutputArrTracks(0), fOutputArrTracksName(""), fGeneratedJetsArr(0), fGeneratedJetsArrName(""), fDistEtaGaussian(0), fDistPhiGaussian(0), fRandom(), fInitialized()
35 {
36 // constructor
37 }
38 
39 //_____________________________________________________________________________________________________
41 {
42 // destructor
44  delete fUEMultDistribution;
45  if(fUEDistribution)
46  delete fUEDistribution;
52  delete fDistPhiGaussian;
54  delete fDistEtaGaussian;
55 }
56 
57 
58 //_____________________________________________________________________________________________________
60 {
61  fRandom = new TRandom3(0);
62 }
63 
64 //_____________________________________________________________________________________________________
66 {
67  // Check if input array can be loaded
68  if(!fInputArrTracksName.IsNull())
69  {
70  fInputArrTracks = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fInputArrTracksName.Data())));
71  if(!fInputArrTracks)
72  AliFatal(Form("Input array '%s' demanded, but not found!", fInputArrTracksName.Data()));
73  if(strcmp(fInputArrTracks->GetClass()->GetName(), "AliAODTrack"))
74  AliFatal(Form("Input array has track type %s. Only AliAODTracks are supported.", fInputArrTracks->GetClass()->GetName()));
75  }
76 
77  // Check if output arrays can be created
78  if((InputEvent()->FindListObject(Form("%s", fOutputArrTracksName.Data()))))
79  AliFatal(Form("Output array '%s' already exists in the event! Rename it.", fOutputArrTracksName.Data()));
80  if((InputEvent()->FindListObject(Form("%s", fGeneratedJetsArrName.Data()))))
81  AliFatal(Form("Output array '%s' already exists in the event! Rename it.", fGeneratedJetsArrName.Data()));
82 
83  // Define functions used in toy model + put the arrays to the event
84  // NOTE: fOutputArrTracks must be added in any case since it contains the tracks of UE and jets
85  fOutputArrTracks = new TClonesArray("AliAODTrack");
86  fOutputArrTracks->SetName(fOutputArrTracksName.Data());
87  InputEvent()->AddObject(fOutputArrTracks);
88 
89  if(fCreateUE)
90  {
91  if(!fUEDistribution)
92  {
93  std::cout << "\n### Distribution for the UE not given -- using default thermal distribution ###\n\n";
94  fUEDistribution = new TF1("fUEDistribution","[0]*exp([1]*x)",0,200.);
95  fUEDistribution->SetParameters(1.0,-1.5);
96  }
97  fUEDistribution->SetNpx(400);
98  }
99 
100  if(fCreateJets)
101  {
102  fGeneratedJetsArr = new TClonesArray("AliEmcalJet");
103  fGeneratedJetsArr->SetName(fGeneratedJetsArrName.Data());
104  InputEvent()->AddObject(fGeneratedJetsArr);
106  {
107  std::cout << "\n### Distribution for the particles in jets not given -- using default power law distribution ###\n\n";
108  fGeneratedJetParticleDistribution = new TF1("fGeneratedJetParticleDistribution","[0]*(x^[1])",4.,120.);
109  fGeneratedJetParticleDistribution->SetParameters(6.0,-2.5);
110  }
112 
114  {
115  std::cout << "\n### Distribution for the jet pt not given -- only min jet pT ###\n\n";
116  }
117  else
118  {
119  Int_t minBin = fGeneratedJetPtDistribution->GetXaxis()->FindBin(fGeneratedJetPtMin);
120  Int_t maxBin = fGeneratedJetPtDistribution->GetXaxis()->FindBin(fGeneratedJetPtMax);
121  for(Int_t i=0; i<=fGeneratedJetPtDistribution->GetNbinsX(); i++)
122  if(i < minBin || i > maxBin)
123  {
124  fGeneratedJetPtDistribution->SetBinContent(i,0.0);
125  fGeneratedJetPtDistribution->SetBinError(i,0.0);
126  }
127 
128  }
129 
130  fDistEtaGaussian = new TF1("fDistEtaGaussian","gaus(0)",-1.,1.);// gaus(0) is [0]*exp(-0.5*((x-[1])/[2])**2)
131  fDistEtaGaussian->SetParameters(1.0,0.,0.5*fGeneratedJetWidthEta);
132  fDistPhiGaussian = new TF1("fDistPhiGaussian","gaus(0)",-1.,1.);
133  fDistPhiGaussian->SetParameters(1.0,0.,0.5*fGeneratedJetWidthPhi);
134  }
135 
136  fInitialized = kTRUE;
137 }
138 
139 //_____________________________________________________________________________________________________
141 {
142  // Run once the exec once (must be here to have the input event ready)
143  if(!fInitialized)
144  ExecOnce();
145 
146  AssembleEvent();
147 // std::cout << fOutputArrTracks->GetName() << ":" << fOutputArrTracks->GetEntries() << std::endl;
148 
149 }
150 
151 //________________________________________________________________________
153 {
154  // Create the event from the several inputs and run the jet finder
155  // Note: those tracks are guaranteed to be AliAODTracks
156  if(fInputArrTracks)
158 
159  // 1. Create a vertex if there is none (needed by correlation task)
160  if(dynamic_cast<AliESDEvent*>(InputEvent()))
161  {
162  if(!(dynamic_cast<AliESDEvent*>(InputEvent()))->GetPrimaryVertexTracks()->GetNContributors())
163  static_cast<AliESDEvent*>(fInputEvent)->SetPrimaryVertexTracks(new AliESDVertex(0.,0., 100));
164  }
165  else if(dynamic_cast<AliAODEvent*>(InputEvent()))
166  {
167  if( (!(dynamic_cast<AliAODEvent*>(InputEvent()))->GetPrimaryVertex()) || (!(dynamic_cast<AliAODEvent*>(InputEvent()))->GetPrimaryVertex()->GetNContributors()) )
168  {
169  Double_t* p = new Double_t[3];
170  p[0] = 0.; p[1] = 0.; p[2] = 0.; // for backwards compatibility
171  AliAODVertex* vertex = new AliAODVertex(p,1.);
172  vertex->SetNContributors(100);
173  vertex->SetName("PrimaryVertex");
174  static_cast<AliAODEvent*>(fInputEvent)->AddVertex(vertex);
175  }
176  }
177 
178  // 2. Create underlying event
179  Int_t UEmultiplicity = fUEMultiplicity;
180  Double_t UEthrownPt = 0.;
181  Double_t etaMin = -0.9;
182  Double_t etaMax = +0.9;
183 
184  if(fCreateUE)
185  {
186  Int_t count = fOutputArrTracks->GetEntries();
187 
189  UEmultiplicity = (Int_t)fUEMultDistribution->GetRandom();
190 
191  for(Int_t i=0;i<UEmultiplicity; i++)
192  {
193  Double_t trackPt = fUEDistribution->GetRandom();
194  Double_t trackEta = etaMin + fRandom->Rndm()*(etaMax-etaMin);
195  Double_t trackTheta = 2.*atan(exp(-trackEta));
196  Double_t trackPhi = fRandom->Rndm()*TMath::TwoPi();
197  Double_t trackCharge = fRandom->Rndm() - 0.5;
198 
199  if(trackCharge>0) trackCharge = 1; else trackCharge = -1;
200 
201  // Add very basic particle to event
202  new ((*fOutputArrTracks)[count]) AliAODTrack();
203  static_cast<AliAODTrack*>(fOutputArrTracks->At(count))->SetPt(trackPt);
204  static_cast<AliAODTrack*>(fOutputArrTracks->At(count))->SetPhi(trackPhi);
205  static_cast<AliAODTrack*>(fOutputArrTracks->At(count))->SetTheta(trackTheta); // AliAODTrack cannot set eta directly
206  static_cast<AliAODTrack*>(fOutputArrTracks->At(count))->SetCharge(trackCharge);
207  count++;
208  UEthrownPt += trackPt;
209  }
210  }
211 
212  // 3. Embed gaussian jets into event
213  Double_t JETthrownPt = 0.;
214  if(fCreateJets)
215  {
216  // Define jets and throw them into the acceptance
217 
218  for(Int_t i=0;i<fGeneratedJetCount; i++)
219  {
221  Double_t jetPhi = fRandom->Rndm()*TMath::TwoPi();
222 
225  jetPt = fGeneratedJetPtDistribution->GetRandom();
226 
227  Int_t count = fOutputArrTracks->GetEntries();
228  Int_t particlesInJet = 0;
229  JETthrownPt = 0;
230  while(JETthrownPt < jetPt)
231  {
232  Double_t trackPt = fGeneratedJetParticleDistribution->GetRandom();
233  Double_t trackEta = jetEta + fDistEtaGaussian->GetRandom();
234  Double_t trackTheta = 2.*atan(exp(-trackEta));
235  Double_t trackPhi = jetPhi + fDistPhiGaussian->GetRandom();
236  Double_t trackCharge = fRandom->Rndm() - 0.5;
237  if(trackCharge>0) trackCharge = 1; else trackCharge = -1;
238 
239  // Add particle to event
240  new ((*fOutputArrTracks)[count]) AliAODTrack();
241  static_cast<AliAODTrack*>(fOutputArrTracks->At(count))->SetPt(trackPt);
242  static_cast<AliAODTrack*>(fOutputArrTracks->At(count))->SetPhi(trackPhi);
243  static_cast<AliAODTrack*>(fOutputArrTracks->At(count))->SetTheta(trackTheta); // AliAODTrack cannot set eta directly
244  static_cast<AliAODTrack*>(fOutputArrTracks->At(count))->SetCharge(trackCharge);
245 
246  count++;
247  particlesInJet++;
248  JETthrownPt += trackPt;
249  }
250 
251  // Save the generated jet for later matching
252  new ((*fGeneratedJetsArr)[i]) AliEmcalJet(JETthrownPt, jetEta, jetPhi, 0);
253  }
254  }
255 
256  std::cout << Form("Event has been generated using %3i particles UE (pT density: %5.2f). ", UEmultiplicity, UEthrownPt/(TMath::TwoPi()*(etaMax-etaMin))) << std::endl;
257  std::cout << Form(" %i embedded gaussian jets (pT density: %5.2f). ", fGeneratedJetCount, JETthrownPt/(TMath::TwoPi()*(etaMax-etaMin))) << std::endl;
258 
259 }
260 
double Double_t
Definition: External.C:58
ClassImp(AliAnalysisTaskChargedJetsHadronToy) AliAnalysisTaskChargedJetsHadronToy
TF1 * fDistPhiGaussian
function for gaussian distribution in toy
TRandom3 * fRandom
function for gaussian distribution in toy
int Int_t
Definition: External.C:63
TString fInputArrTracksName
input array containing tracks from events
TString fOutputArrTracksName
array holding tracks from toy model
TString fGeneratedJetsArrName
array holding generated jets from toy model
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
const char Option_t
Definition: External.C:48