2ad24115f1d8412b78a2e3baf79116e9ee2c3390
[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 // 
17 // Container class for AliGenerator and AfterBurners 
18 // (which are AliGenerators as well) through recursion.
19 // The container is itself an AliGenerator a
20 // what is stored are not the pointers to the generators directly 
21 // but to objects of type
22 // AliGenCocktailAfterBurner entry.   
23 // The class provides also iterator functionality.  
24 // Author: andreas.morsch@cern.ch and piotr.skowronski@cern.ch
25 //
26
27 #include "AliGenCocktailAfterBurner.h"
28 #include "AliGenCocktailEntry.h"
29 #include "AliRun.h"
30 #include "AliStack.h"
31 #include <TObjArray.h>
32 #include <TList.h>
33 #include <TParticle.h>
34 #include <iostream.h>
35
36 static const Bool_t debug = kFALSE;
37 ClassImp(AliGenCocktailAfterBurner)
38
39 AliGenCocktailAfterBurner::AliGenCocktailAfterBurner()
40 {
41 // Constructor
42     if (debug) 
43         cout<<"AliGenCocktailAfterBurner::AliGenCocktailAfterBurner()"<<endl;
44     SetName("AliGenCocktailAfterBurner");
45     SetTitle("AliGenCocktailAfterBurner");
46     fInternalStacks =0;
47     fCurrentEvent =0;
48     fAfterBurnerEntries = new TList();
49     fNAfterBurners = 0;
50     fGenerationDone = kFALSE;
51     
52     fActiveEvent = -1;  
53     fNParticles = -1;
54 }
55
56
57 AliGenCocktailAfterBurner::~AliGenCocktailAfterBurner()
58   {
59     fInternalStacks->SetOwner();
60     delete fInternalStacks;
61     delete fAfterBurnerEntries;
62   }
63
64 void AliGenCocktailAfterBurner::
65 AddAfterBurner(AliGenerator *AfterBurner, char* Name, Float_t RateExp)
66 {
67 //
68 //  Forward parameters to the new AfterBurner
69     
70     if (debug)cout<<"AliGenCocktailAfterBurner::AddAfterBurner  Named "<<Name<<endl;
71
72     if(TestBit(kPtRange)) 
73         AfterBurner->SetPtRange(fPtMin,fPtMax);
74     if(TestBit(kMomentumRange))
75         AfterBurner->SetMomentumRange(fPMin,fPMax);
76     
77     AfterBurner->SetYRange(fYMin,fYMax);
78     AfterBurner->
79         SetPhiRange(fPhiMin*180/TMath::Pi(),fPhiMax*180/TMath::Pi());
80     AfterBurner->
81         SetThetaRange(fThetaMin*180/TMath::Pi(),fThetaMax*180/TMath::Pi());
82     AfterBurner->
83         SetOrigin(fOrigin[0], fOrigin[1], fOrigin[2]);
84     AfterBurner->
85         SetSigma(fOsigma[0], fOsigma[1], fOsigma[2]);
86     AfterBurner->SetVertexSmear(fVertexSmear);
87     AfterBurner->SetTrackingFlag(fTrackIt);    
88 //
89 //  Add AfterBurner to list   
90     
91     AliGenCocktailEntry *entry = 
92         new AliGenCocktailEntry(AfterBurner, Name, RateExp);
93     fAfterBurnerEntries->Add(entry);
94     fNAfterBurners++;
95 }
96
97 void AliGenCocktailAfterBurner::Init()
98 {
99 // Initialisation
100     this->AliGenCocktail::Init(); 
101     
102     if (debug) cout<<"AliGenCocktailAfterBurner::Init"<<endl;
103     TIter next(fAfterBurnerEntries);
104     AliGenCocktailEntry *entry;
105     //
106     // Loop over generators and initialize
107     while((entry = (AliGenCocktailEntry*)next())) {
108         entry->Generator()->Init();
109     }  
110 }
111
112 void AliGenCocktailAfterBurner::Generate()
113 {
114 //
115 // Generate event 
116     if (debug)cout<<"#####################################"<<endl
117                   <<"#AliGenCocktailAfterBurner::Generate#"
118                   <<endl<<"#####################################"<<endl;
119     
120     Int_t i; //iterator
121     AliStack * stack;
122     if (fGenerationDone)
123     { 
124         SetTracks(++fCurrentEvent);
125         cout<<"Returning event "<<fCurrentEvent<<endl;
126         return;  
127     }
128     else
129     { 
130         fCurrentEvent=0;
131         fInternalStacks = new TObjArray(fNumberOfEvents); //Create array of internal stacks
132         for(i=0;i<fNumberOfEvents;i++) 
133         {       
134             stack = new AliStack(10000);
135             stack->Reset();
136             fInternalStacks->Add(stack);
137         }
138 /*********************************************************************/ 
139         TIter next(fEntries);
140         AliGenCocktailEntry *entry;
141         AliGenCocktailEntry *e1;
142         AliGenCocktailEntry *e2;
143         TObjArray *partArray;
144         //
145         // Loop over generators and generate events
146         Int_t igen=0;
147         while((entry = (AliGenCocktailEntry*)next())) 
148         {
149             igen++;
150             cout<<"Generator number "<<igen<<endl;
151             /***********************************************/
152 //First generator for all evenets, than second for all events, etc...
153             for(i=0;i<fNumberOfEvents;i++) 
154             {  
155                 
156                 cout<<"                  EVENT "<<i<<endl;
157                 stack = GetStack(i);
158                 partArray = stack->Particles();
159                 fCurrentGenerator = entry->Generator();
160                 fCurrentGenerator->SetStack(stack);
161                 
162                 if (igen ==1) 
163                 {
164                     entry->SetFirst(0);
165                 } 
166                 else 
167                 {
168                     entry->SetFirst((partArray->GetEntriesFast())+1);
169                 }
170                 fCurrentGenerator->Generate();
171                 entry->SetLast(partArray->GetEntriesFast());
172             }
173             /***********************************************/
174         }
175         next.Reset();
176         while((entry = (AliGenCocktailEntry*)next())) 
177         {
178             entry->PrintInfo();
179         }
180         for ( entry=FirstGenerator();entry;entry=NextGenerator() ) 
181         {
182             entry->PrintInfo();
183         }
184         
185         for (FirstGeneratorPair(e1,e2); (e1&&e2); NextGeneratorPair(e1,e2) )
186         {
187             printf("\n -----------------------------");
188             e1->PrintInfo();
189             e2->PrintInfo();
190         }
191         
192         /***********************************************/
193         /*******After Burners Processing****************/
194         /***********************************************/
195         TIter nextAfterBurner(fAfterBurnerEntries);
196         AliGenCocktailEntry *afterBurnerEntry;
197         Int_t iab =0;
198         cout<<"\n\nRunning After Burners"<<endl;
199         while((afterBurnerEntry = (AliGenCocktailEntry*)nextAfterBurner()))
200         {
201             cout<<"After Burner number "<<iab++<<endl;
202             fCurrentGenerator = afterBurnerEntry->Generator();
203             fCurrentGenerator->Generate();
204         }
205         cout<<endl<<"Finished. Processed "<<iab<<" After Burners"<<endl;
206         
207         /***********************************************/
208         /***********************************************/
209         /***********************************************/       
210         
211         fGenerationDone=kTRUE; 
212         SetTracks(0); //copy event 0 to gAlice stack
213         
214 /*********************************************************************/
215         
216     }//else generated
217 }
218
219
220 AliGenCocktailAfterBurner& AliGenCocktailAfterBurner::operator=(const  AliGenCocktailAfterBurner& rhs)
221 {
222 // Assignment operator
223     return *this;
224 }
225
226
227 AliStack* AliGenCocktailAfterBurner::GetStack(Int_t n)
228 {
229     if( (n<0) || (n>=fNumberOfEvents) )
230     {
231         Fatal("AliGenCocktailAfterBurner::GetStack","Asked for non existing stack (%d)",n);
232         return 0; 
233     }
234     return ((AliStack*) ((*fInternalStacks)[n]) );
235 }
236
237 void AliGenCocktailAfterBurner::SetActiveEventNumber(Int_t actev)
238 {
239     fActiveEvent = actev;
240     fActiveStack = GetStack(actev);
241 }
242
243 void AliGenCocktailAfterBurner::SetTracks(Int_t stackno)
244 {
245     AliStack* instack = GetStack(stackno);
246     Int_t done;
247     Int_t parent; 
248     Int_t pdg;
249     Double_t px, py, pz, e, vx, vy, vz, tof, polx, poly, polz;
250     AliMCProcess mech;
251     Int_t ntr;
252     Float_t weight;
253     TVector3 pol;
254     
255     TParticle * p;
256     Int_t N = instack->GetNtrack();
257     if (debug) 
258     {
259         cout<<"AliGenCocktailAfterBurner::SetTracks("<<stackno<<"). Number of particles is: "<<N<<"\n";
260     }
261     
262     for(Int_t i = 0; i < N; i++)
263     {
264         
265         p = instack->Particle(i);
266         done = !p->TestBit(kDoneBit);
267         if (debug) {cout<<i<<"  "<<done<<"  "; fflush(0);}
268         parent = p->GetMother(0);
269         pdg = p->GetPdgCode();
270         px = p->Px();
271         py = p->Py();
272         pz = p->Pz();
273         e  = p->Energy();
274         vx = p->Vx();
275         vy = p->Vy();
276         vz = p->Vz();
277         tof = p->T();
278         p->GetPolarisation(pol);
279         polx = pol.X();
280         poly = pol.Y();
281         polz = pol.Z();
282         mech = AliGenCocktailAfterBurner::IntToMCProcess(p->GetUniqueID());
283         weight = p->GetWeight();
284         
285         gAlice->SetTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof,
286                          polx, poly, polz, mech, ntr, weight);
287     }
288 }
289
290 AliMCProcess AliGenCocktailAfterBurner::IntToMCProcess(Int_t no)
291 {
292     const AliMCProcess MCprocesses[kMaxMCProcess] = 
293     {kPNoProcess, kPMultipleScattering, kPEnergyLoss, kPMagneticFieldL, 
294      kPDecay, kPPair, kPCompton, kPPhotoelectric, kPBrem, kPDeltaRay,
295      kPAnnihilation, kPHadronic, kPNoProcess, kPEvaporation, kPNuclearFission,
296      kPNuclearAbsorption, kPPbarAnnihilation, kPNCapture, kPHElastic, 
297      kPHInhelastic, kPMuonNuclear, kPTOFlimit,kPPhotoFission, kPNoProcess, 
298      kPRayleigh, kPNoProcess, kPNoProcess, kPNoProcess, kPNull, kPStop};
299     
300     for (Int_t i = 0;i<kMaxMCProcess;i++)
301     {
302         if (MCprocesses[i] == no)
303         {
304             //if (debug) cout<<"IntToMCProcess("<<no<<") returned AliMCProcess Named \""<<AliMCProcessName[MCprocesses[i]]<<"\""<<endl;
305             return MCprocesses[i];
306         }
307     } 
308     return kPNoProcess;
309 }