added method RunSimulationBunchTrain(constInt_t nevents,const Int_t ntracks) to simul...
[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>
609d6d39 11#include <TRandom3.h>
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);
526ddf0e 100
101 Bool_t trackDist = DistortTrack(*tempTrack, time);
102 if(trackDist) retEvent->AddTrack(*tempTrack);
103 delete tempTrack;
104 }
105
106 return retEvent;
107}
108
a1a695e5 109//________________________________________________________________
609d6d39 110void AliToyMCEventGeneratorSimple::RunSimulation(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
a1a695e5 111{
112 //
113 // run simple simulation with equal event spacing
114 //
115
116 if (!ConnectOutputFile()) return;
cd8ed0ac 117 //initialise the space charge. Should be done after the tree was set up
118 InitSpaceCharge();
119
32438f4e 120 fNtracks=ntracks;
a1a695e5 121 Double_t eventTime=0.;
122 const Double_t eventSpacing=1./50e3; //50kHz equally spaced
32438f4e 123 TStopwatch s;
a1a695e5 124 for (Int_t ievent=0; ievent<nevents; ++ievent){
32438f4e 125 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
a1a695e5 126 fEvent = Generate(eventTime);
127 FillTree();
128 delete fEvent;
129 fEvent=0x0;
130 eventTime+=eventSpacing;
131 }
32438f4e 132 s.Stop();
133 s.Print();
134
a1a695e5 135 CloseOutputFile();
136}
526ddf0e 137
f356075c 138//________________________________________________________________
139AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD(AliESDEvent &esdEvent, Double_t time) {
140
141
142 AliToyMCEvent *retEvent = new AliToyMCEvent();
143 retEvent->SetT0(time);
144 retEvent->SetX(esdEvent.GetPrimaryVertex()->GetX());
145 retEvent->SetY(esdEvent.GetPrimaryVertex()->GetY());
146 retEvent->SetZ(esdEvent.GetPrimaryVertex()->GetZ());
147
148
149
150 if(!fNtracks) fNtracks = esdEvent.GetNumberOfTracks();
151 Int_t nUsedTracks = 0;
152 for(Int_t iTrack=0; iTrack<esdEvent.GetNumberOfTracks(); iTrack++){
153
154 AliESDtrack *part = esdEvent.GetTrack(iTrack);
155 if(!part) continue;
156 if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue;
609d6d39 157 if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
f356075c 158 Double_t pxyz[3];
159 pxyz[0]=part->Px();
160 pxyz[1]=part->Py();
161 pxyz[2]=part->Pz();
162 Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
163
164 Double_t cv[21]={0};
165 Int_t sign = part->Charge();
166 AliToyMCTrack *tempTrack = new AliToyMCTrack(vertex,pxyz,cv,sign);
167
168 Bool_t trackDist = DistortTrack(*tempTrack, time);
169 if(trackDist) {
170 retEvent->AddTrack(*tempTrack);
609d6d39 171
f356075c 172 nUsedTracks++;
173 }
174 delete tempTrack;
175
176 if(nUsedTracks >= fNtracks) break;
177 }
178
179 return retEvent;
180}
181//________________________________________________________________
182void AliToyMCEventGeneratorSimple::RunSimulationESD(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
183{
184 //
185 // run simulation using esd input with equal event spacing
186 //
187
188 if (!ConnectOutputFile()) return;
189 TFile f(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"));
190 if(!f.IsOpen()) {
191 std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
192 return;
193 }
194 TTree* esdTree = (TTree*) f.Get("esdTree");
195 if(!esdTree) {
196 std::cout <<"no esdTree in file " << std::endl;
197 return;
198 }
199 InitSpaceCharge();
200 AliESDEvent* esdEvent = new AliESDEvent();
201 fNtracks = ntracks;
202
203 esdEvent->ReadFromTree(esdTree);
204 Double_t eventTime=0.;
205 const Double_t eventSpacing=1./50e3; //50kHz equally spaced
206 TStopwatch s;
207 Int_t nEvents = nevents;
208 fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
209 if(nevents>esdTree->GetEntries()) nEvents = esdTree->GetEntries();
210 Int_t nUsedEvents = 0;
211 for (Int_t ievent=0; ievent<esdTree->GetEntries(); ++ievent){
609d6d39 212 printf("Generating event %3d (%.3g)\n",nUsedEvents,eventTime);
f356075c 213 esdTree->GetEvent(ievent);
214 if(esdEvent->GetNumberOfTracks()==0) {
215 std::cout << " tracks == 0" << std::endl;
216 continue;
217 }
218
219 fEvent = GenerateESD(*esdEvent, eventTime);
609d6d39 220 if(fEvent->GetNumberOfTracks() >=10) {
221 nUsedEvents++;
222 FillTree();
223 eventTime+=eventSpacing;
224 }
f356075c 225 delete fEvent;
226 fEvent=0x0;
609d6d39 227
f356075c 228 if(nUsedEvents>=nevents) break;
229 }
230 s.Stop();
231 s.Print();
232
233 CloseOutputFile();
234}
235
609d6d39 236//________________________________________________________________
237
238void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
239{
240 //
241 // run simple simulation with equal event spacing
242 //
243
244 //Parameters for bunc
245 const Double_t abortGap = 3e-6; //
246 const Double_t collFreq = 50e3;
247 const Double_t bSpacing = 50e-9; //bunch spacing
248 const Int_t nTrainBunches = 48;
249 const Int_t nTrains = 12;
250 const Double_t revFreq = 1.11e4; //revolution frequency
251 const Double_t collProb = collFreq/(nTrainBunches*nTrains*revFreq);
252 const Double_t trainLength = bSpacing*(nTrainBunches-1);
253 const Double_t totTrainLength = nTrains*trainLength;
254 const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1);
255 Bool_t equalSpacing = kFALSE;
256
257 TRandom3 *rand = new TRandom3();
258 //rand->SetSeed();
259
260 if (!ConnectOutputFile()) return;
261 //initialise the space charge. Should be done after the tree was set up
262 InitSpaceCharge();
263
264 fNtracks=ntracks;
265 Double_t eventTime=0.;
266 const Double_t eventSpacing=1./50e3; //50kHz equally spaced
267 TStopwatch s;
268 Int_t nGeneratedEvents = 0;
269 Int_t bunchCounter = 0;
270 Int_t trainCounter = 0;
271
272 //for (Int_t ievent=0; ievent<nevents; ++ievent){
273 while (nGeneratedEvents<nevents){
274 // std::cout <<trainCounter << " " << bunchCounter << " "<< "eventTime " << eventTime << std::endl;
275
276 if(equalSpacing) {
277 printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
278 fEvent = Generate(eventTime);
279 nGeneratedEvents++;
280 FillTree();
281 delete fEvent;
282 fEvent=0x0;
283 eventTime+=1./collFreq;
284 }
285 else{
286 Int_t nCollsInCrossing = rand -> Poisson(collProb);
287 for(Int_t iColl = 0; iColl<nCollsInCrossing; iColl++){
288 printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
289 fEvent = Generate(eventTime);
290 nGeneratedEvents++;
291 FillTree();
292 delete fEvent;
293 fEvent=0x0;
294
295 }
296 bunchCounter++;
297
298 if(bunchCounter>=nTrainBunches){
299
300 trainCounter++;
301 if(trainCounter>=nTrains){
302 eventTime+=abortGap;
303 trainCounter=0;
304 }
305 else eventTime+=trainSpacing;
306
307 bunchCounter=0;
308 }
309 else eventTime+= bSpacing;
310
311 }
312
313
314 }
315 s.Stop();
316 s.Print();
317 delete rand;
318 CloseOutputFile();
319}