fixed bug that set visualization time to wrong value
[u/mrichter/AliRoot.git] / TPC / Upgrade / AliToyMCEventGeneratorSimple.cxx
CommitLineData
526ddf0e 1#include <iostream>
a1a695e5 2
526ddf0e 3#include <TDatabasePDG.h>
4#include <TRandom.h>
5#include <TF1.h>
32438f4e 6#include <TStopwatch.h>
f356075c 7#include <AliESDtrackCuts.h>
8#include <AliESDtrack.h>
9#include <TFile.h>
10#include <TTree.h>
609d6d39 11#include <TRandom3.h>
a1a695e5 12
de0014b7 13#include "AliToyMCEvent.h"
526ddf0e 14
a1a695e5 15#include "AliToyMCEventGeneratorSimple.h"
16
f356075c 17
de0014b7 18ClassImp(AliToyMCEventGeneratorSimple);
526ddf0e 19
20
de0014b7 21AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple()
22 :AliToyMCEventGenerator()
526ddf0e 23 ,fVertexMean(0.)
24 ,fVertexSigma(7.)
32438f4e 25 ,fNtracks(20)
f356075c 26 ,fInputFileNameESD("")
27 ,fESDCuts(0x0)
526ddf0e 28{
29 //
30
31}
32//________________________________________________________________
de0014b7 33AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple(const AliToyMCEventGeneratorSimple &gen)
34 :AliToyMCEventGenerator(gen)
526ddf0e 35 ,fVertexMean(gen.fVertexMean)
36 ,fVertexSigma(gen.fVertexSigma)
32438f4e 37 ,fNtracks(gen.fNtracks)
f356075c 38 ,fInputFileNameESD(gen.fInputFileNameESD)
39 ,fESDCuts(gen.fESDCuts)
526ddf0e 40{
41}
42//________________________________________________________________
de0014b7 43AliToyMCEventGeneratorSimple::~AliToyMCEventGeneratorSimple()
526ddf0e 44{
45}
46 //________________________________________________________________
de0014b7 47AliToyMCEventGeneratorSimple& AliToyMCEventGeneratorSimple::operator = (const AliToyMCEventGeneratorSimple &gen)
526ddf0e 48{
49 //assignment operator
50 if (&gen == this) return *this;
de0014b7 51 new (this) AliToyMCEventGeneratorSimple(gen);
526ddf0e 52
53 return *this;
54}
55//________________________________________________________________
f356075c 56void AliToyMCEventGeneratorSimple::SetParametersSimple(Double_t vertexMean, Double_t vertexSigma) {
526ddf0e 57 fVertexMean = vertexMean;
58 fVertexSigma = vertexSigma;
59}
60//________________________________________________________________
cd8ed0ac 61AliToyMCEvent* AliToyMCEventGeneratorSimple::Generate(Double_t time)
62{
63 //
64 //
65 //
66
de0014b7 67 AliToyMCEvent *retEvent = new AliToyMCEvent();
526ddf0e 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);
32438f4e 78// Int_t nTracks = 400; //TODO: draw from experim dist
79// Int_t nTracks = 20; //TODO: draw from experim dist
526ddf0e 80
32438f4e 81 for(Int_t iTrack=0; iTrack<fNtracks; iTrack++){
526ddf0e 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};
de0014b7 99 AliToyMCTrack *tempTrack = new AliToyMCTrack(vertex,pxyz,cv,sign);
7d26e3ce 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
526ddf0e 104 Bool_t trackDist = DistortTrack(*tempTrack, time);
105 if(trackDist) retEvent->AddTrack(*tempTrack);
106 delete tempTrack;
107 }
108
109 return retEvent;
110}
111
a1a695e5 112//________________________________________________________________
609d6d39 113void AliToyMCEventGeneratorSimple::RunSimulation(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
a1a695e5 114{
115 //
116 // run simple simulation with equal event spacing
117 //
118
119 if (!ConnectOutputFile()) return;
cd8ed0ac 120 //initialise the space charge. Should be done after the tree was set up
121 InitSpaceCharge();
122
32438f4e 123 fNtracks=ntracks;
a1a695e5 124 Double_t eventTime=0.;
125 const Double_t eventSpacing=1./50e3; //50kHz equally spaced
32438f4e 126 TStopwatch s;
a1a695e5 127 for (Int_t ievent=0; ievent<nevents; ++ievent){
32438f4e 128 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
a1a695e5 129 fEvent = Generate(eventTime);
130 FillTree();
131 delete fEvent;
132 fEvent=0x0;
133 eventTime+=eventSpacing;
134 }
32438f4e 135 s.Stop();
136 s.Print();
137
a1a695e5 138 CloseOutputFile();
139}
526ddf0e 140
f356075c 141//________________________________________________________________
142AliToyMCEvent* 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;
609d6d39 160 if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
f356075c 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);
7d26e3ce 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
f356075c 175 Bool_t trackDist = DistortTrack(*tempTrack, time);
176 if(trackDist) {
177 retEvent->AddTrack(*tempTrack);
609d6d39 178
f356075c 179 nUsedTracks++;
180 }
181 delete tempTrack;
182
183 if(nUsedTracks >= fNtracks) break;
184 }
185
186 return retEvent;
187}
188//________________________________________________________________
189void 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){
609d6d39 219 printf("Generating event %3d (%.3g)\n",nUsedEvents,eventTime);
f356075c 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);
609d6d39 227 if(fEvent->GetNumberOfTracks() >=10) {
228 nUsedEvents++;
229 FillTree();
230 eventTime+=eventSpacing;
231 }
f356075c 232 delete fEvent;
233 fEvent=0x0;
609d6d39 234
f356075c 235 if(nUsedEvents>=nevents) break;
236 }
237 s.Stop();
238 s.Print();
239
240 CloseOutputFile();
241}
242
609d6d39 243//________________________________________________________________
244
245void 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}