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