]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EVGEN/AliGenMUONCocktailpp.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Event header added (Bogdan)
[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
00e6c5ee 28// 2009: added possibility to hide x-sections (B. Vulpescu)
29// 2009: added possibility to have the cocktail (fast generator and param.)
30// for pp @ 10 TeV or pp @ 14 TeV (N. Bastid)
31// 2009: added polarization (L. Bianchi)
103ac317 32
33#include <TObjArray.h>
34#include <TParticle.h>
35#include <TF1.h>
aaa95f22 36#include <TVirtualMC.h>
103ac317 37
38#include "AliGenCocktailEventHeader.h"
39
40#include "AliGenCocktailEntry.h"
41#include "AliGenMUONCocktailpp.h"
42#include "AliGenMUONlib.h"
43#include "AliGenParam.h"
44#include "AliMC.h"
45#include "AliRun.h"
46#include "AliStack.h"
aaa95f22 47#include "AliDecayer.h"
103ac317 48#include "AliLog.h"
b44c3901 49#include "AliGenCorrHF.h"
00e6c5ee 50#include "AliDecayerPolarized.h"
103ac317 51
52ClassImp(AliGenMUONCocktailpp)
53
54//________________________________________________________________________
55AliGenMUONCocktailpp::AliGenMUONCocktailpp()
1c56e311 56 :AliGenCocktail(),
aaa95f22 57 fDecayer(0),
58 fDecayModeResonance(kAll),
59 fDecayModePythia(kAll),
1c56e311 60 fMuonMultiplicity(0),
61 fMuonPtCut(0.),
f81a4aca 62 fMuonPCut(0.),
18bc0c8e 63 fMuonThetaMinCut(0.),
f81a4aca 64 fMuonThetaMaxCut(180.),
aaa95f22 65 fMuonOriginCut(-999.),
18bc0c8e 66 fNSucceded(0),
9ff768ee 67 fNGenerated(0),
00e6c5ee 68
69 fJpsiPol(0),
70 fPsiPPol(0),
71 fUpsPol(0),
72 fUpsPPol(0),
73 fUpsPPPol(0),
74 fPolFrame(0),
75
76//x-sections for pp @ 10 TeV
77 fCMSEnergy(10),
78 fSigmaReaction(0.0695),
79 fSigmaJPsi(39.14e-6),
80 fSigmaPsiP(6.039e-6),
81 fSigmaUpsilon(0.658e-6),
82 fSigmaUpsilonP(0.218e-6),
83 fSigmaUpsilonPP(0.122e-6),
84 fSigmaCCbar(8.9e-3),
85 fSigmaBBbar(0.33e-3),
86
87//x-sections for pp @ 14 TeV
88 /*fCMSEnergy(14),
89 fSigmaReaction(0.070),
9ff768ee 90 fSigmaJPsi(49.44e-6),
91 fSigmaPsiP(7.67e-6),
92 fSigmaUpsilon(0.989e-6),
93 fSigmaUpsilonP(0.502e-6),
94 fSigmaUpsilonPP(0.228e-6),
95 fSigmaCCbar(11.2e-3),
00e6c5ee 96 fSigmaBBbar(0.51e-3),*/
97
9ff768ee 98 fSigmaSilent(kFALSE)
103ac317 99{
100// Constructor
103ac317 101}
93a2041b 102
103ac317 103//_________________________________________________________________________
104AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
105// Destructor
106{
107
108}
109
00e6c5ee 110//_________________________________________________________________________
111void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
112 Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
113
114 if (strcmp(PolFrame,"kColSop")==0){
115 fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0;
116 fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0;
117 fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0;
118 fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0;
119 fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
120 fPolFrame = 0;
121 } else if (strcmp(PolFrame,"kHelicity")==0){
122 fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0;
123 fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0;
124 fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0;
125 fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0;
126 fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
127 fPolFrame = 1;
128
129 } else {
130 AliInfo(Form("The polarization frame is not valid"));
131 AliInfo(Form("No polarization will be set"));
132 fJpsiPol=0.;
133 fPsiPPol=0.;
134 fUpsPol=0.;
135 fUpsPPol=0.;
136 fUpsPPPol=0.;
137 }
138}
139
103ac317 140//_________________________________________________________________________
8058ead1 141void AliGenMUONCocktailpp::CreateCocktail()
103ac317 142{
103ac317 143
144// These limits are only used for renormalization of quarkonia cross section
145// Pythia events are generated in 4pi
146 Double_t ptMin = fPtMin;
147 Double_t ptMax = fPtMax;
148 Double_t yMin = fYMin;;
149 Double_t yMax = fYMax;;
150 Double_t phiMin = fPhiMin*180./TMath::Pi();
151 Double_t phiMax = fPhiMax*180./TMath::Pi();
152 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));
153
154// Ratios with respect to the reaction cross-section in the
155// kinematics limit of the MUONCocktail
156 Double_t ratiojpsi;
157 Double_t ratiopsiP;
158 Double_t ratioupsilon;
159 Double_t ratioupsilonP;
160 Double_t ratioupsilonPP;
b44c3901 161 Double_t ratioccbar;
162 Double_t ratiobbbar;
163
00e6c5ee 164// Beam energy
165 Int_t cmsEnergy = fCMSEnergy;
166
18bc0c8e 167// Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and
168// corrected from feed down of higher resonances
103ac317 169
9ff768ee 170 Double_t sigmajpsi = fSigmaJPsi;
171 Double_t sigmapsiP = fSigmaPsiP;
172 Double_t sigmaupsilon = fSigmaUpsilon;
173 Double_t sigmaupsilonP = fSigmaUpsilonP;
174 Double_t sigmaupsilonPP = fSigmaUpsilonPP;
175 Double_t sigmaccbar = fSigmaCCbar;
176 Double_t sigmabbbar = fSigmaBBbar;
3285b419 177
00e6c5ee 178// Cross sections corrected with the BR in mu+mu-
179// (only in case of use of AliDecayerPolarized)
180
181 if(fJpsiPol != 0) {sigmajpsi = fSigmaJPsi*0.0593;}
182 if(fPsiPPol != 0) {sigmapsiP = fSigmaPsiP*0.0075;}
183 if(fUpsPol != 0) {sigmaupsilon = fSigmaUpsilon*0.0248;}
184 if(fUpsPPol != 0) {sigmaupsilonP = fSigmaUpsilonP*0.0193;}
185 if(fUpsPPPol != 0) {sigmaupsilonPP = fSigmaUpsilonPP*0.0218;}
186
187// MinBias NN cross section @ pp 14 TeV -PR Vol II p:473
188 Double_t sigmaReaction = fSigmaReaction;
189
190 Int_t eincStart = 10;
191
aaa95f22 192 AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
193
18bc0c8e 194// Generation using CDF scaled distribution for pp@14TeV (from D. Stocco)
195//----------------------------------------------------------------------
103ac317 196// Jpsi generator
00e6c5ee 197 AliGenParam * genjpsi;
198 if(cmsEnergy == eincStart){
199 genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 10", "Jpsi");
200 } else {
201 genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp ", "Jpsi");
202 }
103ac317 203// first step: generation in 4pi
204 genjpsi->SetPtRange(0.,100.);
205 genjpsi->SetYRange(-8.,8.);
206 genjpsi->SetPhiRange(0.,360.);
aaa95f22 207 genjpsi->SetForceDecay(fDecayModeResonance);
00e6c5ee 208 //if (!gMC) genjpsi->SetDecayer(fDecayer);
aaa95f22 209
103ac317 210 genjpsi->Init(); // generation in 4pi
211 ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
9ff768ee 212 if (!fSigmaSilent) {
00e6c5ee 213 AliInfo(Form("jpsi prod. cross-section in pp %5.3g b",sigmajpsi));
3285b419 214 AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi));
215 }
103ac317 216// second step: generation in selected kinematical range
217 genjpsi->SetPtRange(ptMin, ptMax);
218 genjpsi->SetYRange(yMin, yMax);
219 genjpsi->SetPhiRange(phiMin, phiMax);
220 genjpsi->Init(); // generation in selected kinematic range
00e6c5ee 221
222 TVirtualMCDecayer *Jpsidec = 0;
223 if(fJpsiPol != 0){
224 AliInfo(Form("******Setting polarized decayer for J/psi"));
225 if(fPolFrame==0) {
226 Jpsidec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
227 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fJpsiPol));
228 }
229 if(fPolFrame==1) {
230 Jpsidec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
231 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fJpsiPol));
232 }
233 Jpsidec->SetForceDecay(kAll);
234 Jpsidec->Init();
235 genjpsi->SetDecayer(Jpsidec);
236 }
237
103ac317 238 AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator
103ac317 239//------------------------------------------------------------------
240// Psi prime generator
00e6c5ee 241 AliGenParam * genpsiP;
242 if(cmsEnergy == eincStart){
243 genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp 10", "PsiP");
244 } else {
245 genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
246 }
103ac317 247// first step: generation in 4pi
248 genpsiP->SetPtRange(0.,100.);
249 genpsiP->SetYRange(-8.,8.);
250 genpsiP->SetPhiRange(0.,360.);
aaa95f22 251 genpsiP->SetForceDecay(fDecayModeResonance);
00e6c5ee 252 //if (!gMC) genpsiP->SetDecayer(fDecayer);
aaa95f22 253
103ac317 254 genpsiP->Init(); // generation in 4pi
255 ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
9ff768ee 256 if (!fSigmaSilent) {
00e6c5ee 257 AliInfo(Form("psiP prod. cross-section in pp %5.3g b",sigmapsiP));
3285b419 258 AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
259 }
103ac317 260// second step: generation in selected kinematical range
261 genpsiP->SetPtRange(ptMin, ptMax);
262 genpsiP->SetYRange(yMin, yMax);
263 genpsiP->SetPhiRange(phiMin, phiMax);
264 genpsiP->Init(); // generation in selected kinematic range
00e6c5ee 265
266 TVirtualMCDecayer *PsiPdec = 0;
267 if(fPsiPPol != 0){
268 AliInfo(Form("******Setting polarized decayer for psi'"));
269 if(fPolFrame==0) {
270 PsiPdec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
271 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fPsiPPol));
272 }
273 if(fPolFrame==1) {
274 PsiPdec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
275 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fPsiPPol));
276 }
277 PsiPdec->SetForceDecay(kAll);
278 PsiPdec->Init();
279 genpsiP->SetDecayer(PsiPdec);
280 }
281
103ac317 282 AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
103ac317 283//------------------------------------------------------------------
284// Upsilon 1S generator
00e6c5ee 285 AliGenParam * genupsilon;
286 if(cmsEnergy == eincStart) {
287 genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp 10", "Upsilon");
288 } else {
289 genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
290 }
103ac317 291// first step: generation in 4pi
292 genupsilon->SetPtRange(0.,100.);
293 genupsilon->SetYRange(-8.,8.);
294 genupsilon->SetPhiRange(0.,360.);
aaa95f22 295 genupsilon->SetForceDecay(fDecayModeResonance);
00e6c5ee 296 //if (!gMC) genupsilon->SetDecayer(fDecayer);
103ac317 297 genupsilon->Init(); // generation in 4pi
298 ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
9ff768ee 299 if (!fSigmaSilent) {
00e6c5ee 300 AliInfo(Form("Upsilon 1S prod. cross-section in pp %5.3g b",sigmaupsilon));
3285b419 301 AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
302 }
103ac317 303// second step: generation in selected kinematical range
304 genupsilon->SetPtRange(ptMin, ptMax);
305 genupsilon->SetYRange(yMin, yMax);
306 genupsilon->SetPhiRange(phiMin, phiMax);
307 genupsilon->Init(); // generation in selected kinematical range
00e6c5ee 308
309 TVirtualMCDecayer *Upsdec = 0;
310 if(fUpsPol != 0){
311 AliInfo(Form("******Setting polarized decayer for Upsilon"));
312 if(fPolFrame==0) {
313 Upsdec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
314 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPol));
315 }
316 if(fPolFrame==1) {
317 Upsdec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
318 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPol));
319 }
320 Upsdec->SetForceDecay(kAll);
321 Upsdec->Init();
322 genupsilon->SetDecayer(Upsdec);
323 }
324
103ac317 325 AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
103ac317 326//------------------------------------------------------------------
327// Upsilon 2S generator
00e6c5ee 328 AliGenParam * genupsilonP;
329 if(cmsEnergy == eincStart){
330 genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp 10", "UpsilonP");
331 } else {
332 genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
333 }
334
103ac317 335// first step: generation in 4pi
336 genupsilonP->SetPtRange(0.,100.);
337 genupsilonP->SetYRange(-8.,8.);
338 genupsilonP->SetPhiRange(0.,360.);
aaa95f22 339 genupsilonP->SetForceDecay(fDecayModeResonance);
00e6c5ee 340 //if (!gMC) genupsilonP->SetDecayer(fDecayer);
103ac317 341 genupsilonP->Init(); // generation in 4pi
342 ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
9ff768ee 343 if (!fSigmaSilent) {
00e6c5ee 344 AliInfo(Form("Upsilon 2S prod. cross-section in pp %5.3g b",sigmaupsilonP));
3285b419 345 AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
346 }
103ac317 347// second step: generation in the kinematical range
348 genupsilonP->SetPtRange(ptMin, ptMax);
349 genupsilonP->SetYRange(yMin, yMax);
350 genupsilonP->SetPhiRange(phiMin, phiMax);
351 genupsilonP->Init(); // generation in selected kinematical range
00e6c5ee 352
353 TVirtualMCDecayer *UpsPdec = 0;
354 if(fUpsPPol != 0){
355 AliInfo(Form("******Setting polarized decayer for Upsilon'"));
356 if(fPolFrame==0) {
357 UpsPdec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
358 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPol));
359 }
360 if(fPolFrame==1) {
361 UpsPdec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
362 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPol));
363 }
364 UpsPdec->SetForceDecay(kAll);
365 UpsPdec->Init();
366 genupsilonP->SetDecayer(UpsPdec);
367 }
368
103ac317 369 AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
103ac317 370
371//------------------------------------------------------------------
372// Upsilon 3S generator
00e6c5ee 373 AliGenParam * genupsilonPP;
374 if(cmsEnergy == eincStart){
375 genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp 10", "UpsilonPP");
376 }
377 else {
378 genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");
379 }
380
103ac317 381// first step: generation in 4pi
382 genupsilonPP->SetPtRange(0.,100.);
383 genupsilonPP->SetYRange(-8.,8.);
384 genupsilonPP->SetPhiRange(0.,360.);
aaa95f22 385 genupsilonPP->SetForceDecay(fDecayModeResonance);
00e6c5ee 386 //if (!gMC) genupsilonPP->SetDecayer(fDecayer);
103ac317 387 genupsilonPP->Init(); // generation in 4pi
388 ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
9ff768ee 389 if (!fSigmaSilent) {
00e6c5ee 390 AliInfo(Form("Upsilon 3S prod. cross-section in pp %5.3g b",sigmaupsilonPP));
3285b419 391 AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
392 }
103ac317 393// second step: generation in selected kinematical range
394 genupsilonPP->SetPtRange(ptMin, ptMax);
395 genupsilonPP->SetYRange(yMin, yMax);
396 genupsilonPP->SetPhiRange(phiMin, phiMax);
397 genupsilonPP->Init(); // generation in selected kinematical range
00e6c5ee 398
399 TVirtualMCDecayer *UpsPPdec = 0;
400 if(fUpsPPPol != 0){
401 AliInfo(Form("******Setting polarized decayer for Upsilon''"));
402 if(fPolFrame==0) {
403 UpsPPdec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
404 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPPol));
405 }
406 if(fPolFrame==1) {
407 UpsPPdec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
408 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPPol));
409 }
410 UpsPPdec->SetForceDecay(kAll);
411 UpsPPdec->Init();
412 genupsilonPP->SetDecayer(UpsPPdec);
413 }
414
103ac317 415 AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
b44c3901 416
417//------------------------------------------------------------------
418// Generator of charm
419
00e6c5ee 420 AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
8058ead1 421 gencharm->SetMomentumRange(0,9999);
422 gencharm->SetForceDecay(kAll);
423 ratioccbar = sigmaccbar/sigmaReaction;
00e6c5ee 424 //if (!gMC) gencharm->SetDecayer(fDecayer);
8058ead1 425 gencharm->Init();
9ff768ee 426 if (!fSigmaSilent) {
00e6c5ee 427 AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
3285b419 428 AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
429 }
8058ead1 430 AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
b44c3901 431
103ac317 432//------------------------------------------------------------------
b44c3901 433// Generator of beauty
434
00e6c5ee 435 AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
8058ead1 436 genbeauty->SetMomentumRange(0,9999);
437 genbeauty->SetForceDecay(kAll);
438 ratiobbbar = sigmabbbar/sigmaReaction;
00e6c5ee 439 //if (!gMC) genbeauty->SetDecayer(fDecayer);
8058ead1 440 genbeauty->Init();
9ff768ee 441 if (!fSigmaSilent) {
00e6c5ee 442 AliInfo(Form("b-bbar prod. cross-section in pp %5.3g b",sigmabbbar));
3285b419 443 AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
444 }
8058ead1 445 AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
b44c3901 446
8058ead1 447//-------------------------------------------------------------------
103ac317 448// Pythia generator
8058ead1 449//
450// This has to go into the Config.C
451//
452// AliGenPythia *pythia = new AliGenPythia(1);
453// pythia->SetProcess(kPyMbMSEL1);
454// pythia->SetStrucFunc(kCTEQ5L);
455// pythia->SetEnergyCMS(14000.);
456// AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
457// Decay_t dt = gener->GetDecayModePythia();
458// pythia->SetForceDecay(dt);
459// pythia->SetPtRange(0.,100.);
460// pythia->SetYRange(-8.,8.);
461// pythia->SetPhiRange(0.,360.);
462// pythia->SetPtHard(2.76,-1.0);
463// pythia->SwitchHFOff();
464// pythia->Init();
465// AddGenerator(pythia,"Pythia",1);
466}
103ac317 467
8058ead1 468void AliGenMUONCocktailpp::Init()
469{
470 //
471 // Initialisation
aaa95f22 472 TIter next(fEntries);
473 AliGenCocktailEntry *entry;
474 if (fStack) {
475 while((entry = (AliGenCocktailEntry*)next())) {
476 entry->Generator()->SetStack(fStack);
477 }
478 }
103ac317 479}
480
481//_________________________________________________________________________
482void AliGenMUONCocktailpp::Generate()
483{
484// Generate event
485 TIter next(fEntries);
486 AliGenCocktailEntry *entry = 0;
487 AliGenCocktailEntry *preventry = 0;
488 AliGenerator* gen = 0;
489
01c89ce1 490 if (fHeader) delete fHeader;
491 fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
492
d94af0c1 493 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
103ac317 494
495// Generate the vertex position used by all generators
496 if(fVertexSmear == kPerEvent) Vertex();
497
498// Loop on primordialTrigger:
499// minimum muon multiplicity above a pt cut in a theta acceptance region
500
501 Bool_t primordialTrigger = kFALSE;
502 while(!primordialTrigger) {
503 //Reseting stack
33c3c91a 504 AliRunLoader * runloader = AliRunLoader::Instance();
103ac317 505 if (runloader)
506 if (runloader->Stack())
34da8494 507 runloader->Stack()->Clean();
103ac317 508 // Loop over generators and generate events
509 Int_t igen = 0;
510 Int_t npart = 0;
511 const char* genName = 0;
512 while((entry = (AliGenCocktailEntry*)next())) {
513 gen = entry->Generator();
514 genName = entry->GetName();
515 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
516
517 npart = (strcmp(genName,"Pythia") == 0) ? 1 :
518 gRandom->Poisson(entry->Rate());
519
520 if (npart > 0) {
521 igen++;
522 if (igen == 1) entry->SetFirst(0);
523 else entry->SetFirst((partArray->GetEntriesFast())+1);
524
525 gen->SetNumberParticles(npart);
526 gen->Generate();
527 entry->SetLast(partArray->GetEntriesFast());
528 preventry = entry;
529 }
530 }
531 next.Reset();
532
533// Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
534// in the muon spectrometer acceptance
535 Int_t iPart;
536 fNGenerated++;
b44c3901 537 Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
aaa95f22 538 for(iPart=0; iPart<maxPart; iPart++){
103ac317 539
aaa95f22 540 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
541 if ( TMath::Abs(part->GetPdgCode()) == 13 ){
542 if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
543 (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
544 (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
f81a4aca 545 (part->Pt()>fMuonPtCut) &&
546 (part->P()>fMuonPCut)) {
aaa95f22 547 numberOfMuons++;
103ac317 548 }
aaa95f22 549 }
550 }
01c89ce1 551 if (numberOfMuons >= fMuonMultiplicity) {
552 primordialTrigger = kTRUE;
553 fHeader->SetNProduced(maxPart);
554 }
555
103ac317 556 }
557 fNSucceded++;
558
01c89ce1 559 TArrayF eventVertex;
560 eventVertex.Set(3);
561 for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
562
563 fHeader->SetPrimaryVertex(eventVertex);
564
565 gAlice->SetGenEventHeader(fHeader);
566
aaa95f22 567// AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
103ac317 568 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
569}
570
571