o add STAT in include path
[u/mrichter/AliRoot.git] / TPC / Upgrade / AliToyMCEventGeneratorSimple.cxx
CommitLineData
526ddf0e 1#include <iostream>
a1a695e5 2
61d165ed 3#include <TROOT.h>
526ddf0e 4#include <TDatabasePDG.h>
5#include <TRandom.h>
6#include <TF1.h>
32438f4e 7#include <TStopwatch.h>
f356075c 8#include <TFile.h>
9#include <TTree.h>
609d6d39 10#include <TRandom3.h>
a1a695e5 11
0403120d 12#include <AliESDtrackCuts.h>
13#include <AliESDtrack.h>
14#include <AliTPCLaserTrack.h>
15
de0014b7 16#include "AliToyMCEvent.h"
526ddf0e 17
a1a695e5 18#include "AliToyMCEventGeneratorSimple.h"
19
cdf56507 20/*
21
22
23*/
f356075c 24
de0014b7 25ClassImp(AliToyMCEventGeneratorSimple);
526ddf0e 26
27
de0014b7 28AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple()
29 :AliToyMCEventGenerator()
526ddf0e 30 ,fVertexMean(0.)
31 ,fVertexSigma(7.)
32438f4e 32 ,fNtracks(20)
f356075c 33 ,fInputFileNameESD("")
34 ,fESDCuts(0x0)
29b7f480 35 ,fInputIndex(0)
36 ,fESDEvent(0x0)
37 ,fESDTree(0x0)
38 ,fInputFile(0x0)
1e62e876 39 ,fHPt(0x0)
40 ,fHEta(0x0)
41 ,fHMult(0x0)
42 ,fHistosSet(kFALSE)
43 ,fParamFile(0x0)
526ddf0e 44{
45 //
46
47}
48//________________________________________________________________
de0014b7 49AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple(const AliToyMCEventGeneratorSimple &gen)
50 :AliToyMCEventGenerator(gen)
526ddf0e 51 ,fVertexMean(gen.fVertexMean)
52 ,fVertexSigma(gen.fVertexSigma)
32438f4e 53 ,fNtracks(gen.fNtracks)
f356075c 54 ,fInputFileNameESD(gen.fInputFileNameESD)
55 ,fESDCuts(gen.fESDCuts)
29b7f480 56 ,fInputIndex(gen.fInputIndex)
57 ,fESDEvent(gen.fESDEvent)
58 ,fESDTree(gen.fESDTree)
59 ,fInputFile(gen.fInputFile)
1e62e876 60 ,fHPt(gen.fHPt)
61 ,fHEta(gen.fHEta)
62 ,fHMult(gen.fHMult)
63 ,fHistosSet(gen.fHistosSet)
64 ,fParamFile(0x0)
526ddf0e 65{
66}
67//________________________________________________________________
de0014b7 68AliToyMCEventGeneratorSimple::~AliToyMCEventGeneratorSimple()
526ddf0e 69{
70}
71 //________________________________________________________________
de0014b7 72AliToyMCEventGeneratorSimple& AliToyMCEventGeneratorSimple::operator = (const AliToyMCEventGeneratorSimple &gen)
526ddf0e 73{
74 //assignment operator
75 if (&gen == this) return *this;
de0014b7 76 new (this) AliToyMCEventGeneratorSimple(gen);
526ddf0e 77
78 return *this;
79}
80//________________________________________________________________
1f195c60 81void AliToyMCEventGeneratorSimple::SetParametersToyGen(const Char_t* parfilename/*="files/params.root*/, Double_t vertexMean/*=0*/, Double_t vertexSigma/*=7.*/) {
526ddf0e 82 fVertexMean = vertexMean;
83 fVertexSigma = vertexSigma;
61d165ed 84// fParamFile = new TFile(parfilename, "read");
85 TFile f(parfilename);
86 gROOT->cd();
87 fHPt = (TH1F*) f.Get("hPt");
88 fHEta = (TH1F*) f.Get("hEta");
89 fHMult = (TH1I*) f.Get("hMult") ;
1e62e876 90 fHistosSet = kTRUE;
61d165ed 91 f.Close();
1e62e876 92
93
526ddf0e 94}
95//________________________________________________________________
cd8ed0ac 96AliToyMCEvent* AliToyMCEventGeneratorSimple::Generate(Double_t time)
97{
98 //
99 //
100 //
101
de0014b7 102 AliToyMCEvent *retEvent = new AliToyMCEvent();
526ddf0e 103 retEvent->SetT0(time);
104 retEvent->SetX(0);
105 retEvent->SetX(0);
e83fd282 106 Double_t vertexZ=0;
107 //limit vertex to +- 10cm
108 while ( TMath::Abs(vertexZ=gRandom->Gaus(fVertexMean,fVertexSigma))>10. );
109 retEvent->SetZ(vertexZ);
526ddf0e 110
e83fd282 111 Double_t etaCuts=0.9;
526ddf0e 112 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
e83fd282 113 static TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.3,100);
0403120d 114 if (fpt.GetParameter(0)<1){
115 printf("Set Parameters\n");
116 fpt.SetParameters(7.24,0.120);
117 fpt.SetNpx(200);
118 }
32438f4e 119// Int_t nTracks = 400; //TODO: draw from experim dist
120// Int_t nTracks = 20; //TODO: draw from experim dist
1e62e876 121 Int_t nTracksLocal = fNtracks;
122 if(fHistosSet) {
123 nTracksLocal = fHMult->GetRandom();
124 if(nTracksLocal > fNtracks) nTracksLocal = fNtracks;
125 }
126 std::cout << "Generating " << nTracksLocal << " tracks " << std::endl;
223d9e38 127
1e62e876 128 for(Int_t iTrack=0; iTrack<nTracksLocal; iTrack++){
526ddf0e 129 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
1e62e876 130 Double_t eta = 0.;
131 Double_t pt = 0.;
132
133 if(fHistosSet) {//draw from histograms if set
134 eta = fHEta->GetRandom();
135 pt = fHPt->GetRandom();
136 }
137 else {
138 eta = gRandom->Uniform(-etaCuts, etaCuts);
139 pt = fpt.GetRandom(); // momentum for f1
140 }
526ddf0e 141 Short_t sign=1;
142 if(gRandom->Rndm() < 0.5){
143 sign =1;
144 }else{
145 sign=-1;
146 }
147
148 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
149 Double_t pxyz[3];
150 pxyz[0]=pt*TMath::Cos(phi);
151 pxyz[1]=pt*TMath::Sin(phi);
152 pxyz[2]=pt*TMath::Tan(theta);
153 Double_t vertex[3]={0,0,retEvent->GetZ()};
154 Double_t cv[21]={0};
0403120d 155 AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
7d26e3ce 156 // use unique ID for track number
157 // this will be used in DistortTrack track to set the cluster label
223d9e38 158 // in one simulation the track id should be unique for performance studies
159 tempTrack->SetUniqueID(fCurrentTrack++);
0403120d 160 DistortTrack(*tempTrack, time);
526ddf0e 161 }
162
163 return retEvent;
164}
165
a1a695e5 166//________________________________________________________________
fa97f7bd 167void AliToyMCEventGeneratorSimple::RunSimulation(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/, const Int_t rate/*=50*/)
a1a695e5 168{
169 //
170 // run simple simulation with equal event spacing
fa97f7bd 171 // rate in kHz
a1a695e5 172 //
173
174 if (!ConnectOutputFile()) return;
cd8ed0ac 175 //initialise the space charge. Should be done after the tree was set up
176 InitSpaceCharge();
223d9e38 177
178 // number of tracks to simulate per interaction
32438f4e 179 fNtracks=ntracks;
223d9e38 180 // within one simulation the track count should be unique for effeciency studies
a06336b6 181 // don't use the track ID 0 since the label -0 is not different from 0
182 fCurrentTrack=1;
223d9e38 183
a1a695e5 184 Double_t eventTime=0.;
fa97f7bd 185 const Double_t eventSpacing=1./rate/1e3; //50kHz equally spaced
32438f4e 186 TStopwatch s;
a1a695e5 187 for (Int_t ievent=0; ievent<nevents; ++ievent){
32438f4e 188 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
a1a695e5 189 fEvent = Generate(eventTime);
190 FillTree();
191 delete fEvent;
192 fEvent=0x0;
193 eventTime+=eventSpacing;
194 }
32438f4e 195 s.Stop();
196 s.Print();
197
a1a695e5 198 CloseOutputFile();
199}
526ddf0e 200
f356075c 201//________________________________________________________________
0403120d 202void AliToyMCEventGeneratorSimple::RunSimulationLaser(const Int_t nevents/*=1*/)
203{
204 //
205 // run simple simulation with equal event spacing
206 //
207
208 if (!ConnectOutputFile()) return;
209 //initialise the space charge. Should be done after the tree was set up
210 InitSpaceCharge();
211
212 // within one simulation the track count should be unique for effeciency studies
a06336b6 213 fCurrentTrack=1;
0403120d 214
215 Double_t eventTime=0.;
216 const Double_t eventSpacing=1./10.; //laser is running at 10Hz equally spaced
217 TStopwatch s;
218 for (Int_t ievent=0; ievent<nevents; ++ievent){
219 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
220 fEvent = GenerateLaser(eventTime);
221 FillTree();
222 delete fEvent;
223 fEvent=0x0;
224 eventTime+=eventSpacing;
225 }
226 s.Stop();
227 s.Print();
228
229 CloseOutputFile();
230}
231
232//________________________________________________________________
f356075c 233AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD(AliESDEvent &esdEvent, Double_t time) {
234
235
236 AliToyMCEvent *retEvent = new AliToyMCEvent();
237 retEvent->SetT0(time);
238 retEvent->SetX(esdEvent.GetPrimaryVertex()->GetX());
239 retEvent->SetY(esdEvent.GetPrimaryVertex()->GetY());
240 retEvent->SetZ(esdEvent.GetPrimaryVertex()->GetZ());
241
242
243
244 if(!fNtracks) fNtracks = esdEvent.GetNumberOfTracks();
245 Int_t nUsedTracks = 0;
246 for(Int_t iTrack=0; iTrack<esdEvent.GetNumberOfTracks(); iTrack++){
247
248 AliESDtrack *part = esdEvent.GetTrack(iTrack);
249 if(!part) continue;
250 if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue;
609d6d39 251 if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
f356075c 252 Double_t pxyz[3];
253 pxyz[0]=part->Px();
254 pxyz[1]=part->Py();
255 pxyz[2]=part->Pz();
256 Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
257
258 Double_t cv[21]={0};
259 Int_t sign = part->Charge();
0403120d 260 AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
7d26e3ce 261 // use unique ID for track number
262 // this will be used in DistortTrack track to set the cluster label
223d9e38 263 // in one simulation the track id should be unique for performance studies
264 tempTrack->SetUniqueID(fCurrentTrack++);
0403120d 265 DistortTrack(*tempTrack, time);
266 nUsedTracks++;
f356075c 267
268 if(nUsedTracks >= fNtracks) break;
269 }
270
271 return retEvent;
272}
273//________________________________________________________________
274void AliToyMCEventGeneratorSimple::RunSimulationESD(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
275{
276 //
277 // run simulation using esd input with equal event spacing
278 //
279
280 if (!ConnectOutputFile()) return;
281 TFile f(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"));
282 if(!f.IsOpen()) {
283 std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
284 return;
285 }
286 TTree* esdTree = (TTree*) f.Get("esdTree");
287 if(!esdTree) {
288 std::cout <<"no esdTree in file " << std::endl;
289 return;
290 }
291 InitSpaceCharge();
292 AliESDEvent* esdEvent = new AliESDEvent();
29b7f480 293
f356075c 294
295 esdEvent->ReadFromTree(esdTree);
29b7f480 296
223d9e38 297 fNtracks = ntracks;
298 // within one simulation the track count should be unique for effeciency studies
a06336b6 299 fCurrentTrack=1;
223d9e38 300
f356075c 301 Double_t eventTime=0.;
302 const Double_t eventSpacing=1./50e3; //50kHz equally spaced
303 TStopwatch s;
29b7f480 304
f356075c 305 fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
1e62e876 306 gRandom->SetSeed();
f356075c 307 Int_t nUsedEvents = 0;
308 for (Int_t ievent=0; ievent<esdTree->GetEntries(); ++ievent){
609d6d39 309 printf("Generating event %3d (%.3g)\n",nUsedEvents,eventTime);
f356075c 310 esdTree->GetEvent(ievent);
311 if(esdEvent->GetNumberOfTracks()==0) {
312 std::cout << " tracks == 0" << std::endl;
313 continue;
314 }
315
316 fEvent = GenerateESD(*esdEvent, eventTime);
609d6d39 317 if(fEvent->GetNumberOfTracks() >=10) {
318 nUsedEvents++;
319 FillTree();
320 eventTime+=eventSpacing;
321 }
f356075c 322 delete fEvent;
323 fEvent=0x0;
609d6d39 324
f356075c 325 if(nUsedEvents>=nevents) break;
326 }
327 s.Stop();
328 s.Print();
329
330 CloseOutputFile();
331}
332
609d6d39 333//________________________________________________________________
609d6d39 334void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
335{
336 //
337 // run simple simulation with equal event spacing
338 //
339
340 //Parameters for bunc
341 const Double_t abortGap = 3e-6; //
342 const Double_t collFreq = 50e3;
343 const Double_t bSpacing = 50e-9; //bunch spacing
344 const Int_t nTrainBunches = 48;
345 const Int_t nTrains = 12;
346 const Double_t revFreq = 1.11e4; //revolution frequency
347 const Double_t collProb = collFreq/(nTrainBunches*nTrains*revFreq);
348 const Double_t trainLength = bSpacing*(nTrainBunches-1);
349 const Double_t totTrainLength = nTrains*trainLength;
350 const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1);
351 Bool_t equalSpacing = kFALSE;
352
353 TRandom3 *rand = new TRandom3();
354 //rand->SetSeed();
355
356 if (!ConnectOutputFile()) return;
357 //initialise the space charge. Should be done after the tree was set up
358 InitSpaceCharge();
359
360 fNtracks=ntracks;
223d9e38 361 // within one simulation the track count should be unique for effeciency studies
a06336b6 362 fCurrentTrack=1;
223d9e38 363
609d6d39 364 Double_t eventTime=0.;
29b7f480 365 // const Double_t eventSpacing=1./50e3; //50kHz equally spaced
609d6d39 366 TStopwatch s;
367 Int_t nGeneratedEvents = 0;
368 Int_t bunchCounter = 0;
369 Int_t trainCounter = 0;
370
371 //for (Int_t ievent=0; ievent<nevents; ++ievent){
372 while (nGeneratedEvents<nevents){
373 // std::cout <<trainCounter << " " << bunchCounter << " "<< "eventTime " << eventTime << std::endl;
374
375 if(equalSpacing) {
376 printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
377 fEvent = Generate(eventTime);
378 nGeneratedEvents++;
379 FillTree();
380 delete fEvent;
381 fEvent=0x0;
382 eventTime+=1./collFreq;
383 }
384 else{
385 Int_t nCollsInCrossing = rand -> Poisson(collProb);
386 for(Int_t iColl = 0; iColl<nCollsInCrossing; iColl++){
387 printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
388 fEvent = Generate(eventTime);
389 nGeneratedEvents++;
390 FillTree();
391 delete fEvent;
392 fEvent=0x0;
393
394 }
395 bunchCounter++;
396
397 if(bunchCounter>=nTrainBunches){
398
399 trainCounter++;
400 if(trainCounter>=nTrains){
401 eventTime+=abortGap;
402 trainCounter=0;
403 }
404 else eventTime+=trainSpacing;
405
406 bunchCounter=0;
407 }
408 else eventTime+= bSpacing;
409
410 }
411
412
413 }
414 s.Stop();
415 s.Print();
416 delete rand;
417 CloseOutputFile();
418}
29b7f480 419
420
421
422
423
424//________________________________________________________________
425Int_t AliToyMCEventGeneratorSimple::OpenInputAndGetMaxEvents(const Int_t type, const Int_t nevents) {
426
427
428
429 if(type==0) return nevents;
430
431 if(type==1) {
432
433 fInputFile = new TFile(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"),"read");
434 if(!fInputFile->IsOpen()) {
435 std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
436 return 0;
437 }
438 fESDTree = (TTree*) fInputFile->Get("esdTree");
439 if(!fESDTree) {
440 std::cout <<"no esdTree in file " << std::endl;
441 return 0;
442 }
443 fESDEvent = new AliESDEvent();
444 fESDEvent->ReadFromTree(fESDTree);
445 fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
446
447 fInputIndex = 0;
448
1e62e876 449 gRandom->SetSeed();
36a79b0d 450 return fESDTree->GetEntries();
29b7f480 451 }
452
453
454 std::cout << " no available input type (toymc, esd) specified" << std::endl;
455 return 0;
456
457
458}
459
460
461//________________________________________________________________
462void AliToyMCEventGeneratorSimple::RunSimulation2(const Bool_t equalspacing, const Int_t type, const Int_t nevents, const Int_t ntracks) {
463
464 //type==0 simple toy
465 //type==1 esd input
466
467 Int_t nMaxEvents = OpenInputAndGetMaxEvents(type,nevents);
468 if (!ConnectOutputFile()) return;
469
223d9e38 470 InitSpaceCharge();
29b7f480 471
472
473 fNtracks = ntracks;
223d9e38 474 // within one simulation the track count should be unique for effeciency studies
a06336b6 475 fCurrentTrack=1;
223d9e38 476
29b7f480 477 Double_t eventTime=0.;
478 TStopwatch s;
479 // const Double_t eventSpacing=1./50e3; //50kHz equally spaced
480
481 Int_t nGeneratedEvents = 0;
482
483 for (Int_t ievent=0; ievent<nMaxEvents; ++ievent){
484 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
485
486 Int_t nColls = 0;
487 Double_t spacing = 0;
488 GetNGeneratedEventsAndSpacing(equalspacing, nColls, spacing);
489
490 for(Int_t iColl = 0; iColl<nColls; iColl++){
491 if(type==0) fEvent = Generate(eventTime);
492 else if(type==1) fEvent = GenerateESD2(eventTime);
493 if(fEvent) {
494 FillTree();
495 delete fEvent;
496 fEvent=0x0;
497 nGeneratedEvents++;
498 }
499 }
500 eventTime+=spacing;
501 if(nGeneratedEvents >= nevents) break;
502 }
503 s.Stop();
504 s.Print();
505
506 CloseOutputFile();
507 CloseInputFile();
508
509
510}
511//________________________________________________________________
512AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD2(Double_t time) {
513
514 //test that enough tracks will pass cuts (and that there is tracks at all)
515 Bool_t testEvent = kTRUE;
516 while(fESDTree->GetEvent(fInputIndex) && testEvent) {
517
518 Int_t nPassedCuts = 0;
519 for(Int_t iTrack = 0; iTrack<fESDEvent->GetNumberOfTracks(); iTrack++){
520 if(fESDCuts->AcceptTrack(fESDEvent->GetTrack(iTrack)) && (fESDEvent->GetTrack(iTrack)->Pt() < 0.3) ) nPassedCuts++;
521 }
522 if(nPassedCuts<fNtracks || (fNtracks>100 && nPassedCuts <100) ) fInputIndex++; //100 is a perhaps bit arbitrary, do we want the events were only a small number of tracks are accepted?
523 else testEvent = kFALSE;
524 }
525
526 if(fInputIndex>=fESDTree->GetEntries()) return 0;
527
528
529
530 AliToyMCEvent *retEvent = new AliToyMCEvent();
531 retEvent->SetT0(time);
532 retEvent->SetX(fESDEvent->GetPrimaryVertex()->GetX());
533 retEvent->SetY(fESDEvent->GetPrimaryVertex()->GetY());
534 retEvent->SetZ(fESDEvent->GetPrimaryVertex()->GetZ());
535
536
537
538 if(!fNtracks) fNtracks = fESDEvent->GetNumberOfTracks();
539 Int_t nUsedTracks = 0;
540 for(Int_t iTrack=0; iTrack<fESDEvent->GetNumberOfTracks(); iTrack++){
541
542 AliESDtrack *part = fESDEvent->GetTrack(iTrack);
543 if(!part) continue;
544 if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue;
545 if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
546
547 Double_t pxyz[3];
548 pxyz[0]=part->Px();
549 pxyz[1]=part->Py();
550 pxyz[2]=part->Pz();
551 Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
552 Double_t cv[21]={0};
553 Int_t sign = part->Charge();
554
0403120d 555 AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
29b7f480 556 // use unique ID for track number
557 // this will be used in DistortTrack track to set the cluster label
223d9e38 558 // in one simulation the track id should be unique for performance studies
559 tempTrack->SetUniqueID(fCurrentTrack++);
0403120d 560 DistortTrack(*tempTrack, time);
561
29b7f480 562
563 if(nUsedTracks >= fNtracks) break;
564 }
565 fInputIndex++;
566
567 return retEvent;
568}
569
570//________________________________________________________________
0403120d 571AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateLaser(Double_t time)
572{
573 //
574 // Generate an Event with laser tracks
575 //
576
577 AliToyMCEvent *retEvent = new AliToyMCEvent();
578 retEvent->SetEventType(AliToyMCEvent::kLaser);
579
580 retEvent->SetT0(time);
581 retEvent->SetX(0);
582 retEvent->SetX(0);
583 retEvent->SetZ(0);
584
585 AliTPCLaserTrack::LoadTracks();
586 TObjArray *arr=AliTPCLaserTrack::GetTracks();
587
588 //since we have a laser track force no material budges
589 Bool_t materialBudget=GetUseMaterialBudget();
590 SetUseMaterialBudget(kFALSE);
cedbc66e 591
592 //the laser tracks have partially large inclination angles over the pads
593 // -> relax the propagation contraint and switch off error messages
0403120d 594 SetIsLaser(kTRUE);
cedbc66e 595 Int_t debug=AliLog::GetDebugLevel("","AliToyMCEventGeneratorSimple");
596 AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",-3);
0403120d 597
598 for (Int_t iTrack=0; iTrack<arr->GetEntriesFast(); ++iTrack){
599 AliExternalTrackParam *track=(AliExternalTrackParam*)arr->At(iTrack);
600 AliToyMCTrack *tempTrack = retEvent->AddTrack(AliToyMCTrack(*track));
601 // for laser only TPC clusters exist
602 tempTrack->SetUniqueID(fCurrentTrack++);
603 MakeTPCClusters(*tempTrack, time);
604 }
605
606 SetIsLaser(kFALSE);
cedbc66e 607 AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",debug);
0403120d 608 SetUseMaterialBudget(materialBudget);
609 return retEvent;
610}
611
612//________________________________________________________________
29b7f480 613void AliToyMCEventGeneratorSimple::GetNGeneratedEventsAndSpacing(const Bool_t equalSpacing, Int_t &ngen, Double_t &spacing)
614{
615
616 static Int_t bunchCounter = 0;
617 static Int_t trainCounter = 0;
618
619 if(equalSpacing) {
620 ngen =1;
621 spacing = 1./50e3; //50kHz equally spaced
622 return;
623 }
624 else if(!equalSpacing) {
625 //parameters for bunch train
626 const Double_t abortGap = 3e-6; //
627 const Double_t collFreq = 50e3;
628 const Double_t bSpacing = 50e-9; //bunch spacing
629 const Int_t nTrainBunches = 48;
630 const Int_t nTrains = 12;
631 const Double_t revFreq = 1.11e4; //revolution frequency
632 const Double_t collProb = collFreq/(nTrainBunches*nTrains*revFreq);
633 const Double_t trainLength = bSpacing*(nTrainBunches-1);
634 const Double_t totTrainLength = nTrains*trainLength;
635 const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1);
636 Double_t time = 0;
637 Int_t nCollsInCrossing = 0;
638 while(nCollsInCrossing ==0){
639
640
641 bunchCounter++;
642
643 if(bunchCounter>=nTrainBunches){
644
645 trainCounter++;
646 if(trainCounter>=nTrains){
647 time+=abortGap;
648 trainCounter=0;
649 }
650 else time+=trainSpacing;
651
652 bunchCounter=0;
653 }
654 else time+= bSpacing;
655
656
657 nCollsInCrossing = gRandom -> Poisson(collProb);
658 //std::cout << " nCollsInCrossing " << nCollsInCrossing << std::endl;
659 }
660 ngen = nCollsInCrossing;
661 if(nCollsInCrossing > 1)std::cout << " nCollsInCrossing " << nCollsInCrossing << std::endl;
662 spacing = time;
663 return;
664
665 }
666
667}
668
669//________________________________________________________________
670Bool_t AliToyMCEventGeneratorSimple::CloseInputFile()
671{
672 //
673 // close the input file
674 //
675 if (!fInputFile) return kFALSE;
676 fInputFile->Close();
677 delete fInputFile;
678 fInputFile=0x0;
679
680 return kTRUE;
681}