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