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 ClassImp(AliToyMCEventGeneratorSimple);
  21
  22
  23 AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple()
  24   :AliToyMCEventGenerator()
  25   ,fVertexMean(0.)
  26   ,fVertexSigma(7.)
  27   ,fNtracks(20)
  28   ,fInputFileNameESD("")
  29   ,fESDCuts(0x0)
  30   ,fInputIndex(0)
  31   ,fESDEvent(0x0)
  32   ,fESDTree(0x0)
  33   ,fInputFile(0x0)
  34   ,fHPt(0x0)
  35   ,fHEta(0x0)
  36   ,fHMult(0x0)
  37   ,fHistosSet(kFALSE)
  38   ,fParamFile(0x0)
  39 {
  40   //
  41
  42 }
  43 //________________________________________________________________
  44 AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple(const AliToyMCEventGeneratorSimple &gen)
  45   :AliToyMCEventGenerator(gen)
  46   ,fVertexMean(gen.fVertexMean)
  47   ,fVertexSigma(gen.fVertexSigma)
  48   ,fNtracks(gen.fNtracks)
  49   ,fInputFileNameESD(gen.fInputFileNameESD)
  50   ,fESDCuts(gen.fESDCuts)
  51   ,fInputIndex(gen.fInputIndex)
  52   ,fESDEvent(gen.fESDEvent)
  53   ,fESDTree(gen.fESDTree)
  54   ,fInputFile(gen.fInputFile)
  55   ,fHPt(gen.fHPt)
  56   ,fHEta(gen.fHEta)
  57   ,fHMult(gen.fHMult)
  58   ,fHistosSet(gen.fHistosSet)
  59   ,fParamFile(0x0)
  60 {
  61 }
  62 //________________________________________________________________
  63 AliToyMCEventGeneratorSimple::~AliToyMCEventGeneratorSimple()
  64 {
  65 }
  66  //________________________________________________________________
  67 AliToyMCEventGeneratorSimple& AliToyMCEventGeneratorSimple::operator = (const AliToyMCEventGeneratorSimple &gen)
  68 {
  69   //assignment operator
  70   if (&gen == this) return *this;
  71   new (this) AliToyMCEventGeneratorSimple(gen);
  72
  73   return *this;
  74 }
  75 //________________________________________________________________
  76 void AliToyMCEventGeneratorSimple::SetParametersToyGen(const Char_t* parfilename/*="files/params.root*/, Double_t vertexMean/*=0*/, Double_t vertexSigma/*=7.*/) {
  77   fVertexMean = vertexMean;
  78   fVertexSigma = vertexSigma;
  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
  86 }
  87 //________________________________________________________________
  88 AliToyMCEvent* AliToyMCEventGeneratorSimple::Generate(Double_t time)
  89 {
  90   //
  91   //
  92   //
  93
  94   AliToyMCEvent *retEvent = new AliToyMCEvent();
  95   retEvent->SetT0(time);
  96   retEvent->SetX(0);
  97   retEvent->SetX(0);
  98   Double_t vertexZ=0;
  99   //limit vertex to +- 10cm
 100   while ( TMath::Abs(vertexZ=gRandom->Gaus(fVertexMean,fVertexSigma))>10. );
 101   retEvent->SetZ(vertexZ);
 102
 103   Double_t etaCuts=0.9;
 104   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 105   static TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.3,100);
 106   if (fpt.GetParameter(0)<1){
 107     printf("Set Parameters\n");
 108     fpt.SetParameters(7.24,0.120);
 109     fpt.SetNpx(200);
 110   }
 111 //   Int_t nTracks = 400; //TODO: draw from experim dist
 112 //   Int_t nTracks = 20; //TODO: draw from experim dist
 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;
 119
 120   for(Int_t iTrack=0; iTrack<nTracksLocal; iTrack++){
 121     Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
 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     }
 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};
 147     AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
 148     // use unique ID for track number
 149     // this will be used in DistortTrack track to set the cluster label
 150     // in one simulation the track id should be unique for performance studies
 151     tempTrack->SetUniqueID(fCurrentTrack++);
 152     DistortTrack(*tempTrack, time);
 153   }
 154
 155   return retEvent;
 156 }
 157
 158 //________________________________________________________________
 159 void AliToyMCEventGeneratorSimple::RunSimulation(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/, const Int_t rate/*=50*/)
 160 {
 161   //
 162   // run simple simulation with equal event spacing
 163   // rate in kHz
 164   //
 165
 166   if (!ConnectOutputFile()) return;
 167   //initialise the space charge. Should be done after the tree was set up
 168   InitSpaceCharge();
 169
 170   // number of tracks to simulate per interaction
 171   fNtracks=ntracks;
 172   // within one simulation the track count should be unique for effeciency studies
 173   // don't use the track ID 0 since the label -0 is not different from 0
 174   fCurrentTrack=1;
 175
 176   Double_t eventTime=0.;
 177   const Double_t eventSpacing=1./rate/1e3; //50kHz equally spaced
 178   TStopwatch s;
 179   for (Int_t ievent=0; ievent<nevents; ++ievent){
 180     printf("Generating event %3d (%.3g)\n",ievent,eventTime);
 181     fEvent = Generate(eventTime);
 182     FillTree();
 183     delete fEvent;
 184     fEvent=0x0;
 185     eventTime+=eventSpacing;
 186   }
 187   s.Stop();
 188   s.Print();
 189
 190   CloseOutputFile();
 191 }
 192
 193 //________________________________________________________________
 194 void 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
 205   fCurrentTrack=1;
 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 //________________________________________________________________
 225 AliToyMCEvent* 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;
 243     if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
 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();
 252     AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
 253     // use unique ID for track number
 254     // this will be used in DistortTrack track to set the cluster label
 255     // in one simulation the track id should be unique for performance studies
 256     tempTrack->SetUniqueID(fCurrentTrack++);
 257     DistortTrack(*tempTrack, time);
 258     nUsedTracks++;
 259
 260     if(nUsedTracks >= fNtracks) break;
 261   }
 262
 263   return retEvent;
 264 }
 265 //________________________________________________________________
 266 void 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();
 285
 286
 287   esdEvent->ReadFromTree(esdTree);
 288
 289   fNtracks = ntracks;
 290   // within one simulation the track count should be unique for effeciency studies
 291   fCurrentTrack=1;
 292
 293   Double_t eventTime=0.;
 294   const Double_t eventSpacing=1./50e3; //50kHz equally spaced
 295   TStopwatch s;
 296
 297   fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
 298   gRandom->SetSeed();
 299   Int_t nUsedEvents = 0;
 300   for (Int_t ievent=0; ievent<esdTree->GetEntries(); ++ievent){
 301     printf("Generating event %3d (%.3g)\n",nUsedEvents,eventTime);
 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);
 309     if(fEvent->GetNumberOfTracks() >=10) {
 310       nUsedEvents++;
 311       FillTree();
 312       eventTime+=eventSpacing;
 313     }
 314     delete fEvent;
 315     fEvent=0x0;
 316
 317     if(nUsedEvents>=nevents) break;
 318   }
 319   s.Stop();
 320   s.Print();
 321
 322   CloseOutputFile();
 323 }
 324
 325 //________________________________________________________________
 326 void 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;
 353   // within one simulation the track count should be unique for effeciency studies
 354   fCurrentTrack=1;
 355
 356   Double_t eventTime=0.;
 357   //  const Double_t eventSpacing=1./50e3; //50kHz equally spaced
 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 }
 411
 412
 413
 414
 415
 416 //________________________________________________________________
 417 Int_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();
 442     gRandom->SetSeed();
 443    }
 444
 445
 446   std::cout << " no available input type (toymc, esd) specified" << std::endl;
 447   return 0;
 448
 449
 450 }
 451
 452
 453 //________________________________________________________________
 454 void 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
 462   InitSpaceCharge();
 463
 464
 465   fNtracks = ntracks;
 466   // within one simulation the track count should be unique for effeciency studies
 467   fCurrentTrack=1;
 468
 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 //________________________________________________________________
 504 AliToyMCEvent* 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
 547     AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
 548     // use unique ID for track number
 549     // this will be used in DistortTrack track to set the cluster label
 550     // in one simulation the track id should be unique for performance studies
 551     tempTrack->SetUniqueID(fCurrentTrack++);
 552     DistortTrack(*tempTrack, time);
 553
 554
 555     if(nUsedTracks >= fNtracks) break;
 556   }
 557   fInputIndex++;
 558
 559   return retEvent;
 560 }
 561
 562 //________________________________________________________________
 563 AliToyMCEvent* 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);
 583
 584   //the laser tracks have partially large inclination angles over the pads
 585   // -> relax the propagation contraint and switch off error messages
 586   SetIsLaser(kTRUE);
 587   Int_t debug=AliLog::GetDebugLevel("","AliToyMCEventGeneratorSimple");
 588   AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",-3);
 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);
 599   AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",debug);
 600   SetUseMaterialBudget(materialBudget);
 601   return retEvent;
 602 }
 603
 604 //________________________________________________________________
 605 void 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 //________________________________________________________________
 662 Bool_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 }