Names changed in order to avoid clash with FLUKA
[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 gAlice->GetEventsPerRun()
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 "AliStack.h"
43
44
45 ClassImp(AliGenCocktailAfterBurner)
46 /*********************************************************************/ 
47 /*********************************************************************/ 
48
49 AliGenCocktailAfterBurner::AliGenCocktailAfterBurner()
50 {
51 // Constructor
52     if (gDebug > 0) 
53        cout<<"AliGenCocktailAfterBurner::AliGenCocktailAfterBurner()"<<endl;
54     SetName("AliGenCocktailAfterBurner");
55     SetTitle("AliGenCocktailAfterBurner");
56     fInternalStacks =0;
57     fActiveStack = 0;
58     fCurrentGenerator = 0;
59     fCurrentEvent =0;
60     fAfterBurnerEntries = new TList();
61     fNAfterBurners = 0;
62     fGenerationDone = kFALSE;
63     
64     fActiveEvent = -1;  
65     fNBgEvents = 0;
66 }
67 /*********************************************************************/ 
68 AliGenCocktailAfterBurner::AliGenCocktailAfterBurner(const AliGenCocktailAfterBurner& in)
69 {
70  //cpy ctor
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   }
86 /*********************************************************************/ 
87 /*********************************************************************/ 
88
89 void AliGenCocktailAfterBurner::
90 AddAfterBurner(AliGenerator *AfterBurner, char* Name, Float_t RateExp)
91 {
92 //
93 //  Forward parameters to the new AfterBurner
94     
95     if (gDebug>0)cout<<"AliGenCocktailAfterBurner::AddAfterBurner  Named "<<Name<<endl;
96
97     if(TestBit(kPtRange)) 
98         AfterBurner->SetPtRange(fPtMin,fPtMax);
99     if(TestBit(kMomentumRange))
100         AfterBurner->SetMomentumRange(fPMin,fPMax);
101     
102     AfterBurner->SetYRange(fYMin,fYMax);
103     AfterBurner->SetPhiRange(fPhiMin*180/TMath::Pi(),fPhiMax*180/TMath::Pi());
104     AfterBurner->SetThetaRange(fThetaMin*180/TMath::Pi(),fThetaMax*180/TMath::Pi());
105     AfterBurner->SetOrigin(fOrigin[0], fOrigin[1], fOrigin[2]);
106     AfterBurner->SetSigma(fOsigma[0], fOsigma[1], fOsigma[2]);
107     AfterBurner->SetVertexSmear(fVertexSmear);
108     AfterBurner->SetTrackingFlag(fTrackIt);    
109 //
110 //  Add AfterBurner to list   
111     
112     AliGenCocktailEntry *entry = 
113         new AliGenCocktailEntry(AfterBurner, Name, RateExp);
114     fAfterBurnerEntries->Add(entry);
115     fNAfterBurners++;
116 //
117     
118 }
119 /*********************************************************************/ 
120 /*********************************************************************/ 
121
122 void AliGenCocktailAfterBurner::Init()
123 {
124 // Initialisation
125     fGenerationDone = kFALSE;
126     if (fInternalStacks) //delete stacks
127      { 
128        fInternalStacks->SetOwner();
129        fInternalStacks->Delete(); //clean after previous generation cycle
130      }
131
132     this->AliGenCocktail::Init(); 
133     
134     if (gDebug>0) cout<<"AliGenCocktailAfterBurner::Init"<<endl;
135     TIter next(fAfterBurnerEntries);
136     AliGenCocktailEntry *entry;
137     //
138     // Loop over generators and initialize
139     while((entry = (AliGenCocktailEntry*)next())) {
140         entry->Generator()->Init();
141     }  
142 }
143 /*********************************************************************/ 
144 /*********************************************************************/ 
145
146 void AliGenCocktailAfterBurner::Generate()
147 {
148 //
149 // Generate event
150 //  Firsts runs each generator for all events
151 //  than after burners ones for each event
152 //
153 //  It generates and processes all events during
154 //  first call only.
155 //  In next calls it just returns already generated 
156 //  and processed events to the gAlice
157
158     if (gDebug>0)
159       cout<<"#####################################"<<endl
160           <<"#AliGenCocktailAfterBurner::Generate#"<<endl
161           <<"#####################################"<<endl;
162     
163     Int_t i; //iterator
164     AliStack * stack;
165     
166     if (fGenerationDone)
167     {//if generation is done (in first call) 
168      //just copy particles from the stack to the gAlice
169       SetTracks(++fCurrentEvent);
170       cout<<"Returning event "<<fCurrentEvent<<endl;
171       return;  
172     }
173     else
174     { //Here we are in the first call of the method
175       fCurrentEvent=0;
176       Int_t numberOfEvents = gAlice->GetEventsPerRun();
177       //Create stacks
178       fInternalStacks = new TObjArray(numberOfEvents + fNBgEvents); //Create array of internal stacks
179       for(i=0;i<numberOfEvents + fNBgEvents;i++) 
180        {        
181         stack = new AliStack(10000);
182         stack->Reset();
183         fInternalStacks->Add(stack);
184        }
185 /*********************************************************************/ 
186       TIter next(fEntries);
187       AliGenCocktailEntry *entry;
188       AliGenCocktailEntry *e1;
189       AliGenCocktailEntry *e2;
190       TObjArray *partArray;
191   //
192   // Loop over generators and generate events
193       Int_t igen=0;
194       while((entry = (AliGenCocktailEntry*)next())) 
195       {
196         igen++;
197         cout<<"Generator "<<igen<<"  : "<<entry->GetName()<<endl;
198 /***********************************************/
199 //First generator for all evenets, than second for all events, etc...
200         for(i=0;i<numberOfEvents + fNBgEvents;i++) 
201           {  
202             cout<<"                  EVENT "<<i<<endl;
203             stack = GetStack(i);
204             partArray = stack->Particles();
205             fCurrentGenerator = entry->Generator();
206             fCurrentGenerator->SetStack(stack);
207             if (igen ==1) 
208               {
209                 entry->SetFirst(0);
210               } 
211             else 
212               {
213                 entry->SetFirst((partArray->GetEntriesFast())+1);
214               }
215                 fCurrentGenerator->Generate();
216                 entry->SetLast(partArray->GetEntriesFast());
217            }
218 /***********************************************/
219       }
220       next.Reset();
221       while((entry = (AliGenCocktailEntry*)next())) 
222         {
223           entry->PrintInfo();
224         }
225       for ( entry=FirstGenerator();entry;entry=NextGenerator() ) 
226         {
227           entry->PrintInfo();
228         }
229       for (FirstGeneratorPair(e1,e2); (e1&&e2); NextGeneratorPair(e1,e2) )
230         {
231           printf("\n -----------------------------");
232           e1->PrintInfo();
233           e2->PrintInfo();
234         }
235         
236         
237       /***********************************************/
238       /*******After Burners Processing****************/
239       /***********************************************/
240       TIter nextAfterBurner(fAfterBurnerEntries);
241       AliGenCocktailEntry *afterBurnerEntry;
242       Int_t iab =0; //number of current after burner / counter
243       
244       cout<<"\n\nRunning After Burners"<<endl;
245       while((afterBurnerEntry = (AliGenCocktailEntry*)nextAfterBurner()))
246         {
247           cout<<"After Burner "<<iab++<<"  :"<<afterBurnerEntry->GetName()<<endl;
248           fCurrentGenerator = afterBurnerEntry->Generator();
249           fCurrentGenerator->Generate();
250         }
251       cout<<endl<<"Finished. Processed "<<iab<<" After Burners"<<endl;
252
253       /***********************************************/
254       /***********************************************/
255       /***********************************************/       
256         
257       fGenerationDone=kTRUE; 
258       SetTracks(0); //copy event 0 to gAlice stack
259         
260 /*********************************************************************/
261         
262     }//else generated
263 }
264 /*********************************************************************/
265 /*********************************************************************/ 
266
267 AliGenCocktailAfterBurner& AliGenCocktailAfterBurner::operator=(const  AliGenCocktailAfterBurner& rhs)
268 {
269 // Assignment operator
270     return *this;
271 }
272 /*********************************************************************/
273 /*********************************************************************/ 
274
275 AliStack* AliGenCocktailAfterBurner::GetStack(Int_t n) const
276 {
277 //Returns the pointer to the N'th stack (event)
278   if( ( n<0 ) || ( n>=GetNumberOfEvents() ) )
279     {
280       Fatal("AliGenCocktailAfterBurner::GetStack","Asked for non existing stack (%d)",n);
281       return 0; 
282     }
283     return ((AliStack*) fInternalStacks->At(n) );
284 }
285 /*********************************************************************/ 
286 /*********************************************************************/ 
287
288 void AliGenCocktailAfterBurner::SetActiveEventNumber(Int_t actev)
289 {
290 //Set Active Events Number and Active Stack
291 //There is only one active event number
292 //Made fo convinience of work with AfterBurners (HBT processor)
293
294     fActiveEvent = actev;
295     fActiveStack = GetStack(actev);
296 }
297 /*********************************************************************/ 
298 /*********************************************************************/ 
299
300 void AliGenCocktailAfterBurner::SetTracks(Int_t stackno)
301 {
302 //Method which copies tracks from given stack to the
303 //gAlice's stack
304     AliStack* instack = GetStack(stackno);
305     Int_t done;
306     Int_t parent; 
307     Int_t pdg;
308     Double_t px, py, pz, e, vx, vy, vz, tof, polx, poly, polz;
309     TMCProcess mech;
310     Int_t ntr;
311     Float_t weight;
312     TVector3 pol;
313     
314     TParticle * p;
315     Int_t n = instack->GetNtrack();
316     if (gDebug) 
317     {
318       cout<<"AliGenCocktailAfterBurner::SetTracks("<<stackno<<"). Number of particles is: "<<n<<"\n";
319     }
320     
321     for(Int_t i = 0; i < n; i++)
322     {
323         
324       p = instack->Particle(i);
325       done = !p->TestBit(kDoneBit);
326       parent = p->GetMother(0);
327       pdg = p->GetPdgCode();
328       px = p->Px();
329       py = p->Py();
330       pz = p->Pz();
331       e  = p->Energy();
332       vx = p->Vx();
333       vy = p->Vy();
334       vz = p->Vz();
335       tof = p->T();
336       p->GetPolarisation(pol);
337       polx = pol.X();
338       poly = pol.Y();
339       polz = pol.Z();
340       mech = AliGenCocktailAfterBurner::IntToMCProcess(p->GetUniqueID());
341       weight = p->GetWeight();
342
343       gAlice->SetTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
344                        polx, poly, polz, mech, ntr, weight);
345     }
346 }
347 /*********************************************************************/ 
348 /*********************************************************************/ 
349
350 TMCProcess AliGenCocktailAfterBurner::IntToMCProcess(Int_t no)
351 {
352  //Mothod used to convert uniqueID (integer) to TMCProcess type
353     const TMCProcess kMCprocesses[kMaxMCProcess] = 
354     {
355      kPNoProcess, kPMultipleScattering, kPEnergyLoss, kPMagneticFieldL, 
356      kPDecay, kPPair, kPCompton, kPPhotoelectric, kPBrem, kPDeltaRay,
357      kPAnnihilation, kPHadronic, kPNoProcess, kPEvaporation, kPNuclearFission,
358      kPNuclearAbsorption, kPPbarAnnihilation, kPNCapture, kPHElastic, 
359      kPHInhelastic, kPMuonNuclear, kPTOFlimit,kPPhotoFission, kPNoProcess, 
360      kPRayleigh, kPNoProcess, kPNoProcess, kPNoProcess, kPNull, kPStop
361     };
362     
363     for (Int_t i = 0;i<kMaxMCProcess;i++)
364     {
365       if (kMCprocesses[i] == no)
366         {
367           return kMCprocesses[i];
368         }
369     } 
370     return kPNoProcess;
371 }