TPC/Upgrade/AliToyMCEventGeneratorSimple.cxx

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