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