http://savannah.cern.ch/bugs/?103960
[u/mrichter/AliRoot.git] / EVGEN / AliGenEMCocktail.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: AliGenEMCocktail.cxx 40702 2010-04-26 13:09:52Z morsch $ */
17
18 // Class to create the cocktail for physics with electrons, di-electrons,
19 // and photons from the decay of the following sources:
20 // pizero, eta, rho, omega, etaprime, phi
21 // Kinematic distributions of the sources are taken from AliGenEMlib.
22 // Decay channels can be selected via the method SetDecayMode.
23 // Particles can be generated flat in pT with weights according to the
24 // chosen pT distributions from AliGenEMlib (weighting mode: kNonAnalog),
25 // or they are generated according to the pT distributions themselves
26 // (weighting mode: kAnalog)  
27  
28  
29 #include <TObjArray.h>
30 #include <TParticle.h>
31 #include <TF1.h>
32 #include <TVirtualMC.h>
33 #include <TPDGCode.h>
34 #include "AliGenCocktailEventHeader.h"
35
36 #include "AliGenCocktailEntry.h"
37 #include "AliGenEMCocktail.h"
38 #include "AliGenEMlib.h"
39 #include "AliGenBox.h"
40 #include "AliGenParam.h"
41 #include "AliMC.h"
42 #include "AliRun.h"
43 #include "AliStack.h"
44 #include "AliDecayer.h"
45 #include "AliDecayerPythia.h"
46 #include "AliLog.h"
47 #include "AliGenCorrHF.h"
48
49 ClassImp(AliGenEMCocktail)  
50   
51 //________________________________________________________________________
52 AliGenEMCocktail::AliGenEMCocktail()
53 :AliGenCocktail(),
54    fDecayer(0),
55    fDecayMode(kAll),
56    fWeightingMode(kNonAnalog),
57    fNPart(1000),
58   fYieldArray(),
59   fPtSelect(0),
60   fCentrality(0),
61   fV2Systematic(0),
62   fForceConv(kFALSE),
63   fHeaviestParticle(kGENs)
64 {
65   // Constructor
66
67 }
68
69 //_________________________________________________________________________
70 AliGenEMCocktail::~AliGenEMCocktail()
71 {
72   // Destructor
73
74 }
75
76 //_________________________________________________________________________
77 void AliGenEMCocktail::CreateCocktail()
78 {
79   // create and add sources to the cocktail
80
81   fDecayer->SetForceDecay(fDecayMode);
82   fDecayer->ForceDecay();
83
84   // Set kinematic limits
85   Double_t ptMin  = fPtMin;
86   Double_t ptMax  = fPtMax;
87   Double_t yMin   = fYMin;;
88   Double_t yMax   = fYMax;;
89   Double_t phiMin = fPhiMin*180./TMath::Pi();
90   Double_t phiMax = fPhiMax*180./TMath::Pi();
91   AliInfo(Form("Ranges pT:%4.1f : %4.1f GeV/c, y:%4.2f : %4.2f, Phi:%5.1f : %5.1f degres",ptMin,ptMax,yMin,yMax,phiMin,phiMax));
92   AliInfo(Form("the parametrised sources uses the decay mode %d",fDecayMode));
93
94   //Initialize user selection for Pt Parameterization and centrality:
95   AliGenEMlib::SelectParams(fPtSelect,fCentrality,fV2Systematic);
96
97   // Create and add electron sources to the generator
98
99   // pizero
100   AliGenParam *genpizero=0;
101   Char_t namePizero[10];    
102   snprintf(namePizero,10,"Pizero");    
103   genpizero = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
104   AddSource2Generator(namePizero,genpizero);
105   TF1 *fPtPizero = genpizero->GetPt();
106 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
107   fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
108 #else
109   fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
110 #endif
111
112   // eta  
113   if(fHeaviestParticle<kGenEta)return;
114   AliGenParam *geneta=0;
115   Char_t nameEta[10];    
116   snprintf(nameEta,10,"Eta");    
117   geneta = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
118   AddSource2Generator(nameEta,geneta);
119   TF1 *fPtEta = geneta->GetPt();
120 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
121   fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
122 #else
123   fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
124 #endif
125
126   // rho  
127   if(fHeaviestParticle<kGenRho)return;
128   AliGenParam *genrho=0;
129   Char_t nameRho[10];    
130   snprintf(nameRho,10,"Rho");    
131   genrho = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
132   AddSource2Generator(nameRho,genrho);
133   TF1 *fPtRho = genrho->GetPt();
134 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
135   fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
136 #else
137   fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
138 #endif
139   
140   // omega
141   if(fHeaviestParticle<kGenOmega)return;
142   AliGenParam *genomega=0;
143   Char_t nameOmega[10];    
144   snprintf(nameOmega,10,"Omega");    
145   genomega = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
146   AddSource2Generator(nameOmega,genomega);
147   TF1 *fPtOmega = genomega->GetPt();
148 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
149   fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
150 #else
151   fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
152 #endif
153
154   // etaprime
155   if(fHeaviestParticle<kGenEtaprime)return;
156   AliGenParam *genetaprime=0;
157   Char_t nameEtaprime[10];    
158   snprintf(nameEtaprime,10,"Etaprime");    
159   genetaprime = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
160   AddSource2Generator(nameEtaprime,genetaprime);
161   TF1 *fPtEtaprime = genetaprime->GetPt();
162 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
163   fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
164 #else
165   fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
166 #endif
167
168   // phi  
169   if(fHeaviestParticle<kGenPhi)return;
170   AliGenParam *genphi=0;
171   Char_t namePhi[10];    
172   snprintf(namePhi,10,"Phi");    
173   genphi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
174   AddSource2Generator(namePhi,genphi);
175   TF1 *fPtPhi = genphi->GetPt();
176 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
177   fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
178 #else
179   fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
180 #endif
181
182   // jpsi  
183   if(fHeaviestParticle<kGenJpsi)return;
184   AliGenParam *genjpsi=0;
185   Char_t nameJpsi[10];    
186   snprintf(nameJpsi,10,"Jpsi");    
187   genjpsi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
188   AddSource2Generator(nameJpsi,genjpsi);
189   TF1 *fPtJpsi = genjpsi->GetPt();
190 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
191   fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
192 #else
193   fYieldArray[kGenJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
194 #endif
195
196 }
197
198 //-------------------------------------------------------------------
199 void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource, 
200                                          AliGenParam* const genSource)
201 {
202 // add sources to the cocktail
203   Double_t phiMin = fPhiMin*180./TMath::Pi();
204   Double_t phiMax = fPhiMax*180./TMath::Pi();
205
206   genSource->SetPtRange(fPtMin, fPtMax);  
207   genSource->SetYRange(fYMin, fYMax);
208   genSource->SetPhiRange(phiMin, phiMax);
209   genSource->SetWeighting(fWeightingMode);
210   genSource->SetForceGammaConversion(fForceConv);
211   if (!gMC) genSource->SetDecayer(fDecayer);  
212   genSource->Init();
213     
214   AddGenerator(genSource,nameSource,1.); // Adding Generator    
215 }
216
217 //-------------------------------------------------------------------
218 void AliGenEMCocktail::Init()
219 {
220   // Initialisation
221   TIter next(fEntries);
222   AliGenCocktailEntry *entry;
223   if (fStack) {
224     while((entry = (AliGenCocktailEntry*)next())) {
225       entry->Generator()->SetStack(fStack);
226     }
227   }
228 }
229
230 //_________________________________________________________________________
231 void AliGenEMCocktail::Generate()
232 {
233   // Generate event 
234   TIter next(fEntries);
235   AliGenCocktailEntry *entry = 0;
236   AliGenCocktailEntry *preventry = 0;
237   AliGenerator* gen = 0;
238
239   if (fHeader) delete fHeader;
240   fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
241
242   const TObjArray *partArray = gAlice->GetMCApp()->Particles();
243     
244   // Generate the vertex position used by all generators    
245   if(fVertexSmear == kPerEvent) Vertex();
246
247   //Reseting stack
248   AliRunLoader * runloader = AliRunLoader::Instance();
249   if (runloader)
250     if (runloader->Stack())
251       runloader->Stack()->Clean();
252   
253   // Loop over generators and generate events
254   Int_t igen = 0;
255   Float_t evPlane;
256   Rndm(&evPlane,1);
257   evPlane*=TMath::Pi()*2;
258   while((entry = (AliGenCocktailEntry*)next())) {
259     gen = entry->Generator();
260     gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
261     
262     if (fNPart > 0) {
263       igen++;   
264       if (igen == 1) entry->SetFirst(0);                
265       else  entry->SetFirst((partArray->GetEntriesFast())+1);
266       gen->SetEventPlane(evPlane);
267       gen->SetNumberParticles(fNPart);          
268       gen->Generate();
269       entry->SetLast(partArray->GetEntriesFast());
270       preventry = entry;
271     }
272   }  
273   next.Reset();
274
275   // Setting weights for proper absolute normalization
276   Int_t iPart, iMother;
277   Int_t pdgMother = 0;
278   Double_t weight = 0.;
279   Double_t dNdy = 0.;
280   Int_t maxPart = partArray->GetEntriesFast();
281   for(iPart=0; iPart<maxPart; iPart++){      
282     TParticle *part = gAlice->GetMCApp()->Particle(iPart);
283     iMother = part->GetFirstMother();
284     TParticle *mother = 0;
285     if (iMother>=0){
286       mother = gAlice->GetMCApp()->Particle(iMother);
287       pdgMother = mother->GetPdgCode();
288     }
289     else
290       pdgMother = part->GetPdgCode();
291     switch (pdgMother){
292     case 111:
293       dNdy = fYieldArray[kGenPizero];
294       break;
295     case 221:
296       dNdy = fYieldArray[kGenEta];
297       break;
298     case 113:
299       dNdy = fYieldArray[kGenRho];
300       break;
301     case 223:
302       dNdy = fYieldArray[kGenOmega];
303       break;
304     case 331:
305       dNdy = fYieldArray[kGenEtaprime];
306       break;
307     case 333:
308       dNdy = fYieldArray[kGenPhi];
309       break;
310     case 443:
311       dNdy = fYieldArray[kGenJpsi];
312       break;
313       
314     default:
315       dNdy = 0.;
316     }
317     weight = dNdy*part->GetWeight();
318     part->SetWeight(weight);
319   }     
320   
321   fHeader->SetNProduced(maxPart);
322
323
324   TArrayF eventVertex;
325   eventVertex.Set(3);
326   for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
327   
328   fHeader->SetPrimaryVertex(eventVertex);
329
330   gAlice->SetGenEventHeader(fHeader);
331 }