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