]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - EVGEN/AliGenMUONCocktailpp.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Transition PWG3/base -> PWGDQ/base
[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktailpp.cxx
... / ...
CommitLineData
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// 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
22// The free parameeters are :
23// pp reaction cross-section
24// production cross-sections in pp collisions
25// July 07:added heavy quark production from AliGenCorrHF and heavy quark
26// production switched off in Pythia
27// Aug. 07: added trigger cut on total momentum
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//-----------------
32// 2009: added polarization (L. Bianchi)
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
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
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
57
58#include <TObjArray.h>
59#include <TParticle.h>
60#include <TF1.h>
61#include <TVirtualMC.h>
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"
71#include "AliDecayer.h"
72#include "AliLog.h"
73#include "AliGenCorrHF.h"
74#include "AliDecayerPolarized.h"
75
76ClassImp(AliGenMUONCocktailpp)
77
78//________________________________________________________________________
79AliGenMUONCocktailpp::AliGenMUONCocktailpp()
80 :AliGenCocktail(),
81 fDecayer(0),
82 fDecayModeResonance(kAll),
83 fDecayModePythia(kAll),
84 fMuonMultiplicity(0),
85 fMuonPtCut(0.),
86 fMuonPCut(0.),
87 fMuonThetaMinCut(0.),
88 fMuonThetaMaxCut(180.),
89 fMuonOriginCut(-999.),
90 fNSucceded(0),
91 fNGenerated(0),
92 fCentralityBin(0),
93
94 fJpsiPol(0),
95 fChic1Pol(0),
96 fChic2Pol(0),
97 fPsiPPol(0),
98 fUpsPol(0),
99 fUpsPPol(0),
100 fUpsPPPol(0),
101 fPolFrame(0),
102
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
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
137 fCMSEnergyTeVArray[0] = 7.00;
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;
146 fSigmaCCbarArray[0] = 6.91e-3;
147 fSigmaBBbarArray[0] = 0.232e-3;
148
149//x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers
150// scaled down according to ccbar and bbbar cross-sections
151 fCMSEnergyTeVArray[1] = 10.00;
152 fSigmaReactionArray[1] = 0.070;
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;
162
163//x-sections for pp @ 14 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
164// bottomonium from hep-ph/0311048 Tab.9, page 19 taken into account that
165// feed-down from chib is included
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;
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
219
220}
221
222//_________________________________________________________________________
223AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
224{
225// Destructor
226
227}
228
229//_________________________________________________________________________
230void AliGenMUONCocktailpp::SetCMSEnergy(CMSEnergyCode cmsEnergy)
231{
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];
244}
245
246//_________________________________________________________________________
247void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
248 Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
249// setter for resonances polarization
250 if (strcmp(PolFrame,"kColSop")==0){
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;
257 } else if (strcmp(PolFrame,"kHelicity")==0){
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
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//_________________________________________________________________________
277void AliGenMUONCocktailpp::CreateCocktail()
278{
279// create and add resonances and open HF to the coctail
280 Int_t cmsEnergy = Int_t(fCMSEnergyTeV);
281
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
292// Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and
293// corrected from feed down of higher resonances
294
295 Double_t sigmajpsi = fSigmaJPsi;
296 Double_t sigmachic1 = fSigmaChic1;
297 Double_t sigmachic2 = fSigmaChic2;
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;
304
305// Cross sections corrected with the BR in mu+mu-
306// (only in case of use of AliDecayerPolarized)
307
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;}
315
316 AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
317
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
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");
342
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);
354 } else {
355 AliError("Initialisation failed, wrong cmsEnergy");
356 return;
357 }
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);
407
408 AddReso2Generator(nameUps,genupsilon,bhard[fCentralityBin]*sigmaupsilon,fUpsPol);
409 AddReso2Generator(nameUpsP,genupsilonP,bhard[fCentralityBin]*sigmaupsilonP,fUpsPPol);
410 AddReso2Generator(nameUpsPP,genupsilonPP,bhard[fCentralityBin]*sigmaupsilonPP,fUpsPPPol);
411
412//------------------------------------------------------------------
413// Generator of charm
414 AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
415 gencharm->SetMomentumRange(0,9999);
416 gencharm->SetForceDecay(kAll);
417 Double_t ratioccbar = chard[fCentralityBin]*sigmaccbar/fSigmaReaction;
418 if (!gMC) gencharm->SetDecayer(fDecayer);
419 gencharm->Init();
420 if (!fSigmaSilent) {
421 AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
422 AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
423 }
424 AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
425//------------------------------------------------------------------
426// Generator of beauty
427 AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
428 genbeauty->SetMomentumRange(0,9999);
429 genbeauty->SetForceDecay(kAll);
430 Double_t ratiobbbar = bhard[fCentralityBin]*sigmabbbar/fSigmaReaction;
431 if (!gMC) genbeauty->SetDecayer(fDecayer);
432 genbeauty->Init();
433 if (!fSigmaSilent) {
434 AliInfo(Form("b-bbar prod. cross-section in pp %5.3g b",sigmabbbar));
435 AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
436 }
437 AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
438
439//-------------------------------------------------------------------
440// Pythia generator
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);
458
459}
460
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 }
502 if (decReso) {
503 decReso->SetForceDecay(kAll);
504 decReso->Init();
505 genReso->SetDecayer(decReso);
506 }
507 }
508
509 AddGenerator(genReso,nameReso,ratioReso); // Adding Generator
510}
511
512//-------------------------------------------------------------------
513void AliGenMUONCocktailpp::Init()
514{
515 // Initialisation
516 TIter next(fEntries);
517 AliGenCocktailEntry *entry;
518 if (fStack) {
519 while((entry = (AliGenCocktailEntry*)next())) {
520 entry->Generator()->SetStack(fStack);
521 }
522 }
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
534 if (fHeader) delete fHeader;
535 fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
536
537 const TObjArray *partArray = gAlice->GetMCApp()->Particles();
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
548 AliRunLoader * runloader = AliRunLoader::Instance();
549 if (runloader)
550 if (runloader->Stack())
551 runloader->Stack()->Clean();
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));
560 gen->SetTime(fTime);
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
578
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++;
583 Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
584 for(iPart=0; iPart<maxPart; iPart++){
585
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) &&
591 (part->Pt()>fMuonPtCut) &&
592 (part->P()>fMuonPCut)) {
593 numberOfMuons++;
594 }
595 }
596 }
597 if (numberOfMuons >= fMuonMultiplicity) {
598 primordialTrigger = kTRUE;
599 fHeader->SetNProduced(maxPart);
600 }
601
602 }
603 fNSucceded++;
604
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);
610 fHeader->SetInteractionTime(fTime);
611
612 gAlice->SetGenEventHeader(fHeader);
613
614// AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
615 AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
616}
617
618