AliPhysics  e59a9ba (e59a9ba)
 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 
223  Double_t jetPt = fGeneratedJetPtMin;
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 
ClassImp(AliAnalysisTaskChargedJetsHadronToy) AliAnalysisTaskChargedJetsHadronToy
TF1 * fDistPhiGaussian
function for gaussian distribution in toy
TRandom3 * fRandom
function for gaussian distribution in toy
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