Store total number of produced particles in the header.
[u/mrichter/AliRoot.git] / EVGEN / AliGenBeamGasNew.cxx
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 
6 //  in the same as to AliGenCocktail before calling the Init()
7 //-------------------------------------------------------------
8
9 /* $Id$ */
10
11 #include <TList.h>
12 #include <TObjArray.h>
13 #include "TParticle.h"
14
15 #include "AliGenBeamGasNew.h"
16 #include "AliGenCocktail.h" 
17 #include "AliGenCocktailEntry.h"
18 #include "AliGenCocktailEventHeader.h"
19 #include "../THijing/AliGenHijing.h" 
20 #include "AliCollisionGeometry.h"
21 #include "AliRun.h"
22 #include "AliMC.h"
23 #include "AliStack.h"
24 #include "AliLog.h"
25
26 ClassImp(AliGenBeamGasNew)
27
28 Float_t EtaToTheta(Float_t arg) { return (180./TMath::Pi())*2.*atan(exp(-arg)); }
29
30 AliGenBeamGasNew::AliGenBeamGasNew() : 
31     AliGenCocktail(), 
32     fTwindow(88.e-6)
33 {
34     // Default constructor
35 }
36
37 AliGenBeamGasNew::~AliGenBeamGasNew()
38 {
39     // Destructor
40 }
41
42 AliGenBeamGasNew::AliGenBeamGasNew(const AliGenBeamGasNew& /*rhs*/) : AliGenCocktail()
43 {
44   AliFatal("Copy Constructor not yet implemented!");
45 }
46
47
48 void 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);
77     gen->SetRandomPz(kTRUE);
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())) {
87       if (fStack)  entry->Generator()->SetStack(fStack);
88       entry->Generator()->Init();
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
112 void AliGenBeamGasNew::VertexInternal()
113 {
114   // calculation of the interaction vertex for the beam gas event
115   // both spatial and time coordinate are adjusted.
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));
129 }
130
131 void AliGenBeamGasNew::Generate()
132 {
133   TIter next(fEntries);
134   AliGenCocktailEntry *entry = 0;
135   AliGenCocktailEntry *preventry = 0;
136   AliGenerator* gen = 0;
137   
138   if (fHeader) delete fHeader;
139   fHeader = new AliGenCocktailEventHeader("Beamgas Header");
140   
141   TObjArray *partArray = gAlice->GetMCApp()->Particles();
142   
143   TArrayF eventVertex;
144   int numbg = 0; // variable for counting how many bg interactions happened
145   
146   do {
147     numbg++;
148     Vertex();
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) {
159           entry->SetFirst(0);
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);
193       entry->SetFirst(0);
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     
200     AliStack *stack = gAlice->GetRunLoader()->Stack();
201     TLorentzVector v;
202     for (int k = 0; k < stack->GetNprimary(); k++) {
203       stack->Particle(k)->ProductionVertex(v);
204       v[3] += fItime;
205       stack->Particle(k)->SetProductionVertex(v);
206     }
207     
208     next.Reset();
209   } while (gRandom->Rndm() < 1. / 5.); // must be adjusted to the rate of bg interactions
210   // Event Vertex
211   fHeader->SetPrimaryVertex(eventVertex);
212   // Total number of produced an stored particles
213   fHeader->CalcNProduced();
214
215   gAlice->SetGenEventHeader(fHeader); 
216 }