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 <AliESDtrackCuts.h>
   8 #include <AliESDtrack.h>
   9 #include <TFile.h>
  10 #include <TTree.h>
  11 #include <TRandom3.h>
  12
  13 #include "AliToyMCEvent.h"
  14
  15 #include "AliToyMCEventGeneratorSimple.h"
  16
  17
  18 ClassImp(AliToyMCEventGeneratorSimple);
  19
  20
  21 AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple()
  22   :AliToyMCEventGenerator()
  23   ,fVertexMean(0.)
  24   ,fVertexSigma(7.)
  25   ,fNtracks(20)
  26   ,fInputFileNameESD("")
  27   ,fESDCuts(0x0)
  28 {
  29   //
  30
  31 }
  32 //________________________________________________________________
  33 AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple(const AliToyMCEventGeneratorSimple &gen)
  34   :AliToyMCEventGenerator(gen)
  35   ,fVertexMean(gen.fVertexMean)
  36   ,fVertexSigma(gen.fVertexSigma)
  37   ,fNtracks(gen.fNtracks)
  38   ,fInputFileNameESD(gen.fInputFileNameESD)
  39   ,fESDCuts(gen.fESDCuts)
  40 {
  41 }
  42 //________________________________________________________________
  43 AliToyMCEventGeneratorSimple::~AliToyMCEventGeneratorSimple()
  44 {
  45 }
  46  //________________________________________________________________
  47 AliToyMCEventGeneratorSimple& AliToyMCEventGeneratorSimple::operator = (const AliToyMCEventGeneratorSimple &gen)
  48 {
  49   //assignment operator
  50   if (&gen == this) return *this;
  51   new (this) AliToyMCEventGeneratorSimple(gen);
  52
  53   return *this;
  54 }
  55 //________________________________________________________________
  56 void AliToyMCEventGeneratorSimple::SetParametersSimple(Double_t vertexMean, Double_t vertexSigma) {
  57   fVertexMean = vertexMean;
  58   fVertexSigma = vertexSigma;
  59 }
  60 //________________________________________________________________
  61 AliToyMCEvent* AliToyMCEventGeneratorSimple::Generate(Double_t time)
  62 {
  63   //
  64   //
  65   //
  66
  67   AliToyMCEvent *retEvent = new AliToyMCEvent();
  68   retEvent->SetT0(time);
  69   retEvent->SetX(0);
  70   retEvent->SetX(0);
  71   retEvent->SetZ(gRandom->Gaus(fVertexMean,fVertexSigma));
  72
  73   Double_t etaCuts=.9;
  74   Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
  75   TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
  76   fpt.SetParameters(7.24,0.120);
  77   fpt.SetNpx(10000);
  78 //   Int_t nTracks = 400; //TODO: draw from experim dist
  79 //   Int_t nTracks = 20; //TODO: draw from experim dist
  80
  81   for(Int_t iTrack=0; iTrack<fNtracks; iTrack++){
  82     Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
  83     Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
  84     Double_t pt = fpt.GetRandom(); // momentum for f1
  85     Short_t sign=1;
  86     if(gRandom->Rndm() < 0.5){
  87       sign =1;
  88     }else{
  89       sign=-1;
  90     }
  91
  92     Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
  93     Double_t pxyz[3];
  94     pxyz[0]=pt*TMath::Cos(phi);
  95     pxyz[1]=pt*TMath::Sin(phi);
  96     pxyz[2]=pt*TMath::Tan(theta);
  97     Double_t vertex[3]={0,0,retEvent->GetZ()};
  98     Double_t cv[21]={0};
  99     AliToyMCTrack *tempTrack = new AliToyMCTrack(vertex,pxyz,cv,sign);
 100     // use unique ID for track number
 101     // this will be used in DistortTrack track to set the cluster label
 102     tempTrack->SetUniqueID(iTrack);
 103
 104     Bool_t trackDist = DistortTrack(*tempTrack, time);
 105     if(trackDist) retEvent->AddTrack(*tempTrack);
 106     delete tempTrack;
 107   }
 108
 109   return retEvent;
 110 }
 111
 112 //________________________________________________________________
 113 void AliToyMCEventGeneratorSimple::RunSimulation(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
 114 {
 115   //
 116   // run simple simulation with equal event spacing
 117   //
 118
 119   if (!ConnectOutputFile()) return;
 120   //initialise the space charge. Should be done after the tree was set up
 121   InitSpaceCharge();
 122
 123   fNtracks=ntracks;
 124   Double_t eventTime=0.;
 125   const Double_t eventSpacing=1./50e3; //50kHz equally spaced
 126   TStopwatch s;
 127   for (Int_t ievent=0; ievent<nevents; ++ievent){
 128     printf("Generating event %3d (%.3g)\n",ievent,eventTime);
 129     fEvent = Generate(eventTime);
 130     FillTree();
 131     delete fEvent;
 132     fEvent=0x0;
 133     eventTime+=eventSpacing;
 134   }
 135   s.Stop();
 136   s.Print();
 137
 138   CloseOutputFile();
 139 }
 140
 141 //________________________________________________________________
 142 AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD(AliESDEvent &esdEvent, Double_t time) {
 143
 144
 145   AliToyMCEvent *retEvent = new AliToyMCEvent();
 146   retEvent->SetT0(time);
 147   retEvent->SetX(esdEvent.GetPrimaryVertex()->GetX());
 148   retEvent->SetY(esdEvent.GetPrimaryVertex()->GetY());
 149   retEvent->SetZ(esdEvent.GetPrimaryVertex()->GetZ());
 150
 151
 152
 153   if(!fNtracks) fNtracks =  esdEvent.GetNumberOfTracks();
 154   Int_t nUsedTracks = 0;
 155   for(Int_t iTrack=0; iTrack<esdEvent.GetNumberOfTracks(); iTrack++){
 156
 157     AliESDtrack *part = esdEvent.GetTrack(iTrack);
 158     if(!part) continue;
 159     if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue;
 160     if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
 161     Double_t pxyz[3];
 162     pxyz[0]=part->Px();
 163     pxyz[1]=part->Py();
 164     pxyz[2]=part->Pz();
 165     Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
 166
 167     Double_t cv[21]={0};
 168     Int_t sign = part->Charge();
 169     AliToyMCTrack *tempTrack = new AliToyMCTrack(vertex,pxyz,cv,sign);
 170     // use unique ID for track number
 171     // this will be used in DistortTrack track to set the cluster label
 172     tempTrack->SetUniqueID(iTrack);
 173
 174
 175     Bool_t trackDist = DistortTrack(*tempTrack, time);
 176     if(trackDist) {
 177       retEvent->AddTrack(*tempTrack);
 178
 179       nUsedTracks++;
 180     }
 181     delete tempTrack;
 182
 183     if(nUsedTracks >= fNtracks) break;
 184   }
 185
 186   return retEvent;
 187 }
 188 //________________________________________________________________
 189 void AliToyMCEventGeneratorSimple::RunSimulationESD(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
 190 {
 191   //
 192   // run simulation using esd input with equal event spacing
 193   //
 194
 195   if (!ConnectOutputFile()) return;
 196   TFile f(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"));
 197   if(!f.IsOpen()) {
 198     std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
 199     return;
 200   }
 201   TTree* esdTree = (TTree*) f.Get("esdTree");
 202   if(!esdTree) {
 203     std::cout <<"no esdTree in file " << std::endl;
 204     return;
 205   }
 206   InitSpaceCharge();
 207   AliESDEvent* esdEvent = new AliESDEvent();
 208   fNtracks = ntracks;
 209
 210   esdEvent->ReadFromTree(esdTree);
 211   Double_t eventTime=0.;
 212   const Double_t eventSpacing=1./50e3; //50kHz equally spaced
 213   TStopwatch s;
 214   Int_t nEvents = nevents;
 215   fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
 216   if(nevents>esdTree->GetEntries()) nEvents = esdTree->GetEntries();
 217   Int_t nUsedEvents = 0;
 218   for (Int_t ievent=0; ievent<esdTree->GetEntries(); ++ievent){
 219     printf("Generating event %3d (%.3g)\n",nUsedEvents,eventTime);
 220     esdTree->GetEvent(ievent);
 221     if(esdEvent->GetNumberOfTracks()==0) {
 222       std::cout << " tracks == 0" << std::endl;
 223       continue;
 224     }
 225
 226     fEvent = GenerateESD(*esdEvent, eventTime);
 227     if(fEvent->GetNumberOfTracks() >=10) {
 228       nUsedEvents++;
 229       FillTree();
 230       eventTime+=eventSpacing;
 231     }
 232     delete fEvent;
 233     fEvent=0x0;
 234
 235     if(nUsedEvents>=nevents) break;
 236   }
 237   s.Stop();
 238   s.Print();
 239
 240   CloseOutputFile();
 241 }
 242
 243 //________________________________________________________________
 244
 245 void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
 246 {
 247   //
 248   // run simple simulation with equal event spacing
 249   //
 250
 251   //Parameters for bunc
 252   const Double_t abortGap = 3e-6; //
 253   const Double_t collFreq = 50e3;
 254   const Double_t bSpacing = 50e-9; //bunch spacing
 255   const Int_t nTrainBunches = 48;
 256   const Int_t nTrains = 12;
 257   const Double_t revFreq = 1.11e4; //revolution frequency
 258   const Double_t collProb = collFreq/(nTrainBunches*nTrains*revFreq);
 259   const Double_t trainLength = bSpacing*(nTrainBunches-1);
 260   const Double_t totTrainLength = nTrains*trainLength;
 261   const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1);
 262   Bool_t equalSpacing = kFALSE;
 263
 264   TRandom3 *rand = new TRandom3();
 265   //rand->SetSeed();
 266
 267   if (!ConnectOutputFile()) return;
 268   //initialise the space charge. Should be done after the tree was set up
 269   InitSpaceCharge();
 270
 271   fNtracks=ntracks;
 272   Double_t eventTime=0.;
 273   const Double_t eventSpacing=1./50e3; //50kHz equally spaced
 274   TStopwatch s;
 275   Int_t nGeneratedEvents = 0;
 276   Int_t bunchCounter = 0;
 277   Int_t trainCounter = 0;
 278
 279   //for (Int_t ievent=0; ievent<nevents; ++ievent){
 280   while (nGeneratedEvents<nevents){
 281     //  std::cout <<trainCounter << " " << bunchCounter << " "<< "eventTime " << eventTime << std::endl;
 282
 283     if(equalSpacing)  {
 284       printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
 285       fEvent = Generate(eventTime);
 286       nGeneratedEvents++;
 287       FillTree();
 288       delete fEvent;
 289       fEvent=0x0;
 290       eventTime+=1./collFreq;
 291     }
 292     else{
 293       Int_t nCollsInCrossing = rand -> Poisson(collProb);
 294       for(Int_t iColl = 0; iColl<nCollsInCrossing; iColl++){
 295         printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
 296         fEvent = Generate(eventTime);
 297         nGeneratedEvents++;
 298         FillTree();
 299         delete fEvent;
 300         fEvent=0x0;
 301
 302       }
 303       bunchCounter++;
 304
 305       if(bunchCounter>=nTrainBunches){
 306
 307         trainCounter++;
 308         if(trainCounter>=nTrains){
 309           eventTime+=abortGap;
 310           trainCounter=0;
 311         }
 312         else eventTime+=trainSpacing;
 313
 314         bunchCounter=0;
 315       }
 316       else eventTime+= bSpacing;
 317
 318     }
 319
 320
 321   }
 322   s.Stop();
 323   s.Print();
 324   delete rand;
 325   CloseOutputFile();
 326 }