Store total number of produced particles in the header.
[u/mrichter/AliRoot.git] / EVGEN / AliGenBeamGasNew.cxx
CommitLineData
e3bd319d 1//-------------------------------------------------------------
2// Generate Beam-Gas-Events
3// Default: underlying event: pO @ 7 TeV (fixed target)
4//
5// underlying event can be changed by adding generators
97bd035c 6// in the same as to AliGenCocktail before calling the Init()
e3bd319d 7//-------------------------------------------------------------
8
97bd035c 9/* $Id$ */
10
11#include <TList.h>
12#include <TObjArray.h>
13#include "TParticle.h"
14
e3bd319d 15#include "AliGenBeamGasNew.h"
16#include "AliGenCocktail.h"
17#include "AliGenCocktailEntry.h"
e3bd319d 18#include "AliGenCocktailEventHeader.h"
97bd035c 19#include "../THijing/AliGenHijing.h"
e3bd319d 20#include "AliCollisionGeometry.h"
21#include "AliRun.h"
22#include "AliMC.h"
23#include "AliStack.h"
97bd035c 24#include "AliLog.h"
e3bd319d 25
26ClassImp(AliGenBeamGasNew)
27
28Float_t EtaToTheta(Float_t arg) { return (180./TMath::Pi())*2.*atan(exp(-arg)); }
29
30AliGenBeamGasNew::AliGenBeamGasNew() :
31 AliGenCocktail(),
97bd035c 32 fTwindow(88.e-6)
e3bd319d 33{
97bd035c 34 // Default constructor
e3bd319d 35}
36
37AliGenBeamGasNew::~AliGenBeamGasNew()
38{
97bd035c 39 // Destructor
e3bd319d 40}
41
97bd035c 42AliGenBeamGasNew::AliGenBeamGasNew(const AliGenBeamGasNew& /*rhs*/) : AliGenCocktail()
43{
44 AliFatal("Copy Constructor not yet implemented!");
e3bd319d 45}
46
e3bd319d 47
48void AliGenBeamGasNew::Init()
49{
50 // Initialisation of the class
51 // if no generators were added before calling Init()
52 // AliGenHijing is added configured for pO
53
54 fVertexSmear = kPerEvent;
55 fVertexSource = kInternal;
56 fRandom = kTRUE;
57
58 // Adding default underlying event in case none was specified
59 // p-O-collision at 7 TeV (fixed target)
60 if (!fEntries) {
61 AliGenHijing *gen = new AliGenHijing(-1);
62 gen->SetEnergyCMS(7000);
63 gen->SetReferenceFrame("LAB");
64 gen->SetProjectile("P",1,1);
65 gen->SetTarget("A",16,8);
66 gen->SetPhiRange(0,360);
67 Float_t thmin = EtaToTheta(8);
68 Float_t thmax = EtaToTheta(-8);
69 gen->SetThetaRange(thmin, thmax);
70 gen->SetOrigin(0,0,0);
71 gen->SetSigma(0,0,0);
72 gen->KeepFullEvent();
73 gen->SetShadowing(1);
74 gen->SetDecaysOff(1);
75 gen->SetSpectators(1);
76 gen->SetSelectAll(1);
97bd035c 77 gen->SetRandomPz(kTRUE);
e3bd319d 78 gen->SetPileUpTimeWindow(0.);
79 AddGenerator(gen,"Hijing pO",1);
80 }
81
82 TIter next(fEntries);
83 AliGenCocktailEntry *entry;
84
85 // Loop over generators and initialize
86 while((entry = (AliGenCocktailEntry*)next())) {
97bd035c 87 if (fStack) entry->Generator()->SetStack(fStack);
88 entry->Generator()->Init();
e3bd319d 89 }
90
91 next.Reset();
92
93 if (fRandom) {
94 fProb.Set(fNGenerators);
95 next.Reset();
96 Float_t sum = 0.;
97 while((entry = (AliGenCocktailEntry*)next())) {
98 sum += entry->Rate();
99 }
100
101 next.Reset();
102 Int_t i = 0;
103 Float_t psum = 0.;
104 while((entry = (AliGenCocktailEntry*)next())) {
105 psum += entry->Rate() / sum;
106 fProb[i++] = psum;
107 }
108 }
109 next.Reset();
110}
111
112void AliGenBeamGasNew::VertexInternal()
113{
114 // calculation of the interaction vertex for the beam gas event
115 // both spatial and time coordinate are adjusted.
97bd035c 116
117 Float_t random[2];
118 Float_t nbunch;
119 Int_t bunch;
120
121 Rndm(random,2);
122 fVertex[0] = fVertex[1] = 0.;
123 fVertex[2] = random[1] * 4000. - 2000.;
124 AliInfo(Form("z-Koord.: %13.3f \n", fVertex[2]));
125 nbunch = 2 * fTwindow / 25.e-9;
126 bunch = (Int_t) (nbunch * (random[0] - 0.5));
127 fItime = fVertex[2] / 3.e8 + bunch * 25.e-9;
128 AliInfo(Form("fItime: %13.3e \n", fItime));
e3bd319d 129}
130
131void AliGenBeamGasNew::Generate()
132{
e3bd319d 133 TIter next(fEntries);
134 AliGenCocktailEntry *entry = 0;
135 AliGenCocktailEntry *preventry = 0;
136 AliGenerator* gen = 0;
e3bd319d 137
138 if (fHeader) delete fHeader;
139 fHeader = new AliGenCocktailEventHeader("Beamgas Header");
140
141 TObjArray *partArray = gAlice->GetMCApp()->Particles();
e3bd319d 142
97bd035c 143 TArrayF eventVertex;
144 int numbg = 0; // variable for counting how many bg interactions happened
145
146 do {
147 numbg++;
e3bd319d 148 Vertex();
e3bd319d 149 eventVertex.Set(3);
150 for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
151
152 if (!fRandom) {
153 Int_t igen=0;
154 while((entry = (AliGenCocktailEntry*)next())) {
155 if (fUsePerEventRate && (gRandom->Rndm() > entry->Rate())) continue;
156
157 igen++;
158 if (igen ==1) {
97bd035c 159 entry->SetFirst(0);
e3bd319d 160 } else {
161 entry->SetFirst((partArray->GetEntriesFast())+1);
162 }
163 //
164 // Handle case in which current generator needs collision geometry from previous generator
165 //
166 gen = entry->Generator();
167 if (gen->NeedsCollisionGeometry())
168 {
169 if (preventry && preventry->Generator()->ProvidesCollisionGeometry())
170 {
171 gen->SetCollisionGeometry(preventry->Generator()->CollisionGeometry());
172 } else {
173 Fatal("Generate()", "No Collision Geometry Provided");
174 }
175 }
176 entry->Generator()->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
177 entry->Generator()->Generate();
178 entry->SetLast(partArray->GetEntriesFast());
179 preventry = entry;
180 }
181 } else {
182 //
183 // Select a generator randomly
184 //
185 Int_t i;
186 Float_t p0 = gRandom->Rndm();
187
188 for (i = 0; i < fNGenerators; i++) {
189 if (p0 < fProb[i]) break;
190 }
191
192 entry = (AliGenCocktailEntry*) fEntries->At(i);
97bd035c 193 entry->SetFirst(0);
e3bd319d 194 gen = entry->Generator();
195 entry->Generator()->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
196 entry->Generator()->Generate();
197 entry->SetLast(partArray->GetEntriesFast());
198 }
199
97bd035c 200 AliStack *stack = gAlice->GetRunLoader()->Stack();
201 TLorentzVector v;
202 for (int k = 0; k < stack->GetNprimary(); k++) {
e3bd319d 203 stack->Particle(k)->ProductionVertex(v);
97bd035c 204 v[3] += fItime;
e3bd319d 205 stack->Particle(k)->SetProductionVertex(v);
e3bd319d 206 }
97bd035c 207
e3bd319d 208 next.Reset();
97bd035c 209 } while (gRandom->Rndm() < 1. / 5.); // must be adjusted to the rate of bg interactions
e3bd319d 210 // Event Vertex
211 fHeader->SetPrimaryVertex(eventVertex);
97bd035c 212 // Total number of produced an stored particles
213 fHeader->CalcNProduced();
214
215 gAlice->SetGenEventHeader(fHeader);
e3bd319d 216}