]>
Commit | Line | Data |
---|---|---|
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 <TFile.h> |
8 | #include <TTree.h> | |
609d6d39 | 9 | #include <TRandom3.h> |
a1a695e5 | 10 | |
0403120d | 11 | #include <AliESDtrackCuts.h> |
12 | #include <AliESDtrack.h> | |
13 | #include <AliTPCLaserTrack.h> | |
14 | ||
de0014b7 | 15 | #include "AliToyMCEvent.h" |
526ddf0e | 16 | |
a1a695e5 | 17 | #include "AliToyMCEventGeneratorSimple.h" |
18 | ||
cdf56507 | 19 | /* |
20 | ||
21 | ||
22 | */ | |
f356075c | 23 | |
de0014b7 | 24 | ClassImp(AliToyMCEventGeneratorSimple); |
526ddf0e | 25 | |
26 | ||
de0014b7 | 27 | AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple() |
28 | :AliToyMCEventGenerator() | |
526ddf0e | 29 | ,fVertexMean(0.) |
30 | ,fVertexSigma(7.) | |
32438f4e | 31 | ,fNtracks(20) |
f356075c | 32 | ,fInputFileNameESD("") |
33 | ,fESDCuts(0x0) | |
29b7f480 | 34 | ,fInputIndex(0) |
35 | ,fESDEvent(0x0) | |
36 | ,fESDTree(0x0) | |
37 | ,fInputFile(0x0) | |
1e62e876 | 38 | ,fHPt(0x0) |
39 | ,fHEta(0x0) | |
40 | ,fHMult(0x0) | |
41 | ,fHistosSet(kFALSE) | |
42 | ,fParamFile(0x0) | |
526ddf0e | 43 | { |
44 | // | |
45 | ||
46 | } | |
47 | //________________________________________________________________ | |
de0014b7 | 48 | AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple(const AliToyMCEventGeneratorSimple &gen) |
49 | :AliToyMCEventGenerator(gen) | |
526ddf0e | 50 | ,fVertexMean(gen.fVertexMean) |
51 | ,fVertexSigma(gen.fVertexSigma) | |
32438f4e | 52 | ,fNtracks(gen.fNtracks) |
f356075c | 53 | ,fInputFileNameESD(gen.fInputFileNameESD) |
54 | ,fESDCuts(gen.fESDCuts) | |
29b7f480 | 55 | ,fInputIndex(gen.fInputIndex) |
56 | ,fESDEvent(gen.fESDEvent) | |
57 | ,fESDTree(gen.fESDTree) | |
58 | ,fInputFile(gen.fInputFile) | |
1e62e876 | 59 | ,fHPt(gen.fHPt) |
60 | ,fHEta(gen.fHEta) | |
61 | ,fHMult(gen.fHMult) | |
62 | ,fHistosSet(gen.fHistosSet) | |
63 | ,fParamFile(0x0) | |
526ddf0e | 64 | { |
65 | } | |
66 | //________________________________________________________________ | |
de0014b7 | 67 | AliToyMCEventGeneratorSimple::~AliToyMCEventGeneratorSimple() |
526ddf0e | 68 | { |
69 | } | |
70 | //________________________________________________________________ | |
de0014b7 | 71 | AliToyMCEventGeneratorSimple& AliToyMCEventGeneratorSimple::operator = (const AliToyMCEventGeneratorSimple &gen) |
526ddf0e | 72 | { |
73 | //assignment operator | |
74 | if (&gen == this) return *this; | |
de0014b7 | 75 | new (this) AliToyMCEventGeneratorSimple(gen); |
526ddf0e | 76 | |
77 | return *this; | |
78 | } | |
79 | //________________________________________________________________ | |
1f195c60 | 80 | void AliToyMCEventGeneratorSimple::SetParametersToyGen(const Char_t* parfilename/*="files/params.root*/, Double_t vertexMean/*=0*/, Double_t vertexSigma/*=7.*/) { |
526ddf0e | 81 | fVertexMean = vertexMean; |
82 | fVertexSigma = vertexSigma; | |
1e62e876 | 83 | fParamFile = new TFile(parfilename, "read"); |
84 | fHPt = (TH1F*) fParamFile->Get("hPt"); | |
85 | fHEta = (TH1F*) fParamFile->Get("hEta"); | |
86 | fHMult = (TH1I*) fParamFile->Get("hMult") ; | |
87 | fHistosSet = kTRUE; | |
88 | ||
89 | ||
526ddf0e | 90 | } |
91 | //________________________________________________________________ | |
cd8ed0ac | 92 | AliToyMCEvent* AliToyMCEventGeneratorSimple::Generate(Double_t time) |
93 | { | |
94 | // | |
95 | // | |
96 | // | |
97 | ||
de0014b7 | 98 | AliToyMCEvent *retEvent = new AliToyMCEvent(); |
526ddf0e | 99 | retEvent->SetT0(time); |
100 | retEvent->SetX(0); | |
101 | retEvent->SetX(0); | |
e83fd282 | 102 | Double_t vertexZ=0; |
103 | //limit vertex to +- 10cm | |
104 | while ( TMath::Abs(vertexZ=gRandom->Gaus(fVertexMean,fVertexSigma))>10. ); | |
105 | retEvent->SetZ(vertexZ); | |
526ddf0e | 106 | |
e83fd282 | 107 | Double_t etaCuts=0.9; |
526ddf0e | 108 | Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass(); |
e83fd282 | 109 | static TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.3,100); |
0403120d | 110 | if (fpt.GetParameter(0)<1){ |
111 | printf("Set Parameters\n"); | |
112 | fpt.SetParameters(7.24,0.120); | |
113 | fpt.SetNpx(200); | |
114 | } | |
32438f4e | 115 | // Int_t nTracks = 400; //TODO: draw from experim dist |
116 | // Int_t nTracks = 20; //TODO: draw from experim dist | |
1e62e876 | 117 | Int_t nTracksLocal = fNtracks; |
118 | if(fHistosSet) { | |
119 | nTracksLocal = fHMult->GetRandom(); | |
120 | if(nTracksLocal > fNtracks) nTracksLocal = fNtracks; | |
121 | } | |
122 | std::cout << "Generating " << nTracksLocal << " tracks " << std::endl; | |
223d9e38 | 123 | |
1e62e876 | 124 | for(Int_t iTrack=0; iTrack<nTracksLocal; iTrack++){ |
526ddf0e | 125 | Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi()); |
1e62e876 | 126 | Double_t eta = 0.; |
127 | Double_t pt = 0.; | |
128 | ||
129 | if(fHistosSet) {//draw from histograms if set | |
130 | eta = fHEta->GetRandom(); | |
131 | pt = fHPt->GetRandom(); | |
132 | } | |
133 | else { | |
134 | eta = gRandom->Uniform(-etaCuts, etaCuts); | |
135 | pt = fpt.GetRandom(); // momentum for f1 | |
136 | } | |
526ddf0e | 137 | Short_t sign=1; |
138 | if(gRandom->Rndm() < 0.5){ | |
139 | sign =1; | |
140 | }else{ | |
141 | sign=-1; | |
142 | } | |
143 | ||
144 | Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.; | |
145 | Double_t pxyz[3]; | |
146 | pxyz[0]=pt*TMath::Cos(phi); | |
147 | pxyz[1]=pt*TMath::Sin(phi); | |
148 | pxyz[2]=pt*TMath::Tan(theta); | |
149 | Double_t vertex[3]={0,0,retEvent->GetZ()}; | |
150 | Double_t cv[21]={0}; | |
0403120d | 151 | AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign); |
7d26e3ce | 152 | // use unique ID for track number |
153 | // this will be used in DistortTrack track to set the cluster label | |
223d9e38 | 154 | // in one simulation the track id should be unique for performance studies |
155 | tempTrack->SetUniqueID(fCurrentTrack++); | |
0403120d | 156 | DistortTrack(*tempTrack, time); |
526ddf0e | 157 | } |
158 | ||
159 | return retEvent; | |
160 | } | |
161 | ||
a1a695e5 | 162 | //________________________________________________________________ |
fa97f7bd | 163 | void AliToyMCEventGeneratorSimple::RunSimulation(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/, const Int_t rate/*=50*/) |
a1a695e5 | 164 | { |
165 | // | |
166 | // run simple simulation with equal event spacing | |
fa97f7bd | 167 | // rate in kHz |
a1a695e5 | 168 | // |
169 | ||
170 | if (!ConnectOutputFile()) return; | |
cd8ed0ac | 171 | //initialise the space charge. Should be done after the tree was set up |
172 | InitSpaceCharge(); | |
223d9e38 | 173 | |
174 | // number of tracks to simulate per interaction | |
32438f4e | 175 | fNtracks=ntracks; |
223d9e38 | 176 | // within one simulation the track count should be unique for effeciency studies |
a06336b6 | 177 | // don't use the track ID 0 since the label -0 is not different from 0 |
178 | fCurrentTrack=1; | |
223d9e38 | 179 | |
a1a695e5 | 180 | Double_t eventTime=0.; |
fa97f7bd | 181 | const Double_t eventSpacing=1./rate/1e3; //50kHz equally spaced |
32438f4e | 182 | TStopwatch s; |
a1a695e5 | 183 | for (Int_t ievent=0; ievent<nevents; ++ievent){ |
32438f4e | 184 | printf("Generating event %3d (%.3g)\n",ievent,eventTime); |
a1a695e5 | 185 | fEvent = Generate(eventTime); |
186 | FillTree(); | |
187 | delete fEvent; | |
188 | fEvent=0x0; | |
189 | eventTime+=eventSpacing; | |
190 | } | |
32438f4e | 191 | s.Stop(); |
192 | s.Print(); | |
193 | ||
a1a695e5 | 194 | CloseOutputFile(); |
195 | } | |
526ddf0e | 196 | |
0403120d | 197 | //________________________________________________________________ |
198 | void AliToyMCEventGeneratorSimple::RunSimulationLaser(const Int_t nevents/*=1*/) | |
199 | { | |
200 | // | |
201 | // run simple simulation with equal event spacing | |
202 | // | |
203 | ||
204 | if (!ConnectOutputFile()) return; | |
205 | //initialise the space charge. Should be done after the tree was set up | |
206 | InitSpaceCharge(); | |
207 | ||
208 | // within one simulation the track count should be unique for effeciency studies | |
a06336b6 | 209 | fCurrentTrack=1; |
0403120d | 210 | |
211 | Double_t eventTime=0.; | |
212 | const Double_t eventSpacing=1./10.; //laser is running at 10Hz equally spaced | |
213 | TStopwatch s; | |
214 | for (Int_t ievent=0; ievent<nevents; ++ievent){ | |
215 | printf("Generating event %3d (%.3g)\n",ievent,eventTime); | |
216 | fEvent = GenerateLaser(eventTime); | |
217 | FillTree(); | |
218 | delete fEvent; | |
219 | fEvent=0x0; | |
220 | eventTime+=eventSpacing; | |
221 | } | |
222 | s.Stop(); | |
223 | s.Print(); | |
224 | ||
225 | CloseOutputFile(); | |
226 | } | |
227 | ||
f356075c | 228 | //________________________________________________________________ |
229 | AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD(AliESDEvent &esdEvent, Double_t time) { | |
230 | ||
231 | ||
232 | AliToyMCEvent *retEvent = new AliToyMCEvent(); | |
233 | retEvent->SetT0(time); | |
234 | retEvent->SetX(esdEvent.GetPrimaryVertex()->GetX()); | |
235 | retEvent->SetY(esdEvent.GetPrimaryVertex()->GetY()); | |
236 | retEvent->SetZ(esdEvent.GetPrimaryVertex()->GetZ()); | |
237 | ||
238 | ||
239 | ||
240 | if(!fNtracks) fNtracks = esdEvent.GetNumberOfTracks(); | |
241 | Int_t nUsedTracks = 0; | |
242 | for(Int_t iTrack=0; iTrack<esdEvent.GetNumberOfTracks(); iTrack++){ | |
243 | ||
244 | AliESDtrack *part = esdEvent.GetTrack(iTrack); | |
245 | if(!part) continue; | |
246 | if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue; | |
609d6d39 | 247 | if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume |
f356075c | 248 | Double_t pxyz[3]; |
249 | pxyz[0]=part->Px(); | |
250 | pxyz[1]=part->Py(); | |
251 | pxyz[2]=part->Pz(); | |
252 | Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()}; | |
253 | ||
254 | Double_t cv[21]={0}; | |
255 | Int_t sign = part->Charge(); | |
0403120d | 256 | AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign); |
7d26e3ce | 257 | // use unique ID for track number |
258 | // this will be used in DistortTrack track to set the cluster label | |
223d9e38 | 259 | // in one simulation the track id should be unique for performance studies |
260 | tempTrack->SetUniqueID(fCurrentTrack++); | |
0403120d | 261 | DistortTrack(*tempTrack, time); |
262 | nUsedTracks++; | |
f356075c | 263 | |
264 | if(nUsedTracks >= fNtracks) break; | |
265 | } | |
266 | ||
267 | return retEvent; | |
268 | } | |
269 | //________________________________________________________________ | |
270 | void AliToyMCEventGeneratorSimple::RunSimulationESD(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/) | |
271 | { | |
272 | // | |
273 | // run simulation using esd input with equal event spacing | |
274 | // | |
275 | ||
276 | if (!ConnectOutputFile()) return; | |
277 | TFile f(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root")); | |
278 | if(!f.IsOpen()) { | |
279 | std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl; | |
280 | return; | |
281 | } | |
282 | TTree* esdTree = (TTree*) f.Get("esdTree"); | |
283 | if(!esdTree) { | |
284 | std::cout <<"no esdTree in file " << std::endl; | |
285 | return; | |
286 | } | |
287 | InitSpaceCharge(); | |
288 | AliESDEvent* esdEvent = new AliESDEvent(); | |
29b7f480 | 289 | |
f356075c | 290 | |
291 | esdEvent->ReadFromTree(esdTree); | |
29b7f480 | 292 | |
223d9e38 | 293 | fNtracks = ntracks; |
294 | // within one simulation the track count should be unique for effeciency studies | |
a06336b6 | 295 | fCurrentTrack=1; |
223d9e38 | 296 | |
f356075c | 297 | Double_t eventTime=0.; |
298 | const Double_t eventSpacing=1./50e3; //50kHz equally spaced | |
299 | TStopwatch s; | |
29b7f480 | 300 | |
f356075c | 301 | fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1); |
1e62e876 | 302 | gRandom->SetSeed(); |
f356075c | 303 | Int_t nUsedEvents = 0; |
304 | for (Int_t ievent=0; ievent<esdTree->GetEntries(); ++ievent){ | |
609d6d39 | 305 | printf("Generating event %3d (%.3g)\n",nUsedEvents,eventTime); |
f356075c | 306 | esdTree->GetEvent(ievent); |
307 | if(esdEvent->GetNumberOfTracks()==0) { | |
308 | std::cout << " tracks == 0" << std::endl; | |
309 | continue; | |
310 | } | |
311 | ||
312 | fEvent = GenerateESD(*esdEvent, eventTime); | |
609d6d39 | 313 | if(fEvent->GetNumberOfTracks() >=10) { |
314 | nUsedEvents++; | |
315 | FillTree(); | |
316 | eventTime+=eventSpacing; | |
317 | } | |
f356075c | 318 | delete fEvent; |
319 | fEvent=0x0; | |
609d6d39 | 320 | |
f356075c | 321 | if(nUsedEvents>=nevents) break; |
322 | } | |
323 | s.Stop(); | |
324 | s.Print(); | |
325 | ||
326 | CloseOutputFile(); | |
327 | } | |
328 | ||
609d6d39 | 329 | //________________________________________________________________ |
609d6d39 | 330 | void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/) |
331 | { | |
332 | // | |
333 | // run simple simulation with equal event spacing | |
334 | // | |
335 | ||
336 | //Parameters for bunc | |
337 | const Double_t abortGap = 3e-6; // | |
338 | const Double_t collFreq = 50e3; | |
339 | const Double_t bSpacing = 50e-9; //bunch spacing | |
340 | const Int_t nTrainBunches = 48; | |
341 | const Int_t nTrains = 12; | |
342 | const Double_t revFreq = 1.11e4; //revolution frequency | |
343 | const Double_t collProb = collFreq/(nTrainBunches*nTrains*revFreq); | |
344 | const Double_t trainLength = bSpacing*(nTrainBunches-1); | |
345 | const Double_t totTrainLength = nTrains*trainLength; | |
346 | const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1); | |
347 | Bool_t equalSpacing = kFALSE; | |
348 | ||
349 | TRandom3 *rand = new TRandom3(); | |
350 | //rand->SetSeed(); | |
351 | ||
352 | if (!ConnectOutputFile()) return; | |
353 | //initialise the space charge. Should be done after the tree was set up | |
354 | InitSpaceCharge(); | |
355 | ||
356 | fNtracks=ntracks; | |
223d9e38 | 357 | // within one simulation the track count should be unique for effeciency studies |
a06336b6 | 358 | fCurrentTrack=1; |
223d9e38 | 359 | |
609d6d39 | 360 | Double_t eventTime=0.; |
29b7f480 | 361 | // const Double_t eventSpacing=1./50e3; //50kHz equally spaced |
609d6d39 | 362 | TStopwatch s; |
363 | Int_t nGeneratedEvents = 0; | |
364 | Int_t bunchCounter = 0; | |
365 | Int_t trainCounter = 0; | |
366 | ||
367 | //for (Int_t ievent=0; ievent<nevents; ++ievent){ | |
368 | while (nGeneratedEvents<nevents){ | |
369 | // std::cout <<trainCounter << " " << bunchCounter << " "<< "eventTime " << eventTime << std::endl; | |
370 | ||
371 | if(equalSpacing) { | |
372 | printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime); | |
373 | fEvent = Generate(eventTime); | |
374 | nGeneratedEvents++; | |
375 | FillTree(); | |
376 | delete fEvent; | |
377 | fEvent=0x0; | |
378 | eventTime+=1./collFreq; | |
379 | } | |
380 | else{ | |
381 | Int_t nCollsInCrossing = rand -> Poisson(collProb); | |
382 | for(Int_t iColl = 0; iColl<nCollsInCrossing; iColl++){ | |
383 | printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime); | |
384 | fEvent = Generate(eventTime); | |
385 | nGeneratedEvents++; | |
386 | FillTree(); | |
387 | delete fEvent; | |
388 | fEvent=0x0; | |
389 | ||
390 | } | |
391 | bunchCounter++; | |
392 | ||
393 | if(bunchCounter>=nTrainBunches){ | |
394 | ||
395 | trainCounter++; | |
396 | if(trainCounter>=nTrains){ | |
397 | eventTime+=abortGap; | |
398 | trainCounter=0; | |
399 | } | |
400 | else eventTime+=trainSpacing; | |
401 | ||
402 | bunchCounter=0; | |
403 | } | |
404 | else eventTime+= bSpacing; | |
405 | ||
406 | } | |
407 | ||
408 | ||
409 | } | |
410 | s.Stop(); | |
411 | s.Print(); | |
412 | delete rand; | |
413 | CloseOutputFile(); | |
414 | } | |
29b7f480 | 415 | |
416 | ||
417 | ||
418 | ||
419 | ||
420 | //________________________________________________________________ | |
421 | Int_t AliToyMCEventGeneratorSimple::OpenInputAndGetMaxEvents(const Int_t type, const Int_t nevents) { | |
422 | ||
423 | ||
424 | ||
425 | if(type==0) return nevents; | |
426 | ||
427 | if(type==1) { | |
428 | ||
429 | fInputFile = new TFile(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"),"read"); | |
430 | if(!fInputFile->IsOpen()) { | |
431 | std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl; | |
432 | return 0; | |
433 | } | |
434 | fESDTree = (TTree*) fInputFile->Get("esdTree"); | |
435 | if(!fESDTree) { | |
436 | std::cout <<"no esdTree in file " << std::endl; | |
437 | return 0; | |
438 | } | |
439 | fESDEvent = new AliESDEvent(); | |
440 | fESDEvent->ReadFromTree(fESDTree); | |
441 | fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1); | |
442 | ||
443 | fInputIndex = 0; | |
444 | ||
445 | return fESDTree->GetEntries(); | |
1e62e876 | 446 | gRandom->SetSeed(); |
29b7f480 | 447 | } |
448 | ||
449 | ||
450 | std::cout << " no available input type (toymc, esd) specified" << std::endl; | |
451 | return 0; | |
452 | ||
453 | ||
454 | } | |
455 | ||
456 | ||
457 | //________________________________________________________________ | |
458 | void AliToyMCEventGeneratorSimple::RunSimulation2(const Bool_t equalspacing, const Int_t type, const Int_t nevents, const Int_t ntracks) { | |
459 | ||
460 | //type==0 simple toy | |
461 | //type==1 esd input | |
462 | ||
463 | Int_t nMaxEvents = OpenInputAndGetMaxEvents(type,nevents); | |
464 | if (!ConnectOutputFile()) return; | |
465 | ||
223d9e38 | 466 | InitSpaceCharge(); |
29b7f480 | 467 | |
468 | ||
469 | fNtracks = ntracks; | |
223d9e38 | 470 | // within one simulation the track count should be unique for effeciency studies |
a06336b6 | 471 | fCurrentTrack=1; |
223d9e38 | 472 | |
29b7f480 | 473 | Double_t eventTime=0.; |
474 | TStopwatch s; | |
475 | // const Double_t eventSpacing=1./50e3; //50kHz equally spaced | |
476 | ||
477 | Int_t nGeneratedEvents = 0; | |
478 | ||
479 | for (Int_t ievent=0; ievent<nMaxEvents; ++ievent){ | |
480 | printf("Generating event %3d (%.3g)\n",ievent,eventTime); | |
481 | ||
482 | Int_t nColls = 0; | |
483 | Double_t spacing = 0; | |
484 | GetNGeneratedEventsAndSpacing(equalspacing, nColls, spacing); | |
485 | ||
486 | for(Int_t iColl = 0; iColl<nColls; iColl++){ | |
487 | if(type==0) fEvent = Generate(eventTime); | |
488 | else if(type==1) fEvent = GenerateESD2(eventTime); | |
489 | if(fEvent) { | |
490 | FillTree(); | |
491 | delete fEvent; | |
492 | fEvent=0x0; | |
493 | nGeneratedEvents++; | |
494 | } | |
495 | } | |
496 | eventTime+=spacing; | |
497 | if(nGeneratedEvents >= nevents) break; | |
498 | } | |
499 | s.Stop(); | |
500 | s.Print(); | |
501 | ||
502 | CloseOutputFile(); | |
503 | CloseInputFile(); | |
504 | ||
505 | ||
506 | } | |
507 | //________________________________________________________________ | |
508 | AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD2(Double_t time) { | |
509 | ||
510 | //test that enough tracks will pass cuts (and that there is tracks at all) | |
511 | Bool_t testEvent = kTRUE; | |
512 | while(fESDTree->GetEvent(fInputIndex) && testEvent) { | |
513 | ||
514 | Int_t nPassedCuts = 0; | |
515 | for(Int_t iTrack = 0; iTrack<fESDEvent->GetNumberOfTracks(); iTrack++){ | |
516 | if(fESDCuts->AcceptTrack(fESDEvent->GetTrack(iTrack)) && (fESDEvent->GetTrack(iTrack)->Pt() < 0.3) ) nPassedCuts++; | |
517 | } | |
518 | 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? | |
519 | else testEvent = kFALSE; | |
520 | } | |
521 | ||
522 | if(fInputIndex>=fESDTree->GetEntries()) return 0; | |
523 | ||
524 | ||
525 | ||
526 | AliToyMCEvent *retEvent = new AliToyMCEvent(); | |
527 | retEvent->SetT0(time); | |
528 | retEvent->SetX(fESDEvent->GetPrimaryVertex()->GetX()); | |
529 | retEvent->SetY(fESDEvent->GetPrimaryVertex()->GetY()); | |
530 | retEvent->SetZ(fESDEvent->GetPrimaryVertex()->GetZ()); | |
531 | ||
532 | ||
533 | ||
534 | if(!fNtracks) fNtracks = fESDEvent->GetNumberOfTracks(); | |
535 | Int_t nUsedTracks = 0; | |
536 | for(Int_t iTrack=0; iTrack<fESDEvent->GetNumberOfTracks(); iTrack++){ | |
537 | ||
538 | AliESDtrack *part = fESDEvent->GetTrack(iTrack); | |
539 | if(!part) continue; | |
540 | if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue; | |
541 | if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume | |
542 | ||
543 | Double_t pxyz[3]; | |
544 | pxyz[0]=part->Px(); | |
545 | pxyz[1]=part->Py(); | |
546 | pxyz[2]=part->Pz(); | |
547 | Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()}; | |
548 | Double_t cv[21]={0}; | |
549 | Int_t sign = part->Charge(); | |
550 | ||
0403120d | 551 | AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign); |
29b7f480 | 552 | // use unique ID for track number |
553 | // this will be used in DistortTrack track to set the cluster label | |
223d9e38 | 554 | // in one simulation the track id should be unique for performance studies |
555 | tempTrack->SetUniqueID(fCurrentTrack++); | |
0403120d | 556 | DistortTrack(*tempTrack, time); |
557 | ||
29b7f480 | 558 | |
559 | if(nUsedTracks >= fNtracks) break; | |
560 | } | |
561 | fInputIndex++; | |
562 | ||
563 | return retEvent; | |
564 | } | |
565 | ||
0403120d | 566 | //________________________________________________________________ |
567 | AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateLaser(Double_t time) | |
568 | { | |
569 | // | |
570 | // Generate an Event with laser tracks | |
571 | // | |
572 | ||
573 | AliToyMCEvent *retEvent = new AliToyMCEvent(); | |
574 | retEvent->SetEventType(AliToyMCEvent::kLaser); | |
575 | ||
576 | retEvent->SetT0(time); | |
577 | retEvent->SetX(0); | |
578 | retEvent->SetX(0); | |
579 | retEvent->SetZ(0); | |
580 | ||
581 | AliTPCLaserTrack::LoadTracks(); | |
582 | TObjArray *arr=AliTPCLaserTrack::GetTracks(); | |
583 | ||
584 | //since we have a laser track force no material budges | |
585 | Bool_t materialBudget=GetUseMaterialBudget(); | |
586 | SetUseMaterialBudget(kFALSE); | |
cedbc66e | 587 | |
588 | //the laser tracks have partially large inclination angles over the pads | |
589 | // -> relax the propagation contraint and switch off error messages | |
0403120d | 590 | SetIsLaser(kTRUE); |
cedbc66e | 591 | Int_t debug=AliLog::GetDebugLevel("","AliToyMCEventGeneratorSimple"); |
592 | AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",-3); | |
0403120d | 593 | |
594 | for (Int_t iTrack=0; iTrack<arr->GetEntriesFast(); ++iTrack){ | |
595 | AliExternalTrackParam *track=(AliExternalTrackParam*)arr->At(iTrack); | |
596 | AliToyMCTrack *tempTrack = retEvent->AddTrack(AliToyMCTrack(*track)); | |
597 | // for laser only TPC clusters exist | |
598 | tempTrack->SetUniqueID(fCurrentTrack++); | |
599 | MakeTPCClusters(*tempTrack, time); | |
600 | } | |
601 | ||
602 | SetIsLaser(kFALSE); | |
cedbc66e | 603 | AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",debug); |
0403120d | 604 | SetUseMaterialBudget(materialBudget); |
605 | return retEvent; | |
606 | } | |
607 | ||
29b7f480 | 608 | //________________________________________________________________ |
609 | void AliToyMCEventGeneratorSimple::GetNGeneratedEventsAndSpacing(const Bool_t equalSpacing, Int_t &ngen, Double_t &spacing) | |
610 | { | |
611 | ||
612 | static Int_t bunchCounter = 0; | |
613 | static Int_t trainCounter = 0; | |
614 | ||
615 | if(equalSpacing) { | |
616 | ngen =1; | |
617 | spacing = 1./50e3; //50kHz equally spaced | |
618 | return; | |
619 | } | |
620 | else if(!equalSpacing) { | |
621 | //parameters for bunch train | |
622 | const Double_t abortGap = 3e-6; // | |
623 | const Double_t collFreq = 50e3; | |
624 | const Double_t bSpacing = 50e-9; //bunch spacing | |
625 | const Int_t nTrainBunches = 48; | |
626 | const Int_t nTrains = 12; | |
627 | const Double_t revFreq = 1.11e4; //revolution frequency | |
628 | const Double_t collProb = collFreq/(nTrainBunches*nTrains*revFreq); | |
629 | const Double_t trainLength = bSpacing*(nTrainBunches-1); | |
630 | const Double_t totTrainLength = nTrains*trainLength; | |
631 | const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1); | |
632 | Double_t time = 0; | |
633 | Int_t nCollsInCrossing = 0; | |
634 | while(nCollsInCrossing ==0){ | |
635 | ||
636 | ||
637 | bunchCounter++; | |
638 | ||
639 | if(bunchCounter>=nTrainBunches){ | |
640 | ||
641 | trainCounter++; | |
642 | if(trainCounter>=nTrains){ | |
643 | time+=abortGap; | |
644 | trainCounter=0; | |
645 | } | |
646 | else time+=trainSpacing; | |
647 | ||
648 | bunchCounter=0; | |
649 | } | |
650 | else time+= bSpacing; | |
651 | ||
652 | ||
653 | nCollsInCrossing = gRandom -> Poisson(collProb); | |
654 | //std::cout << " nCollsInCrossing " << nCollsInCrossing << std::endl; | |
655 | } | |
656 | ngen = nCollsInCrossing; | |
657 | if(nCollsInCrossing > 1)std::cout << " nCollsInCrossing " << nCollsInCrossing << std::endl; | |
658 | spacing = time; | |
659 | return; | |
660 | ||
661 | } | |
662 | ||
663 | } | |
664 | ||
665 | //________________________________________________________________ | |
666 | Bool_t AliToyMCEventGeneratorSimple::CloseInputFile() | |
667 | { | |
668 | // | |
669 | // close the input file | |
670 | // | |
671 | if (!fInputFile) return kFALSE; | |
672 | fInputFile->Close(); | |
673 | delete fInputFile; | |
674 | fInputFile=0x0; | |
675 | ||
676 | return kTRUE; | |
677 | } |