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