Coding rule violations corrected.
[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)
ba8bf3f5 31//-----------------
00e6c5ee 32// 2009: added polarization (L. Bianchi)
ba8bf3f5 33//------------------
34// 11/2009: added chi_c1 & chi_c2 (P.Crochet & N.Bastid).
35// Cross-sections for charmonia are now directly taken from the Yellow Report
36// (hep-ph/0311048) Tab.9, page 19. See below for details w.r.t. beam energy.
37// usage: see example of Config in $ALICE_ROOT/prod/LHC09a10/Config.C
103ac317 38
39#include <TObjArray.h>
40#include <TParticle.h>
41#include <TF1.h>
aaa95f22 42#include <TVirtualMC.h>
103ac317 43
44#include "AliGenCocktailEventHeader.h"
45
46#include "AliGenCocktailEntry.h"
47#include "AliGenMUONCocktailpp.h"
48#include "AliGenMUONlib.h"
49#include "AliGenParam.h"
50#include "AliMC.h"
51#include "AliRun.h"
52#include "AliStack.h"
aaa95f22 53#include "AliDecayer.h"
103ac317 54#include "AliLog.h"
b44c3901 55#include "AliGenCorrHF.h"
00e6c5ee 56#include "AliDecayerPolarized.h"
103ac317 57
58ClassImp(AliGenMUONCocktailpp)
59
60//________________________________________________________________________
61AliGenMUONCocktailpp::AliGenMUONCocktailpp()
1c56e311 62 :AliGenCocktail(),
aaa95f22 63 fDecayer(0),
64 fDecayModeResonance(kAll),
65 fDecayModePythia(kAll),
1c56e311 66 fMuonMultiplicity(0),
67 fMuonPtCut(0.),
f81a4aca 68 fMuonPCut(0.),
18bc0c8e 69 fMuonThetaMinCut(0.),
f81a4aca 70 fMuonThetaMaxCut(180.),
aaa95f22 71 fMuonOriginCut(-999.),
18bc0c8e 72 fNSucceded(0),
9ff768ee 73 fNGenerated(0),
00e6c5ee 74
75 fJpsiPol(0),
ba8bf3f5 76 fChic1Pol(0),
77 fChic2Pol(0),
00e6c5ee 78 fPsiPPol(0),
79 fUpsPol(0),
80 fUpsPPol(0),
81 fUpsPPPol(0),
82 fPolFrame(0),
83
ba8bf3f5 84// x-sections for pp @ 7 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
85// bottomnium as for 10 TeV
86/* fCMSEnergy(7),
87 fSigmaReaction(0.0695),
88 fSigmaJPsi(21.8e-6),
89 fSigmaChic1(21.1e-6),
90 fSigmaChic2(34.9e-6),
91 fSigmaPsiP(4.93e-6),
92 fSigmaUpsilon(0.463e-6),
93 fSigmaUpsilonP(0.154e-6),
94 fSigmaUpsilonPP(0.0886e-6),
95 fSigmaCCbar(6.91e-3),
96 fSigmaBBbar(0.232e-3),
97*/
98
99//x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers
100// scaled down according to ccbar and bbar cross-sections
00e6c5ee 101 fCMSEnergy(10),
102 fSigmaReaction(0.0695),
ba8bf3f5 103 fSigmaJPsi(26.06e-6),
104 fSigmaChic1(25.18e-6),
105 fSigmaChic2(41.58e-6),
106 fSigmaPsiP(5.88e-6),
00e6c5ee 107 fSigmaUpsilon(0.658e-6),
108 fSigmaUpsilonP(0.218e-6),
109 fSigmaUpsilonPP(0.122e-6),
110 fSigmaCCbar(8.9e-3),
111 fSigmaBBbar(0.33e-3),
112
ba8bf3f5 113
114//x-sections for pp @ 14 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
115// bottomonium from hep-ph/0311048 Tab.9, page 19 taken inton account that
116// feed-down from chib is included
117/* fCMSEnergy(14),
00e6c5ee 118 fSigmaReaction(0.070),
ba8bf3f5 119 fSigmaJPsi(32.9e-6),
120 fSigmaChic1(31.8e-6),
121 fSigmaChic2(52.5e-6),
122 fSigmaPsiP(7.43e-6),
9ff768ee 123 fSigmaUpsilon(0.989e-6),
124 fSigmaUpsilonP(0.502e-6),
125 fSigmaUpsilonPP(0.228e-6),
126 fSigmaCCbar(11.2e-3),
ba8bf3f5 127 fSigmaBBbar(0.51e-3),
128*/
9ff768ee 129 fSigmaSilent(kFALSE)
103ac317 130{
131// Constructor
103ac317 132}
93a2041b 133
103ac317 134//_________________________________________________________________________
135AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
136// Destructor
137{
138
139}
140
141//_________________________________________________________________________
00e6c5ee 142void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
143 Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
144
145 if (strcmp(PolFrame,"kColSop")==0){
146 fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0;
147 fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0;
148 fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0;
149 fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0;
150 fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
151 fPolFrame = 0;
152 } else if (strcmp(PolFrame,"kHelicity")==0){
153 fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0;
154 fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0;
155 fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0;
156 fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0;
157 fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
158 fPolFrame = 1;
159
160 } else {
161 AliInfo(Form("The polarization frame is not valid"));
162 AliInfo(Form("No polarization will be set"));
163 fJpsiPol=0.;
164 fPsiPPol=0.;
165 fUpsPol=0.;
166 fUpsPPol=0.;
167 fUpsPPPol=0.;
168 }
169}
170
171//_________________________________________________________________________
8058ead1 172void AliGenMUONCocktailpp::CreateCocktail()
103ac317 173{
103ac317 174
175// These limits are only used for renormalization of quarkonia cross section
176// Pythia events are generated in 4pi
177 Double_t ptMin = fPtMin;
178 Double_t ptMax = fPtMax;
179 Double_t yMin = fYMin;;
180 Double_t yMax = fYMax;;
181 Double_t phiMin = fPhiMin*180./TMath::Pi();
182 Double_t phiMax = fPhiMax*180./TMath::Pi();
183 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));
184
185// Ratios with respect to the reaction cross-section in the
186// kinematics limit of the MUONCocktail
187 Double_t ratiojpsi;
ba8bf3f5 188 Double_t ratiochic1;
189 Double_t ratiochic2;
103ac317 190 Double_t ratiopsiP;
191 Double_t ratioupsilon;
192 Double_t ratioupsilonP;
193 Double_t ratioupsilonPP;
b44c3901 194 Double_t ratioccbar;
195 Double_t ratiobbbar;
196
00e6c5ee 197// Beam energy
198 Int_t cmsEnergy = fCMSEnergy;
199
18bc0c8e 200// Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and
201// corrected from feed down of higher resonances
103ac317 202
9ff768ee 203 Double_t sigmajpsi = fSigmaJPsi;
ba8bf3f5 204 Double_t sigmachic1 = fSigmaChic1;
205 Double_t sigmachic2 = fSigmaChic2;
9ff768ee 206 Double_t sigmapsiP = fSigmaPsiP;
207 Double_t sigmaupsilon = fSigmaUpsilon;
208 Double_t sigmaupsilonP = fSigmaUpsilonP;
209 Double_t sigmaupsilonPP = fSigmaUpsilonPP;
210 Double_t sigmaccbar = fSigmaCCbar;
211 Double_t sigmabbbar = fSigmaBBbar;
3285b419 212
00e6c5ee 213// Cross sections corrected with the BR in mu+mu-
214// (only in case of use of AliDecayerPolarized)
215
e5ab62dd 216 if(TMath::Abs(fJpsiPol) > 1.e-30) {sigmajpsi = fSigmaJPsi*0.0593;}
217 if(TMath::Abs(fChic1Pol) > 1.e-30) {sigmachic1 = fSigmaChic1*0.;} // tb consistent
218 if(TMath::Abs(fChic2Pol) > 1.e-30) {sigmachic2 = fSigmaChic2*0.;} // tb consistent
219 if(TMath::Abs(fPsiPPol) > 1.e-30) {sigmapsiP = fSigmaPsiP*0.0075;}
220 if(TMath::Abs(fUpsPol) > 1.e-30) {sigmaupsilon = fSigmaUpsilon*0.0248;}
221 if(TMath::Abs(fUpsPPol) > 1.e-30) {sigmaupsilonP = fSigmaUpsilonP*0.0193;}
222 if(TMath::Abs(fUpsPPPol) > 1.e-30) {sigmaupsilonPP = fSigmaUpsilonPP*0.0218;}
00e6c5ee 223
224// MinBias NN cross section @ pp 14 TeV -PR Vol II p:473
225 Double_t sigmaReaction = fSigmaReaction;
226
227 Int_t eincStart = 10;
228
aaa95f22 229 AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
230
18bc0c8e 231// Generation using CDF scaled distribution for pp@14TeV (from D. Stocco)
232//----------------------------------------------------------------------
103ac317 233// Jpsi generator
00e6c5ee 234 AliGenParam * genjpsi;
235 if(cmsEnergy == eincStart){
236 genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp 10", "Jpsi");
237 } else {
238 genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp ", "Jpsi");
239 }
103ac317 240// first step: generation in 4pi
241 genjpsi->SetPtRange(0.,100.);
242 genjpsi->SetYRange(-8.,8.);
243 genjpsi->SetPhiRange(0.,360.);
aaa95f22 244 genjpsi->SetForceDecay(fDecayModeResonance);
00e6c5ee 245 //if (!gMC) genjpsi->SetDecayer(fDecayer);
aaa95f22 246
103ac317 247 genjpsi->Init(); // generation in 4pi
248 ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
9ff768ee 249 if (!fSigmaSilent) {
00e6c5ee 250 AliInfo(Form("jpsi prod. cross-section in pp %5.3g b",sigmajpsi));
3285b419 251 AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi));
252 }
103ac317 253// second step: generation in selected kinematical range
254 genjpsi->SetPtRange(ptMin, ptMax);
255 genjpsi->SetYRange(yMin, yMax);
256 genjpsi->SetPhiRange(phiMin, phiMax);
257 genjpsi->Init(); // generation in selected kinematic range
00e6c5ee 258
e5ab62dd 259 TVirtualMCDecayer *jpsiDec = 0;
260 if(TMath::Abs(fJpsiPol) > 1.e-30){
00e6c5ee 261 AliInfo(Form("******Setting polarized decayer for J/psi"));
262 if(fPolFrame==0) {
e5ab62dd 263 jpsiDec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
00e6c5ee 264 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fJpsiPol));
265 }
266 if(fPolFrame==1) {
e5ab62dd 267 jpsiDec = new AliDecayerPolarized(fJpsiPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
00e6c5ee 268 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fJpsiPol));
269 }
e5ab62dd 270 jpsiDec->SetForceDecay(kAll);
271 jpsiDec->Init();
272 genjpsi->SetDecayer(jpsiDec);
00e6c5ee 273 }
274
103ac317 275 AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator
ba8bf3f5 276
277
278//----------------------------------------------------------------------
279// Chic1 generator
280 AliGenParam * genchic1;
281 if(cmsEnergy == eincStart){
282 genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp 10", "Chic1");
283 } else {
284 genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, "CDF pp ", "Chic1");
285 }
286// first step: generation in 4pi
287 genchic1->SetPtRange(0.,100.);
288 genchic1->SetYRange(-8.,8.);
289 genchic1->SetPhiRange(0.,360.);
290 genchic1->SetForceDecay(fDecayModeResonance);
291 //if (!gMC) genchic1->SetDecayer(fDecayer);
292
293 genchic1->Init(); // generation in 4pi
294 ratiochic1 = sigmachic1 / sigmaReaction * genchic1->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
295 if (!fSigmaSilent) {
296 AliInfo(Form("chic1 prod. cross-section in pp %5.3g b",sigmachic1));
297 AliInfo(Form("chic1 prod. probability per collision in acceptance %5.3g",ratiochic1));
298 }
299// second step: generation in selected kinematical range
300 genchic1->SetPtRange(ptMin, ptMax);
301 genchic1->SetYRange(yMin, yMax);
302 genchic1->SetPhiRange(phiMin, phiMax);
303 genchic1->Init(); // generation in selected kinematic range
304
e5ab62dd 305 TVirtualMCDecayer *chic1Dec = 0;
ba8bf3f5 306 if(fChic1Pol != 0){
307 AliInfo(Form("******Setting polarized decayer for Chic1"));
308 if(fPolFrame==0) {
e5ab62dd 309 chic1Dec = new AliDecayerPolarized(fChic1Pol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
ba8bf3f5 310 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fChic1Pol));
311 }
312 if(fPolFrame==1) {
e5ab62dd 313 chic1Dec = new AliDecayerPolarized(fChic1Pol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
ba8bf3f5 314 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fChic1Pol));
315 }
e5ab62dd 316 chic1Dec->SetForceDecay(kAll);
317 chic1Dec->Init();
318 genchic1->SetDecayer(chic1Dec);
ba8bf3f5 319 }
320
321 AddGenerator(genchic1, "Chic1", ratiochic1); // Adding Generator
322
323//----------------------------------------------------------------------
324// Chic2 generator
325 AliGenParam * genchic2;
326 if(cmsEnergy == eincStart){
327 genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp 10", "Chic2");
328 } else {
329 genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, "CDF pp ", "Chic2");
330 }
331// first step: generation in 4pi
332 genchic2->SetPtRange(0.,100.);
333 genchic2->SetYRange(-8.,8.);
334 genchic2->SetPhiRange(0.,360.);
335 genchic2->SetForceDecay(fDecayModeResonance);
336 //if (!gMC) genchic1->SetDecayer(fDecayer);
337
338 genchic2->Init(); // generation in 4pi
339 ratiochic2 = sigmachic2 / sigmaReaction * genchic2->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
340 if (!fSigmaSilent) {
341 AliInfo(Form("chic2 prod. cross-section in pp %5.3g b",sigmachic2));
342 AliInfo(Form("chic2 prod. probability per collision in acceptance %5.3g",ratiochic2));
343 }
344// second step: generation in selected kinematical range
345 genchic2->SetPtRange(ptMin, ptMax);
346 genchic2->SetYRange(yMin, yMax);
347 genchic2->SetPhiRange(phiMin, phiMax);
348 genchic2->Init(); // generation in selected kinematic range
349
e5ab62dd 350 TVirtualMCDecayer *chic2Dec = 0;
ba8bf3f5 351 if(fChic2Pol != 0){
352 AliInfo(Form("******Setting polarized decayer for Chic2"));
353 if(fPolFrame==0) {
e5ab62dd 354 chic2Dec = new AliDecayerPolarized(fChic2Pol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
ba8bf3f5 355 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fChic2Pol));
356 }
357 if(fPolFrame==1) {
e5ab62dd 358 chic2Dec = new AliDecayerPolarized(fChic2Pol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
ba8bf3f5 359 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fChic2Pol));
360 }
e5ab62dd 361 chic2Dec->SetForceDecay(kAll);
362 chic2Dec->Init();
363 genchic2->SetDecayer(chic2Dec);
ba8bf3f5 364 }
365
366 AddGenerator(genchic2, "Chic2", ratiochic2); // Adding Generator
367
103ac317 368//------------------------------------------------------------------
369// Psi prime generator
00e6c5ee 370 AliGenParam * genpsiP;
371 if(cmsEnergy == eincStart){
372 genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp 10", "PsiP");
373 } else {
374 genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
375 }
103ac317 376// first step: generation in 4pi
377 genpsiP->SetPtRange(0.,100.);
378 genpsiP->SetYRange(-8.,8.);
379 genpsiP->SetPhiRange(0.,360.);
aaa95f22 380 genpsiP->SetForceDecay(fDecayModeResonance);
00e6c5ee 381 //if (!gMC) genpsiP->SetDecayer(fDecayer);
aaa95f22 382
103ac317 383 genpsiP->Init(); // generation in 4pi
384 ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
9ff768ee 385 if (!fSigmaSilent) {
00e6c5ee 386 AliInfo(Form("psiP prod. cross-section in pp %5.3g b",sigmapsiP));
3285b419 387 AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
388 }
103ac317 389// second step: generation in selected kinematical range
390 genpsiP->SetPtRange(ptMin, ptMax);
391 genpsiP->SetYRange(yMin, yMax);
392 genpsiP->SetPhiRange(phiMin, phiMax);
393 genpsiP->Init(); // generation in selected kinematic range
00e6c5ee 394
e5ab62dd 395 TVirtualMCDecayer *psipDec = 0;
396 if(TMath::Abs(fPsiPPol) > 1.e-30){
00e6c5ee 397 AliInfo(Form("******Setting polarized decayer for psi'"));
398 if(fPolFrame==0) {
e5ab62dd 399 psipDec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
00e6c5ee 400 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fPsiPPol));
401 }
402 if(fPolFrame==1) {
e5ab62dd 403 psipDec = new AliDecayerPolarized(fPsiPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
00e6c5ee 404 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fPsiPPol));
405 }
e5ab62dd 406 psipDec->SetForceDecay(kAll);
407 psipDec->Init();
408 genpsiP->SetDecayer(psipDec);
00e6c5ee 409 }
410
103ac317 411 AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
ba8bf3f5 412
103ac317 413//------------------------------------------------------------------
414// Upsilon 1S generator
00e6c5ee 415 AliGenParam * genupsilon;
416 if(cmsEnergy == eincStart) {
417 genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp 10", "Upsilon");
418 } else {
419 genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
420 }
103ac317 421// first step: generation in 4pi
422 genupsilon->SetPtRange(0.,100.);
423 genupsilon->SetYRange(-8.,8.);
424 genupsilon->SetPhiRange(0.,360.);
aaa95f22 425 genupsilon->SetForceDecay(fDecayModeResonance);
00e6c5ee 426 //if (!gMC) genupsilon->SetDecayer(fDecayer);
103ac317 427 genupsilon->Init(); // generation in 4pi
428 ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
9ff768ee 429 if (!fSigmaSilent) {
00e6c5ee 430 AliInfo(Form("Upsilon 1S prod. cross-section in pp %5.3g b",sigmaupsilon));
3285b419 431 AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
432 }
103ac317 433// second step: generation in selected kinematical range
434 genupsilon->SetPtRange(ptMin, ptMax);
435 genupsilon->SetYRange(yMin, yMax);
436 genupsilon->SetPhiRange(phiMin, phiMax);
437 genupsilon->Init(); // generation in selected kinematical range
00e6c5ee 438
e5ab62dd 439 TVirtualMCDecayer *upsDec = 0;
440 if(TMath::Abs(fUpsPol) > 1.e-30){
00e6c5ee 441 AliInfo(Form("******Setting polarized decayer for Upsilon"));
442 if(fPolFrame==0) {
e5ab62dd 443 upsDec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
00e6c5ee 444 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPol));
445 }
446 if(fPolFrame==1) {
e5ab62dd 447 upsDec = new AliDecayerPolarized(fUpsPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
00e6c5ee 448 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPol));
449 }
e5ab62dd 450 upsDec->SetForceDecay(kAll);
451 upsDec->Init();
452 genupsilon->SetDecayer(upsDec);
00e6c5ee 453 }
454
103ac317 455 AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
103ac317 456//------------------------------------------------------------------
457// Upsilon 2S generator
00e6c5ee 458 AliGenParam * genupsilonP;
459 if(cmsEnergy == eincStart){
460 genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp 10", "UpsilonP");
461 } else {
462 genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
463 }
464
103ac317 465// first step: generation in 4pi
466 genupsilonP->SetPtRange(0.,100.);
467 genupsilonP->SetYRange(-8.,8.);
468 genupsilonP->SetPhiRange(0.,360.);
aaa95f22 469 genupsilonP->SetForceDecay(fDecayModeResonance);
00e6c5ee 470 //if (!gMC) genupsilonP->SetDecayer(fDecayer);
103ac317 471 genupsilonP->Init(); // generation in 4pi
472 ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
9ff768ee 473 if (!fSigmaSilent) {
00e6c5ee 474 AliInfo(Form("Upsilon 2S prod. cross-section in pp %5.3g b",sigmaupsilonP));
3285b419 475 AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
476 }
103ac317 477// second step: generation in the kinematical range
478 genupsilonP->SetPtRange(ptMin, ptMax);
479 genupsilonP->SetYRange(yMin, yMax);
480 genupsilonP->SetPhiRange(phiMin, phiMax);
481 genupsilonP->Init(); // generation in selected kinematical range
00e6c5ee 482
e5ab62dd 483 TVirtualMCDecayer *upspDec = 0;
484 if(TMath::Abs(fUpsPPol) > 1.e-30){
00e6c5ee 485 AliInfo(Form("******Setting polarized decayer for Upsilon'"));
486 if(fPolFrame==0) {
e5ab62dd 487 upspDec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
00e6c5ee 488 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPol));
489 }
490 if(fPolFrame==1) {
e5ab62dd 491 upspDec = new AliDecayerPolarized(fUpsPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
00e6c5ee 492 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPol));
493 }
e5ab62dd 494 upspDec->SetForceDecay(kAll);
495 upspDec->Init();
496 genupsilonP->SetDecayer(upspDec);
00e6c5ee 497 }
498
103ac317 499 AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
103ac317 500
501//------------------------------------------------------------------
502// Upsilon 3S generator
00e6c5ee 503 AliGenParam * genupsilonPP;
504 if(cmsEnergy == eincStart){
505 genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp 10", "UpsilonPP");
506 }
507 else {
508 genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");
509 }
510
103ac317 511// first step: generation in 4pi
512 genupsilonPP->SetPtRange(0.,100.);
513 genupsilonPP->SetYRange(-8.,8.);
514 genupsilonPP->SetPhiRange(0.,360.);
aaa95f22 515 genupsilonPP->SetForceDecay(fDecayModeResonance);
00e6c5ee 516 //if (!gMC) genupsilonPP->SetDecayer(fDecayer);
103ac317 517 genupsilonPP->Init(); // generation in 4pi
518 ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
9ff768ee 519 if (!fSigmaSilent) {
00e6c5ee 520 AliInfo(Form("Upsilon 3S prod. cross-section in pp %5.3g b",sigmaupsilonPP));
3285b419 521 AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
522 }
103ac317 523// second step: generation in selected kinematical range
524 genupsilonPP->SetPtRange(ptMin, ptMax);
525 genupsilonPP->SetYRange(yMin, yMax);
526 genupsilonPP->SetPhiRange(phiMin, phiMax);
527 genupsilonPP->Init(); // generation in selected kinematical range
00e6c5ee 528
e5ab62dd 529 TVirtualMCDecayer *upsppDec = 0;
00e6c5ee 530 if(fUpsPPPol != 0){
531 AliInfo(Form("******Setting polarized decayer for Upsilon''"));
532 if(fPolFrame==0) {
e5ab62dd 533 upsppDec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
00e6c5ee 534 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",fUpsPPPol));
535 }
536 if(fPolFrame==1) {
e5ab62dd 537 upsppDec = new AliDecayerPolarized(fUpsPPPol,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
00e6c5ee 538 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",fUpsPPPol));
539 }
e5ab62dd 540 upsppDec->SetForceDecay(kAll);
541 upsppDec->Init();
542 genupsilonPP->SetDecayer(upsppDec);
00e6c5ee 543 }
544
103ac317 545 AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
b44c3901 546
547//------------------------------------------------------------------
548// Generator of charm
549
00e6c5ee 550 AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
8058ead1 551 gencharm->SetMomentumRange(0,9999);
552 gencharm->SetForceDecay(kAll);
553 ratioccbar = sigmaccbar/sigmaReaction;
00e6c5ee 554 //if (!gMC) gencharm->SetDecayer(fDecayer);
8058ead1 555 gencharm->Init();
9ff768ee 556 if (!fSigmaSilent) {
00e6c5ee 557 AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
3285b419 558 AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
559 }
8058ead1 560 AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
b44c3901 561
103ac317 562//------------------------------------------------------------------
b44c3901 563// Generator of beauty
564
00e6c5ee 565 AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
8058ead1 566 genbeauty->SetMomentumRange(0,9999);
567 genbeauty->SetForceDecay(kAll);
568 ratiobbbar = sigmabbbar/sigmaReaction;
00e6c5ee 569 //if (!gMC) genbeauty->SetDecayer(fDecayer);
8058ead1 570 genbeauty->Init();
9ff768ee 571 if (!fSigmaSilent) {
00e6c5ee 572 AliInfo(Form("b-bbar prod. cross-section in pp %5.3g b",sigmabbbar));
3285b419 573 AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
574 }
8058ead1 575 AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
b44c3901 576
8058ead1 577//-------------------------------------------------------------------
103ac317 578// Pythia generator
8058ead1 579//
580// This has to go into the Config.C
581//
582// AliGenPythia *pythia = new AliGenPythia(1);
583// pythia->SetProcess(kPyMbMSEL1);
584// pythia->SetStrucFunc(kCTEQ5L);
585// pythia->SetEnergyCMS(14000.);
586// AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
587// Decay_t dt = gener->GetDecayModePythia();
588// pythia->SetForceDecay(dt);
589// pythia->SetPtRange(0.,100.);
590// pythia->SetYRange(-8.,8.);
591// pythia->SetPhiRange(0.,360.);
592// pythia->SetPtHard(2.76,-1.0);
593// pythia->SwitchHFOff();
594// pythia->Init();
595// AddGenerator(pythia,"Pythia",1);
ba8bf3f5 596
8058ead1 597}
103ac317 598
8058ead1 599void AliGenMUONCocktailpp::Init()
600{
601 //
602 // Initialisation
aaa95f22 603 TIter next(fEntries);
604 AliGenCocktailEntry *entry;
605 if (fStack) {
606 while((entry = (AliGenCocktailEntry*)next())) {
607 entry->Generator()->SetStack(fStack);
608 }
609 }
103ac317 610}
611
612//_________________________________________________________________________
613void AliGenMUONCocktailpp::Generate()
614{
615// Generate event
616 TIter next(fEntries);
617 AliGenCocktailEntry *entry = 0;
618 AliGenCocktailEntry *preventry = 0;
619 AliGenerator* gen = 0;
620
01c89ce1 621 if (fHeader) delete fHeader;
622 fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
623
d94af0c1 624 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
103ac317 625
626// Generate the vertex position used by all generators
627 if(fVertexSmear == kPerEvent) Vertex();
628
629// Loop on primordialTrigger:
630// minimum muon multiplicity above a pt cut in a theta acceptance region
631
632 Bool_t primordialTrigger = kFALSE;
633 while(!primordialTrigger) {
634 //Reseting stack
33c3c91a 635 AliRunLoader * runloader = AliRunLoader::Instance();
103ac317 636 if (runloader)
637 if (runloader->Stack())
34da8494 638 runloader->Stack()->Clean();
103ac317 639 // Loop over generators and generate events
640 Int_t igen = 0;
641 Int_t npart = 0;
642 const char* genName = 0;
643 while((entry = (AliGenCocktailEntry*)next())) {
644 gen = entry->Generator();
645 genName = entry->GetName();
646 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
647
648 npart = (strcmp(genName,"Pythia") == 0) ? 1 :
649 gRandom->Poisson(entry->Rate());
650
651 if (npart > 0) {
652 igen++;
653 if (igen == 1) entry->SetFirst(0);
654 else entry->SetFirst((partArray->GetEntriesFast())+1);
655
656 gen->SetNumberParticles(npart);
657 gen->Generate();
658 entry->SetLast(partArray->GetEntriesFast());
659 preventry = entry;
660 }
661 }
662 next.Reset();
663
664// Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
665// in the muon spectrometer acceptance
666 Int_t iPart;
667 fNGenerated++;
b44c3901 668 Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
aaa95f22 669 for(iPart=0; iPart<maxPart; iPart++){
103ac317 670
aaa95f22 671 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
672 if ( TMath::Abs(part->GetPdgCode()) == 13 ){
673 if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
674 (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
675 (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
f81a4aca 676 (part->Pt()>fMuonPtCut) &&
677 (part->P()>fMuonPCut)) {
aaa95f22 678 numberOfMuons++;
103ac317 679 }
aaa95f22 680 }
681 }
01c89ce1 682 if (numberOfMuons >= fMuonMultiplicity) {
683 primordialTrigger = kTRUE;
684 fHeader->SetNProduced(maxPart);
685 }
686
103ac317 687 }
688 fNSucceded++;
689
01c89ce1 690 TArrayF eventVertex;
691 eventVertex.Set(3);
692 for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
693
694 fHeader->SetPrimaryVertex(eventVertex);
695
696 gAlice->SetGenEventHeader(fHeader);
697
aaa95f22 698// AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
103ac317 699 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
700}
701
702