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