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