SetSeed implementation
[u/mrichter/AliRoot.git] / EVGEN / AliGenCocktail.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 // Container class for AliGenerator through recursion.
19 // Container is itself an AliGenerator.
20 // What is stored are not the pointers to the generators directly but to objects of type
21 // AliGenCocktail entry.   
22 // The class provides also iterator functionality.  
23 // Author: andreas.morsch@cern.ch 
24 //
25
26 #include <TList.h>
27 #include <TObjArray.h>
28 #include <TFormula.h>
29
30 #include "AliGenCocktail.h"
31 #include "AliGenCocktailEntry.h"
32 #include "AliCollisionGeometry.h"
33 #include "AliRun.h"
34 #include "AliLog.h"
35 #include "AliMC.h"
36 #include "AliGenCocktailEventHeader.h"
37
38 ClassImp(AliGenCocktail)
39
40 AliGenCocktail::AliGenCocktail()
41     :AliGenerator(), 
42      fNGenerators(0),
43      fTotalRate(0.),
44      fSRandom(kFALSE),
45      fUsePerEventRate(kFALSE),
46      fProb(0),
47      fEntries(0),
48      flnk1(0),
49      flnk2(0), 
50      fHeader(0)
51 {
52 // Constructor
53     fName = "Cocktail";
54     fTitle= "Particle Generator using cocktail of generators";
55 }
56
57 AliGenCocktail::~AliGenCocktail()
58 {
59 // Destructor
60     delete fEntries;
61     fEntries = 0;
62     //    delete fHeader; // It is removed in AliRunLoader
63     fHeader = 0;
64 }
65
66 void AliGenCocktail::
67 AddGenerator(AliGenerator *Generator, const char* Name, Float_t RateExp, TFormula* formula, Int_t ntimes)
68 {
69 //
70 // Add a generator to the list 
71 // First check that list exists
72     if (!fEntries) fEntries = new TList();
73     fTotalRate += RateExp;
74 //
75 //  Forward parameters to the new generator
76     if(TestBit(kPtRange) && !(Generator->TestBit(kPtRange)) && !(Generator->TestBit(kMomentumRange))) 
77         Generator->SetPtRange(fPtMin,fPtMax);
78     if(TestBit(kMomentumRange) && !(Generator->TestBit(kPtRange)) && !(Generator->TestBit(kMomentumRange)))
79         Generator->SetMomentumRange(fPMin,fPMax);
80     
81     if (TestBit(kYRange) && !(Generator->TestBit(kYRange)))
82         Generator->SetYRange(fYMin,fYMax);
83     if (TestBit(kPhiRange) && !(Generator->TestBit(kPhiRange)))
84         Generator->SetPhiRange(fPhiMin*180/TMath::Pi(),fPhiMax*180/TMath::Pi());
85     if (TestBit(kThetaRange) && !(Generator->TestBit(kThetaRange)) && !(Generator->TestBit(kEtaRange)))
86         Generator->SetThetaRange(fThetaMin*180/TMath::Pi(),fThetaMax*180/TMath::Pi());
87     if (!(Generator->TestBit(kVertexRange))) {
88         Generator->SetOrigin(fOrigin[0], fOrigin[1], fOrigin[2]);
89         Generator->SetSigma(fOsigma[0], fOsigma[1], fOsigma[2]);
90         Generator->SetVertexSmear(fVertexSmear);
91         Generator->SetVertexSource(kContainer);
92     }
93     Generator->SetTrackingFlag(fTrackIt);
94     Generator->SetContainer(this);
95
96         
97 //
98 //  Add generator to list   
99     char theName[256];
100     snprintf(theName, 256, "%s_%d",Name, fNGenerators);
101     Generator->SetName(theName);
102
103     AliGenCocktailEntry *entry = 
104         new AliGenCocktailEntry(Generator, Name, RateExp);
105     if (formula) entry->SetFormula(formula);  
106     entry->SetNTimes(ntimes);
107      fEntries->Add(entry);
108      fNGenerators++;
109      flnk1 = 0;
110      flnk2 = 0;
111      fSRandom  = kFALSE;
112      fHeader  = 0;
113 }
114
115   void AliGenCocktail::Init()
116 {
117 // Initialisation
118     TIter next(fEntries);
119     AliGenCocktailEntry *entry;
120     //
121     // Loop over generators and initialize
122     while((entry = (AliGenCocktailEntry*)next())) {
123         if (fStack)  entry->Generator()->SetStack(fStack);
124         entry->Generator()->SetSeed(fSeed);
125         entry->Generator()->Init();
126     }  
127
128     next.Reset();
129
130     if (fSRandom) {
131         fProb.Set(fNGenerators);
132         next.Reset();
133         Float_t sum = 0.;
134         while((entry = (AliGenCocktailEntry*)next())) {
135             sum += entry->Rate();
136         } 
137
138         next.Reset();
139         Int_t i = 0;
140         Float_t psum = 0.;
141         while((entry = (AliGenCocktailEntry*)next())) {
142             psum +=  entry->Rate() / sum;
143             fProb[i++] = psum;
144         }
145     }
146         next.Reset();
147 }
148
149   void AliGenCocktail::FinishRun()
150 {
151 // Initialisation
152     TIter next(fEntries);
153     AliGenCocktailEntry *entry;
154     //
155     // Loop over generators and initialize
156     while((entry = (AliGenCocktailEntry*)next())) {
157         entry->Generator()->FinishRun();
158     }  
159 }
160
161  void AliGenCocktail::Generate()
162 {
163 //
164 // Generate event 
165     TIter next(fEntries);
166     AliGenCocktailEntry *entry = 0;
167     AliGenCocktailEntry *preventry = 0;
168     AliGenCocktailEntry *collentry = 0;
169     AliGenerator* gen = 0;
170     if (fHeader) delete fHeader;
171
172     
173     fHeader = new AliGenCocktailEventHeader("Cocktail Header");
174
175     const TObjArray *partArray = gAlice->GetMCApp()->Particles();
176
177 //
178 //  Generate the vertex position used by all generators
179 //    
180     if(fVertexSmear == kPerEvent) Vertex();
181
182     TArrayF eventVertex;
183     eventVertex.Set(3);
184     for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
185
186     if (!fSRandom) {
187         //
188         // Loop over generators and generate events
189         Int_t igen   = 0;
190         while((entry = (AliGenCocktailEntry*)next())) {
191           Int_t ntimes = entry->NTimes();
192           if (fUsePerEventRate && (gRandom->Rndm() > entry->Rate())) continue;
193           
194           igen++;
195           if (igen ==1) {
196             entry->SetFirst(0);
197           } else {
198             entry->SetFirst((partArray->GetEntriesFast())+1);
199           }
200           gen = entry->Generator();
201           if (gen->ProvidesCollisionGeometry()) collentry = entry; 
202           //
203           //      Handle case in which current generator needs collision geometry from previous generator
204           //
205           if (gen->NeedsCollisionGeometry() && (entry->Formula() == 0))
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           //
215           //      Number of signals is calculated from Collision Geometry
216           //      and entry with given centrality bin is selected
217           //
218           if (entry->Formula() != 0)
219             {
220               if (!collentry) {
221                 Fatal("Generate()", "No Collision Geometry Provided");
222                 return;
223               }
224               AliCollisionGeometry* coll = (collentry->Generator())->CollisionGeometry();
225               Float_t b  = coll->ImpactParameter();
226               Int_t nsig = Int_t(entry->Formula()->Eval(b));
227               Int_t bin = entry->Bin() - 100;
228               if (bin > 0) {
229                 if (bin != nsig) continue;
230               } else {
231                 if (nsig < 1) nsig = 1;
232                 AliInfo(Form("Signal Events %13.3f %5d %5d\n", b, coll->HardScatters(), nsig));
233                 ntimes = nsig;
234               }
235             }
236           gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), fTime);
237           
238           gen->GenerateN(ntimes);
239           entry->SetLast(partArray->GetEntriesFast());
240           preventry = entry;
241         }
242     } else if (fSRandom) {
243         //
244         // Select a generator randomly
245         //
246         Int_t i;
247         Float_t p0 =  gRandom->Rndm();
248
249         for (i = 0; i < fNGenerators; i++) {
250             if (p0 < fProb[i]) break;
251         }
252
253         entry = (AliGenCocktailEntry*) fEntries->At(i);
254         entry->SetFirst(0);
255         gen = entry->Generator();
256         gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), fTime);
257         gen->Generate();
258         entry->SetLast(partArray->GetEntriesFast());
259     } 
260     
261     next.Reset();
262
263     // Event Vertex
264     fHeader->SetPrimaryVertex(eventVertex);
265     fHeader->CalcNProduced();
266     if (fContainer) {
267       fHeader->SetName(fName);
268       fContainer->AddHeader(fHeader);
269     } else {
270       gAlice->SetGenEventHeader(fHeader);       
271     }
272 }
273
274 void AliGenCocktail::SetVertexSmear(VertexSmear_t smear)
275 {
276 // Set vertex smearing and propagate it to the generators
277
278   AliGenerator::SetVertexSmear(smear);
279   TIter next(fEntries);
280   while (AliGenCocktailEntry* entry = (AliGenCocktailEntry*)next()) {
281     entry->Generator()->SetVertexSmear(smear);
282   }
283 }
284
285 AliGenCocktailEntry *  AliGenCocktail::FirstGenerator()
286 {
287 // Iterator over generators: Initialisation
288     flnk1 = fEntries->FirstLink();
289     if (flnk1) {
290         return (AliGenCocktailEntry*) (flnk1->GetObject());
291     } else {
292         return 0;
293     }
294 }
295
296 AliGenCocktailEntry*  AliGenCocktail::NextGenerator()
297 {
298 // Iterator over generators: Increment
299     flnk1 = flnk1->Next();
300     if (flnk1) {
301         return (AliGenCocktailEntry*) (flnk1->GetObject());
302     } else {
303         return 0;
304     }
305 }
306
307 void AliGenCocktail::
308 FirstGeneratorPair(AliGenCocktailEntry*& e1, AliGenCocktailEntry*& e2)
309 {
310 // Iterator over generator pairs: Initialisation
311     flnk2 = flnk1 = fEntries->FirstLink();
312     if (flnk1) {
313         e2 = e1 = (AliGenCocktailEntry*) (flnk1->GetObject());
314     } else {
315         e2= e1 = 0;
316     }
317 }
318
319 void AliGenCocktail::
320 NextGeneratorPair(AliGenCocktailEntry*& e1, AliGenCocktailEntry*& e2)
321 {
322 // Iterator over generators: Increment
323     flnk2 = flnk2->Next();
324     if (flnk2) {
325         e1 = (AliGenCocktailEntry*) (flnk1->GetObject());
326         e2 = (AliGenCocktailEntry*) (flnk2->GetObject());       
327     } else {
328         flnk2 = flnk1 = flnk1->Next();
329         if (flnk1) {
330             e1 = (AliGenCocktailEntry*) (flnk1->GetObject());
331             e2 = (AliGenCocktailEntry*) (flnk2->GetObject());
332         } else {
333             e1=0;
334             e2=0;
335         }
336     }
337 }
338
339 void AliGenCocktail::AddHeader(AliGenEventHeader* header)
340 {
341 // Add a header to the list 
342     if (fHeader) fHeader->AddHeader(header);
343 }
344
345
346
347