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