corrected TRD/TOF MV position
[u/mrichter/AliRoot.git] / THijing / 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 <TList.h>
31 #include <TObjArray.h>
32 #include <TParticle.h>
33
34 #include "AliGenBeamGasNew.h"
35 #include "AliGenCocktail.h" 
36 #include "AliGenCocktailEntry.h"
37 #include "AliGenCocktailEventHeader.h"
38 #include "AliGenHijing.h" 
39 #include "AliCollisionGeometry.h"
40 #include "AliRun.h"
41 #include "AliMC.h"
42 #include "AliStack.h"
43 #include "AliLog.h"
44
45 #include "AliHeader.h"
46 #include "AliGenEventHeader.h"
47
48 ClassImp(AliGenBeamGasNew)
49
50 Float_t EtaToTheta(Float_t arg) { return (180./TMath::Pi())*2.*atan(exp(-arg)); }
51
52 AliGenBeamGasNew::AliGenBeamGasNew() : 
53     AliGenCocktail(), 
54     fItime(0), 
55     fTwindow(88e-6), 
56     fZwindow(2000), 
57     fRate(1e3)
58 {
59   // Default constructor
60 }
61
62 AliGenBeamGasNew::~AliGenBeamGasNew()
63 {
64 }
65
66 void AliGenBeamGasNew::SetTimeWindow(Float_t twindow) { fTwindow = twindow; }
67
68 bool AliGenBeamGasNew::SetRate(Float_t rate) {
69 // Set the Rate
70   if (rate >= 0) {  
71     fRate = rate; 
72     return true;
73   } else
74     return false;
75 }
76
77 void AliGenBeamGasNew::SetZWindow(Float_t zwindow) { fZwindow = zwindow; }
78
79 void 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;
87   fSRandom = kTRUE;
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);
93     gen->SetReferenceFrame("LAB");
94     gen->SetEnergyCMS(7000);
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);
108     gen->SetRandomPz(kFALSE);
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())) {
118     if (fStack)  entry->Generator()->SetStack(fStack);
119     entry->Generator()->Init();
120   }  
121   
122   next.Reset();
123   
124   if (fSRandom) {
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
143 void AliGenBeamGasNew::VertexInternal()
144 {
145   // calculation of the interaction vertex for the beam gas event
146   // both spatial and time coordinate are adjusted.
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));
159 }
160
161 void AliGenBeamGasNew::Generate()
162 {
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   
169   TIter next(fEntries);
170   AliGenCocktailEntry *entry = 0;
171   AliGenCocktailEntry *preventry = 0;
172   AliGenerator* gen = 0;
173   TLorentzVector v;
174   TArrayF eventVertex;
175   Int_t lastpart=0;
176   Float_t sign;
177   
178   if (fHeader) delete fHeader;
179   fHeader = new AliGenCocktailEventHeader("Beamgas Header");
180   
181   const TObjArray *partArray = gAlice->GetMCApp()->Particles();
182   AliStack *stack = AliRunLoader::Instance()->Stack();
183   
184   for (Int_t l = 0; l < numbg; l++) {
185     Vertex();
186     sign = (gRandom->Rndm() < 0.5)? -1. : 1.;
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) {
197           entry->SetFirst(lastpart);
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       //
223       Int_t i = 0;
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);
231       entry->SetFirst(lastpart);
232       gen = entry->Generator();
233       gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
234       gen->Generate();
235       entry->SetLast(partArray->GetEntriesFast());
236     } 
237     
238     for (Int_t k = lastpart; k < stack->GetNprimary(); k++) {
239       stack->Particle(k)->ProductionVertex(v);
240       v[2] *= sign;
241       eventVertex[2] *= sign;
242       v[3] = fItime;  // ??? +=
243       stack->Particle(k)->SetProductionVertex(v);
244       stack->Particle(k)->Momentum(v);
245       v[2] *= sign;
246       stack->Particle(k)->SetMomentum(v);
247     }
248
249     ((AliGenEventHeader*) fHeader->GetHeaders()->Last())->SetPrimaryVertex(eventVertex);
250
251     lastpart = stack->GetNprimary();
252     next.Reset();
253   } 
254
255   fHeader->CalcNProduced();
256
257   if (fContainer) {
258     fContainer->AddHeader(fHeader);
259   } else {
260     gAlice->SetGenEventHeader(fHeader); 
261   }
262 }