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     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(const 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(const Int_t nevents/*=10*/, const 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(const Int_t nevents/*=10*/, const 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(const Int_t type, const 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(const Bool_t equalspacing, const Int_t type, const Int_t nevents, const 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(const 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 }