Fixing coding violations and removing warnings from compilation.
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktailpp.cxx
CommitLineData
103ac317 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// Classe to create the coktail for physics with muons for pp at 14 TeV
20// The followoing sources:
21// jpsi,psiP, upsilon, upsilonP, upsilonPP are added to Pythia
22// The free parameeters are :
23// pp reaction cross-section
24// production cross-sections in pp collisions
25//
26
27#include <TObjArray.h>
28#include <TParticle.h>
29#include <TF1.h>
aaa95f22 30#include <TVirtualMC.h>
103ac317 31
32#include "AliGenCocktailEventHeader.h"
33
34#include "AliGenCocktailEntry.h"
35#include "AliGenMUONCocktailpp.h"
36#include "AliGenMUONlib.h"
37#include "AliGenParam.h"
38#include "AliMC.h"
39#include "AliRun.h"
40#include "AliStack.h"
aaa95f22 41#include "AliDecayer.h"
103ac317 42#include "AliLog.h"
43#include "AliGenPythia.h"
44
45
46ClassImp(AliGenMUONCocktailpp)
47
48//________________________________________________________________________
49AliGenMUONCocktailpp::AliGenMUONCocktailpp()
1c56e311 50 :AliGenCocktail(),
aaa95f22 51 fDecayer(0),
52 fDecayModeResonance(kAll),
53 fDecayModePythia(kAll),
1c56e311 54 fTotalRate(0),
55 fMuonMultiplicity(0),
56 fMuonPtCut(0.),
18bc0c8e 57 fMuonThetaMinCut(0.),
1c56e311 58 fMuonThetaMaxCut(0.),
aaa95f22 59 fMuonOriginCut(-999.),
18bc0c8e 60 fNSucceded(0),
1c56e311 61 fNGenerated(0)
103ac317 62{
63// Constructor
103ac317 64}
93a2041b 65
103ac317 66//_________________________________________________________________________
67AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
68// Destructor
69{
70
71}
72
73//_________________________________________________________________________
74void AliGenMUONCocktailpp::Init()
75{
76// MinBias NN cross section @ pp 14 TeV -PR Vol II p:473
77 Double_t sigmaReaction = 0.070;
78
79// These limits are only used for renormalization of quarkonia cross section
80// Pythia events are generated in 4pi
81 Double_t ptMin = fPtMin;
82 Double_t ptMax = fPtMax;
83 Double_t yMin = fYMin;;
84 Double_t yMax = fYMax;;
85 Double_t phiMin = fPhiMin*180./TMath::Pi();
86 Double_t phiMax = fPhiMax*180./TMath::Pi();
87 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));
88
89// Ratios with respect to the reaction cross-section in the
90// kinematics limit of the MUONCocktail
91 Double_t ratiojpsi;
92 Double_t ratiopsiP;
93 Double_t ratioupsilon;
94 Double_t ratioupsilonP;
95 Double_t ratioupsilonPP;
96
18bc0c8e 97// Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and
98// corrected from feed down of higher resonances
103ac317 99
18bc0c8e 100 Double_t sigmajpsi = 49.44e-6;
103ac317 101 Double_t sigmapsiP = 7.67e-6;
18bc0c8e 102 Double_t sigmaupsilon = 0.989e-6;
103 Double_t sigmaupsilonP = 0.502e-6;
103ac317 104 Double_t sigmaupsilonPP = 0.228e-6;
105
aaa95f22 106 AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
107
18bc0c8e 108// Generation using CDF scaled distribution for pp@14TeV (from D. Stocco)
109//----------------------------------------------------------------------
103ac317 110// Jpsi generator
18bc0c8e 111 AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp", "Jpsi");
103ac317 112// first step: generation in 4pi
113 genjpsi->SetPtRange(0.,100.);
114 genjpsi->SetYRange(-8.,8.);
115 genjpsi->SetPhiRange(0.,360.);
aaa95f22 116 genjpsi->SetForceDecay(fDecayModeResonance);
117 if (!gMC) genjpsi->SetDecayer(fDecayer);
118
103ac317 119 genjpsi->Init(); // generation in 4pi
120 ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
121 AliInfo(Form("jpsi prod. cross-section in pp(14 TeV) %5.3g b",sigmajpsi));
122 AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi));
123// second step: generation in selected kinematical range
124 genjpsi->SetPtRange(ptMin, ptMax);
125 genjpsi->SetYRange(yMin, yMax);
126 genjpsi->SetPhiRange(phiMin, phiMax);
127 genjpsi->Init(); // generation in selected kinematic range
128 AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator
129 fTotalRate+=ratiojpsi;
130
131//------------------------------------------------------------------
132// Psi prime generator
18bc0c8e 133 AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
103ac317 134// first step: generation in 4pi
135 genpsiP->SetPtRange(0.,100.);
136 genpsiP->SetYRange(-8.,8.);
137 genpsiP->SetPhiRange(0.,360.);
aaa95f22 138 genpsiP->SetForceDecay(fDecayModeResonance);
139 if (!gMC) genpsiP->SetDecayer(fDecayer);
140
103ac317 141 genpsiP->Init(); // generation in 4pi
142 ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
143 AliInfo(Form("psiP prod. cross-section in pp(14 TeV) %5.3g b",sigmapsiP));
144 AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
145// second step: generation in selected kinematical range
146 genpsiP->SetPtRange(ptMin, ptMax);
147 genpsiP->SetYRange(yMin, yMax);
148 genpsiP->SetPhiRange(phiMin, phiMax);
149 genpsiP->Init(); // generation in selected kinematic range
150 AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
151 fTotalRate+=ratiopsiP;
152
153//------------------------------------------------------------------
154// Upsilon 1S generator
18bc0c8e 155 AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
103ac317 156// first step: generation in 4pi
157 genupsilon->SetPtRange(0.,100.);
158 genupsilon->SetYRange(-8.,8.);
159 genupsilon->SetPhiRange(0.,360.);
aaa95f22 160 genupsilon->SetForceDecay(fDecayModeResonance);
161 if (!gMC) genupsilon->SetDecayer(fDecayer);
103ac317 162 genupsilon->Init(); // generation in 4pi
163 ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
164 AliInfo(Form("Upsilon 1S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilon));
165 AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
166// second step: generation in selected kinematical range
167 genupsilon->SetPtRange(ptMin, ptMax);
168 genupsilon->SetYRange(yMin, yMax);
169 genupsilon->SetPhiRange(phiMin, phiMax);
170 genupsilon->Init(); // generation in selected kinematical range
171 AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
172 fTotalRate+=ratioupsilon;
173
174//------------------------------------------------------------------
175// Upsilon 2S generator
18bc0c8e 176 AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
103ac317 177// first step: generation in 4pi
178 genupsilonP->SetPtRange(0.,100.);
179 genupsilonP->SetYRange(-8.,8.);
180 genupsilonP->SetPhiRange(0.,360.);
aaa95f22 181 genupsilonP->SetForceDecay(fDecayModeResonance);
182 if (!gMC) genupsilonP->SetDecayer(fDecayer);
103ac317 183 genupsilonP->Init(); // generation in 4pi
184 ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
185 AliInfo(Form("Upsilon 2S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonP));
186 AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
187// second step: generation in the kinematical range
188 genupsilonP->SetPtRange(ptMin, ptMax);
189 genupsilonP->SetYRange(yMin, yMax);
190 genupsilonP->SetPhiRange(phiMin, phiMax);
191 genupsilonP->Init(); // generation in selected kinematical range
192 AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
193 fTotalRate+=ratioupsilonP;
194
195//------------------------------------------------------------------
196// Upsilon 3S generator
18bc0c8e 197 AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");
103ac317 198// first step: generation in 4pi
199 genupsilonPP->SetPtRange(0.,100.);
200 genupsilonPP->SetYRange(-8.,8.);
201 genupsilonPP->SetPhiRange(0.,360.);
aaa95f22 202 genupsilonPP->SetForceDecay(fDecayModeResonance);
203 if (!gMC) genupsilonPP->SetDecayer(fDecayer);
103ac317 204 genupsilonPP->Init(); // generation in 4pi
205 ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
206 AliInfo(Form("Upsilon 3S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonPP));
207 AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
208// second step: generation in selected kinematical range
209 genupsilonPP->SetPtRange(ptMin, ptMax);
210 genupsilonPP->SetYRange(yMin, yMax);
211 genupsilonPP->SetPhiRange(phiMin, phiMax);
212 genupsilonPP->Init(); // generation in selected kinematical range
213 AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
214 fTotalRate+=ratioupsilonPP;
103ac317 215//------------------------------------------------------------------
216// Pythia generator
217 AliGenPythia *pythia = new AliGenPythia(1);
218 pythia->SetProcess(kPyMbMSEL1);
219 pythia->SetStrucFunc(kCTEQ5L);
220 pythia->SetEnergyCMS(14000.);
aaa95f22 221 AliInfo(Form("\n\npythia uses the decay mode %d",fDecayModePythia));
222 pythia->SetForceDecay(fDecayModePythia);
103ac317 223 pythia->SetPtRange(0.,100.);
224 pythia->SetYRange(-8.,8.);
225 pythia->SetPhiRange(0.,360.);
f58d19bb 226 pythia->SetPtHard(2.76,-1.0);
103ac317 227 pythia->Init();
228 AddGenerator(pythia,"Pythia",1);
229 fTotalRate+=1.;
230
aaa95f22 231 TIter next(fEntries);
232 AliGenCocktailEntry *entry;
233 if (fStack) {
234 while((entry = (AliGenCocktailEntry*)next())) {
235 entry->Generator()->SetStack(fStack);
236 }
237 }
103ac317 238}
239
240//_________________________________________________________________________
241void AliGenMUONCocktailpp::Generate()
242{
243// Generate event
244 TIter next(fEntries);
245 AliGenCocktailEntry *entry = 0;
246 AliGenCocktailEntry *preventry = 0;
247 AliGenerator* gen = 0;
248
249 TObjArray *partArray = gAlice->GetMCApp()->Particles();
250
251// Generate the vertex position used by all generators
252 if(fVertexSmear == kPerEvent) Vertex();
253
254// Loop on primordialTrigger:
255// minimum muon multiplicity above a pt cut in a theta acceptance region
256
257 Bool_t primordialTrigger = kFALSE;
258 while(!primordialTrigger) {
259 //Reseting stack
260 AliRunLoader * runloader = gAlice->GetRunLoader();
261 if (runloader)
262 if (runloader->Stack())
263 runloader->Stack()->Reset();
264 // Loop over generators and generate events
265 Int_t igen = 0;
266 Int_t npart = 0;
267 const char* genName = 0;
268 while((entry = (AliGenCocktailEntry*)next())) {
269 gen = entry->Generator();
270 genName = entry->GetName();
271 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
272
273 npart = (strcmp(genName,"Pythia") == 0) ? 1 :
274 gRandom->Poisson(entry->Rate());
275
276 if (npart > 0) {
277 igen++;
278 if (igen == 1) entry->SetFirst(0);
279 else entry->SetFirst((partArray->GetEntriesFast())+1);
280
281 gen->SetNumberParticles(npart);
282 gen->Generate();
283 entry->SetLast(partArray->GetEntriesFast());
284 preventry = entry;
285 }
286 }
287 next.Reset();
288
289// Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
290// in the muon spectrometer acceptance
291 Int_t iPart;
292 fNGenerated++;
293 Int_t numberOfMuons=0;
aaa95f22 294
295 Int_t maxPart = partArray->GetEntriesFast();
296 for(iPart=0; iPart<maxPart; iPart++){
103ac317 297
aaa95f22 298 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
299 if ( TMath::Abs(part->GetPdgCode()) == 13 ){
300 if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
301 (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
302 (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
303 (part->Pt()>fMuonPtCut)) {
304 numberOfMuons++;
103ac317 305 }
aaa95f22 306 }
307 }
103ac317 308 if (numberOfMuons >= fMuonMultiplicity) primordialTrigger = kTRUE;
309 }
310 fNSucceded++;
311
aaa95f22 312// AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
103ac317 313 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
314}
315
316