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