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