restore threshold getters after parameter dynamics update (fw v. >= A012)
[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//
ba762ab4 19// Class to create the coktail for physics with muons for pp collisions
20// using the The followoing sources:
21// jpsi, psiP, upsilon, upsilonP, upsilonPP, open charm and open beauty
103ac317 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
ba762ab4 38//------------------------
39// 04/2010:
40// - CMS energy passed via parameter
41// i.e. gener->SetCMSEnergy(AliGenMUONCocktailpp::kCMS07TeV) in Config.C
42// - resonances now added to the cocktail via AddReso2Generator
43// - cleaning
7455632e 44// B.Vulpescu & P.Crochet
45//-----------------------
46// 10/2011:
47// - added the cocktail for p-Pb & Pb-p @ 8.8 TeV with 4 centrality bins and
48// for Pb-Pb @ 2.76 TeV with 11 centrality bins. Bins should be defined also
49// in the Config.C with one AliGenMUONCocktailpp per bin. These generators
50// included in a AliGenCocktail together with an event generator (e.g. Hijing)
51// providing the underlying event and collision centrality. The bin number n
52// passed via AliGenMUONCocktailpp::SetCentralityBin(n).
53// See details in my presentation at the PWG3-Muon meeting (05.10.2011):
54// https://indico.cern.ch/conferenceDisplay.py?confId=157367
55// - simplifications and bug fix in CreateCocktail()
56// S. Grigoryan
ba762ab4 57
103ac317 58#include <TObjArray.h>
59#include <TParticle.h>
60#include <TF1.h>
aaa95f22 61#include <TVirtualMC.h>
103ac317 62#include "AliGenCocktailEventHeader.h"
63
64#include "AliGenCocktailEntry.h"
65#include "AliGenMUONCocktailpp.h"
66#include "AliGenMUONlib.h"
67#include "AliGenParam.h"
68#include "AliMC.h"
69#include "AliRun.h"
70#include "AliStack.h"
aaa95f22 71#include "AliDecayer.h"
103ac317 72#include "AliLog.h"
b44c3901 73#include "AliGenCorrHF.h"
00e6c5ee 74#include "AliDecayerPolarized.h"
103ac317 75
76ClassImp(AliGenMUONCocktailpp)
77
78//________________________________________________________________________
79AliGenMUONCocktailpp::AliGenMUONCocktailpp()
1c56e311 80 :AliGenCocktail(),
aaa95f22 81 fDecayer(0),
82 fDecayModeResonance(kAll),
83 fDecayModePythia(kAll),
1c56e311 84 fMuonMultiplicity(0),
85 fMuonPtCut(0.),
f81a4aca 86 fMuonPCut(0.),
18bc0c8e 87 fMuonThetaMinCut(0.),
f81a4aca 88 fMuonThetaMaxCut(180.),
aaa95f22 89 fMuonOriginCut(-999.),
18bc0c8e 90 fNSucceded(0),
9ff768ee 91 fNGenerated(0),
7455632e 92 fCentralityBin(0),
00e6c5ee 93
94 fJpsiPol(0),
ba8bf3f5 95 fChic1Pol(0),
96 fChic2Pol(0),
00e6c5ee 97 fPsiPPol(0),
98 fUpsPol(0),
99 fUpsPPol(0),
100 fUpsPPPol(0),
101 fPolFrame(0),
102
ba762ab4 103 fCMSEnergyTeV(0),
104 fCMSEnergyTeVArray(),
105 fSigmaReaction(0),
106 fSigmaReactionArray(),
107 fSigmaJPsi(0),
108 fSigmaJPsiArray(),
109 fSigmaChic1(0),
110 fSigmaChic1Array(),
111 fSigmaChic2(0),
112 fSigmaChic2Array(),
113 fSigmaPsiP(0),
114 fSigmaPsiPArray(),
115 fSigmaUpsilon(0),
116 fSigmaUpsilonArray(),
117 fSigmaUpsilonP(0),
118 fSigmaUpsilonPArray(),
119 fSigmaUpsilonPP(0),
120 fSigmaUpsilonPPArray(),
121 fSigmaCCbar(0),
122 fSigmaCCbarArray(),
123 fSigmaBBbar(0),
124 fSigmaBBbarArray(),
125
126 fSigmaSilent(kFALSE)
127{
128// Constructor
129
7455632e 130// x-sections for pp @ 7 TeV:
131// -charmonia: 4pi integral of fit function for inclusive J/psi dsigma/dy LHC data
132// gives 60 mub; so sigma_prompt = 54 mub, while Ref = R.Vogt_arXiv:1003.3497 (Table 2)
133// gives 35 mub. Below we use sigma_direct from the Ref scaled by the factor 54/35.
134// -bottomonia: 4pi integral of fit function for inclusive Upsilon1S dsigma/dy LHC data
135// gives 0.56 mub, sigmas for 2S & 3S obtained using CMS data for ratios 2S/1S & 3S/1S
136// -ccbar & bbbar: NLO pQCD computations - http://www-alice.gsi.de/ana/MNR/results.html
ba762ab4 137 fCMSEnergyTeVArray[0] = 7.00;
7455632e 138 fSigmaReactionArray[0] = 0.070;
139 fSigmaJPsiArray[0] = 33.6e-6;
140 fSigmaChic1Array[0] = 32.6e-6;
141 fSigmaChic2Array[0] = 53.8e-6;
142 fSigmaPsiPArray[0] = 7.6e-6;
143 fSigmaUpsilonArray[0] = 0.56e-6;
144 fSigmaUpsilonPArray[0] = 0.19e-6;
145 fSigmaUpsilonPPArray[0] = 0.09e-6;
ba762ab4 146 fSigmaCCbarArray[0] = 6.91e-3;
147 fSigmaBBbarArray[0] = 0.232e-3;
148
ba8bf3f5 149//x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers
7455632e 150// scaled down according to ccbar and bbbar cross-sections
ba762ab4 151 fCMSEnergyTeVArray[1] = 10.00;
7455632e 152 fSigmaReactionArray[1] = 0.070;
ba762ab4 153 fSigmaJPsiArray[1] = 26.06e-6;
154 fSigmaChic1Array[1] = 25.18e-6;
155 fSigmaChic2Array[1] = 41.58e-6;
156 fSigmaPsiPArray[1] = 5.88e-6;
157 fSigmaUpsilonArray[1] = 0.658e-6;
158 fSigmaUpsilonPArray[1] = 0.218e-6;
159 fSigmaUpsilonPPArray[1] = 0.122e-6;
160 fSigmaCCbarArray[1] = 8.9e-3;
161 fSigmaBBbarArray[1] = 0.33e-3;
7455632e 162
ba8bf3f5 163//x-sections for pp @ 14 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
7455632e 164// bottomonium from hep-ph/0311048 Tab.9, page 19 taken into account that
ba8bf3f5 165// feed-down from chib is included
ba762ab4 166 fCMSEnergyTeVArray[2] = 14.00;
167 fSigmaReactionArray[2] = 0.070;
168 fSigmaJPsiArray[2] = 32.9e-6;
169 fSigmaChic1Array[2] = 31.8e-6;
170 fSigmaChic2Array[2] = 52.5e-6;
171 fSigmaPsiPArray[2] = 7.43e-6;
172 fSigmaUpsilonArray[2] = 0.989e-6;
173 fSigmaUpsilonPArray[2] = 0.502e-6;
174 fSigmaUpsilonPPArray[2] = 0.228e-6;
175 fSigmaCCbarArray[2] = 11.2e-3;
7455632e 176 fSigmaBBbarArray[2] = 0.445e-3;
177
178// x-sections for Min. Bias p-Pb & Pb-p @ 8.8 TeV: charmonia and bottomonia
179// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
180// and with Glauber scaling
181 fCMSEnergyTeVArray[3] = 9.00; // for 8.8 TeV
182 fSigmaReactionArray[3] = 2.10;
183 fSigmaJPsiArray[3] = 8.19e-3; // 208*1.172*33.6e-6
184 fSigmaChic1Array[3] = 7.95e-3;
185 fSigmaChic2Array[3] = 13.1e-3;
186 fSigmaPsiPArray[3] = 1.85e-3;
187 fSigmaUpsilonArray[3] = 0.146e-3; // 208*1.25*0.56e-6
188 fSigmaUpsilonPArray[3] = 0.049e-3;
189 fSigmaUpsilonPPArray[3] = 0.023e-3;
190 fSigmaCCbarArray[3] = 1.68; // 208*8.1e-3
191 fSigmaBBbarArray[3] = 0.061; // 208*0.29e-3
192
193 fCMSEnergyTeVArray[4] = -fCMSEnergyTeVArray[3];
194 fSigmaReactionArray[4] = fSigmaReactionArray[3];
195 fSigmaJPsiArray[4] = fSigmaJPsiArray[3];
196 fSigmaChic1Array[4] = fSigmaChic1Array[3];
197 fSigmaChic2Array[4] = fSigmaChic2Array[3];
198 fSigmaPsiPArray[4] = fSigmaPsiPArray[3];
199 fSigmaUpsilonArray[4] = fSigmaUpsilonArray[3];
200 fSigmaUpsilonPArray[4] = fSigmaUpsilonPArray[3];
201 fSigmaUpsilonPPArray[4] = fSigmaUpsilonPPArray[3];
202 fSigmaCCbarArray[4] = fSigmaCCbarArray[3];
203 fSigmaBBbarArray[4] = fSigmaBBbarArray[3];
204
205// x-sections for Min. Bias Pb-Pb @ 2.76 TeV: charmonia and bottomonia
206// from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
207// and with Glauber scaling
208 fCMSEnergyTeVArray[5] = 3.00; // for 2.76 TeV
209 fSigmaReactionArray[5] = 7.65;
210 fSigmaJPsiArray[5] = 0.734; // 208*208*0.505*33.6e-6
211 fSigmaChic1Array[5] = 0.712;
212 fSigmaChic2Array[5] = 1.175;
213 fSigmaPsiPArray[5] = 0.166;
214 fSigmaUpsilonArray[5] = 0.0092; // 208*208*0.379*0.56e-6
215 fSigmaUpsilonPArray[5] = 0.0031;
216 fSigmaUpsilonPPArray[5] = 0.0015;
217 fSigmaCCbarArray[5] = 151.; // 208*208*3.49e-3
218 fSigmaBBbarArray[5] = 3.8; // 208*208*0.088e-3
ba762ab4 219
103ac317 220}
93a2041b 221
103ac317 222//_________________________________________________________________________
223AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
ba762ab4 224{
103ac317 225// Destructor
ba762ab4 226
227}
228
229//_________________________________________________________________________
230void AliGenMUONCocktailpp::SetCMSEnergy(CMSEnergyCode cmsEnergy)
103ac317 231{
ba762ab4 232// setter for CMSEnergy and corresponding cross-sections
233 fCMSEnergyTeV = fCMSEnergyTeVArray[cmsEnergy];
234 fSigmaReaction = fSigmaReactionArray[cmsEnergy];
235 fSigmaJPsi = fSigmaJPsiArray[cmsEnergy];
236 fSigmaChic1 = fSigmaChic1Array[cmsEnergy];
237 fSigmaChic2 = fSigmaChic2Array[cmsEnergy];
238 fSigmaPsiP = fSigmaPsiPArray[cmsEnergy];
239 fSigmaUpsilon = fSigmaUpsilonArray[cmsEnergy];
240 fSigmaUpsilonP = fSigmaUpsilonPArray[cmsEnergy];
241 fSigmaUpsilonPP = fSigmaUpsilonPPArray[cmsEnergy];
242 fSigmaCCbar = fSigmaCCbarArray[cmsEnergy];
243 fSigmaBBbar = fSigmaBBbarArray[cmsEnergy];
103ac317 244}
245
246//_________________________________________________________________________
00e6c5ee 247void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
248 Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
ba762ab4 249// setter for resonances polarization
00e6c5ee 250 if (strcmp(PolFrame,"kColSop")==0){
ba762ab4 251 fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0;
252 fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0;
253 fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0;
254 fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0;
255 fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
256 fPolFrame = 0;
00e6c5ee 257 } else if (strcmp(PolFrame,"kHelicity")==0){
ba762ab4 258 fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0;
259 fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0;
260 fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0;
261 fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0;
262 fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
263 fPolFrame = 1;
264
00e6c5ee 265 } else {
266 AliInfo(Form("The polarization frame is not valid"));
267 AliInfo(Form("No polarization will be set"));
268 fJpsiPol=0.;
269 fPsiPPol=0.;
270 fUpsPol=0.;
271 fUpsPPol=0.;
272 fUpsPPPol=0.;
273 }
274}
275
276//_________________________________________________________________________
8058ead1 277void AliGenMUONCocktailpp::CreateCocktail()
103ac317 278{
ba762ab4 279// create and add resonances and open HF to the coctail
280 Int_t cmsEnergy = Int_t(fCMSEnergyTeV);
281
103ac317 282// These limits are only used for renormalization of quarkonia cross section
283// Pythia events are generated in 4pi
284 Double_t ptMin = fPtMin;
285 Double_t ptMax = fPtMax;
286 Double_t yMin = fYMin;;
287 Double_t yMax = fYMax;;
288 Double_t phiMin = fPhiMin*180./TMath::Pi();
289 Double_t phiMax = fPhiMax*180./TMath::Pi();
290 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));
291
18bc0c8e 292// Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and
293// corrected from feed down of higher resonances
103ac317 294
9ff768ee 295 Double_t sigmajpsi = fSigmaJPsi;
ba8bf3f5 296 Double_t sigmachic1 = fSigmaChic1;
297 Double_t sigmachic2 = fSigmaChic2;
9ff768ee 298 Double_t sigmapsiP = fSigmaPsiP;
299 Double_t sigmaupsilon = fSigmaUpsilon;
300 Double_t sigmaupsilonP = fSigmaUpsilonP;
301 Double_t sigmaupsilonPP = fSigmaUpsilonPP;
302 Double_t sigmaccbar = fSigmaCCbar;
303 Double_t sigmabbbar = fSigmaBBbar;
3285b419 304
00e6c5ee 305// Cross sections corrected with the BR in mu+mu-
306// (only in case of use of AliDecayerPolarized)
307
ba762ab4 308 if(TMath::Abs(fJpsiPol) > 1.e-30) {sigmajpsi = fSigmaJPsi*0.0593;}
309 if(TMath::Abs(fChic1Pol) > 1.e-30) {sigmachic1 = fSigmaChic1*0.;} // tb consistent
310 if(TMath::Abs(fChic2Pol) > 1.e-30) {sigmachic2 = fSigmaChic2*0.;} // tb consistent
311 if(TMath::Abs(fPsiPPol) > 1.e-30) {sigmapsiP = fSigmaPsiP*0.0075;}
312 if(TMath::Abs(fUpsPol) > 1.e-30) {sigmaupsilon = fSigmaUpsilon*0.0248;}
313 if(TMath::Abs(fUpsPPol) > 1.e-30) {sigmaupsilonP = fSigmaUpsilonP*0.0193;}
314 if(TMath::Abs(fUpsPPPol) > 1.e-30) {sigmaupsilonPP = fSigmaUpsilonPP*0.0218;}
00e6c5ee 315
aaa95f22 316 AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
317
ba762ab4 318// Create and add resonances to the generator
319 AliGenParam * genjpsi=0;
320 AliGenParam * genchic1=0;
321 AliGenParam * genchic2=0;
322 AliGenParam * genpsiP=0;
323 AliGenParam * genupsilon=0;
324 AliGenParam * genupsilonP=0;
325 AliGenParam * genupsilonPP=0;
326
327 Char_t nameJpsi[10];
328 Char_t nameChic1[10];
329 Char_t nameChic2[10];
330 Char_t namePsiP[10];
331 Char_t nameUps[10];
332 Char_t nameUpsP[10];
333 Char_t nameUpsPP[10];
334
c478a271 335 snprintf(nameJpsi,10, "Jpsi");
336 snprintf(nameChic1,10, "Chic1");
337 snprintf(nameChic2,10, "Chic2");
338 snprintf(namePsiP,10, "PsiP");
339 snprintf(nameUps,10, "Ups");
340 snprintf(nameUpsP,10, "UpsP");
341 snprintf(nameUpsPP,10, "UpsPP");
ba762ab4 342
7455632e 343 Char_t tname[40] = "";
344 if(cmsEnergy == 10) {snprintf(tname, 40, "CDF pp 10");
345 } else if (cmsEnergy == 14){snprintf(tname, 40, "CDF pp");
346 } else if (cmsEnergy == 7) {snprintf(tname, 40, "pp 7");
347 // } else if (cmsEnergy == 2) {snprintf(tname, 40, "pp 2.76");
348 } else if (cmsEnergy == 9) {snprintf(tname, 40, "pPb 8.8");
349 if (fCentralityBin > 0) snprintf(tname, 40, "pPb 8.8c%d",fCentralityBin);
350 } else if (cmsEnergy == -9){snprintf(tname, 40, "Pbp 8.8");
351 if (fCentralityBin > 0) snprintf(tname, 40, "Pbp 8.8c%d",fCentralityBin);
352 } else if (cmsEnergy == 3) {snprintf(tname, 40, "PbPb 2.76");
353 if (fCentralityBin > 0) snprintf(tname, 40, "PbPb 2.76c%d",fCentralityBin);
418e091f 354 } else {
7455632e 355 AliError("Initialisation failed, wrong cmsEnergy");
418e091f 356 return;
00e6c5ee 357 }
7455632e 358 genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, tname, "Jpsi");
359 genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, tname, "Chic1");
360 genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, tname, "Chic2");
361 genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, tname, "PsiP");
362 genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, tname, "Upsilon");
363 genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, tname, "UpsilonP");
364 genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, tname, "UpsilonPP");
365
366// Hard process yield per pA or AA collision for i-th centrality bin is R*r[i]*shad[i]
367// where R is the ratio of hard and geometrical x-sections, r[i] is the ratio of these
368// x-section fractions for given centrality and shad[i] is the shadowing factor (in 4pi).
369// The latter is assumed to be the same for HF-hadrons & quarkonia of the same flavour.
370 Int_t i = 0;
371 Double_t chard[20] = {0}; // charm & beauty shadowing factors are different
372 Double_t bhard[20] = {0};
373 chard[0] = 1; // 1st element for pp and min. bias (MB) collisions
374 bhard[0] = 1;
375
376// 4 centrality bins for p-Pb & Pb-p: 0-20-40-60-100 %
377 if (cmsEnergy == 9 || cmsEnergy == -9) {
378 const Int_t n9 = 5; // 1st element for MB collisions
379 Double_t r9[n9] = {1, 1.936, 1.473, 0.914, 0.333}; // ratio of hard-over-geo fractions
380 Double_t cshad9[n9] = {0.785, 0.715, 0.775, 0.856, 0.951};// EKS98 shadowing factors
381 Double_t bshad9[n9] = {0.889, 0.853, 0.884, 0.926, 0.975};
382 for(i=0; i<n9; i++) {
383 chard[i] = cshad9[i]*r9[i];
384 bhard[i] = bshad9[i]*r9[i];
385 }
386 }
387
388// 11 centrality bins for Pb-Pb: 0-5-10-20-30-40-50-60-70-80-90-100 %
389 if (cmsEnergy == 3) {
390 const Int_t n3 = 12; // 1st element for MB collisions
391 Double_t r3[n3] = {1, 4.661, 3.647, 2.551, 1.544, 0.887, 0.474,
392 0.235, 0.106, 0.044, 0.017, 0.007}; // ratio of hard-over-geo fractions
393 Double_t cshad3[n3] = {0.662, 0.622, 0.631, 0.650, 0.681, 0.718,
394 0.760, 0.805, 0.849, 0.888, 0.918, 0.944};// EKS98 shadowing factors
395 Double_t bshad3[n3] = {0.874, 0.856, 0.861, 0.869, 0.883, 0.898,
396 0.915, 0.932, 0.948, 0.962, 0.972, 0.981};
397 for(i=0; i<n3; i++) {
398 chard[i] = cshad3[i]*r3[i];
399 bhard[i] = bshad3[i]*r3[i];
400 }
401 }
402
403 AddReso2Generator(nameJpsi,genjpsi,chard[fCentralityBin]*sigmajpsi,fJpsiPol);
404 AddReso2Generator(nameChic1,genchic1,chard[fCentralityBin]*sigmachic1,fChic1Pol);
405 AddReso2Generator(nameChic2,genchic2,chard[fCentralityBin]*sigmachic2,fChic2Pol);
406 AddReso2Generator(namePsiP,genpsiP,chard[fCentralityBin]*sigmapsiP,fPsiPPol);
00e6c5ee 407
7455632e 408 AddReso2Generator(nameUps,genupsilon,bhard[fCentralityBin]*sigmaupsilon,fUpsPol);
409 AddReso2Generator(nameUpsP,genupsilonP,bhard[fCentralityBin]*sigmaupsilonP,fUpsPPol);
410 AddReso2Generator(nameUpsPP,genupsilonPP,bhard[fCentralityBin]*sigmaupsilonPP,fUpsPPPol);
b44c3901 411
412//------------------------------------------------------------------
413// Generator of charm
ba762ab4 414 AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
8058ead1 415 gencharm->SetMomentumRange(0,9999);
416 gencharm->SetForceDecay(kAll);
7455632e 417 Double_t ratioccbar = chard[fCentralityBin]*sigmaccbar/fSigmaReaction;
ba762ab4 418 if (!gMC) gencharm->SetDecayer(fDecayer);
8058ead1 419 gencharm->Init();
9ff768ee 420 if (!fSigmaSilent) {
00e6c5ee 421 AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
3285b419 422 AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
423 }
8058ead1 424 AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
103ac317 425//------------------------------------------------------------------
b44c3901 426// Generator of beauty
ba762ab4 427 AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
8058ead1 428 genbeauty->SetMomentumRange(0,9999);
429 genbeauty->SetForceDecay(kAll);
7455632e 430 Double_t ratiobbbar = bhard[fCentralityBin]*sigmabbbar/fSigmaReaction;
ba762ab4 431 if (!gMC) genbeauty->SetDecayer(fDecayer);
8058ead1 432 genbeauty->Init();
9ff768ee 433 if (!fSigmaSilent) {
00e6c5ee 434 AliInfo(Form("b-bbar prod. cross-section in pp %5.3g b",sigmabbbar));
3285b419 435 AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
436 }
ba762ab4 437 AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
b44c3901 438
8058ead1 439//-------------------------------------------------------------------
103ac317 440// Pythia generator
8058ead1 441//
442// This has to go into the Config.C
443//
444// AliGenPythia *pythia = new AliGenPythia(1);
445// pythia->SetProcess(kPyMbMSEL1);
446// pythia->SetStrucFunc(kCTEQ5L);
447// pythia->SetEnergyCMS(14000.);
448// AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
449// Decay_t dt = gener->GetDecayModePythia();
450// pythia->SetForceDecay(dt);
451// pythia->SetPtRange(0.,100.);
452// pythia->SetYRange(-8.,8.);
453// pythia->SetPhiRange(0.,360.);
454// pythia->SetPtHard(2.76,-1.0);
455// pythia->SwitchHFOff();
456// pythia->Init();
457// AddGenerator(pythia,"Pythia",1);
ba8bf3f5 458
8058ead1 459}
103ac317 460
ba762ab4 461//-------------------------------------------------------------------
462void AliGenMUONCocktailpp::AddReso2Generator(Char_t* nameReso,
463 AliGenParam* const genReso,
464 Double_t sigmaReso,
465 Double_t polReso)
466{
467// add resonances to the cocktail
468 Double_t phiMin = fPhiMin*180./TMath::Pi();
469 Double_t phiMax = fPhiMax*180./TMath::Pi();
470
471 // first step: generation in 4pi
472 genReso->SetPtRange(0.,100.);
473 genReso->SetYRange(-8.,8.);
474 genReso->SetPhiRange(0.,360.);
475 genReso->SetForceDecay(fDecayModeResonance);
476 if (!gMC) genReso->SetDecayer(fDecayer);
477 genReso->Init(); // generation in 4pi
478// Ratios with respect to the reaction cross-section in the
479// kinematics limit of the MUONCocktail
480 Double_t ratioReso = sigmaReso / fSigmaReaction * genReso->GetRelativeArea(fPtMin,fPtMax,fYMin,fYMax,phiMin,phiMax);
481 if (!fSigmaSilent) {
482 AliInfo(Form("%s prod. cross-section in pp %5.3g b",nameReso,sigmaReso));
483 AliInfo(Form("%s prod. probability per collision in acceptance %5.3g",nameReso,ratioReso));
484 }
485// second step: generation in selected kinematical range
486 genReso->SetPtRange(fPtMin, fPtMax);
487 genReso->SetYRange(fYMin, fYMax);
488 genReso->SetPhiRange(phiMin, phiMax);
489 genReso->Init(); // generation in selected kinematical range
490
491 TVirtualMCDecayer *decReso = 0;
492 if(TMath::Abs(polReso) > 1.e-30){
493 AliInfo(Form("******Setting polarized decayer for %s''",nameReso));
494 if(fPolFrame==0) {
495 decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
496 AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",polReso));
497 }
498 if(fPolFrame==1) {
499 decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
500 AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",polReso));
501 }
418e091f 502 if (decReso) {
503 decReso->SetForceDecay(kAll);
504 decReso->Init();
505 genReso->SetDecayer(decReso);
506 }
ba762ab4 507 }
508
509 AddGenerator(genReso,nameReso,ratioReso); // Adding Generator
510}
511
512//-------------------------------------------------------------------
8058ead1 513void AliGenMUONCocktailpp::Init()
514{
8058ead1 515 // Initialisation
aaa95f22 516 TIter next(fEntries);
517 AliGenCocktailEntry *entry;
518 if (fStack) {
519 while((entry = (AliGenCocktailEntry*)next())) {
520 entry->Generator()->SetStack(fStack);
521 }
522 }
103ac317 523}
524
525//_________________________________________________________________________
526void AliGenMUONCocktailpp::Generate()
527{
528// Generate event
529 TIter next(fEntries);
530 AliGenCocktailEntry *entry = 0;
531 AliGenCocktailEntry *preventry = 0;
532 AliGenerator* gen = 0;
533
01c89ce1 534 if (fHeader) delete fHeader;
535 fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
536
d94af0c1 537 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
103ac317 538
539// Generate the vertex position used by all generators
540 if(fVertexSmear == kPerEvent) Vertex();
541
542// Loop on primordialTrigger:
543// minimum muon multiplicity above a pt cut in a theta acceptance region
544
545 Bool_t primordialTrigger = kFALSE;
546 while(!primordialTrigger) {
547 //Reseting stack
33c3c91a 548 AliRunLoader * runloader = AliRunLoader::Instance();
103ac317 549 if (runloader)
550 if (runloader->Stack())
34da8494 551 runloader->Stack()->Clean();
103ac317 552 // Loop over generators and generate events
553 Int_t igen = 0;
554 Int_t npart = 0;
555 const char* genName = 0;
556 while((entry = (AliGenCocktailEntry*)next())) {
557 gen = entry->Generator();
558 genName = entry->GetName();
559 gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
21391258 560 gen->SetTime(fTime);
103ac317 561
562 npart = (strcmp(genName,"Pythia") == 0) ? 1 :
563 gRandom->Poisson(entry->Rate());
564
565 if (npart > 0) {
566 igen++;
567 if (igen == 1) entry->SetFirst(0);
568 else entry->SetFirst((partArray->GetEntriesFast())+1);
569
570 gen->SetNumberParticles(npart);
571 gen->Generate();
572 entry->SetLast(partArray->GetEntriesFast());
573 preventry = entry;
574 }
575 }
576 next.Reset();
577
ba762ab4 578
103ac317 579// Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
580// in the muon spectrometer acceptance
581 Int_t iPart;
582 fNGenerated++;
b44c3901 583 Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
aaa95f22 584 for(iPart=0; iPart<maxPart; iPart++){
103ac317 585
aaa95f22 586 TParticle *part = gAlice->GetMCApp()->Particle(iPart);
587 if ( TMath::Abs(part->GetPdgCode()) == 13 ){
588 if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
589 (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
590 (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
f81a4aca 591 (part->Pt()>fMuonPtCut) &&
592 (part->P()>fMuonPCut)) {
aaa95f22 593 numberOfMuons++;
103ac317 594 }
aaa95f22 595 }
596 }
01c89ce1 597 if (numberOfMuons >= fMuonMultiplicity) {
598 primordialTrigger = kTRUE;
599 fHeader->SetNProduced(maxPart);
600 }
601
103ac317 602 }
603 fNSucceded++;
604
01c89ce1 605 TArrayF eventVertex;
606 eventVertex.Set(3);
607 for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
608
609 fHeader->SetPrimaryVertex(eventVertex);
21391258 610 fHeader->SetInteractionTime(fTime);
01c89ce1 611
612 gAlice->SetGenEventHeader(fHeader);
613
aaa95f22 614// AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
103ac317 615 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
616}
617
618