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
 101     Bool_t trackDist = DistortTrack(*tempTrack, time);
 102     if(trackDist) retEvent->AddTrack(*tempTrack);
 103     delete tempTrack;
 104   }
 105
 106   return retEvent;
 107 }
 108
 109 //________________________________________________________________
 110 void AliToyMCEventGeneratorSimple::RunSimulation(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
 111 {
 112   //
 113   // run simple simulation with equal event spacing
 114   //
 115
 116   if (!ConnectOutputFile()) return;
 117   //initialise the space charge. Should be done after the tree was set up
 118   InitSpaceCharge();
 119
 120   fNtracks=ntracks;
 121   Double_t eventTime=0.;
 122   const Double_t eventSpacing=1./50e3; //50kHz equally spaced
 123   TStopwatch s;
 124   for (Int_t ievent=0; ievent<nevents; ++ievent){
 125     printf("Generating event %3d (%.3g)\n",ievent,eventTime);
 126     fEvent = Generate(eventTime);
 127     FillTree();
 128     delete fEvent;
 129     fEvent=0x0;
 130     eventTime+=eventSpacing;
 131   }
 132   s.Stop();
 133   s.Print();
 134
 135   CloseOutputFile();
 136 }
 137
 138 //________________________________________________________________
 139 AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD(AliESDEvent &esdEvent, Double_t time) {
 140
 141
 142   AliToyMCEvent *retEvent = new AliToyMCEvent();
 143   retEvent->SetT0(time);
 144   retEvent->SetX(esdEvent.GetPrimaryVertex()->GetX());
 145   retEvent->SetY(esdEvent.GetPrimaryVertex()->GetY());
 146   retEvent->SetZ(esdEvent.GetPrimaryVertex()->GetZ());
 147
 148
 149
 150   if(!fNtracks) fNtracks =  esdEvent.GetNumberOfTracks();
 151   Int_t nUsedTracks = 0;
 152   for(Int_t iTrack=0; iTrack<esdEvent.GetNumberOfTracks(); iTrack++){
 153
 154     AliESDtrack *part = esdEvent.GetTrack(iTrack);
 155     if(!part) continue;
 156     if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue;
 157     if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
 158     Double_t pxyz[3];
 159     pxyz[0]=part->Px();
 160     pxyz[1]=part->Py();
 161     pxyz[2]=part->Pz();
 162     Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
 163
 164     Double_t cv[21]={0};
 165     Int_t sign = part->Charge();
 166     AliToyMCTrack *tempTrack = new AliToyMCTrack(vertex,pxyz,cv,sign);
 167
 168     Bool_t trackDist = DistortTrack(*tempTrack, time);
 169     if(trackDist) {
 170       retEvent->AddTrack(*tempTrack);
 171
 172       nUsedTracks++;
 173     }
 174     delete tempTrack;
 175
 176     if(nUsedTracks >= fNtracks) break;
 177   }
 178
 179   return retEvent;
 180 }
 181 //________________________________________________________________
 182 void AliToyMCEventGeneratorSimple::RunSimulationESD(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
 183 {
 184   //
 185   // run simulation using esd input with equal event spacing
 186   //
 187
 188   if (!ConnectOutputFile()) return;
 189   TFile f(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"));
 190   if(!f.IsOpen()) {
 191     std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
 192     return;
 193   }
 194   TTree* esdTree = (TTree*) f.Get("esdTree");
 195   if(!esdTree) {
 196     std::cout <<"no esdTree in file " << std::endl;
 197     return;
 198   }
 199   InitSpaceCharge();
 200   AliESDEvent* esdEvent = new AliESDEvent();
 201   fNtracks = ntracks;
 202
 203   esdEvent->ReadFromTree(esdTree);
 204   Double_t eventTime=0.;
 205   const Double_t eventSpacing=1./50e3; //50kHz equally spaced
 206   TStopwatch s;
 207   Int_t nEvents = nevents;
 208   fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
 209   if(nevents>esdTree->GetEntries()) nEvents = esdTree->GetEntries();
 210   Int_t nUsedEvents = 0;
 211   for (Int_t ievent=0; ievent<esdTree->GetEntries(); ++ievent){
 212     printf("Generating event %3d (%.3g)\n",nUsedEvents,eventTime);
 213     esdTree->GetEvent(ievent);
 214     if(esdEvent->GetNumberOfTracks()==0) {
 215       std::cout << " tracks == 0" << std::endl;
 216       continue;
 217     }
 218
 219     fEvent = GenerateESD(*esdEvent, eventTime);
 220     if(fEvent->GetNumberOfTracks() >=10) {
 221       nUsedEvents++;
 222       FillTree();
 223       eventTime+=eventSpacing;
 224     }
 225     delete fEvent;
 226     fEvent=0x0;
 227
 228     if(nUsedEvents>=nevents) break;
 229   }
 230   s.Stop();
 231   s.Print();
 232
 233   CloseOutputFile();
 234 }
 235
 236 //________________________________________________________________
 237
 238 void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
 239 {
 240   //
 241   // run simple simulation with equal event spacing
 242   //
 243
 244   //Parameters for bunc
 245   const Double_t abortGap = 3e-6; //
 246   const Double_t collFreq = 50e3;
 247   const Double_t bSpacing = 50e-9; //bunch spacing
 248   const Int_t nTrainBunches = 48;
 249   const Int_t nTrains = 12;
 250   const Double_t revFreq = 1.11e4; //revolution frequency
 251   const Double_t collProb = collFreq/(nTrainBunches*nTrains*revFreq);
 252   const Double_t trainLength = bSpacing*(nTrainBunches-1);
 253   const Double_t totTrainLength = nTrains*trainLength;
 254   const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1);
 255   Bool_t equalSpacing = kFALSE;
 256
 257   TRandom3 *rand = new TRandom3();
 258   //rand->SetSeed();
 259
 260   if (!ConnectOutputFile()) return;
 261   //initialise the space charge. Should be done after the tree was set up
 262   InitSpaceCharge();
 263
 264   fNtracks=ntracks;
 265   Double_t eventTime=0.;
 266   const Double_t eventSpacing=1./50e3; //50kHz equally spaced
 267   TStopwatch s;
 268   Int_t nGeneratedEvents = 0;
 269   Int_t bunchCounter = 0;
 270   Int_t trainCounter = 0;
 271
 272   //for (Int_t ievent=0; ievent<nevents; ++ievent){
 273   while (nGeneratedEvents<nevents){
 274     //  std::cout <<trainCounter << " " << bunchCounter << " "<< "eventTime " << eventTime << std::endl;
 275
 276     if(equalSpacing)  {
 277       printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
 278       fEvent = Generate(eventTime);
 279       nGeneratedEvents++;
 280       FillTree();
 281       delete fEvent;
 282       fEvent=0x0;
 283       eventTime+=1./collFreq;
 284     }
 285     else{
 286       Int_t nCollsInCrossing = rand -> Poisson(collProb);
 287       for(Int_t iColl = 0; iColl<nCollsInCrossing; iColl++){
 288         printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
 289         fEvent = Generate(eventTime);
 290         nGeneratedEvents++;
 291         FillTree();
 292         delete fEvent;
 293         fEvent=0x0;
 294
 295       }
 296       bunchCounter++;
 297
 298       if(bunchCounter>=nTrainBunches){
 299
 300         trainCounter++;
 301         if(trainCounter>=nTrains){
 302           eventTime+=abortGap;
 303           trainCounter=0;
 304         }
 305         else eventTime+=trainSpacing;
 306
 307         bunchCounter=0;
 308       }
 309       else eventTime+= bSpacing;
 310
 311     }
 312
 313
 314   }
 315   s.Stop();
 316   s.Print();
 317   delete rand;
 318   CloseOutputFile();
 319 }