bugfix: use time, when looking for next cluster, exception: kIdealTracking
[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)
29b7f480 28 ,fInputIndex(0)
29 ,fESDEvent(0x0)
30 ,fESDTree(0x0)
31 ,fInputFile(0x0)
32
526ddf0e 33{
34 //
35
36}
37//________________________________________________________________
de0014b7 38AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple(const AliToyMCEventGeneratorSimple &gen)
39 :AliToyMCEventGenerator(gen)
526ddf0e 40 ,fVertexMean(gen.fVertexMean)
41 ,fVertexSigma(gen.fVertexSigma)
32438f4e 42 ,fNtracks(gen.fNtracks)
f356075c 43 ,fInputFileNameESD(gen.fInputFileNameESD)
44 ,fESDCuts(gen.fESDCuts)
29b7f480 45 ,fInputIndex(gen.fInputIndex)
46 ,fESDEvent(gen.fESDEvent)
47 ,fESDTree(gen.fESDTree)
48 ,fInputFile(gen.fInputFile)
49
526ddf0e 50{
51}
52//________________________________________________________________
de0014b7 53AliToyMCEventGeneratorSimple::~AliToyMCEventGeneratorSimple()
526ddf0e 54{
55}
56 //________________________________________________________________
de0014b7 57AliToyMCEventGeneratorSimple& AliToyMCEventGeneratorSimple::operator = (const AliToyMCEventGeneratorSimple &gen)
526ddf0e 58{
59 //assignment operator
60 if (&gen == this) return *this;
de0014b7 61 new (this) AliToyMCEventGeneratorSimple(gen);
526ddf0e 62
63 return *this;
64}
65//________________________________________________________________
f356075c 66void AliToyMCEventGeneratorSimple::SetParametersSimple(Double_t vertexMean, Double_t vertexSigma) {
526ddf0e 67 fVertexMean = vertexMean;
68 fVertexSigma = vertexSigma;
69}
70//________________________________________________________________
cd8ed0ac 71AliToyMCEvent* AliToyMCEventGeneratorSimple::Generate(Double_t time)
72{
73 //
74 //
75 //
76
de0014b7 77 AliToyMCEvent *retEvent = new AliToyMCEvent();
526ddf0e 78 retEvent->SetT0(time);
79 retEvent->SetX(0);
80 retEvent->SetX(0);
81 retEvent->SetZ(gRandom->Gaus(fVertexMean,fVertexSigma));
82
83 Double_t etaCuts=.9;
84 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
85 TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
86 fpt.SetParameters(7.24,0.120);
87 fpt.SetNpx(10000);
32438f4e 88// Int_t nTracks = 400; //TODO: draw from experim dist
89// Int_t nTracks = 20; //TODO: draw from experim dist
526ddf0e 90
32438f4e 91 for(Int_t iTrack=0; iTrack<fNtracks; iTrack++){
526ddf0e 92 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
93 Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
94 Double_t pt = fpt.GetRandom(); // momentum for f1
95 Short_t sign=1;
96 if(gRandom->Rndm() < 0.5){
97 sign =1;
98 }else{
99 sign=-1;
100 }
101
102 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
103 Double_t pxyz[3];
104 pxyz[0]=pt*TMath::Cos(phi);
105 pxyz[1]=pt*TMath::Sin(phi);
106 pxyz[2]=pt*TMath::Tan(theta);
107 Double_t vertex[3]={0,0,retEvent->GetZ()};
108 Double_t cv[21]={0};
de0014b7 109 AliToyMCTrack *tempTrack = new AliToyMCTrack(vertex,pxyz,cv,sign);
7d26e3ce 110 // use unique ID for track number
111 // this will be used in DistortTrack track to set the cluster label
112 tempTrack->SetUniqueID(iTrack);
113
526ddf0e 114 Bool_t trackDist = DistortTrack(*tempTrack, time);
115 if(trackDist) retEvent->AddTrack(*tempTrack);
116 delete tempTrack;
117 }
118
119 return retEvent;
120}
121
a1a695e5 122//________________________________________________________________
609d6d39 123void AliToyMCEventGeneratorSimple::RunSimulation(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
a1a695e5 124{
125 //
126 // run simple simulation with equal event spacing
127 //
128
129 if (!ConnectOutputFile()) return;
cd8ed0ac 130 //initialise the space charge. Should be done after the tree was set up
131 InitSpaceCharge();
132
32438f4e 133 fNtracks=ntracks;
a1a695e5 134 Double_t eventTime=0.;
135 const Double_t eventSpacing=1./50e3; //50kHz equally spaced
32438f4e 136 TStopwatch s;
a1a695e5 137 for (Int_t ievent=0; ievent<nevents; ++ievent){
32438f4e 138 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
a1a695e5 139 fEvent = Generate(eventTime);
140 FillTree();
141 delete fEvent;
142 fEvent=0x0;
143 eventTime+=eventSpacing;
144 }
32438f4e 145 s.Stop();
146 s.Print();
147
a1a695e5 148 CloseOutputFile();
149}
526ddf0e 150
f356075c 151//________________________________________________________________
152AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD(AliESDEvent &esdEvent, Double_t time) {
153
154
155 AliToyMCEvent *retEvent = new AliToyMCEvent();
156 retEvent->SetT0(time);
157 retEvent->SetX(esdEvent.GetPrimaryVertex()->GetX());
158 retEvent->SetY(esdEvent.GetPrimaryVertex()->GetY());
159 retEvent->SetZ(esdEvent.GetPrimaryVertex()->GetZ());
160
161
162
163 if(!fNtracks) fNtracks = esdEvent.GetNumberOfTracks();
164 Int_t nUsedTracks = 0;
165 for(Int_t iTrack=0; iTrack<esdEvent.GetNumberOfTracks(); iTrack++){
166
167 AliESDtrack *part = esdEvent.GetTrack(iTrack);
168 if(!part) continue;
169 if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue;
609d6d39 170 if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
f356075c 171 Double_t pxyz[3];
172 pxyz[0]=part->Px();
173 pxyz[1]=part->Py();
174 pxyz[2]=part->Pz();
175 Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
176
177 Double_t cv[21]={0};
178 Int_t sign = part->Charge();
179 AliToyMCTrack *tempTrack = new AliToyMCTrack(vertex,pxyz,cv,sign);
7d26e3ce 180 // use unique ID for track number
181 // this will be used in DistortTrack track to set the cluster label
182 tempTrack->SetUniqueID(iTrack);
183
184
f356075c 185 Bool_t trackDist = DistortTrack(*tempTrack, time);
186 if(trackDist) {
187 retEvent->AddTrack(*tempTrack);
609d6d39 188
f356075c 189 nUsedTracks++;
190 }
191 delete tempTrack;
192
193 if(nUsedTracks >= fNtracks) break;
194 }
195
196 return retEvent;
197}
198//________________________________________________________________
199void AliToyMCEventGeneratorSimple::RunSimulationESD(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
200{
201 //
202 // run simulation using esd input with equal event spacing
203 //
204
205 if (!ConnectOutputFile()) return;
206 TFile f(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"));
207 if(!f.IsOpen()) {
208 std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
209 return;
210 }
211 TTree* esdTree = (TTree*) f.Get("esdTree");
212 if(!esdTree) {
213 std::cout <<"no esdTree in file " << std::endl;
214 return;
215 }
216 InitSpaceCharge();
217 AliESDEvent* esdEvent = new AliESDEvent();
29b7f480 218
f356075c 219
220 esdEvent->ReadFromTree(esdTree);
29b7f480 221
222 fNtracks = ntracks;
f356075c 223 Double_t eventTime=0.;
224 const Double_t eventSpacing=1./50e3; //50kHz equally spaced
225 TStopwatch s;
29b7f480 226
f356075c 227 fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
29b7f480 228
f356075c 229 Int_t nUsedEvents = 0;
230 for (Int_t ievent=0; ievent<esdTree->GetEntries(); ++ievent){
609d6d39 231 printf("Generating event %3d (%.3g)\n",nUsedEvents,eventTime);
f356075c 232 esdTree->GetEvent(ievent);
233 if(esdEvent->GetNumberOfTracks()==0) {
234 std::cout << " tracks == 0" << std::endl;
235 continue;
236 }
237
238 fEvent = GenerateESD(*esdEvent, eventTime);
609d6d39 239 if(fEvent->GetNumberOfTracks() >=10) {
240 nUsedEvents++;
241 FillTree();
242 eventTime+=eventSpacing;
243 }
f356075c 244 delete fEvent;
245 fEvent=0x0;
609d6d39 246
f356075c 247 if(nUsedEvents>=nevents) break;
248 }
249 s.Stop();
250 s.Print();
251
252 CloseOutputFile();
253}
254
609d6d39 255//________________________________________________________________
256
257void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
258{
259 //
260 // run simple simulation with equal event spacing
261 //
262
263 //Parameters for bunc
264 const Double_t abortGap = 3e-6; //
265 const Double_t collFreq = 50e3;
266 const Double_t bSpacing = 50e-9; //bunch spacing
267 const Int_t nTrainBunches = 48;
268 const Int_t nTrains = 12;
269 const Double_t revFreq = 1.11e4; //revolution frequency
270 const Double_t collProb = collFreq/(nTrainBunches*nTrains*revFreq);
271 const Double_t trainLength = bSpacing*(nTrainBunches-1);
272 const Double_t totTrainLength = nTrains*trainLength;
273 const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1);
274 Bool_t equalSpacing = kFALSE;
275
276 TRandom3 *rand = new TRandom3();
277 //rand->SetSeed();
278
279 if (!ConnectOutputFile()) return;
280 //initialise the space charge. Should be done after the tree was set up
281 InitSpaceCharge();
282
283 fNtracks=ntracks;
284 Double_t eventTime=0.;
29b7f480 285 // const Double_t eventSpacing=1./50e3; //50kHz equally spaced
609d6d39 286 TStopwatch s;
287 Int_t nGeneratedEvents = 0;
288 Int_t bunchCounter = 0;
289 Int_t trainCounter = 0;
290
291 //for (Int_t ievent=0; ievent<nevents; ++ievent){
292 while (nGeneratedEvents<nevents){
293 // std::cout <<trainCounter << " " << bunchCounter << " "<< "eventTime " << eventTime << std::endl;
294
295 if(equalSpacing) {
296 printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
297 fEvent = Generate(eventTime);
298 nGeneratedEvents++;
299 FillTree();
300 delete fEvent;
301 fEvent=0x0;
302 eventTime+=1./collFreq;
303 }
304 else{
305 Int_t nCollsInCrossing = rand -> Poisson(collProb);
306 for(Int_t iColl = 0; iColl<nCollsInCrossing; iColl++){
307 printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
308 fEvent = Generate(eventTime);
309 nGeneratedEvents++;
310 FillTree();
311 delete fEvent;
312 fEvent=0x0;
313
314 }
315 bunchCounter++;
316
317 if(bunchCounter>=nTrainBunches){
318
319 trainCounter++;
320 if(trainCounter>=nTrains){
321 eventTime+=abortGap;
322 trainCounter=0;
323 }
324 else eventTime+=trainSpacing;
325
326 bunchCounter=0;
327 }
328 else eventTime+= bSpacing;
329
330 }
331
332
333 }
334 s.Stop();
335 s.Print();
336 delete rand;
337 CloseOutputFile();
338}
29b7f480 339
340
341
342
343
344//________________________________________________________________
345Int_t AliToyMCEventGeneratorSimple::OpenInputAndGetMaxEvents(const Int_t type, const Int_t nevents) {
346
347
348
349 if(type==0) return nevents;
350
351 if(type==1) {
352
353 fInputFile = new TFile(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"),"read");
354 if(!fInputFile->IsOpen()) {
355 std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
356 return 0;
357 }
358 fESDTree = (TTree*) fInputFile->Get("esdTree");
359 if(!fESDTree) {
360 std::cout <<"no esdTree in file " << std::endl;
361 return 0;
362 }
363 fESDEvent = new AliESDEvent();
364 fESDEvent->ReadFromTree(fESDTree);
365 fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
366
367 fInputIndex = 0;
368
369 return fESDTree->GetEntries();
370 }
371
372
373 std::cout << " no available input type (toymc, esd) specified" << std::endl;
374 return 0;
375
376
377}
378
379
380//________________________________________________________________
381void AliToyMCEventGeneratorSimple::RunSimulation2(const Bool_t equalspacing, const Int_t type, const Int_t nevents, const Int_t ntracks) {
382
383 //type==0 simple toy
384 //type==1 esd input
385
386 Int_t nMaxEvents = OpenInputAndGetMaxEvents(type,nevents);
387 if (!ConnectOutputFile()) return;
388
389 InitSpaceCharge();
390
391
392
393
394 fNtracks = ntracks;
395
396 Double_t eventTime=0.;
397 TStopwatch s;
398 // const Double_t eventSpacing=1./50e3; //50kHz equally spaced
399
400 Int_t nGeneratedEvents = 0;
401
402 for (Int_t ievent=0; ievent<nMaxEvents; ++ievent){
403 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
404
405 Int_t nColls = 0;
406 Double_t spacing = 0;
407 GetNGeneratedEventsAndSpacing(equalspacing, nColls, spacing);
408
409 for(Int_t iColl = 0; iColl<nColls; iColl++){
410 if(type==0) fEvent = Generate(eventTime);
411 else if(type==1) fEvent = GenerateESD2(eventTime);
412 if(fEvent) {
413 FillTree();
414 delete fEvent;
415 fEvent=0x0;
416 nGeneratedEvents++;
417 }
418 }
419 eventTime+=spacing;
420 if(nGeneratedEvents >= nevents) break;
421 }
422 s.Stop();
423 s.Print();
424
425 CloseOutputFile();
426 CloseInputFile();
427
428
429}
430//________________________________________________________________
431AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD2(Double_t time) {
432
433 //test that enough tracks will pass cuts (and that there is tracks at all)
434 Bool_t testEvent = kTRUE;
435 while(fESDTree->GetEvent(fInputIndex) && testEvent) {
436
437 Int_t nPassedCuts = 0;
438 for(Int_t iTrack = 0; iTrack<fESDEvent->GetNumberOfTracks(); iTrack++){
439 if(fESDCuts->AcceptTrack(fESDEvent->GetTrack(iTrack)) && (fESDEvent->GetTrack(iTrack)->Pt() < 0.3) ) nPassedCuts++;
440 }
441 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?
442 else testEvent = kFALSE;
443 }
444
445 if(fInputIndex>=fESDTree->GetEntries()) return 0;
446
447
448
449 AliToyMCEvent *retEvent = new AliToyMCEvent();
450 retEvent->SetT0(time);
451 retEvent->SetX(fESDEvent->GetPrimaryVertex()->GetX());
452 retEvent->SetY(fESDEvent->GetPrimaryVertex()->GetY());
453 retEvent->SetZ(fESDEvent->GetPrimaryVertex()->GetZ());
454
455
456
457 if(!fNtracks) fNtracks = fESDEvent->GetNumberOfTracks();
458 Int_t nUsedTracks = 0;
459 for(Int_t iTrack=0; iTrack<fESDEvent->GetNumberOfTracks(); iTrack++){
460
461 AliESDtrack *part = fESDEvent->GetTrack(iTrack);
462 if(!part) continue;
463 if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue;
464 if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
465
466 Double_t pxyz[3];
467 pxyz[0]=part->Px();
468 pxyz[1]=part->Py();
469 pxyz[2]=part->Pz();
470 Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
471 Double_t cv[21]={0};
472 Int_t sign = part->Charge();
473
474 AliToyMCTrack *tempTrack = new AliToyMCTrack(vertex,pxyz,cv,sign);
475 // use unique ID for track number
476 // this will be used in DistortTrack track to set the cluster label
477 tempTrack->SetUniqueID(iTrack);
478
479
480 Bool_t trackDist = DistortTrack(*tempTrack, time);
481 if(trackDist) {
482 retEvent->AddTrack(*tempTrack);
483
484 nUsedTracks++;
485 }
486 delete tempTrack;
487
488 if(nUsedTracks >= fNtracks) break;
489 }
490 fInputIndex++;
491
492 return retEvent;
493}
494
495//________________________________________________________________
496void AliToyMCEventGeneratorSimple::GetNGeneratedEventsAndSpacing(const Bool_t equalSpacing, Int_t &ngen, Double_t &spacing)
497{
498
499 static Int_t bunchCounter = 0;
500 static Int_t trainCounter = 0;
501
502 if(equalSpacing) {
503 ngen =1;
504 spacing = 1./50e3; //50kHz equally spaced
505 return;
506 }
507 else if(!equalSpacing) {
508 //parameters for bunch train
509 const Double_t abortGap = 3e-6; //
510 const Double_t collFreq = 50e3;
511 const Double_t bSpacing = 50e-9; //bunch spacing
512 const Int_t nTrainBunches = 48;
513 const Int_t nTrains = 12;
514 const Double_t revFreq = 1.11e4; //revolution frequency
515 const Double_t collProb = collFreq/(nTrainBunches*nTrains*revFreq);
516 const Double_t trainLength = bSpacing*(nTrainBunches-1);
517 const Double_t totTrainLength = nTrains*trainLength;
518 const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1);
519 Double_t time = 0;
520 Int_t nCollsInCrossing = 0;
521 while(nCollsInCrossing ==0){
522
523
524 bunchCounter++;
525
526 if(bunchCounter>=nTrainBunches){
527
528 trainCounter++;
529 if(trainCounter>=nTrains){
530 time+=abortGap;
531 trainCounter=0;
532 }
533 else time+=trainSpacing;
534
535 bunchCounter=0;
536 }
537 else time+= bSpacing;
538
539
540 nCollsInCrossing = gRandom -> Poisson(collProb);
541 //std::cout << " nCollsInCrossing " << nCollsInCrossing << std::endl;
542 }
543 ngen = nCollsInCrossing;
544 if(nCollsInCrossing > 1)std::cout << " nCollsInCrossing " << nCollsInCrossing << std::endl;
545 spacing = time;
546 return;
547
548 }
549
550}
551
552//________________________________________________________________
553Bool_t AliToyMCEventGeneratorSimple::CloseInputFile()
554{
555 //
556 // close the input file
557 //
558 if (!fInputFile) return kFALSE;
559 fInputFile->Close();
560 delete fInputFile;
561 fInputFile=0x0;
562
563 return kTRUE;
564}