Fixed ClassDef() accidentally removed
[u/mrichter/AliRoot.git] / TPC / TPCupgrade / 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(Int_t nevents/*=10*/, Int_t ntracks/*=400*/, 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     SetSCScalingFactor();
194     FillTree();
195     delete fEvent;
196     fEvent=0x0;
197     eventTime+=eventSpacing;
198   }
199   s.Stop();
200   s.Print();
201   
202   CloseOutputFile();
203 }
204
205 //________________________________________________________________
206 void AliToyMCEventGeneratorSimple::RunSimulationLaser(Int_t nevents/*=1*/)
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
217   fCurrentTrack=1;
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 //________________________________________________________________
237 AliToyMCEvent* 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;
255     if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
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();
264     AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
265     // use unique ID for track number
266     // this will be used in DistortTrack track to set the cluster label
267     // in one simulation the track id should be unique for performance studies
268     tempTrack->SetUniqueID(fCurrentTrack++);
269     DistortTrack(*tempTrack, time);
270     nUsedTracks++;
271     
272     if(nUsedTracks >= fNtracks) break;
273   }
274
275   return retEvent;
276 }
277 //________________________________________________________________
278 void AliToyMCEventGeneratorSimple::RunSimulationESD(Int_t nevents/*=10*/, Int_t ntracks/*=400*/)
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();
297  
298  
299   esdEvent->ReadFromTree(esdTree);
300
301   fNtracks = ntracks;
302   // within one simulation the track count should be unique for effeciency studies
303   fCurrentTrack=1;
304   
305   Double_t eventTime=0.;
306   const Double_t eventSpacing=1./50e3; //50kHz equally spaced
307   TStopwatch s;
308   
309   fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
310   gRandom->SetSeed();
311   Int_t nUsedEvents = 0;
312   for (Int_t ievent=0; ievent<esdTree->GetEntries(); ++ievent){
313     printf("Generating event %3d (%.3g)\n",nUsedEvents,eventTime);
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);
321     if(fEvent->GetNumberOfTracks() >=10) {
322       nUsedEvents++;
323       FillTree();
324       eventTime+=eventSpacing;
325     }
326     delete fEvent;
327     fEvent=0x0;
328     
329     if(nUsedEvents>=nevents) break;
330   }
331   s.Stop();
332   s.Print();
333   
334   CloseOutputFile();
335 }
336
337 //________________________________________________________________
338 void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(Int_t nevents/*=10*/, Int_t ntracks/*=400*/)
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;
365   // within one simulation the track count should be unique for effeciency studies
366   fCurrentTrack=1;
367   
368   Double_t eventTime=0.;
369   //  const Double_t eventSpacing=1./50e3; //50kHz equally spaced
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);
382       SetSCScalingFactor();
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);
394         SetSCScalingFactor();
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 }
425
426
427
428
429
430 //________________________________________________________________
431 Int_t AliToyMCEventGeneratorSimple::OpenInputAndGetMaxEvents(Int_t type, Int_t nevents) {
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
455     gRandom->SetSeed();
456     return fESDTree->GetEntries();
457    }
458
459  
460   std::cout << " no available input type (toymc, esd) specified" << std::endl;
461   return 0;
462  
463
464 }
465
466
467 //________________________________________________________________
468 void AliToyMCEventGeneratorSimple::RunSimulation2(Bool_t equalspacing, Int_t type, Int_t nevents, Int_t ntracks) {
469
470   //type==0 simple toy
471   //type==1 esd input
472
473   Int_t nMaxEvents = OpenInputAndGetMaxEvents(type,nevents);
474   if (!ConnectOutputFile()) return;
475   
476   InitSpaceCharge();
477
478
479   fNtracks = ntracks;
480   // within one simulation the track count should be unique for effeciency studies
481   fCurrentTrack=1;
482   
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 //________________________________________________________________
518 AliToyMCEvent* 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;
522
523   // iterate over space charge maps in case they are set
524   IterateSC();
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
564     AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
565     // use unique ID for track number
566     // this will be used in DistortTrack track to set the cluster label
567     // in one simulation the track id should be unique for performance studies
568     tempTrack->SetUniqueID(fCurrentTrack++);
569     DistortTrack(*tempTrack, time);
570     
571     
572     if(nUsedTracks >= fNtracks) break;
573   }
574   fInputIndex++;
575  
576   return retEvent;
577 }
578
579 //________________________________________________________________
580 AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateLaser(Double_t time)
581 {
582   //
583   // Generate an Event with laser tracks
584   //
585
586   // iterate over space charge maps in case they are set
587   IterateSC();
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);
603
604   //the laser tracks have partially large inclination angles over the pads
605   // -> relax the propagation contraint and switch off error messages
606   SetIsLaser(kTRUE);
607   Int_t debug=AliLog::GetDebugLevel("","AliToyMCEventGeneratorSimple");
608   AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",-3);
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);
619   AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",debug);
620   SetUseMaterialBudget(materialBudget);
621   return retEvent;
622 }
623
624 //________________________________________________________________
625 void AliToyMCEventGeneratorSimple::GetNGeneratedEventsAndSpacing(Bool_t equalSpacing, Int_t &ngen, Double_t &spacing)
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 //________________________________________________________________
682 Bool_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 }