#include <iostream>
+#include <TROOT.h>
#include <TDatabasePDG.h>
#include <TRandom.h>
#include <TF1.h>
#include <TStopwatch.h>
-#include <AliESDtrackCuts.h>
-#include <AliESDtrack.h>
#include <TFile.h>
#include <TTree.h>
#include <TRandom3.h>
+#include <AliESDtrackCuts.h>
+#include <AliESDtrack.h>
+#include <AliTPCLaserTrack.h>
+
#include "AliToyMCEvent.h"
#include "AliToyMCEventGeneratorSimple.h"
+/*
+
+
+*/
ClassImp(AliToyMCEventGeneratorSimple);
,fESDEvent(0x0)
,fESDTree(0x0)
,fInputFile(0x0)
-
+ ,fHPt(0x0)
+ ,fHEta(0x0)
+ ,fHMult(0x0)
+ ,fHistosSet(kFALSE)
+ ,fParamFile(0x0)
{
//
,fESDEvent(gen.fESDEvent)
,fESDTree(gen.fESDTree)
,fInputFile(gen.fInputFile)
-
+ ,fHPt(gen.fHPt)
+ ,fHEta(gen.fHEta)
+ ,fHMult(gen.fHMult)
+ ,fHistosSet(gen.fHistosSet)
+ ,fParamFile(0x0)
{
}
//________________________________________________________________
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; iTrack<fNtracks; iTrack++){
+ Int_t nTracksLocal = fNtracks;
+ if(fHistosSet) {
+ nTracksLocal = fHMult->GetRandom();
+ if(nTracksLocal > fNtracks) nTracksLocal = fNtracks;
+ }
+ std::cout << "Generating " << nTracksLocal << " tracks " << std::endl;
+
+ for(Int_t iTrack=0; iTrack<nTracksLocal; iTrack++){
Double_t phi = gRandom->Uniform(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;
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; ievent<nevents; ++ievent){
printf("Generating event %3d (%.3g)\n",ievent,eventTime);
fEvent = Generate(eventTime);
+ SetSCScalingFactor();
+ FillTree();
+ delete fEvent;
+ fEvent=0x0;
+ eventTime+=eventSpacing;
+ }
+ s.Stop();
+ s.Print();
+
+ CloseOutputFile();
+}
+
+//________________________________________________________________
+void AliToyMCEventGeneratorSimple::RunSimulationLaser(const Int_t nevents/*=1*/)
+{
+ //
+ // run simple simulation with equal event spacing
+ //
+
+ if (!ConnectOutputFile()) return;
+ //initialise the space charge. Should be done after the tree was set up
+ InitSpaceCharge();
+
+ // within one simulation the track count should be unique for effeciency studies
+ fCurrentTrack=1;
+
+ Double_t eventTime=0.;
+ const Double_t eventSpacing=1./10.; //laser is running at 10Hz equally spaced
+ TStopwatch s;
+ for (Int_t ievent=0; ievent<nevents; ++ievent){
+ printf("Generating event %3d (%.3g)\n",ievent,eventTime);
+ fEvent = GenerateLaser(eventTime);
FillTree();
delete fEvent;
fEvent=0x0;
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);
+ nUsedTracks++;
if(nUsedTracks >= fNtracks) break;
}
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; ievent<esdTree->GetEntries(); ++ievent){
printf("Generating event %3d (%.3g)\n",nUsedEvents,eventTime);
}
//________________________________________________________________
-
void AliToyMCEventGeneratorSimple::RunSimulationBunchTrain(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
{
//
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;
if(equalSpacing) {
printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
fEvent = Generate(eventTime);
+ SetSCScalingFactor();
nGeneratedEvents++;
FillTree();
delete fEvent;
for(Int_t iColl = 0; iColl<nCollsInCrossing; iColl++){
printf("Generating event %3d (%.3g)\n",nGeneratedEvents,eventTime);
fEvent = Generate(eventTime);
+ SetSCScalingFactor();
nGeneratedEvents++;
FillTree();
delete fEvent;
fInputIndex = 0;
+ gRandom->SetSeed();
return fESDTree->GetEntries();
}
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
//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;
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;
}
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; iTrack<arr->GetEntriesFast(); ++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)
{
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff