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"
24 ClassImp(AliToyMCEventGeneratorSimple);
27 AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple()
28 :AliToyMCEventGenerator()
32 ,fInputFileNameESD("")
47 //________________________________________________________________
48 AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple(const AliToyMCEventGeneratorSimple &gen)
49 :AliToyMCEventGenerator(gen)
50 ,fVertexMean(gen.fVertexMean)
51 ,fVertexSigma(gen.fVertexSigma)
52 ,fNtracks(gen.fNtracks)
53 ,fInputFileNameESD(gen.fInputFileNameESD)
54 ,fESDCuts(gen.fESDCuts)
55 ,fInputIndex(gen.fInputIndex)
56 ,fESDEvent(gen.fESDEvent)
57 ,fESDTree(gen.fESDTree)
58 ,fInputFile(gen.fInputFile)
62 ,fHistosSet(gen.fHistosSet)
66 //________________________________________________________________
67 AliToyMCEventGeneratorSimple::~AliToyMCEventGeneratorSimple()
70 //________________________________________________________________
71 AliToyMCEventGeneratorSimple& AliToyMCEventGeneratorSimple::operator = (const AliToyMCEventGeneratorSimple &gen)
74 if (&gen == this) return *this;
75 new (this) AliToyMCEventGeneratorSimple(gen);
79 //________________________________________________________________
80 void AliToyMCEventGeneratorSimple::SetParametersToyGen(const Char_t* parfilename/*="files/params.root*/, Double_t vertexMean/*=0*/, Double_t vertexSigma/*=7.*/) {
81 fVertexMean = vertexMean;
82 fVertexSigma = vertexSigma;
83 fParamFile = new TFile(parfilename, "read");
84 fHPt = (TH1F*) fParamFile->Get("hPt");
85 fHEta = (TH1F*) fParamFile->Get("hEta");
86 fHMult = (TH1I*) fParamFile->Get("hMult") ;
91 //________________________________________________________________
92 AliToyMCEvent* AliToyMCEventGeneratorSimple::Generate(Double_t time)
98 AliToyMCEvent *retEvent = new AliToyMCEvent();
99 retEvent->SetT0(time);
103 //limit vertex to +- 10cm
104 while ( TMath::Abs(vertexZ=gRandom->Gaus(fVertexMean,fVertexSigma))>10. );
105 retEvent->SetZ(vertexZ);
107 Double_t etaCuts=0.9;
108 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
109 static TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.3,100);
110 if (fpt.GetParameter(0)<1){
111 printf("Set Parameters\n");
112 fpt.SetParameters(7.24,0.120);
115 // Int_t nTracks = 400; //TODO: draw from experim dist
116 // Int_t nTracks = 20; //TODO: draw from experim dist
117 Int_t nTracksLocal = fNtracks;
119 nTracksLocal = fHMult->GetRandom();
120 if(nTracksLocal > fNtracks) nTracksLocal = fNtracks;
122 std::cout << "Generating " << nTracksLocal << " tracks " << std::endl;
124 for(Int_t iTrack=0; iTrack<nTracksLocal; iTrack++){
125 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
129 if(fHistosSet) {//draw from histograms if set
130 eta = fHEta->GetRandom();
131 pt = fHPt->GetRandom();
134 eta = gRandom->Uniform(-etaCuts, etaCuts);
135 pt = fpt.GetRandom(); // momentum for f1
138 if(gRandom->Rndm() < 0.5){
144 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
146 pxyz[0]=pt*TMath::Cos(phi);
147 pxyz[1]=pt*TMath::Sin(phi);
148 pxyz[2]=pt*TMath::Tan(theta);
149 Double_t vertex[3]={0,0,retEvent->GetZ()};
151 AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
152 // use unique ID for track number
153 // this will be used in DistortTrack track to set the cluster label
154 // in one simulation the track id should be unique for performance studies
155 tempTrack->SetUniqueID(fCurrentTrack++);
156 DistortTrack(*tempTrack, time);
162 //________________________________________________________________
163 void AliToyMCEventGeneratorSimple::RunSimulation(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/, const Int_t rate/*=50*/)
166 // run simple simulation with equal event spacing
170 if (!ConnectOutputFile()) return;
171 //initialise the space charge. Should be done after the tree was set up
174 // number of tracks to simulate per interaction
176 // within one simulation the track count should be unique for effeciency studies
177 // don't use the track ID 0 since the label -0 is not different from 0
180 Double_t eventTime=0.;
181 const Double_t eventSpacing=1./rate/1e3; //50kHz equally spaced
183 for (Int_t ievent=0; ievent<nevents; ++ievent){
184 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
185 fEvent = Generate(eventTime);
189 eventTime+=eventSpacing;
197 //________________________________________________________________
198 void AliToyMCEventGeneratorSimple::RunSimulationLaser(const Int_t nevents/*=1*/)
201 // run simple simulation with equal event spacing
204 if (!ConnectOutputFile()) return;
205 //initialise the space charge. Should be done after the tree was set up
208 // within one simulation the track count should be unique for effeciency studies
211 Double_t eventTime=0.;
212 const Double_t eventSpacing=1./10.; //laser is running at 10Hz equally spaced
214 for (Int_t ievent=0; ievent<nevents; ++ievent){
215 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
216 fEvent = GenerateLaser(eventTime);
220 eventTime+=eventSpacing;
228 //________________________________________________________________
229 AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD(AliESDEvent &esdEvent, Double_t time) {
232 AliToyMCEvent *retEvent = new AliToyMCEvent();
233 retEvent->SetT0(time);
234 retEvent->SetX(esdEvent.GetPrimaryVertex()->GetX());
235 retEvent->SetY(esdEvent.GetPrimaryVertex()->GetY());
236 retEvent->SetZ(esdEvent.GetPrimaryVertex()->GetZ());
240 if(!fNtracks) fNtracks = esdEvent.GetNumberOfTracks();
241 Int_t nUsedTracks = 0;
242 for(Int_t iTrack=0; iTrack<esdEvent.GetNumberOfTracks(); iTrack++){
244 AliESDtrack *part = esdEvent.GetTrack(iTrack);
246 if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue;
247 if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
252 Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
255 Int_t sign = part->Charge();
256 AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
257 // use unique ID for track number
258 // this will be used in DistortTrack track to set the cluster label
259 // in one simulation the track id should be unique for performance studies
260 tempTrack->SetUniqueID(fCurrentTrack++);
261 DistortTrack(*tempTrack, time);
264 if(nUsedTracks >= fNtracks) break;
269 //________________________________________________________________
270 void AliToyMCEventGeneratorSimple::RunSimulationESD(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
273 // run simulation using esd input with equal event spacing
276 if (!ConnectOutputFile()) return;
277 TFile f(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"));
279 std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
282 TTree* esdTree = (TTree*) f.Get("esdTree");
284 std::cout <<"no esdTree in file " << std::endl;
288 AliESDEvent* esdEvent = new AliESDEvent();
291 esdEvent->ReadFromTree(esdTree);
294 // within one simulation the track count should be unique for effeciency studies
297 Double_t eventTime=0.;
298 const Double_t eventSpacing=1./50e3; //50kHz equally spaced
301 fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
303 Int_t nUsedEvents = 0;
304 for (Int_t ievent=0; ievent<esdTree->GetEntries(); ++ievent){
305 printf("Generating event %3d (%.3g)\n",nUsedEvents,eventTime);
306 esdTree->GetEvent(ievent);
307 if(esdEvent->GetNumberOfTracks()==0) {
308 std::cout << " tracks == 0" << std::endl;
312 fEvent = GenerateESD(*esdEvent, eventTime);
313 if(fEvent->GetNumberOfTracks() >=10) {
316 eventTime+=eventSpacing;
321 if(nUsedEvents>=nevents) break;
329 //________________________________________________________________
330 void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
333 // run simple simulation with equal event spacing
336 //Parameters for bunc
337 const Double_t abortGap = 3e-6; //
338 const Double_t collFreq = 50e3;
339 const Double_t bSpacing = 50e-9; //bunch spacing
340 const Int_t nTrainBunches = 48;
341 const Int_t nTrains = 12;
342 const Double_t revFreq = 1.11e4; //revolution frequency
343 const Double_t collProb = collFreq/(nTrainBunches*nTrains*revFreq);
344 const Double_t trainLength = bSpacing*(nTrainBunches-1);
345 const Double_t totTrainLength = nTrains*trainLength;
346 const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1);
347 Bool_t equalSpacing = kFALSE;
349 TRandom3 *rand = new TRandom3();
352 if (!ConnectOutputFile()) return;
353 //initialise the space charge. Should be done after the tree was set up
357 // within one simulation the track count should be unique for effeciency studies
360 Double_t eventTime=0.;
361 // const Double_t eventSpacing=1./50e3; //50kHz equally spaced
363 Int_t nGeneratedEvents = 0;
364 Int_t bunchCounter = 0;
365 Int_t trainCounter = 0;
367 //for (Int_t ievent=0; ievent<nevents; ++ievent){
368 while (nGeneratedEvents<nevents){
369 // std::cout <<trainCounter << " " << bunchCounter << " "<< "eventTime " << eventTime << std::endl;
372 printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
373 fEvent = Generate(eventTime);
378 eventTime+=1./collFreq;
381 Int_t nCollsInCrossing = rand -> Poisson(collProb);
382 for(Int_t iColl = 0; iColl<nCollsInCrossing; iColl++){
383 printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
384 fEvent = Generate(eventTime);
393 if(bunchCounter>=nTrainBunches){
396 if(trainCounter>=nTrains){
400 else eventTime+=trainSpacing;
404 else eventTime+= bSpacing;
420 //________________________________________________________________
421 Int_t AliToyMCEventGeneratorSimple::OpenInputAndGetMaxEvents(const Int_t type, const Int_t nevents) {
425 if(type==0) return nevents;
429 fInputFile = new TFile(Form("%s%s",fInputFileNameESD.Data(),"#AliESDs.root"),"read");
430 if(!fInputFile->IsOpen()) {
431 std::cout << "file "<<fInputFileNameESD.Data() <<" could not be opened" << std::endl;
434 fESDTree = (TTree*) fInputFile->Get("esdTree");
436 std::cout <<"no esdTree in file " << std::endl;
439 fESDEvent = new AliESDEvent();
440 fESDEvent->ReadFromTree(fESDTree);
441 fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1);
445 return fESDTree->GetEntries();
450 std::cout << " no available input type (toymc, esd) specified" << std::endl;
457 //________________________________________________________________
458 void AliToyMCEventGeneratorSimple::RunSimulation2(const Bool_t equalspacing, const Int_t type, const Int_t nevents, const Int_t ntracks) {
463 Int_t nMaxEvents = OpenInputAndGetMaxEvents(type,nevents);
464 if (!ConnectOutputFile()) return;
470 // within one simulation the track count should be unique for effeciency studies
473 Double_t eventTime=0.;
475 // const Double_t eventSpacing=1./50e3; //50kHz equally spaced
477 Int_t nGeneratedEvents = 0;
479 for (Int_t ievent=0; ievent<nMaxEvents; ++ievent){
480 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
483 Double_t spacing = 0;
484 GetNGeneratedEventsAndSpacing(equalspacing, nColls, spacing);
486 for(Int_t iColl = 0; iColl<nColls; iColl++){
487 if(type==0) fEvent = Generate(eventTime);
488 else if(type==1) fEvent = GenerateESD2(eventTime);
497 if(nGeneratedEvents >= nevents) break;
507 //________________________________________________________________
508 AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD2(Double_t time) {
510 //test that enough tracks will pass cuts (and that there is tracks at all)
511 Bool_t testEvent = kTRUE;
512 while(fESDTree->GetEvent(fInputIndex) && testEvent) {
514 Int_t nPassedCuts = 0;
515 for(Int_t iTrack = 0; iTrack<fESDEvent->GetNumberOfTracks(); iTrack++){
516 if(fESDCuts->AcceptTrack(fESDEvent->GetTrack(iTrack)) && (fESDEvent->GetTrack(iTrack)->Pt() < 0.3) ) nPassedCuts++;
518 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?
519 else testEvent = kFALSE;
522 if(fInputIndex>=fESDTree->GetEntries()) return 0;
526 AliToyMCEvent *retEvent = new AliToyMCEvent();
527 retEvent->SetT0(time);
528 retEvent->SetX(fESDEvent->GetPrimaryVertex()->GetX());
529 retEvent->SetY(fESDEvent->GetPrimaryVertex()->GetY());
530 retEvent->SetZ(fESDEvent->GetPrimaryVertex()->GetZ());
534 if(!fNtracks) fNtracks = fESDEvent->GetNumberOfTracks();
535 Int_t nUsedTracks = 0;
536 for(Int_t iTrack=0; iTrack<fESDEvent->GetNumberOfTracks(); iTrack++){
538 AliESDtrack *part = fESDEvent->GetTrack(iTrack);
540 if (!fESDCuts->AcceptTrack(/*(AliESDtrack*)*/part))continue;
541 if(part->Pt() < 0.3) continue; //avoid tracks that fail to propagate through entire volume
547 Double_t vertex[3]={retEvent->GetX(),retEvent->GetY(),retEvent->GetZ()};
549 Int_t sign = part->Charge();
551 AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign);
552 // use unique ID for track number
553 // this will be used in DistortTrack track to set the cluster label
554 // in one simulation the track id should be unique for performance studies
555 tempTrack->SetUniqueID(fCurrentTrack++);
556 DistortTrack(*tempTrack, time);
559 if(nUsedTracks >= fNtracks) break;
566 //________________________________________________________________
567 AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateLaser(Double_t time)
570 // Generate an Event with laser tracks
573 AliToyMCEvent *retEvent = new AliToyMCEvent();
574 retEvent->SetEventType(AliToyMCEvent::kLaser);
576 retEvent->SetT0(time);
581 AliTPCLaserTrack::LoadTracks();
582 TObjArray *arr=AliTPCLaserTrack::GetTracks();
584 //since we have a laser track force no material budges
585 Bool_t materialBudget=GetUseMaterialBudget();
586 SetUseMaterialBudget(kFALSE);
588 //the laser tracks have partially large inclination angles over the pads
589 // -> relax the propagation contraint and switch off error messages
591 Int_t debug=AliLog::GetDebugLevel("","AliToyMCEventGeneratorSimple");
592 AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",-3);
594 for (Int_t iTrack=0; iTrack<arr->GetEntriesFast(); ++iTrack){
595 AliExternalTrackParam *track=(AliExternalTrackParam*)arr->At(iTrack);
596 AliToyMCTrack *tempTrack = retEvent->AddTrack(AliToyMCTrack(*track));
597 // for laser only TPC clusters exist
598 tempTrack->SetUniqueID(fCurrentTrack++);
599 MakeTPCClusters(*tempTrack, time);
603 AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",debug);
604 SetUseMaterialBudget(materialBudget);
608 //________________________________________________________________
609 void AliToyMCEventGeneratorSimple::GetNGeneratedEventsAndSpacing(const Bool_t equalSpacing, Int_t &ngen, Double_t &spacing)
612 static Int_t bunchCounter = 0;
613 static Int_t trainCounter = 0;
617 spacing = 1./50e3; //50kHz equally spaced
620 else if(!equalSpacing) {
621 //parameters for bunch train
622 const Double_t abortGap = 3e-6; //
623 const Double_t collFreq = 50e3;
624 const Double_t bSpacing = 50e-9; //bunch spacing
625 const Int_t nTrainBunches = 48;
626 const Int_t nTrains = 12;
627 const Double_t revFreq = 1.11e4; //revolution frequency
628 const Double_t collProb = collFreq/(nTrainBunches*nTrains*revFreq);
629 const Double_t trainLength = bSpacing*(nTrainBunches-1);
630 const Double_t totTrainLength = nTrains*trainLength;
631 const Double_t trainSpacing = (1./revFreq - abortGap - totTrainLength)/(nTrains-1);
633 Int_t nCollsInCrossing = 0;
634 while(nCollsInCrossing ==0){
639 if(bunchCounter>=nTrainBunches){
642 if(trainCounter>=nTrains){
646 else time+=trainSpacing;
650 else time+= bSpacing;
653 nCollsInCrossing = gRandom -> Poisson(collProb);
654 //std::cout << " nCollsInCrossing " << nCollsInCrossing << std::endl;
656 ngen = nCollsInCrossing;
657 if(nCollsInCrossing > 1)std::cout << " nCollsInCrossing " << nCollsInCrossing << std::endl;
665 //________________________________________________________________
666 Bool_t AliToyMCEventGeneratorSimple::CloseInputFile()
669 // close the input file
671 if (!fInputFile) return kFALSE;