3 #include <TDatabasePDG.h>
6 #include <TStopwatch.h>
11 #include <AliESDtrackCuts.h>
12 #include <AliESDtrack.h>
13 #include <AliTPCLaserTrack.h>
15 #include "AliToyMCEvent.h"
17 #include "AliToyMCEventGeneratorSimple.h"
20 ClassImp(AliToyMCEventGeneratorSimple);
23 AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple()
24 :AliToyMCEventGenerator()
28 ,fInputFileNameESD("")
43 //________________________________________________________________
44 AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple(const AliToyMCEventGeneratorSimple &gen)
45 :AliToyMCEventGenerator(gen)
46 ,fVertexMean(gen.fVertexMean)
47 ,fVertexSigma(gen.fVertexSigma)
48 ,fNtracks(gen.fNtracks)
49 ,fInputFileNameESD(gen.fInputFileNameESD)
50 ,fESDCuts(gen.fESDCuts)
51 ,fInputIndex(gen.fInputIndex)
52 ,fESDEvent(gen.fESDEvent)
53 ,fESDTree(gen.fESDTree)
54 ,fInputFile(gen.fInputFile)
58 ,fHistosSet(gen.fHistosSet)
62 //________________________________________________________________
63 AliToyMCEventGeneratorSimple::~AliToyMCEventGeneratorSimple()
66 //________________________________________________________________
67 AliToyMCEventGeneratorSimple& AliToyMCEventGeneratorSimple::operator = (const AliToyMCEventGeneratorSimple &gen)
70 if (&gen == this) return *this;
71 new (this) AliToyMCEventGeneratorSimple(gen);
75 //________________________________________________________________
76 void AliToyMCEventGeneratorSimple::SetParametersToyGen(const Char_t* parfilename/*="files/params.root*/, Double_t vertexMean/*=0*/, Double_t vertexSigma/*=7.*/) {
77 fVertexMean = vertexMean;
78 fVertexSigma = vertexSigma;
79 fParamFile = new TFile(parfilename, "read");
80 fHPt = (TH1F*) fParamFile->Get("hPt");
81 fHEta = (TH1F*) fParamFile->Get("hEta");
82 fHMult = (TH1I*) fParamFile->Get("hMult") ;
87 //________________________________________________________________
88 AliToyMCEvent* AliToyMCEventGeneratorSimple::Generate(Double_t time)
94 AliToyMCEvent *retEvent = new AliToyMCEvent();
95 retEvent->SetT0(time);
99 //limit vertex to +- 10cm
100 while ( TMath::Abs(vertexZ=gRandom->Gaus(fVertexMean,fVertexSigma))>10. );
101 retEvent->SetZ(vertexZ);
103 Double_t etaCuts=0.9;
104 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
105 static TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.3,100);
106 if (fpt.GetParameter(0)<1){
107 printf("Set Parameters\n");
108 fpt.SetParameters(7.24,0.120);
111 // Int_t nTracks = 400; //TODO: draw from experim dist
112 // Int_t nTracks = 20; //TODO: draw from experim dist
113 Int_t nTracksLocal = fNtracks;
115 nTracksLocal = fHMult->GetRandom();
116 if(nTracksLocal > fNtracks) nTracksLocal = fNtracks;
118 std::cout << "Generating " << nTracksLocal << " tracks " << std::endl;
120 for(Int_t iTrack=0; iTrack<nTracksLocal; iTrack++){
121 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
125 if(fHistosSet) {//draw from histograms if set
126 eta = fHEta->GetRandom();
127 pt = fHPt->GetRandom();
130 eta = gRandom->Uniform(-etaCuts, etaCuts);
131 pt = fpt.GetRandom(); // momentum for f1
134 if(gRandom->Rndm() < 0.5){
140 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
142 pxyz[0]=pt*TMath::Cos(phi);
143 pxyz[1]=pt*TMath::Sin(phi);
144 pxyz[2]=pt*TMath::Tan(theta);
145 Double_t vertex[3]={0,0,retEvent->GetZ()};
147 AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
148 // use unique ID for track number
149 // this will be used in DistortTrack track to set the cluster label
150 // in one simulation the track id should be unique for performance studies
151 tempTrack->SetUniqueID(fCurrentTrack++);
152 DistortTrack(*tempTrack, time);
158 //________________________________________________________________
159 void AliToyMCEventGeneratorSimple::RunSimulation(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/, const Int_t rate/*=50*/)
162 // run simple simulation with equal event spacing
166 if (!ConnectOutputFile()) return;
167 //initialise the space charge. Should be done after the tree was set up
170 // number of tracks to simulate per interaction
172 // within one simulation the track count should be unique for effeciency studies
173 // don't use the track ID 0 since the label -0 is not different from 0
176 Double_t eventTime=0.;
177 const Double_t eventSpacing=1./rate/1e3; //50kHz equally spaced
179 for (Int_t ievent=0; ievent<nevents; ++ievent){
180 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
181 fEvent = Generate(eventTime);
185 eventTime+=eventSpacing;
193 //________________________________________________________________
194 void AliToyMCEventGeneratorSimple::RunSimulationLaser(const Int_t nevents/*=1*/)
197 // run simple simulation with equal event spacing
200 if (!ConnectOutputFile()) return;
201 //initialise the space charge. Should be done after the tree was set up
204 // within one simulation the track count should be unique for effeciency studies
207 Double_t eventTime=0.;
208 const Double_t eventSpacing=1./10.; //laser is running at 10Hz equally spaced
210 for (Int_t ievent=0; ievent<nevents; ++ievent){
211 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
212 fEvent = GenerateLaser(eventTime);
216 eventTime+=eventSpacing;
224 //________________________________________________________________
225 AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD(AliESDEvent &esdEvent, Double_t time) {
228 AliToyMCEvent *retEvent = new AliToyMCEvent();
229 retEvent->SetT0(time);
230 retEvent->SetX(esdEvent.GetPrimaryVertex()->GetX());
231 retEvent->SetY(esdEvent.GetPrimaryVertex()->GetY());
232 retEvent->SetZ(esdEvent.GetPrimaryVertex()->GetZ());
236 if(!fNtracks) fNtracks = esdEvent.GetNumberOfTracks();
237 Int_t nUsedTracks = 0;
238 for(Int_t iTrack=0; iTrack<esdEvent.GetNumberOfTracks(); iTrack++){
240 AliESDtrack *part = esdEvent.GetTrack(iTrack);
242 if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue;
243 if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
248 Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
251 Int_t sign = part->Charge();
252 AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
253 // use unique ID for track number
254 // this will be used in DistortTrack track to set the cluster label
255 // in one simulation the track id should be unique for performance studies
256 tempTrack->SetUniqueID(fCurrentTrack++);
257 DistortTrack(*tempTrack, time);
260 if(nUsedTracks >= fNtracks) break;
265 //________________________________________________________________
266 void AliToyMCEventGeneratorSimple::RunSimulationESD(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
269 // run simulation using esd input with equal event spacing
272 if (!ConnectOutputFile()) return;
273 TFile f(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"));
275 std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
278 TTree* esdTree = (TTree*) f.Get("esdTree");
280 std::cout <<"no esdTree in file " << std::endl;
284 AliESDEvent* esdEvent = new AliESDEvent();
287 esdEvent->ReadFromTree(esdTree);
290 // within one simulation the track count should be unique for effeciency studies
293 Double_t eventTime=0.;
294 const Double_t eventSpacing=1./50e3; //50kHz equally spaced
297 fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
299 Int_t nUsedEvents = 0;
300 for (Int_t ievent=0; ievent<esdTree->GetEntries(); ++ievent){
301 printf("Generating event %3d (%.3g)\n",nUsedEvents,eventTime);
302 esdTree->GetEvent(ievent);
303 if(esdEvent->GetNumberOfTracks()==0) {
304 std::cout << " tracks == 0" << std::endl;
308 fEvent = GenerateESD(*esdEvent, eventTime);
309 if(fEvent->GetNumberOfTracks() >=10) {
312 eventTime+=eventSpacing;
317 if(nUsedEvents>=nevents) break;
325 //________________________________________________________________
326 void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
329 // run simple simulation with equal event spacing
332 //Parameters for bunc
333 const Double_t abortGap = 3e-6; //
334 const Double_t collFreq = 50e3;
335 const Double_t bSpacing = 50e-9; //bunch spacing
336 const Int_t nTrainBunches = 48;
337 const Int_t nTrains = 12;
338 const Double_t revFreq = 1.11e4; //revolution frequency
339 const Double_t collProb = collFreq/(nTrainBunches*nTrains*revFreq);
340 const Double_t trainLength = bSpacing*(nTrainBunches-1);
341 const Double_t totTrainLength = nTrains*trainLength;
342 const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1);
343 Bool_t equalSpacing = kFALSE;
345 TRandom3 *rand = new TRandom3();
348 if (!ConnectOutputFile()) return;
349 //initialise the space charge. Should be done after the tree was set up
353 // within one simulation the track count should be unique for effeciency studies
356 Double_t eventTime=0.;
357 // const Double_t eventSpacing=1./50e3; //50kHz equally spaced
359 Int_t nGeneratedEvents = 0;
360 Int_t bunchCounter = 0;
361 Int_t trainCounter = 0;
363 //for (Int_t ievent=0; ievent<nevents; ++ievent){
364 while (nGeneratedEvents<nevents){
365 // std::cout <<trainCounter << " " << bunchCounter << " "<< "eventTime " << eventTime << std::endl;
368 printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
369 fEvent = Generate(eventTime);
374 eventTime+=1./collFreq;
377 Int_t nCollsInCrossing = rand -> Poisson(collProb);
378 for(Int_t iColl = 0; iColl<nCollsInCrossing; iColl++){
379 printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
380 fEvent = Generate(eventTime);
389 if(bunchCounter>=nTrainBunches){
392 if(trainCounter>=nTrains){
396 else eventTime+=trainSpacing;
400 else eventTime+= bSpacing;
416 //________________________________________________________________
417 Int_t AliToyMCEventGeneratorSimple::OpenInputAndGetMaxEvents(const Int_t type, const Int_t nevents) {
421 if(type==0) return nevents;
425 fInputFile = new TFile(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"),"read");
426 if(!fInputFile->IsOpen()) {
427 std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
430 fESDTree = (TTree*) fInputFile->Get("esdTree");
432 std::cout <<"no esdTree in file " << std::endl;
435 fESDEvent = new AliESDEvent();
436 fESDEvent->ReadFromTree(fESDTree);
437 fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
441 return fESDTree->GetEntries();
446 std::cout << " no available input type (toymc, esd) specified" << std::endl;
453 //________________________________________________________________
454 void AliToyMCEventGeneratorSimple::RunSimulation2(const Bool_t equalspacing, const Int_t type, const Int_t nevents, const Int_t ntracks) {
459 Int_t nMaxEvents = OpenInputAndGetMaxEvents(type,nevents);
460 if (!ConnectOutputFile()) return;
466 // within one simulation the track count should be unique for effeciency studies
469 Double_t eventTime=0.;
471 // const Double_t eventSpacing=1./50e3; //50kHz equally spaced
473 Int_t nGeneratedEvents = 0;
475 for (Int_t ievent=0; ievent<nMaxEvents; ++ievent){
476 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
479 Double_t spacing = 0;
480 GetNGeneratedEventsAndSpacing(equalspacing, nColls, spacing);
482 for(Int_t iColl = 0; iColl<nColls; iColl++){
483 if(type==0) fEvent = Generate(eventTime);
484 else if(type==1) fEvent = GenerateESD2(eventTime);
493 if(nGeneratedEvents >= nevents) break;
503 //________________________________________________________________
504 AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD2(Double_t time) {
506 //test that enough tracks will pass cuts (and that there is tracks at all)
507 Bool_t testEvent = kTRUE;
508 while(fESDTree->GetEvent(fInputIndex) && testEvent) {
510 Int_t nPassedCuts = 0;
511 for(Int_t iTrack = 0; iTrack<fESDEvent->GetNumberOfTracks(); iTrack++){
512 if(fESDCuts->AcceptTrack(fESDEvent->GetTrack(iTrack)) && (fESDEvent->GetTrack(iTrack)->Pt() < 0.3) ) nPassedCuts++;
514 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?
515 else testEvent = kFALSE;
518 if(fInputIndex>=fESDTree->GetEntries()) return 0;
522 AliToyMCEvent *retEvent = new AliToyMCEvent();
523 retEvent->SetT0(time);
524 retEvent->SetX(fESDEvent->GetPrimaryVertex()->GetX());
525 retEvent->SetY(fESDEvent->GetPrimaryVertex()->GetY());
526 retEvent->SetZ(fESDEvent->GetPrimaryVertex()->GetZ());
530 if(!fNtracks) fNtracks = fESDEvent->GetNumberOfTracks();
531 Int_t nUsedTracks = 0;
532 for(Int_t iTrack=0; iTrack<fESDEvent->GetNumberOfTracks(); iTrack++){
534 AliESDtrack *part = fESDEvent->GetTrack(iTrack);
536 if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue;
537 if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
543 Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
545 Int_t sign = part->Charge();
547 AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
548 // use unique ID for track number
549 // this will be used in DistortTrack track to set the cluster label
550 // in one simulation the track id should be unique for performance studies
551 tempTrack->SetUniqueID(fCurrentTrack++);
552 DistortTrack(*tempTrack, time);
555 if(nUsedTracks >= fNtracks) break;
562 //________________________________________________________________
563 AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateLaser(Double_t time)
566 // Generate an Event with laser tracks
569 AliToyMCEvent *retEvent = new AliToyMCEvent();
570 retEvent->SetEventType(AliToyMCEvent::kLaser);
572 retEvent->SetT0(time);
577 AliTPCLaserTrack::LoadTracks();
578 TObjArray *arr=AliTPCLaserTrack::GetTracks();
580 //since we have a laser track force no material budges
581 Bool_t materialBudget=GetUseMaterialBudget();
582 SetUseMaterialBudget(kFALSE);
584 //the laser tracks have partially large inclination angles over the pads
585 // -> relax the propagation contraint and switch off error messages
587 Int_t debug=AliLog::GetDebugLevel("","AliToyMCEventGeneratorSimple");
588 AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",-3);
590 for (Int_t iTrack=0; iTrack<arr->GetEntriesFast(); ++iTrack){
591 AliExternalTrackParam *track=(AliExternalTrackParam*)arr->At(iTrack);
592 AliToyMCTrack *tempTrack = retEvent->AddTrack(AliToyMCTrack(*track));
593 // for laser only TPC clusters exist
594 tempTrack->SetUniqueID(fCurrentTrack++);
595 MakeTPCClusters(*tempTrack, time);
599 AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",debug);
600 SetUseMaterialBudget(materialBudget);
604 //________________________________________________________________
605 void AliToyMCEventGeneratorSimple::GetNGeneratedEventsAndSpacing(const Bool_t equalSpacing, Int_t &ngen, Double_t &spacing)
608 static Int_t bunchCounter = 0;
609 static Int_t trainCounter = 0;
613 spacing = 1./50e3; //50kHz equally spaced
616 else if(!equalSpacing) {
617 //parameters for bunch train
618 const Double_t abortGap = 3e-6; //
619 const Double_t collFreq = 50e3;
620 const Double_t bSpacing = 50e-9; //bunch spacing
621 const Int_t nTrainBunches = 48;
622 const Int_t nTrains = 12;
623 const Double_t revFreq = 1.11e4; //revolution frequency
624 const Double_t collProb = collFreq/(nTrainBunches*nTrains*revFreq);
625 const Double_t trainLength = bSpacing*(nTrainBunches-1);
626 const Double_t totTrainLength = nTrains*trainLength;
627 const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1);
629 Int_t nCollsInCrossing = 0;
630 while(nCollsInCrossing ==0){
635 if(bunchCounter>=nTrainBunches){
638 if(trainCounter>=nTrains){
642 else time+=trainSpacing;
646 else time+= bSpacing;
649 nCollsInCrossing = gRandom -> Poisson(collProb);
650 //std::cout << " nCollsInCrossing " << nCollsInCrossing << std::endl;
652 ngen = nCollsInCrossing;
653 if(nCollsInCrossing > 1)std::cout << " nCollsInCrossing " << nCollsInCrossing << std::endl;
661 //________________________________________________________________
662 Bool_t AliToyMCEventGeneratorSimple::CloseInputFile()
665 // close the input file
667 if (!fInputFile) return kFALSE;