2259e952a810758a32677c9a5553d7bb7ffb6fd8
[u/mrichter/AliRoot.git] / EVGEN / AliGenBeamGasNew.cxx
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
45 ClassImp(AliGenBeamGasNew)
46
47 Float_t EtaToTheta(Float_t arg) { return (180./TMath::Pi())*2.*atan(exp(-arg)); }
48
49 AliGenBeamGasNew::AliGenBeamGasNew() : 
50     AliGenCocktail(), 
51     fItime(0), 
52     fTwindow(88e-6), 
53     fZwindow(2000), 
54     fRate(1e3)
55 {
56   printf("AliGenBeamGasNew instatiated!\n");
57 }
58
59 AliGenBeamGasNew::~AliGenBeamGasNew()
60 {
61   delete fEntries;
62   fEntries = 0;
63   fHeader = 0;
64 }
65
66 void AliGenBeamGasNew::SetTimeWindow(Float_t twindow) { fTwindow = twindow; }
67
68 bool AliGenBeamGasNew::SetRate(Float_t rate) {
69   if (rate >= 0) {  
70     fRate = rate; 
71     return true;
72   } else
73     return false;
74 }
75
76 void AliGenBeamGasNew::SetZWindow(Float_t zwindow) { fZwindow = zwindow; }
77
78 void 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
142 void 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
159 void 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 }