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