Do not include from subdirectories
[u/mrichter/AliRoot.git] / THijing / AliGenBeamGasNew.cxx
CommitLineData
520b9406 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
e3bd319d 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
520b9406 25// in the same way as to AliGenCocktail before calling
26// the Init()
27// Author: Jochen Klein
e3bd319d 28//-------------------------------------------------------------
29
97bd035c 30#include <TList.h>
31#include <TObjArray.h>
520b9406 32#include <TParticle.h>
97bd035c 33
e3bd319d 34#include "AliGenBeamGasNew.h"
35#include "AliGenCocktail.h"
36#include "AliGenCocktailEntry.h"
e3bd319d 37#include "AliGenCocktailEventHeader.h"
c93255fe 38#include "AliGenHijing.h"
e3bd319d 39#include "AliCollisionGeometry.h"
40#include "AliRun.h"
41#include "AliMC.h"
42#include "AliStack.h"
97bd035c 43#include "AliLog.h"
e3bd319d 44
520b9406 45#include "AliHeader.h"
46#include "AliGenEventHeader.h"
47
e3bd319d 48ClassImp(AliGenBeamGasNew)
49
50Float_t EtaToTheta(Float_t arg) { return (180./TMath::Pi())*2.*atan(exp(-arg)); }
51
52AliGenBeamGasNew::AliGenBeamGasNew() :
53 AliGenCocktail(),
520b9406 54 fItime(0),
55 fTwindow(88e-6),
56 fZwindow(2000),
57 fRate(1e3)
e3bd319d 58{
520b9406 59 // Default constructor
e3bd319d 60}
61
62AliGenBeamGasNew::~AliGenBeamGasNew()
63{
e3bd319d 64}
65
520b9406 66void AliGenBeamGasNew::SetTimeWindow(Float_t twindow) { fTwindow = twindow; }
67
68bool AliGenBeamGasNew::SetRate(Float_t rate) {
904a753d 69// Set the Rate
520b9406 70 if (rate >= 0) {
71 fRate = rate;
72 return true;
73 } else
74 return false;
e3bd319d 75}
76
520b9406 77void AliGenBeamGasNew::SetZWindow(Float_t zwindow) { fZwindow = zwindow; }
e3bd319d 78
79void AliGenBeamGasNew::Init()
80{
81 // Initialisation of the class
82 // if no generators were added before calling Init()
83 // AliGenHijing is added configured for pO
84
85 fVertexSmear = kPerEvent;
86 fVertexSource = kInternal;
f0e207ce 87 fSRandom = kTRUE;
e3bd319d 88
89 // Adding default underlying event in case none was specified
90 // p-O-collision at 7 TeV (fixed target)
91 if (!fEntries) {
92 AliGenHijing *gen = new AliGenHijing(-1);
e3bd319d 93 gen->SetReferenceFrame("LAB");
520b9406 94 gen->SetEnergyCMS(7000);
e3bd319d 95 gen->SetProjectile("P",1,1);
96 gen->SetTarget("A",16,8);
97 gen->SetPhiRange(0,360);
98 Float_t thmin = EtaToTheta(8);
99 Float_t thmax = EtaToTheta(-8);
100 gen->SetThetaRange(thmin, thmax);
101 gen->SetOrigin(0,0,0);
102 gen->SetSigma(0,0,0);
103 gen->KeepFullEvent();
104 gen->SetShadowing(1);
105 gen->SetDecaysOff(1);
106 gen->SetSpectators(1);
107 gen->SetSelectAll(1);
520b9406 108 gen->SetRandomPz(kFALSE);
e3bd319d 109 gen->SetPileUpTimeWindow(0.);
110 AddGenerator(gen,"Hijing pO",1);
111 }
112
113 TIter next(fEntries);
114 AliGenCocktailEntry *entry;
115
116 // Loop over generators and initialize
117 while((entry = (AliGenCocktailEntry*)next())) {
520b9406 118 if (fStack) entry->Generator()->SetStack(fStack);
119 entry->Generator()->Init();
e3bd319d 120 }
121
122 next.Reset();
123
f0e207ce 124 if (fSRandom) {
e3bd319d 125 fProb.Set(fNGenerators);
126 next.Reset();
127 Float_t sum = 0.;
128 while((entry = (AliGenCocktailEntry*)next())) {
129 sum += entry->Rate();
130 }
131
132 next.Reset();
133 Int_t i = 0;
134 Float_t psum = 0.;
135 while((entry = (AliGenCocktailEntry*)next())) {
136 psum += entry->Rate() / sum;
137 fProb[i++] = psum;
138 }
139 }
140 next.Reset();
141}
142
143void AliGenBeamGasNew::VertexInternal()
144{
145 // calculation of the interaction vertex for the beam gas event
146 // both spatial and time coordinate are adjusted.
520b9406 147
148 Float_t random[2];
149 Float_t nbunch;
150 Float_t bunch;
151
152 gRandom->RndmArray(2,random);
153 fVertex[0] = fVertex[1] = 0;
154 fVertex[2] = random[0] * 2 * fZwindow - fZwindow;
155 nbunch = fTwindow/25e-9;
156 bunch = floor(2*nbunch * random[1] - nbunch);
157 fItime = fVertex[2] / 100 / 3e8 + bunch * 25e-9;
158 AliInfo(Form("fItime: %13.3e \n", fItime));
e3bd319d 159}
160
161void AliGenBeamGasNew::Generate()
162{
520b9406 163//
164// Generate the collisions for one event
165//
166 Int_t numbg = gRandom->Poisson(fRate*fZwindow/100*2*fTwindow);
167 if (numbg == 0) return;
168
e3bd319d 169 TIter next(fEntries);
170 AliGenCocktailEntry *entry = 0;
171 AliGenCocktailEntry *preventry = 0;
172 AliGenerator* gen = 0;
520b9406 173 TLorentzVector v;
174 TArrayF eventVertex;
175 Int_t lastpart=0;
176 Float_t sign;
e3bd319d 177
178 if (fHeader) delete fHeader;
179 fHeader = new AliGenCocktailEventHeader("Beamgas Header");
180
d94af0c1 181 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
33c3c91a 182 AliStack *stack = AliRunLoader::Instance()->Stack();
e3bd319d 183
520b9406 184 for (Int_t l = 0; l < numbg; l++) {
e3bd319d 185 Vertex();
520b9406 186 sign = (gRandom->Rndm() < 0.5)? -1. : 1.;
e3bd319d 187 eventVertex.Set(3);
188 for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
189
190 if (!fRandom) {
191 Int_t igen=0;
192 while((entry = (AliGenCocktailEntry*)next())) {
193 if (fUsePerEventRate && (gRandom->Rndm() > entry->Rate())) continue;
194
195 igen++;
196 if (igen ==1) {
520b9406 197 entry->SetFirst(lastpart);
e3bd319d 198 } else {
199 entry->SetFirst((partArray->GetEntriesFast())+1);
200 }
201 //
202 // Handle case in which current generator needs collision geometry from previous generator
203 //
204 gen = entry->Generator();
205 if (gen->NeedsCollisionGeometry())
206 {
207 if (preventry && preventry->Generator()->ProvidesCollisionGeometry())
208 {
209 gen->SetCollisionGeometry(preventry->Generator()->CollisionGeometry());
210 } else {
211 Fatal("Generate()", "No Collision Geometry Provided");
212 }
213 }
214 entry->Generator()->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
215 entry->Generator()->Generate();
216 entry->SetLast(partArray->GetEntriesFast());
217 preventry = entry;
218 }
219 } else {
220 //
221 // Select a generator randomly
222 //
9bce4ac1 223 Int_t i = 0;
e3bd319d 224 Float_t p0 = gRandom->Rndm();
225
226 for (i = 0; i < fNGenerators; i++) {
227 if (p0 < fProb[i]) break;
228 }
229
230 entry = (AliGenCocktailEntry*) fEntries->At(i);
520b9406 231 entry->SetFirst(lastpart);
e3bd319d 232 gen = entry->Generator();
0dca4638 233 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
234 gen->Generate();
e3bd319d 235 entry->SetLast(partArray->GetEntriesFast());
236 }
237
9bce4ac1 238 for (Int_t k = lastpart; k < stack->GetNprimary(); k++) {
e3bd319d 239 stack->Particle(k)->ProductionVertex(v);
520b9406 240 v[2] *= sign;
9bce4ac1 241 eventVertex[2] *= sign;
520b9406 242 v[3] = fItime; // ??? +=
e3bd319d 243 stack->Particle(k)->SetProductionVertex(v);
520b9406 244 stack->Particle(k)->Momentum(v);
245 v[2] *= sign;
246 stack->Particle(k)->SetMomentum(v);
e3bd319d 247 }
520b9406 248
9bce4ac1 249 ((AliGenEventHeader*) fHeader->GetHeaders()->Last())->SetPrimaryVertex(eventVertex);
250
520b9406 251 lastpart = stack->GetNprimary();
e3bd319d 252 next.Reset();
520b9406 253 }
9bce4ac1 254
97bd035c 255 fHeader->CalcNProduced();
256
520b9406 257 if (fContainer) {
258 fContainer->AddHeader(fHeader);
259 } else {
260 gAlice->SetGenEventHeader(fHeader);
261 }
e3bd319d 262}