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