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 }