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