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