AliToyMCEventGenerator would crash if a sector only gets one space point (TSpline3...
[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>
a1a695e5 7
de0014b7 8#include "AliToyMCEvent.h"
526ddf0e 9
a1a695e5 10#include "AliToyMCEventGeneratorSimple.h"
11
de0014b7 12ClassImp(AliToyMCEventGeneratorSimple);
526ddf0e 13
14
de0014b7 15AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple()
16 :AliToyMCEventGenerator()
526ddf0e 17 ,fVertexMean(0.)
18 ,fVertexSigma(7.)
32438f4e 19 ,fNtracks(20)
526ddf0e 20{
21 //
22
23}
24//________________________________________________________________
de0014b7 25AliToyMCEventGeneratorSimple::AliToyMCEventGeneratorSimple(const AliToyMCEventGeneratorSimple &gen)
26 :AliToyMCEventGenerator(gen)
526ddf0e 27 ,fVertexMean(gen.fVertexMean)
28 ,fVertexSigma(gen.fVertexSigma)
32438f4e 29 ,fNtracks(gen.fNtracks)
526ddf0e 30{
31}
32//________________________________________________________________
de0014b7 33AliToyMCEventGeneratorSimple::~AliToyMCEventGeneratorSimple()
526ddf0e 34{
35}
36 //________________________________________________________________
de0014b7 37AliToyMCEventGeneratorSimple& AliToyMCEventGeneratorSimple::operator = (const AliToyMCEventGeneratorSimple &gen)
526ddf0e 38{
39 //assignment operator
40 if (&gen == this) return *this;
de0014b7 41 new (this) AliToyMCEventGeneratorSimple(gen);
526ddf0e 42
43 return *this;
44}
45//________________________________________________________________
de0014b7 46void AliToyMCEventGeneratorSimple::SetParameters(Double_t vertexMean, Double_t vertexSigma) {
526ddf0e 47 fVertexMean = vertexMean;
48 fVertexSigma = vertexSigma;
49}
50//________________________________________________________________
cd8ed0ac 51AliToyMCEvent* AliToyMCEventGeneratorSimple::Generate(Double_t time)
52{
53 //
54 //
55 //
56
de0014b7 57 AliToyMCEvent *retEvent = new AliToyMCEvent();
526ddf0e 58 retEvent->SetT0(time);
59 retEvent->SetX(0);
60 retEvent->SetX(0);
61 retEvent->SetZ(gRandom->Gaus(fVertexMean,fVertexSigma));
62
63 Double_t etaCuts=.9;
64 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
65 TF1 fpt("fpt",Form("x*(1+(sqrt(x*x+%f^2)-%f)/([0]*[1]))^(-[0])",mass,mass),0.4,10);
66 fpt.SetParameters(7.24,0.120);
67 fpt.SetNpx(10000);
32438f4e 68// Int_t nTracks = 400; //TODO: draw from experim dist
69// Int_t nTracks = 20; //TODO: draw from experim dist
526ddf0e 70
32438f4e 71 for(Int_t iTrack=0; iTrack<fNtracks; iTrack++){
526ddf0e 72 Double_t phi = gRandom->Uniform(0.0, 2*TMath::Pi());
73 Double_t eta = gRandom->Uniform(-etaCuts, etaCuts);
74 Double_t pt = fpt.GetRandom(); // momentum for f1
75 Short_t sign=1;
76 if(gRandom->Rndm() < 0.5){
77 sign =1;
78 }else{
79 sign=-1;
80 }
81
82 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta))-TMath::Pi()/2.;
83 Double_t pxyz[3];
84 pxyz[0]=pt*TMath::Cos(phi);
85 pxyz[1]=pt*TMath::Sin(phi);
86 pxyz[2]=pt*TMath::Tan(theta);
87 Double_t vertex[3]={0,0,retEvent->GetZ()};
88 Double_t cv[21]={0};
de0014b7 89 AliToyMCTrack *tempTrack = new AliToyMCTrack(vertex,pxyz,cv,sign);
526ddf0e 90
91 Bool_t trackDist = DistortTrack(*tempTrack, time);
92 if(trackDist) retEvent->AddTrack(*tempTrack);
93 delete tempTrack;
94 }
95
96 return retEvent;
97}
98
a1a695e5 99//________________________________________________________________
32438f4e 100void AliToyMCEventGeneratorSimple::RunSimulation(const Int_t nevents/*=10*/, const Int_t ntracks/*=400*/)
a1a695e5 101{
102 //
103 // run simple simulation with equal event spacing
104 //
105
106 if (!ConnectOutputFile()) return;
cd8ed0ac 107 //initialise the space charge. Should be done after the tree was set up
108 InitSpaceCharge();
109
32438f4e 110 fNtracks=ntracks;
a1a695e5 111 Double_t eventTime=0.;
112 const Double_t eventSpacing=1./50e3; //50kHz equally spaced
32438f4e 113 TStopwatch s;
a1a695e5 114 for (Int_t ievent=0; ievent<nevents; ++ievent){
32438f4e 115 printf("Generating event %3d (%.3g)\n",ievent,eventTime);
a1a695e5 116 fEvent = Generate(eventTime);
117 FillTree();
118 delete fEvent;
119 fEvent=0x0;
120 eventTime+=eventSpacing;
121 }
32438f4e 122 s.Stop();
123 s.Print();
124
a1a695e5 125 CloseOutputFile();
126}
526ddf0e 127