X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FUpgrade%2FAliToyMCEventGeneratorSimple.cxx;h=8123181b39dffb431461a902d0783c3a823ba097;hb=70cc54298863351c61ccbcce95252920485b9afa;hp=19588671ea97bc63fe6b8a8e0bbfc1c01bcd65f9;hpb=29b7f480eb4d33bdf00817a2e23b55d01ad99ea0;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/Upgrade/AliToyMCEventGeneratorSimple.cxx b/TPC/Upgrade/AliToyMCEventGeneratorSimple.cxx index 19588671ea9..8123181b39d 100644 --- a/TPC/Upgrade/AliToyMCEventGeneratorSimple.cxx +++ b/TPC/Upgrade/AliToyMCEventGeneratorSimple.cxx @@ -1,19 +1,26 @@ #include +#include #include #include #include #include -#include -#include #include #include #include +#include +#include +#include + #include "AliToyMCEvent.h" #include "AliToyMCEventGeneratorSimple.h" +/* + + +*/ ClassImp(AliToyMCEventGeneratorSimple); @@ -29,7 +36,11 @@ AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple() ,fESDEvent(0x0) ,fESDTree(0x0) ,fInputFile(0x0) - + ,fHPt(0x0) + ,fHEta(0x0) + ,fHMult(0x0) + ,fHistosSet(kFALSE) + ,fParamFile(0x0) { // @@ -46,7 +57,11 @@ AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple(const AliToyMCEventGe ,fESDEvent(gen.fESDEvent) ,fESDTree(gen.fESDTree) ,fInputFile(gen.fInputFile) - + ,fHPt(gen.fHPt) + ,fHEta(gen.fHEta) + ,fHMult(gen.fHMult) + ,fHistosSet(gen.fHistosSet) + ,fParamFile(0x0) { } //________________________________________________________________ @@ -63,35 +78,69 @@ AliToyMCEventGeneratorSimple& AliToyMCEventGeneratorSimple::operator = (const Al return *this; } //________________________________________________________________ -void AliToyMCEventGeneratorSimple::SetParametersSimple(Double_t vertexMean, Double_t vertexSigma) { +void AliToyMCEventGeneratorSimple::SetParametersToyGen(const Char_t* parfilename/*="files/params.root*/, Double_t vertexMean/*=0*/, Double_t vertexSigma/*=7.*/) { fVertexMean = vertexMean; fVertexSigma = vertexSigma; +// fParamFile = new TFile(parfilename, "read"); + TFile f(parfilename); + gROOT->cd(); + fHPt = (TH1F*) f.Get("hPt"); + fHEta = (TH1F*) f.Get("hEta"); + fHMult = (TH1I*) f.Get("hMult") ; + fHistosSet = kTRUE; + f.Close(); + + } //________________________________________________________________ AliToyMCEvent* AliToyMCEventGeneratorSimple::Generate(Double_t time) { // + // Generate an event at 'time' // - // + + // iterate over space charge maps in case they are set + IterateSC(); AliToyMCEvent *retEvent = new AliToyMCEvent(); retEvent->SetT0(time); retEvent->SetX(0); retEvent->SetX(0); - retEvent->SetZ(gRandom->Gaus(fVertexMean,fVertexSigma)); + Double_t vertexZ=0; + //limit vertex to +- 10cm + while ( TMath::Abs(vertexZ=gRandom->Gaus(fVertexMean,fVertexSigma))>10. ); + retEvent->SetZ(vertexZ); - Double_t etaCuts=.9; + Double_t etaCuts=0.9; Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass(); - TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10); - fpt.SetParameters(7.24,0.120); - fpt.SetNpx(10000); + static TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.3,100); + if (fpt.GetParameter(0)<1){ + printf("Set Parameters\n"); + fpt.SetParameters(7.24,0.120); + fpt.SetNpx(200); + } // Int_t nTracks = 400; //TODO: draw from experim dist // Int_t nTracks = 20; //TODO: draw from experim dist - - for(Int_t iTrack=0; iTrackGetRandom(); + if(nTracksLocal > fNtracks) nTracksLocal = fNtracks; + } + std::cout << "Generating " << nTracksLocal << " tracks " << std::endl; + + for(Int_t iTrack=0; iTrackUniform(0.0, 2*TMath::Pi()); - Double_t eta = gRandom->Uniform(-etaCuts, etaCuts); - Double_t pt = fpt.GetRandom(); // momentum for f1 + Double_t eta = 0.; + Double_t pt = 0.; + + if(fHistosSet) {//draw from histograms if set + eta = fHEta->GetRandom(); + pt = fHPt->GetRandom(); + } + else { + eta = gRandom->Uniform(-etaCuts, etaCuts); + pt = fpt.GetRandom(); // momentum for f1 + } Short_t sign=1; if(gRandom->Rndm() < 0.5){ sign =1; @@ -106,37 +155,73 @@ AliToyMCEvent* AliToyMCEventGeneratorSimple::Generate(Double_t time) pxyz[2]=pt*TMath::Tan(theta); Double_t vertex[3]={0,0,retEvent->GetZ()}; Double_t cv[21]={0}; - AliToyMCTrack *tempTrack = new AliToyMCTrack(vertex,pxyz,cv,sign); + AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign); // use unique ID for track number // this will be used in DistortTrack track to set the cluster label - tempTrack->SetUniqueID(iTrack); - - Bool_t trackDist = DistortTrack(*tempTrack, time); - if(trackDist) retEvent->AddTrack(*tempTrack); - delete tempTrack; + // in one simulation the track id should be unique for performance studies + tempTrack->SetUniqueID(fCurrentTrack++); + DistortTrack(*tempTrack, time); } return retEvent; } //________________________________________________________________ -void AliToyMCEventGeneratorSimple::RunSimulation(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/) +void AliToyMCEventGeneratorSimple::RunSimulation(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/, const Int_t rate/*=50*/) { // // run simple simulation with equal event spacing + // rate in kHz // if (!ConnectOutputFile()) return; //initialise the space charge. Should be done after the tree was set up InitSpaceCharge(); - + + // number of tracks to simulate per interaction fNtracks=ntracks; + // within one simulation the track count should be unique for effeciency studies + // don't use the track ID 0 since the label -0 is not different from 0 + fCurrentTrack=1; + Double_t eventTime=0.; - const Double_t eventSpacing=1./50e3; //50kHz equally spaced + const Double_t eventSpacing=1./rate/1e3; //50kHz equally spaced TStopwatch s; for (Int_t ievent=0; ieventCharge(); - AliToyMCTrack *tempTrack = new AliToyMCTrack(vertex,pxyz,cv,sign); + AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign); // use unique ID for track number // this will be used in DistortTrack track to set the cluster label - tempTrack->SetUniqueID(iTrack); - - - Bool_t trackDist = DistortTrack(*tempTrack, time); - if(trackDist) { - retEvent->AddTrack(*tempTrack); - - nUsedTracks++; - } - delete tempTrack; + // in one simulation the track id should be unique for performance studies + tempTrack->SetUniqueID(fCurrentTrack++); + DistortTrack(*tempTrack, time); + nUsedTracks++; if(nUsedTracks >= fNtracks) break; } @@ -219,13 +298,16 @@ void AliToyMCEventGeneratorSimple::RunSimulationESD(const Int_t nevents/*=10*/, esdEvent->ReadFromTree(esdTree); - fNtracks = ntracks; + fNtracks = ntracks; + // within one simulation the track count should be unique for effeciency studies + fCurrentTrack=1; + Double_t eventTime=0.; const Double_t eventSpacing=1./50e3; //50kHz equally spaced TStopwatch s; fESDCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kTRUE,1); - + gRandom->SetSeed(); Int_t nUsedEvents = 0; for (Int_t ievent=0; ieventGetEntries(); ++ievent){ printf("Generating event %3d (%.3g)\n",nUsedEvents,eventTime); @@ -253,7 +335,6 @@ void AliToyMCEventGeneratorSimple::RunSimulationESD(const Int_t nevents/*=10*/, } //________________________________________________________________ - void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/) { // @@ -281,6 +362,9 @@ void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(const Int_t nevents/* InitSpaceCharge(); fNtracks=ntracks; + // within one simulation the track count should be unique for effeciency studies + fCurrentTrack=1; + Double_t eventTime=0.; // const Double_t eventSpacing=1./50e3; //50kHz equally spaced TStopwatch s; @@ -295,6 +379,7 @@ void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(const Int_t nevents/* if(equalSpacing) { printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime); fEvent = Generate(eventTime); + SetSCScalingFactor(); nGeneratedEvents++; FillTree(); delete fEvent; @@ -306,6 +391,7 @@ void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(const Int_t nevents/* for(Int_t iColl = 0; iCollSetSeed(); return fESDTree->GetEntries(); } @@ -386,13 +473,13 @@ void AliToyMCEventGeneratorSimple::RunSimulation2(const Bool_t equalspacing, con Int_t nMaxEvents = OpenInputAndGetMaxEvents(type,nevents); if (!ConnectOutputFile()) return; - InitSpaceCharge(); - - + InitSpaceCharge(); fNtracks = ntracks; - + // within one simulation the track count should be unique for effeciency studies + fCurrentTrack=1; + Double_t eventTime=0.; TStopwatch s; // const Double_t eventSpacing=1./50e3; //50kHz equally spaced @@ -432,6 +519,9 @@ AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD2(Double_t time) { //test that enough tracks will pass cuts (and that there is tracks at all) Bool_t testEvent = kTRUE; + + // iterate over space charge maps in case they are set + IterateSC(); while(fESDTree->GetEvent(fInputIndex) && testEvent) { Int_t nPassedCuts = 0; @@ -471,19 +561,13 @@ AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD2(Double_t time) { Double_t cv[21]={0}; Int_t sign = part->Charge(); - AliToyMCTrack *tempTrack = new AliToyMCTrack(vertex,pxyz,cv,sign); + AliToyMCTrack *tempTrack = retEvent->AddTrack(vertex,pxyz,cv,sign); // use unique ID for track number // this will be used in DistortTrack track to set the cluster label - tempTrack->SetUniqueID(iTrack); - - - Bool_t trackDist = DistortTrack(*tempTrack, time); - if(trackDist) { - retEvent->AddTrack(*tempTrack); - - nUsedTracks++; - } - delete tempTrack; + // in one simulation the track id should be unique for performance studies + tempTrack->SetUniqueID(fCurrentTrack++); + DistortTrack(*tempTrack, time); + if(nUsedTracks >= fNtracks) break; } @@ -492,6 +576,51 @@ AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateESD2(Double_t time) { return retEvent; } +//________________________________________________________________ +AliToyMCEvent* AliToyMCEventGeneratorSimple::GenerateLaser(Double_t time) +{ + // + // Generate an Event with laser tracks + // + + // iterate over space charge maps in case they are set + IterateSC(); + + AliToyMCEvent *retEvent = new AliToyMCEvent(); + retEvent->SetEventType(AliToyMCEvent::kLaser); + + retEvent->SetT0(time); + retEvent->SetX(0); + retEvent->SetX(0); + retEvent->SetZ(0); + + AliTPCLaserTrack::LoadTracks(); + TObjArray *arr=AliTPCLaserTrack::GetTracks(); + + //since we have a laser track force no material budges + Bool_t materialBudget=GetUseMaterialBudget(); + SetUseMaterialBudget(kFALSE); + + //the laser tracks have partially large inclination angles over the pads + // -> relax the propagation contraint and switch off error messages + SetIsLaser(kTRUE); + Int_t debug=AliLog::GetDebugLevel("","AliToyMCEventGeneratorSimple"); + AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",-3); + + for (Int_t iTrack=0; iTrackGetEntriesFast(); ++iTrack){ + AliExternalTrackParam *track=(AliExternalTrackParam*)arr->At(iTrack); + AliToyMCTrack *tempTrack = retEvent->AddTrack(AliToyMCTrack(*track)); + // for laser only TPC clusters exist + tempTrack->SetUniqueID(fCurrentTrack++); + MakeTPCClusters(*tempTrack, time); + } + + SetIsLaser(kFALSE); + AliLog::SetClassDebugLevel("AliToyMCEventGeneratorSimple",debug); + SetUseMaterialBudget(materialBudget); + return retEvent; +} + //________________________________________________________________ void AliToyMCEventGeneratorSimple::GetNGeneratedEventsAndSpacing(const Bool_t equalSpacing, Int_t &ngen, Double_t &spacing) {