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