Visualization class. To use
[u/mrichter/AliRoot.git] / TPC / Upgrade / AliToyMCEventGeneratorSimple.cxx
CommitLineData
526ddf0e 1#include <iostream>
a1a695e5 2
526ddf0e 3#include <TDatabasePDG.h>
4#include <TRandom.h>
5#include <TF1.h>
32438f4e 6#include <TStopwatch.h>
f356075c 7#include <AliESDtrackCuts.h>
8#include <AliESDtrack.h>
9#include <TFile.h>
10#include <TTree.h>
11
a1a695e5 12
de0014b7 13#include "AliToyMCEvent.h"
526ddf0e 14
a1a695e5 15#include "AliToyMCEventGeneratorSimple.h"
16
f356075c 17
de0014b7 18ClassImp(AliToyMCEventGeneratorSimple);
526ddf0e 19
20
de0014b7 21AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple()
22 :AliToyMCEventGenerator()
526ddf0e 23 ,fVertexMean(0.)
24 ,fVertexSigma(7.)
32438f4e 25 ,fNtracks(20)
f356075c 26 ,fInputFileNameESD("")
27 ,fESDCuts(0x0)
526ddf0e 28{
29 //
30
31}
32//________________________________________________________________
de0014b7 33AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple(const AliToyMCEventGeneratorSimple &gen)
34 :AliToyMCEventGenerator(gen)
526ddf0e 35 ,fVertexMean(gen.fVertexMean)
36 ,fVertexSigma(gen.fVertexSigma)
32438f4e 37 ,fNtracks(gen.fNtracks)
f356075c 38 ,fInputFileNameESD(gen.fInputFileNameESD)
39 ,fESDCuts(gen.fESDCuts)
526ddf0e 40{
41}
42//________________________________________________________________
de0014b7 43AliToyMCEventGeneratorSimple::~AliToyMCEventGeneratorSimple()
526ddf0e 44{
45}
46 //________________________________________________________________
de0014b7 47AliToyMCEventGeneratorSimple& AliToyMCEventGeneratorSimple::operator = (const AliToyMCEventGeneratorSimple &gen)
526ddf0e 48{
49 //assignment operator
50 if (&gen == this) return *this;
de0014b7 51 new (this) AliToyMCEventGeneratorSimple(gen);
526ddf0e 52
53 return *this;
54}
55//________________________________________________________________
f356075c 56void AliToyMCEventGeneratorSimple::SetParametersSimple(Double_t vertexMean, Double_t vertexSigma) {
526ddf0e 57 fVertexMean = vertexMean;
58 fVertexSigma = vertexSigma;
59}
60//________________________________________________________________
cd8ed0ac 61AliToyMCEvent* AliToyMCEventGeneratorSimple::Generate(Double_t time)
62{
63 //
64 //
65 //
66
de0014b7 67 AliToyMCEvent *retEvent = new AliToyMCEvent();
526ddf0e 68 retEvent->SetT0(time);
69 retEvent->SetX(0);
70 retEvent->SetX(0);
71 retEvent->SetZ(gRandom->Gaus(fVertexMean,fVertexSigma));
72
73 Double_t etaCuts=.9;
74 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
75 TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
76 fpt.SetParameters(7.24,0.120);
77 fpt.SetNpx(10000);
32438f4e 78// Int_t nTracks = 400; //TODO: draw from experim dist
79// Int_t nTracks = 20; //TODO: draw from experim dist
526ddf0e 80
32438f4e 81 for(Int_t iTrack=0; iTrack<fNtracks; iTrack++){
526ddf0e 82 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
83 Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
84 Double_t pt = fpt.GetRandom(); // momentum for f1
85 Short_t sign=1;
86 if(gRandom->Rndm() < 0.5){
87 sign =1;
88 }else{
89 sign=-1;
90 }
91
92 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
93 Double_t pxyz[3];
94 pxyz[0]=pt*TMath::Cos(phi);
95 pxyz[1]=pt*TMath::Sin(phi);
96 pxyz[2]=pt*TMath::Tan(theta);
97 Double_t vertex[3]={0,0,retEvent->GetZ()};
98 Double_t cv[21]={0};
de0014b7 99 AliToyMCTrack *tempTrack = new AliToyMCTrack(vertex,pxyz,cv,sign);
d1cf83f5 100 // use unique ID for track number
101 // this will be used in DistortTrack track to set the cluster label
102 tempTrack->SetUniqueID(iTrack);
526ddf0e 103
104 Bool_t trackDist = DistortTrack(*tempTrack, time);
105 if(trackDist) retEvent->AddTrack(*tempTrack);
106 delete tempTrack;
107 }
108
109 return retEvent;
110}
111
a1a695e5 112//________________________________________________________________
f356075c 113void AliToyMCEventGeneratorSimple::RunSimulationSimple(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
a1a695e5 114{
115 //
116 // run simple simulation with equal event spacing
117 //
118
119 if (!ConnectOutputFile()) return;
cd8ed0ac 120 //initialise the space charge. Should be done after the tree was set up
121 InitSpaceCharge();
122
32438f4e 123 fNtracks=ntracks;
a1a695e5 124 Double_t eventTime=0.;
125 const Double_t eventSpacing=1./50e3; //50kHz equally spaced
32438f4e 126 TStopwatch s;
a1a695e5 127 for (Int_t ievent=0; ievent<nevents; ++ievent){
32438f4e 128 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
a1a695e5 129 fEvent = Generate(eventTime);
130 FillTree();
131 delete fEvent;
132 fEvent=0x0;
133 eventTime+=eventSpacing;
134 }
32438f4e 135 s.Stop();
136 s.Print();
137
a1a695e5 138 CloseOutputFile();
139}
526ddf0e 140
f356075c 141
142
143//________________________________________________________________
144AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD(AliESDEvent &esdEvent, Double_t time) {
145
146
147 AliToyMCEvent *retEvent = new AliToyMCEvent();
148 retEvent->SetT0(time);
149 retEvent->SetX(esdEvent.GetPrimaryVertex()->GetX());
150 retEvent->SetY(esdEvent.GetPrimaryVertex()->GetY());
151 retEvent->SetZ(esdEvent.GetPrimaryVertex()->GetZ());
152
153
154
155 if(!fNtracks) fNtracks = esdEvent.GetNumberOfTracks();
156 Int_t nUsedTracks = 0;
157 for(Int_t iTrack=0; iTrack<esdEvent.GetNumberOfTracks(); iTrack++){
158
159 AliESDtrack *part = esdEvent.GetTrack(iTrack);
160 if(!part) continue;
161 if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue;
162
163 Double_t pxyz[3];
164 pxyz[0]=part->Px();
165 pxyz[1]=part->Py();
166 pxyz[2]=part->Pz();
167 Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
168
169 Double_t cv[21]={0};
170 Int_t sign = part->Charge();
171 AliToyMCTrack *tempTrack = new AliToyMCTrack(vertex,pxyz,cv,sign);
172
173 Bool_t trackDist = DistortTrack(*tempTrack, time);
174 if(trackDist) {
175 retEvent->AddTrack(*tempTrack);
176 nUsedTracks++;
177 }
178 delete tempTrack;
179
180 if(nUsedTracks >= fNtracks) break;
181 }
182
183 return retEvent;
184}
185//________________________________________________________________
186void AliToyMCEventGeneratorSimple::RunSimulationESD(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
187{
188 //
189 // run simulation using esd input with equal event spacing
190 //
191
192 if (!ConnectOutputFile()) return;
193 TFile f(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"));
194 if(!f.IsOpen()) {
195 std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
196 return;
197 }
198 TTree* esdTree = (TTree*) f.Get("esdTree");
199 if(!esdTree) {
200 std::cout <<"no esdTree in file " << std::endl;
201 return;
202 }
203 InitSpaceCharge();
204 AliESDEvent* esdEvent = new AliESDEvent();
205 fNtracks = ntracks;
206
207 esdEvent->ReadFromTree(esdTree);
208 Double_t eventTime=0.;
209 const Double_t eventSpacing=1./50e3; //50kHz equally spaced
210 TStopwatch s;
211 Int_t nEvents = nevents;
212 fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
213 if(nevents>esdTree->GetEntries()) nEvents = esdTree->GetEntries();
214 Int_t nUsedEvents = 0;
215 for (Int_t ievent=0; ievent<esdTree->GetEntries(); ++ievent){
216 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
217 esdTree->GetEvent(ievent);
218 if(esdEvent->GetNumberOfTracks()==0) {
219 std::cout << " tracks == 0" << std::endl;
220 continue;
221 }
222
223 fEvent = GenerateESD(*esdEvent, eventTime);
224 nUsedEvents++;
225 FillTree();
226 delete fEvent;
227 fEvent=0x0;
228 eventTime+=eventSpacing;
229 if(nUsedEvents>=nevents) break;
230 }
231 s.Stop();
232 s.Print();
233
234 CloseOutputFile();
235}
236