[bug #84578] Request to extend AliGenBox for using YRange
[u/mrichter/AliRoot.git] / EVGEN / AliGenCocktailAfterBurner.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 // Container class for AliGenerator and AfterBurners 
20 // (which are AliGenerators as well) through recursion.
21 // The container is itself an AliGenerator a
22 // what is stored are not the pointers to the generators directly 
23 // but to objects of type
24 // AliGenCocktailAfterBurner entry.   
25 // The class provides also iterator functionality.  
26 // Author: andreas.morsch@cern.ch and piotr.skowronski@cern.ch
27 //
28 // 24.09.2001  Piotr Skowronski
29 //             debug -> gDebug,
30 //             fNEvents replaced with AliRunLoader::GetNumberOfEvents()
31 //
32
33
34 #include <Riostream.h>
35
36 #include <TList.h>
37 #include <TObjArray.h>
38 #include <TParticle.h>
39
40 #include "AliGenCocktailAfterBurner.h"
41 #include "AliGenCocktailEntry.h"
42 #include "AliGenCocktailEventHeader.h"
43 #include "AliCollisionGeometry.h"
44 #include "AliStack.h"
45 #include "AliMC.h"
46 #include "AliRun.h"
47
48
49 ClassImp(AliGenCocktailAfterBurner)
50 /*********************************************************************/ 
51 /*********************************************************************/ 
52
53     AliGenCocktailAfterBurner::AliGenCocktailAfterBurner():
54         fNAfterBurners(0),
55         fAfterBurnerEntries(0),
56         fGenerationDone(kFALSE),
57         fInternalStacks(0),
58         fCollisionGeometries(0),
59         fHeaders(0),
60         fCurrentEvent(0),
61         fActiveStack(0),
62         fActiveEvent(-1),
63         fCurrentGenerator(0),
64         fNBgEvents(0)
65 {
66 // Constructor
67     if (gDebug > 0) 
68         cout<<"AliGenCocktailAfterBurner::AliGenCocktailAfterBurner()"<<endl;
69     SetName("AliGenCocktailAfterBurner");
70     SetTitle("AliGenCocktailAfterBurner");
71 }
72
73 /*********************************************************************/ 
74
75 AliGenCocktailAfterBurner::~AliGenCocktailAfterBurner()
76   {
77 //destructor
78
79     if (fInternalStacks) //delete stacks
80      { 
81        fInternalStacks->SetOwner();
82        delete fInternalStacks;
83     }
84     if (fAfterBurnerEntries) delete fAfterBurnerEntries; //delete entries
85     Int_t numberOfEvents = AliRunLoader::Instance()->GetNumberOfEventsPerRun();
86     if (fCollisionGeometries) {
87       for (Int_t i = 0; i < (numberOfEvents + fNBgEvents); i++)
88         if (fCollisionGeometries[i]) delete fCollisionGeometries[i];
89       delete[] fCollisionGeometries;
90     }
91     if (fHeaders) {
92       for (Int_t i = 0; i < (numberOfEvents + fNBgEvents); i++)
93         if (fHeaders[i]) delete fHeaders[i];
94       delete[] fHeaders;
95     }
96   }
97 /*********************************************************************/ 
98 /*********************************************************************/ 
99
100 void AliGenCocktailAfterBurner::
101 AddAfterBurner(AliGenerator *AfterBurner, char* Name, Float_t RateExp)
102 {
103 //
104 //  Forward parameters to the new AfterBurner
105     
106     if (gDebug>0)cout<<"AliGenCocktailAfterBurner::AddAfterBurner  Named "<<Name<<endl;
107
108     if(TestBit(kPtRange) && !(AfterBurner->TestBit(kPtRange)) && !(AfterBurner->TestBit(kMomentumRange))) 
109         AfterBurner->SetPtRange(fPtMin,fPtMax);
110     if(TestBit(kMomentumRange) && !(AfterBurner->TestBit(kPtRange)) && !(AfterBurner->TestBit(kMomentumRange)))
111         AfterBurner->SetMomentumRange(fPMin,fPMax);
112     
113     if (TestBit(kYRange) && !(AfterBurner->TestBit(kYRange)))
114         AfterBurner->SetYRange(fYMin,fYMax);
115     if (TestBit(kPhiRange) && !(AfterBurner->TestBit(kPhiRange)))
116         AfterBurner->SetPhiRange(fPhiMin*180/TMath::Pi(),fPhiMax*180/TMath::Pi());
117     if (TestBit(kThetaRange) && !(AfterBurner->TestBit(kThetaRange)) && !(AfterBurner->TestBit(kEtaRange)))
118         AfterBurner->SetThetaRange(fThetaMin*180/TMath::Pi(),fThetaMax*180/TMath::Pi());
119     if (TestBit(kVertexRange) && !(AfterBurner->TestBit(kVertexRange))) {
120         AfterBurner->SetOrigin(fOrigin[0], fOrigin[1], fOrigin[2]);
121         AfterBurner->SetSigma(fOsigma[0], fOsigma[1], fOsigma[2]);
122         AfterBurner->SetVertexSmear(fVertexSmear);
123         AfterBurner->SetVertexSource(kContainer);
124     }
125     AfterBurner->SetTrackingFlag(fTrackIt);
126     //AfterBurner->SetContainer(this);
127
128 //
129 //  Add AfterBurner to list   
130     
131     AliGenCocktailEntry *entry = 
132         new AliGenCocktailEntry(AfterBurner, Name, RateExp);
133     if (!fAfterBurnerEntries) fAfterBurnerEntries = new TList();
134     
135     fAfterBurnerEntries->Add(entry);
136     fNAfterBurners++;
137 //
138     
139 }
140 /*********************************************************************/ 
141 /*********************************************************************/ 
142
143 void AliGenCocktailAfterBurner::Init()
144 {
145 // Initialisation
146   Int_t numberOfEvents = AliRunLoader::Instance()->GetNumberOfEventsPerRun();
147     fGenerationDone = kFALSE;
148     if (fInternalStacks) //delete stacks
149      { 
150        fInternalStacks->SetOwner();
151        fInternalStacks->Delete(); //clean after previous generation cycle
152      }
153
154     if (fCollisionGeometries) {
155       for (Int_t i = 0; i < (numberOfEvents + fNBgEvents); i++)
156         if (fCollisionGeometries[i]) delete fCollisionGeometries[i];
157       delete[] fCollisionGeometries;
158     }
159     if (fHeaders) {
160       for (Int_t i = 0; i < (numberOfEvents + fNBgEvents); i++)
161         if (fHeaders[i]) delete fHeaders[i];
162       delete[] fHeaders;
163     }
164
165     this->AliGenCocktail::Init(); 
166     
167     if (gDebug>0) cout<<"AliGenCocktailAfterBurner::Init"<<endl;
168     TIter next(fAfterBurnerEntries);
169     AliGenCocktailEntry *entry;
170     //
171     // Loop over generators and initialize
172     while((entry = (AliGenCocktailEntry*)next())) {
173         entry->Generator()->Init();
174     }  
175 }
176 /*********************************************************************/ 
177 /*********************************************************************/ 
178
179 void AliGenCocktailAfterBurner::Generate()
180 {
181 //
182 // Generate event
183 //  Firsts runs each generator for all events
184 //  than after burners ones for each event
185 //
186 //  It generates and processes all events during
187 //  first call only.
188 //  In next calls it just returns already generated 
189 //  and processed events to the gAlice
190
191     if (gDebug>0)
192       cout<<"#####################################"<<endl
193           <<"#AliGenCocktailAfterBurner::Generate#"<<endl
194           <<"#####################################"<<endl;
195     // Initialize header
196     //
197     Int_t i; //iterator
198     AliStack * stack;
199     
200     if (fGenerationDone)
201     {//if generation is done (in first call) 
202      //just copy particles from the stack to the gAlice
203       SetTracks(++fCurrentEvent);
204       fHeader = fHeaders[fCurrentEvent];
205       gAlice->SetGenEventHeader(fHeader); 
206       cout<<"Returning event " << fCurrentEvent<<endl;
207       return;  
208     }
209     else
210     { //Here we are in the first call of the method
211         Int_t numberOfEvents = AliRunLoader::Instance()->GetNumberOfEventsPerRun();
212         cout << "Number of events per run" <<  numberOfEvents << endl;
213         TArrayF eventVertex;
214         eventVertex.Set(3 * (numberOfEvents + fNBgEvents));
215         fCurrentEvent=0;
216       //Create stacks
217         fInternalStacks      = new TObjArray(numberOfEvents + fNBgEvents); //Create array of internal stacks
218         fCollisionGeometries = new AliCollisionGeometry*[numberOfEvents + fNBgEvents]; //Create array of collision geometries
219         fHeaders             = new AliGenCocktailEventHeader*[numberOfEvents + fNBgEvents]; //Create array of headers   
220
221         for(i = 0; i < numberOfEvents + fNBgEvents; i++) 
222        {        
223            stack = new AliStack(10000);
224            stack->Reset();
225            fInternalStacks->Add(stack);
226            Vertex();
227            for (Int_t j = 0; j < 3; j++) eventVertex[3 * i +  j] = fVertex[j];
228            fHeaders[i] = new AliGenCocktailEventHeader();
229            fCollisionGeometries[i] = 0;
230        }
231 /*********************************************************************/ 
232       TIter next(fEntries);
233       AliGenCocktailEntry *entry;
234       AliGenCocktailEntry *e1;
235       AliGenCocktailEntry *e2;
236       const TObjArray *partArray;
237   //
238   // Loop over generators and generate events
239       Int_t igen=0;
240       while((entry = (AliGenCocktailEntry*)next())) 
241       {
242         igen++;
243         cout<<"Generator "<<igen<<"  : "<<entry->GetName()<<endl;
244 /***********************************************/
245 //First generator for all evenets, than second for all events, etc...
246         for(i = 0; i < numberOfEvents + fNBgEvents; i++) 
247         {  
248             cout<<"                  EVENT "<<i << endl;
249             stack = GetStack(i);
250             partArray = stack->Particles();
251             fCurrentGenerator = entry->Generator();
252             fCurrentGenerator->SetStack(stack);
253             if (igen ==1) 
254             {
255                 entry->SetFirst(0);
256             } 
257             else 
258             {
259                 entry->SetFirst((partArray->GetEntriesFast())+1);
260             }
261             // Set the vertex for the generator
262             Int_t ioff = 3 * i;
263             fCurrentGenerator->SetVertex(eventVertex.At(ioff), eventVertex.At(ioff + 1), eventVertex.At(ioff + 2));
264             fHeader = fHeaders[i];
265             // Set the vertex for the cocktail
266             TArrayF v(3);
267             for (Int_t j=0; j<3; j++) v[j] = eventVertex.At(ioff + j);
268             fHeader->SetPrimaryVertex(v);
269             // Generate event
270             fCurrentGenerator->Generate();
271             //
272             entry->SetLast(partArray->GetEntriesFast());
273             
274             if (fCurrentGenerator->ProvidesCollisionGeometry())  
275               fCollisionGeometries[i] = 
276                 new AliCollisionGeometry(*(fCurrentGenerator->CollisionGeometry()));
277         } // event loop
278 /***********************************************/
279       } // generator loop
280       next.Reset();
281       while((entry = (AliGenCocktailEntry*)next())) 
282         {
283           entry->PrintInfo();
284         }
285       for ( entry=FirstGenerator();entry;entry=NextGenerator() ) 
286         {
287           entry->PrintInfo();
288         }
289       for (FirstGeneratorPair(e1,e2); (e1&&e2); NextGeneratorPair(e1,e2) )
290         {
291           printf("\n -----------------------------");
292           e1->PrintInfo();
293           e2->PrintInfo();
294         }
295         
296         
297       /***********************************************/
298       /*******After Burners Processing****************/
299       /***********************************************/
300       TIter nextAfterBurner(fAfterBurnerEntries);
301       AliGenCocktailEntry *afterBurnerEntry;
302       Int_t iab =0; //number of current after burner / counter
303       
304       cout<<"\n\nRunning After Burners"<<endl;
305       while((afterBurnerEntry = (AliGenCocktailEntry*)nextAfterBurner()))
306         {
307           cout<<"After Burner "<<iab++<<"  :"<<afterBurnerEntry->GetName()<<endl;
308           fCurrentGenerator = afterBurnerEntry->Generator();
309           fCurrentGenerator->Generate();
310         }
311       cout<<endl<<"Finished. Processed "<<iab<<" After Burners"<<endl;
312
313       /***********************************************/
314       /***********************************************/
315       /***********************************************/       
316         
317       fGenerationDone=kTRUE; 
318       SetTracks(0); //copy event 0 to gAlice stack
319         
320 /*********************************************************************/
321       // Pass the header to gAlice
322       fHeader = fHeaders[0];
323       gAlice->SetGenEventHeader(fHeader); 
324     } //else generated
325 }
326 /*********************************************************************/
327 /*********************************************************************/ 
328
329 AliStack* AliGenCocktailAfterBurner::GetStack(Int_t n) const
330 {
331 //Returns the pointer to the N'th stack (event)
332   if( ( n<0 ) || ( n >= (GetNumberOfEvents()) ) )
333     {
334       Fatal("AliGenCocktailAfterBurner::GetStack","Asked for non existing stack (%d)",n);
335       return 0; 
336     }
337     return ((AliStack*) fInternalStacks->At(n) );
338 }
339
340 /*********************************************************************/ 
341 /*********************************************************************/ 
342
343 AliCollisionGeometry* AliGenCocktailAfterBurner::GetCollisionGeometry(Int_t n) const
344 {
345 //Returns the pointer to the N'th stack (event)
346   if( ( n<0 ) || ( n>=GetNumberOfEvents() ) )
347     {
348       Fatal("AliGenCocktailAfterBurner::GetCollisionGeometry","Asked for non existing stack (%d)",n);
349       return 0; 
350     }
351     return fCollisionGeometries[n];
352 }
353
354 /*********************************************************************/ 
355 /*********************************************************************/ 
356
357 void AliGenCocktailAfterBurner::SetActiveEventNumber(Int_t actev)
358 {
359 //Set Active Events Number and Active Stack
360 //There is only one active event number
361 //Made fo convinience of work with AfterBurners (HBT processor)
362
363     fActiveEvent = actev;
364     fActiveStack = GetStack(actev);
365 }
366 /*********************************************************************/ 
367 /*********************************************************************/ 
368
369 void AliGenCocktailAfterBurner::SetTracks(Int_t stackno)
370 {
371 //Method which copies tracks from given stack to the
372 //gAlice's stack
373     AliStack* instack = GetStack(stackno);
374     Int_t done;
375     Int_t parent; 
376     Int_t pdg;
377     Double_t px, py, pz, e, vx, vy, vz, tof, polx, poly, polz;
378     TMCProcess mech;
379     Int_t ntr;
380     Float_t weight;
381     TVector3 pol;
382     
383     TParticle * p;
384     Int_t n = instack->GetNtrack();
385     if (gDebug) 
386     {
387       cout<<"AliGenCocktailAfterBurner::SetTracks("<<stackno<<"). Number of particles is: "<<n<<"\n";
388     }
389     
390     for(Int_t i = 0; i < n; i++)
391     {
392         
393       p = instack->Particle(i);
394       done = !p->TestBit(kDoneBit);
395       parent = p->GetMother(0);
396       pdg = p->GetPdgCode();
397       px = p->Px();
398       py = p->Py();
399       pz = p->Pz();
400       e  = p->Energy();
401       vx = p->Vx();
402       vy = p->Vy();
403       vz = p->Vz();
404       tof = p->T();
405       p->GetPolarisation(pol);
406       polx = pol.X();
407       poly = pol.Y();
408       polz = pol.Z();
409       mech = AliGenCocktailAfterBurner::IntToMCProcess(p->GetUniqueID());
410       weight = p->GetWeight();
411
412       gAlice->GetMCApp()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,polx, poly, polz, mech, ntr, weight);
413
414       SetHighWaterMark(ntr) ; 
415
416     }
417 }
418 /*********************************************************************/ 
419 /*********************************************************************/ 
420
421 TMCProcess AliGenCocktailAfterBurner::IntToMCProcess(Int_t no)
422 {
423  //Mothod used to convert uniqueID (integer) to TMCProcess type
424     const TMCProcess kMCprocesses[kMaxMCProcess] = 
425     {
426      kPNoProcess, kPMultipleScattering, kPEnergyLoss, kPMagneticFieldL, 
427      kPDecay, kPPair, kPCompton, kPPhotoelectric, kPBrem, kPDeltaRay,
428      kPAnnihilation, kPHadronic, kPNoProcess, kPEvaporation, kPNuclearFission,
429      kPNuclearAbsorption, kPPbarAnnihilation, kPNCapture, kPHElastic, 
430      kPHInhelastic, kPMuonNuclear, kPTOFlimit,kPPhotoFission, kPNoProcess, 
431      kPRayleigh, kPNoProcess, kPNoProcess, kPNoProcess, kPNull, kPStop
432     };
433     
434     for (Int_t i = 0;i<kMaxMCProcess;i++)
435     {
436       if (kMCprocesses[i] == no)
437         {
438           return kMCprocesses[i];
439         }
440     } 
441     return kPNoProcess;
442 }
443