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