newly implemented model for slow nucleon production
[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 cocktails 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 #include <TObjArray.h>
29 #include <TParticle.h>
30 #include <TF1.h>
31 #include <TVirtualMC.h>
32 #include <TPDGCode.h>
33 #include "AliGenCocktailEventHeader.h"
34
35 #include "AliGenCocktailEntry.h"
36 #include "AliGenEMCocktail.h"
37 #include "AliGenEMlib.h"
38 #include "AliGenBox.h"
39 #include "AliGenParam.h"
40 #include "AliMC.h"
41 #include "AliRun.h"
42 #include "AliStack.h"
43 #include "AliDecayer.h"
44 #include "AliDecayerPythia.h"
45 #include "AliLog.h"
46 #include "AliGenCorrHF.h"
47
48 ClassImp(AliGenEMCocktail)  
49   
50 //________________________________________________________________________
51 AliGenEMCocktail::AliGenEMCocktail()
52   :AliGenCocktail(),
53    fDecayer(0),
54    fDecayMode(kAll),
55    fWeightingMode(kNonAnalog),
56    fNPart(1000),
57    fYieldArray()
58 {
59 // Constructor
60
61 }
62
63 //_________________________________________________________________________
64 AliGenEMCocktail::~AliGenEMCocktail()
65 {
66 // Destructor
67
68 }
69
70 //_________________________________________________________________________
71 void AliGenEMCocktail::CreateCocktail()
72 {
73 // create and add sources to the cocktail
74
75   fDecayer->SetForceDecay(fDecayMode);
76   fDecayer->ForceDecay();
77
78 // Set kinematic limits
79   Double_t ptMin  = fPtMin;
80   Double_t ptMax  = fPtMax;
81   Double_t yMin   = fYMin;;
82   Double_t yMax   = fYMax;;
83   Double_t phiMin = fPhiMin*180./TMath::Pi();
84   Double_t phiMax = fPhiMax*180./TMath::Pi();
85   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));
86   AliInfo(Form("the parametrised sources uses the decay mode %d",fDecayMode));
87
88 // Create and add electron sources to the generator
89
90 // pizero
91   AliGenParam * genpizero=0;
92   Char_t namePizero[10];    
93   snprintf(namePizero,10, "Pizero");    
94   genpizero = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPizero, "DUMMY");
95   AddSource2Generator(namePizero,genpizero);
96   TF1 *fPtPizero = genpizero->GetPt();
97   fYieldArray[kGenPizero] = fPtPizero->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
98 // eta  
99   AliGenParam * geneta=0;
100   Char_t nameEta[10];    
101   snprintf(nameEta,10, "Eta");    
102   geneta = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEta, "DUMMY");
103   AddSource2Generator(nameEta,geneta);
104   TF1 *fPtEta = geneta->GetPt();
105   fYieldArray[kGenEta] = fPtEta->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
106 // rho  
107   AliGenParam * genrho=0;
108   Char_t nameRho[10];    
109   snprintf(nameRho,10, "Rho");    
110   genrho = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kRho, "DUMMY");
111   AddSource2Generator(nameRho,genrho);
112   TF1 *fPtRho = genrho->GetPt();
113   fYieldArray[kGenRho] = fPtRho->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
114 // omega
115   AliGenParam * genomega=0;
116   Char_t nameOmega[10];    
117   snprintf(nameOmega,10, "Omega");    
118   genomega = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kOmega, "DUMMY");
119   AddSource2Generator(nameOmega,genomega);
120   TF1 *fPtOmega = genomega->GetPt();
121   fYieldArray[kGenOmega] = fPtOmega->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
122 // etaprime
123   AliGenParam * genetaprime=0;
124   Char_t nameEtaprime[10];    
125   snprintf(nameEtaprime,10, "Etaprime");    
126   genetaprime = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kEtaprime, "DUMMY");
127   AddSource2Generator(nameEtaprime,genetaprime);
128   TF1 *fPtEtaprime = genetaprime->GetPt();
129   fYieldArray[kGenEtaprime] = fPtEtaprime->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
130 // phi  
131   AliGenParam * genphi=0;
132   Char_t namePhi[10];    
133   snprintf(namePhi, 10, "Phi");    
134   genphi = new AliGenParam(fNPart, new AliGenEMlib(), AliGenEMlib::kPhi, "DUMMY");
135   AddSource2Generator(namePhi,genphi);
136   TF1 *fPtPhi = genphi->GetPt();
137   fYieldArray[kGenPhi] = fPtPhi->Integral(fPtMin,fPtMax,(Double_t *) 0x0,1.e-6);
138 }
139
140 //-------------------------------------------------------------------
141 void AliGenEMCocktail::AddSource2Generator(Char_t* nameSource, 
142                                          AliGenParam* const genSource)
143 {
144 // add sources to the cocktail
145   Double_t phiMin = fPhiMin*180./TMath::Pi();
146   Double_t phiMax = fPhiMax*180./TMath::Pi();
147
148   genSource->SetPtRange(fPtMin, fPtMax);  
149   genSource->SetYRange(fYMin, fYMax);
150   genSource->SetPhiRange(phiMin, phiMax);
151   genSource->SetWeighting(fWeightingMode);
152   if (!gMC) genSource->SetDecayer(fDecayer);  
153   genSource->Init();
154     
155   AddGenerator(genSource,nameSource,1.); // Adding Generator    
156 }
157
158 //-------------------------------------------------------------------
159 void AliGenEMCocktail::Init()
160 {
161 // Initialisation
162   TIter next(fEntries);
163   AliGenCocktailEntry *entry;
164   if (fStack) {
165     while((entry = (AliGenCocktailEntry*)next())) {
166       entry->Generator()->SetStack(fStack);
167     }
168   }
169 }
170
171 //_________________________________________________________________________
172 void AliGenEMCocktail::Generate()
173 {
174 // Generate event 
175   TIter next(fEntries);
176   AliGenCocktailEntry *entry = 0;
177   AliGenCocktailEntry *preventry = 0;
178   AliGenerator* gen = 0;
179
180   if (fHeader) delete fHeader;
181   fHeader = new AliGenCocktailEventHeader("Electromagnetic Cocktail Header");
182
183   const TObjArray *partArray = gAlice->GetMCApp()->Particles();
184     
185 // Generate the vertex position used by all generators    
186   if(fVertexSmear == kPerEvent) Vertex();
187
188 //Reseting stack
189   AliRunLoader * runloader = AliRunLoader::Instance();
190   if (runloader)
191     if (runloader->Stack())
192       runloader->Stack()->Clean();
193   
194 // Loop over generators and generate events
195   Int_t igen = 0;
196   Float_t evPlane;
197   Rndm(&evPlane,1);
198   evPlane*=TMath::Pi()*2;
199   while((entry = (AliGenCocktailEntry*)next())) {
200     gen = entry->Generator();
201     gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
202     gen->SetTime(fTime);
203     
204     if (fNPart > 0) {
205       igen++;   
206       if (igen == 1) entry->SetFirst(0);                
207       else  entry->SetFirst((partArray->GetEntriesFast())+1);
208       gen->SetEventPlane(evPlane);
209       gen->SetNumberParticles(fNPart);          
210       gen->Generate();
211       entry->SetLast(partArray->GetEntriesFast());
212       preventry = entry;
213     }
214   }  
215   next.Reset();
216
217 // Setting weights for proper absolute normalization
218   Int_t iPart, iMother;
219   Int_t pdgMother = 0;
220   Double_t weight = 0.;
221   Double_t dNdy = 0.;
222   Int_t maxPart = partArray->GetEntriesFast();
223   for(iPart=0; iPart<maxPart; iPart++){      
224     TParticle *part = gAlice->GetMCApp()->Particle(iPart);
225     iMother = part->GetFirstMother();
226     TParticle *mother = 0;
227     if (iMother>=0){
228       mother = gAlice->GetMCApp()->Particle(iMother);
229       pdgMother = mother->GetPdgCode();
230     }
231     else
232       pdgMother = part->GetPdgCode();
233     switch (pdgMother){
234     case 111:
235       dNdy = fYieldArray[kGenPizero];
236       break;
237     case 221:
238       dNdy = fYieldArray[kGenEta];
239       break;
240     case 113:
241       dNdy = fYieldArray[kGenRho];
242       break;
243     case 223:
244       dNdy = fYieldArray[kGenOmega];
245       break;
246     case 331:
247       dNdy = fYieldArray[kGenEtaprime];
248       break;
249     case 333:
250       dNdy = fYieldArray[kGenPhi];
251       break;
252       
253     default:
254       dNdy = 0.;
255     }
256     
257     weight = dNdy*part->GetWeight();
258     part->SetWeight(weight);
259   }     
260   
261   fHeader->SetNProduced(maxPart);
262
263
264   TArrayF eventVertex;
265   eventVertex.Set(3);
266   for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
267   
268   fHeader->SetPrimaryVertex(eventVertex);
269   fHeader->SetInteractionTime(fTime);
270
271   gAlice->SetGenEventHeader(fHeader);
272 }
273
274