Generator for beam gas collisions. (Jochen Klein)
[u/mrichter/AliRoot.git] / EVGEN / AliGenBeamGasNew.cxx
CommitLineData
e3bd319d 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18
19
20//-------------------------------------------------------------
21// Generate Beam-Gas-Events
22// Default: underlying event: pO @ 7 TeV (fixed target)
23//
24// underlying event can be changed by adding generators
25// in the same way as to AliGenCocktail before calling
26// the Init()
27// Author: Jochen Klein
28//-------------------------------------------------------------
29
30#include "AliGenBeamGasNew.h"
31#include "AliGenCocktail.h"
32#include "AliGenCocktailEntry.h"
33#include "AliGenEventHeader.h"
34#include "AliGenCocktailEventHeader.h"
35#include "AliGenHijing.h"
36#include <TList.h>
37#include <TObjArray.h>
38#include "AliCollisionGeometry.h"
39#include "AliRun.h"
40#include "AliMC.h"
41#include "AliStack.h"
42#include "AliHeader.h"
43#include "TParticle.h"
44
45ClassImp(AliGenBeamGasNew)
46
47Float_t EtaToTheta(Float_t arg) { return (180./TMath::Pi())*2.*atan(exp(-arg)); }
48
49AliGenBeamGasNew::AliGenBeamGasNew() :
50 AliGenCocktail(),
51 fItime(0),
52 fTwindow(88e-6),
53 fZwindow(2000),
54 fRate(1e3)
55{
56 printf("AliGenBeamGasNew instatiated!\n");
57}
58
59AliGenBeamGasNew::~AliGenBeamGasNew()
60{
61 delete fEntries;
62 fEntries = 0;
63 fHeader = 0;
64}
65
66void AliGenBeamGasNew::SetTimeWindow(Float_t twindow) { fTwindow = twindow; }
67
68bool AliGenBeamGasNew::SetRate(Float_t rate) {
69 if (rate >= 0) {
70 fRate = rate;
71 return true;
72 } else
73 return false;
74}
75
76void AliGenBeamGasNew::SetZWindow(Float_t zwindow) { fZwindow = zwindow; }
77
78void AliGenBeamGasNew::Init()
79{
80 // Initialisation of the class
81 // if no generators were added before calling Init()
82 // AliGenHijing is added configured for pO
83
84 fVertexSmear = kPerEvent;
85 fVertexSource = kInternal;
86 fRandom = kTRUE;
87
88 // Adding default underlying event in case none was specified
89 // p-O-collision at 7 TeV (fixed target)
90 if (!fEntries) {
91 AliGenHijing *gen = new AliGenHijing(-1);
92 gen->SetEnergyCMS(7000);
93 gen->SetReferenceFrame("LAB");
94 gen->SetProjectile("P",1,1);
95 gen->SetTarget("A",16,8);
96 gen->SetPhiRange(0,360);
97 Float_t thmin = EtaToTheta(8);
98 Float_t thmax = EtaToTheta(-8);
99 gen->SetThetaRange(thmin, thmax);
100 gen->SetOrigin(0,0,0);
101 gen->SetSigma(0,0,0);
102 gen->KeepFullEvent();
103 gen->SetShadowing(1);
104 gen->SetDecaysOff(1);
105 gen->SetSpectators(1);
106 gen->SetSelectAll(1);
107 gen->SetRandomPz(kFALSE);
108 gen->SetPileUpTimeWindow(0.);
109 AddGenerator(gen,"Hijing pO",1);
110 }
111
112 TIter next(fEntries);
113 AliGenCocktailEntry *entry;
114
115 // Loop over generators and initialize
116 while((entry = (AliGenCocktailEntry*)next())) {
117 if (fStack) entry->Generator()->SetStack(fStack);
118 entry->Generator()->Init();
119 }
120
121 next.Reset();
122
123 if (fRandom) {
124 fProb.Set(fNGenerators);
125 next.Reset();
126 Float_t sum = 0.;
127 while((entry = (AliGenCocktailEntry*)next())) {
128 sum += entry->Rate();
129 }
130
131 next.Reset();
132 Int_t i = 0;
133 Float_t psum = 0.;
134 while((entry = (AliGenCocktailEntry*)next())) {
135 psum += entry->Rate() / sum;
136 fProb[i++] = psum;
137 }
138 }
139 next.Reset();
140}
141
142void AliGenBeamGasNew::VertexInternal()
143{
144 // calculation of the interaction vertex for the beam gas event
145 // both spatial and time coordinate are adjusted.
146
147 Float_t random[2];
148 Float_t nbunch;
149 Float_t bunch;
150
151 gRandom->RndmArray(2,random);
152 fVertex[0] = fVertex[1] = 0;
153 fVertex[2] = random[0] * 4000 - 2000;
154 nbunch = fTwindow/25e-9;
155 bunch = floor(2*nbunch * random[1] - nbunch);
156 fItime = fVertex[2] / 100 / 3e8 + bunch * 25e-9;
157}
158
159void AliGenBeamGasNew::Generate()
160{
161//
162// Generate the collisions for one event
163//
164 Int_t numbg = gRandom->Poisson(fRate*fZwindow/100*2*fTwindow);
165 if (numbg == 0) return;
166
167 TIter next(fEntries);
168 AliGenCocktailEntry *entry = 0;
169 AliGenCocktailEntry *preventry = 0;
170 AliGenerator* gen = 0;
171 TLorentzVector v;
172 TArrayF eventVertex;
173 Int_t lastpart=0;
174 Float_t sign;
175
176 if (fHeader) delete fHeader;
177 fHeader = new AliGenCocktailEventHeader("Beamgas Header");
178
179 TObjArray *partArray = gAlice->GetMCApp()->Particles();
180 AliStack *stack = gAlice->GetRunLoader()->Stack();
181
182 for (Int_t l = 0; l < numbg; l++) {
183 Vertex();
184 sign = (gRandom->Rndm() < 0.5)? -1. : 1.;
185 eventVertex.Set(3);
186 for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
187
188 if (!fRandom) {
189 Int_t igen=0;
190 while((entry = (AliGenCocktailEntry*)next())) {
191 if (fUsePerEventRate && (gRandom->Rndm() > entry->Rate())) continue;
192
193 igen++;
194 if (igen ==1) {
195 entry->SetFirst(lastpart);
196 } else {
197 entry->SetFirst((partArray->GetEntriesFast())+1);
198 }
199 //
200 // Handle case in which current generator needs collision geometry from previous generator
201 //
202 gen = entry->Generator();
203 if (gen->NeedsCollisionGeometry())
204 {
205 if (preventry && preventry->Generator()->ProvidesCollisionGeometry())
206 {
207 gen->SetCollisionGeometry(preventry->Generator()->CollisionGeometry());
208 } else {
209 Fatal("Generate()", "No Collision Geometry Provided");
210 }
211 }
212 entry->Generator()->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
213 entry->Generator()->Generate();
214 entry->SetLast(partArray->GetEntriesFast());
215 preventry = entry;
216 }
217 } else {
218 //
219 // Select a generator randomly
220 //
221 Int_t i;
222 Float_t p0 = gRandom->Rndm();
223
224 for (i = 0; i < fNGenerators; i++) {
225 if (p0 < fProb[i]) break;
226 }
227
228 entry = (AliGenCocktailEntry*) fEntries->At(i);
229 entry->SetFirst(lastpart);
230 gen = entry->Generator();
231 entry->Generator()->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
232 entry->Generator()->Generate();
233 entry->SetLast(partArray->GetEntriesFast());
234 }
235
236 for (int k = lastpart; k < stack->GetNprimary(); k++) {
237 stack->Particle(k)->ProductionVertex(v);
238 v[2] *= sign;
239 v[3] = fItime; // ??? +=
240 stack->Particle(k)->SetProductionVertex(v);
241 stack->Particle(k)->Momentum(v);
242 v[2] *= sign;
243 stack->Particle(k)->SetMomentum(v);
244 }
245
246 lastpart = stack->GetNprimary();
247 next.Reset();
248 }
249 // Event Vertex
250 fHeader->SetPrimaryVertex(eventVertex);
251 // fHeader->SetNvertex(numbg);
252 if (fContainer) {
253 fContainer->AddHeader(fHeader);
254 } else {
255 gAlice->SetGenEventHeader(fHeader);
256 }
257}