Changes to compile with Root6 on macosx64
[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 <TDatabasePDG.h>
35 #include "AliGenCocktailEventHeader.h"
36
37 #include "AliGenCocktailEntry.h"
38 #include "AliGenEMCocktail.h"
39 #include "AliGenEMlib.h"
40 #include "AliGenBox.h"
41 #include "AliGenParam.h"
42 #include "AliMC.h"
43 #include "AliRun.h"
44 #include "AliStack.h"
45 #include "AliDecayer.h"
46 #include "AliDecayerPythia.h"
47 #include "AliLog.h"
48 #include "AliGenCorrHF.h"
49
50 ClassImp(AliGenEMCocktail)  
51   
52 //________________________________________________________________________
53 AliGenEMCocktail::AliGenEMCocktail()
54 :AliGenCocktail(),
55    fDecayer(0),
56    fDecayMode(kAll),
57    fWeightingMode(kNonAnalog),
58    fNPart(1000),
59   fYieldArray(),
60   fPtSelect(0),
61   fCentrality(0),
62   fV2Systematic(0),
63   fForceConv(kFALSE),
64   fSelectedParticles(kGenHadrons)
65 {
66   // Constructor
67
68 }
69
70 //_________________________________________________________________________
71 AliGenEMCocktail::~AliGenEMCocktail()
72 {
73   // Destructor
74
75 }
76
77 //_________________________________________________________________________
78 void AliGenEMCocktail::CreateCocktail()
79 {
80   // create and add sources to the cocktail
81
82   fDecayer->SetForceDecay(fDecayMode);
83   fDecayer->ForceDecay();
84
85   // Set kinematic limits
86   Double_t ptMin  = fPtMin;
87   Double_t ptMax  = fPtMax;
88   Double_t yMin   = fYMin;;
89   Double_t yMax   = fYMax;;
90   Double_t phiMin = fPhiMin*180./TMath::Pi();
91   Double_t phiMax = fPhiMax*180./TMath::Pi();
92   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));
93   AliInfo(Form("the parametrised sources uses the decay mode %d",fDecayMode));
94
95   //Initialize user selection for Pt Parameterization and centrality:
96   AliGenEMlib::SelectParams(fPtSelect,fCentrality,fV2Systematic);
97
98   // Create and add electron sources to the generator
99   // pizero
100   if(fSelectedParticles&kGenPizero){
101   AliGenParam *genpizero=0;
102   Char_t namePizero[10];    
103   snprintf(namePizero,10,"Pizero");    
104   genpizero = new AliGenParam(fNPart/0.925, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
105   genpizero->SetYRange(fYMin/0.925, fYMax/0.925);
106   AddSource2Generator(namePizero,genpizero);
107   TF1 *fPtPizero = genpizero->GetPt();
108 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
109     fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,1.e-6);
110 #else
111     fYieldArray[kPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
112 #endif
113   }
114
115   // eta  
116   if(fSelectedParticles&kGenEta){
117   AliGenParam *geneta=0;
118   Char_t nameEta[10];    
119   snprintf(nameEta,10,"Eta");    
120   geneta = new AliGenParam(fNPart/0.825, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
121   geneta->SetYRange(fYMin/0.825, fYMax/0.825);
122   AddSource2Generator(nameEta,geneta);
123   TF1 *fPtEta = geneta->GetPt();
124 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
125     fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,1.e-6);
126 #else
127     fYieldArray[kEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
128 #endif
129   }
130
131   // rho  
132   if(fSelectedParticles&kGenRho){
133   AliGenParam *genrho=0;
134   Char_t nameRho[10];    
135   snprintf(nameRho,10,"Rho");    
136   genrho = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
137   genrho->SetYRange(fYMin/0.775, fYMax/0.775);
138   AddSource2Generator(nameRho,genrho);
139   TF1 *fPtRho = genrho->GetPt();
140 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
141     fYieldArray[kRho] = fPtRho->Integral(fPtMin,fPtMax,1.e-6);
142 #else
143     fYieldArray[kRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
144 #endif
145   }
146   
147   // omega
148   if(fSelectedParticles&kGenOmega){
149   AliGenParam *genomega=0;
150   Char_t nameOmega[10];    
151   snprintf(nameOmega,10,"Omega");    
152   genomega = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
153   genomega->SetYRange(fYMin/0.775, fYMax/0.775);
154   AddSource2Generator(nameOmega,genomega);
155   TF1 *fPtOmega = genomega->GetPt();
156 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
157     fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,1.e-6);
158 #else
159     fYieldArray[kOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
160 #endif
161   }
162
163   // etaprime
164   if(fSelectedParticles&kGenEtaprime){
165   AliGenParam *genetaprime=0;
166   Char_t nameEtaprime[10];    
167   snprintf(nameEtaprime,10,"Etaprime");    
168   genetaprime = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
169   genetaprime->SetYRange(fYMin/0.725, fYMax/0.725);
170   AddSource2Generator(nameEtaprime,genetaprime);
171   TF1 *fPtEtaprime = genetaprime->GetPt();
172 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
173     fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,1.e-6);
174 #else
175     fYieldArray[kEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
176 #endif
177   }
178
179   // phi  
180   if(fSelectedParticles&kGenPhi){
181   AliGenParam *genphi=0;
182   Char_t namePhi[10];    
183   snprintf(namePhi,10,"Phi");    
184   genphi = new AliGenParam(fNPart/0.725, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
185   genphi->SetYRange(fYMin/0.725, fYMax/0.725);
186   AddSource2Generator(namePhi,genphi);
187   TF1 *fPtPhi = genphi->GetPt();
188 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
189     fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,1.e-6);
190 #else
191     fYieldArray[kPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
192 #endif
193   }
194
195   // jpsi  
196   if(fSelectedParticles&kGenJpsi){
197   AliGenParam *genjpsi=0;
198   Char_t nameJpsi[10];    
199   snprintf(nameJpsi,10,"Jpsi");    
200   genjpsi = new AliGenParam(fNPart/0.525, new AliGenEMlib(), AliGenEMlib::kJpsi, "DUMMY");
201   genjpsi->SetYRange(fYMin/0.525, fYMax/0.525);
202   AddSource2Generator(nameJpsi,genjpsi);
203   TF1 *fPtJpsi = genjpsi->GetPt();
204 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
205     fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,1.e-6);
206 #else
207     fYieldArray[kJpsi] = fPtJpsi->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
208 #endif
209   }
210
211   // direct gamma
212   if(fDecayMode!=kGammaEM) return;
213   if(fSelectedParticles&kGenDirectRealGamma){
214     AliGenParam *genDirectRealG=0;
215     Char_t nameDirectRealG[10];    
216     snprintf(nameDirectRealG,10,"DirectRealGamma");    
217     genDirectRealG = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kDirectRealGamma, "DUMMY");
218     genDirectRealG->SetYRange(fYMin, fYMax);
219     AddSource2Generator(nameDirectRealG,genDirectRealG);
220     TF1 *fPtDirectRealG = genDirectRealG->GetPt();
221 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
222     fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,1.e-6);
223 #else
224     fYieldArray[kDirectRealGamma] = fPtDirectRealG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
225 #endif
226   }
227
228   if(fSelectedParticles&kGenDirectVirtGamma){
229     TDatabasePDG::Instance()->AddParticle("DirectVirtGamma","DirectVirtGamma",0,true,0,0,"GaugeBoson",220000);
230     AliGenParam *genDirectVirtG=0;
231     Char_t nameDirectVirtG[10];    
232     snprintf(nameDirectVirtG,10,"DirectVirtGamma");    
233     genDirectVirtG = new AliGenParam(fNPart/0.775, new AliGenEMlib(), AliGenEMlib::kDirectVirtGamma, "DUMMY");
234     genDirectVirtG->SetYRange(fYMin/0.775, fYMax/0.775);
235     AddSource2Generator(nameDirectVirtG,genDirectVirtG);
236     TF1 *fPtDirectVirtG = genDirectVirtG->GetPt();
237 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,99,0)
238     fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,1.e-6);
239 #else
240     fYieldArray[kDirectVirtGamma] = fPtDirectVirtG->Integral(fPtMin,fPtMax,(Double_t *)0,1.e-6);
241 #endif
242   }
243 }
244
245 //-------------------------------------------------------------------
246 void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource, 
247                                          AliGenParam* const genSource)
248 {
249 // add sources to the cocktail
250   Double_t phiMin = fPhiMin*180./TMath::Pi();
251   Double_t phiMax = fPhiMax*180./TMath::Pi();
252
253   genSource->SetPtRange(fPtMin, fPtMax);  
254   genSource->SetPhiRange(phiMin, phiMax);
255   genSource->SetWeighting(fWeightingMode);
256   genSource->SetForceGammaConversion(fForceConv);
257   if (!TVirtualMC::GetMC()) genSource->SetDecayer(fDecayer);  
258   genSource->Init();
259     
260   AddGenerator(genSource,nameSource,1.); // Adding Generator    
261 }
262
263 //-------------------------------------------------------------------
264 void AliGenEMCocktail::Init()
265 {
266   // Initialisation
267   TIter next(fEntries);
268   AliGenCocktailEntry *entry;
269   if (fStack) {
270     while((entry = (AliGenCocktailEntry*)next())) {
271       entry->Generator()->SetStack(fStack);
272     }
273   }
274 }
275
276 //_________________________________________________________________________
277 void AliGenEMCocktail::Generate()
278 {
279   // Generate event 
280   TIter next(fEntries);
281   AliGenCocktailEntry *entry = 0;
282   AliGenCocktailEntry *preventry = 0;
283   AliGenerator* gen = 0;
284
285   if (fHeader) delete fHeader;
286   fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
287
288   const TObjArray *partArray = gAlice->GetMCApp()->Particles();
289     
290   // Generate the vertex position used by all generators    
291   if(fVertexSmear == kPerEvent) Vertex();
292
293   //Reseting stack
294   AliRunLoader * runloader = AliRunLoader::Instance();
295   if (runloader)
296     if (runloader->Stack())
297       runloader->Stack()->Clean();
298   
299   // Loop over generators and generate events
300   Int_t igen = 0;
301   Float_t evPlane;
302   Rndm(&evPlane,1);
303   evPlane*=TMath::Pi()*2;
304   while((entry = (AliGenCocktailEntry*)next())) {
305     gen = entry->Generator();
306     gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
307     
308     if (fNPart > 0) {
309       igen++;   
310       if (igen == 1) entry->SetFirst(0);                
311       else  entry->SetFirst((partArray->GetEntriesFast())+1);
312       gen->SetEventPlane(evPlane);
313       gen->Generate();
314       entry->SetLast(partArray->GetEntriesFast());
315       preventry = entry;
316     }
317   }  
318   next.Reset();
319
320   // Setting weights for proper absolute normalization
321   Int_t iPart, iMother;
322   Int_t pdgMother = 0;
323   Double_t weight = 0.;
324   Double_t dNdy = 0.;
325   Int_t maxPart = partArray->GetEntriesFast();
326   for(iPart=0; iPart<maxPart; iPart++){      
327     TParticle *part = gAlice->GetMCApp()->Particle(iPart);
328     iMother = part->GetFirstMother();
329     TParticle *mother = 0;
330     if (iMother>=0){
331       mother = gAlice->GetMCApp()->Particle(iMother);
332       pdgMother = mother->GetPdgCode();
333     }
334     else
335       pdgMother = part->GetPdgCode();
336
337     switch (pdgMother){
338     case 111:
339       dNdy = fYieldArray[kPizero];
340       break;
341     case 221:
342       dNdy = fYieldArray[kEta];
343       break;
344     case 113:
345       dNdy = fYieldArray[kRho];
346       break;
347     case 223:
348       dNdy = fYieldArray[kOmega];
349       break;
350     case 331:
351       dNdy = fYieldArray[kEtaprime];
352       break;
353     case 333:
354       dNdy = fYieldArray[kPhi];
355       break;
356     case 443:
357       dNdy = fYieldArray[kJpsi];
358       break;
359     case 22:
360       dNdy = fYieldArray[kDirectRealGamma];
361       break;
362     case 220000:
363       dNdy = fYieldArray[kDirectVirtGamma];
364       break;
365       
366     default:
367       dNdy = 0.;
368     }
369
370     weight = dNdy*part->GetWeight();
371     part->SetWeight(weight);
372   }     
373   
374   fHeader->SetNProduced(maxPart);
375
376
377   TArrayF eventVertex;
378   eventVertex.Set(3);
379   for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
380   
381   fHeader->SetPrimaryVertex(eventVertex);
382
383   gAlice->SetGenEventHeader(fHeader);
384 }