Line data Source code
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 : // 06/2012:
59 : // - added the cocktail for p-Pb & Pb-p @ 2.76, 4.4 & 5.03 TeV with 4 centrality
60 : // bins, using the EPS09-LO shadowing computed for 5.03 TeV. Energies are set by
61 : // AliGenMUONCocktailpp::SetCMSEnergy(int CMSEnergyCode), CMSEnergy codes are
62 : // defined in AliGenMUONCocktailpp.h.
63 : // - added functions to scale x-section of JPsi, Charmonia, Bottomonia, CCbar & BBbar
64 : // in Config.C to manage the statistics. Example of usage (in a cocktail with Hijing):
65 : /*
66 : AliGenCocktail *cocktail = new AliGenCocktail();
67 : cocktail->AddGenerator(hijing,"hijing",1);
68 : TFormula *form[nb]; // nb - number of centrality bins with impact params bBins[i]
69 : AliGenMUONCocktailpp *gen[nb];
70 : for (Int_t i=0; i<nb; i++) {
71 : form[i] = new TFormula(Form("f%d",i),"[0]+(x-[1])/([2]-[1])");
72 : form[i]->SetParameters(i+1, bBins[i], bBins[i+1]);
73 : gen[i] = MuonCocktail();
74 : gen[i]->SetCentralityBin(i+1);
75 : gen[i]->ScaleJPsi(100.);
76 : gen[i]->CreateCocktail();
77 : cocktail->AddGenerator(gen[i], Form("g%d",i), 101+i, form[i]);
78 : }
79 : AliGenMUONCocktailpp* MuonCocktail() {
80 : AliGenMUONCocktailpp *mc = new AliGenMUONCocktailpp();
81 : ....................................
82 : return mc;
83 : }
84 : */
85 : // - a bug fixed in the function Generate()
86 : // S. Grigoryan
87 :
88 : #include <TObjArray.h>
89 : #include <TParticle.h>
90 : #include <TF1.h>
91 : #include <TVirtualMC.h>
92 : #include "AliGenCocktailEventHeader.h"
93 :
94 : #include "AliGenCocktailEntry.h"
95 : #include "AliGenMUONCocktailpp.h"
96 : #include "AliGenMUONlib.h"
97 : #include "AliGenParam.h"
98 : #include "AliMC.h"
99 : #include "AliRun.h"
100 : #include "AliStack.h"
101 : #include "AliDecayer.h"
102 : #include "AliLog.h"
103 : #include "AliGenCorrHF.h"
104 : #include "AliDecayerPolarized.h"
105 :
106 6 : ClassImp(AliGenMUONCocktailpp)
107 :
108 : //________________________________________________________________________
109 : AliGenMUONCocktailpp::AliGenMUONCocktailpp()
110 0 : :AliGenCocktail(),
111 0 : fDecayer(0),
112 0 : fDecayModeResonance(kAll),
113 0 : fDecayModePythia(kAll),
114 0 : fMuonMultiplicity(0),
115 0 : fMuonPtCut(0.),
116 0 : fMuonPCut(0.),
117 0 : fMuonThetaMinCut(0.),
118 0 : fMuonThetaMaxCut(180.),
119 0 : fMuonOriginCut(-999.),
120 0 : fNSucceded(0),
121 0 : fNGenerated(0),
122 0 : fCentralityBin(0),
123 0 : fScaleJPsi(1),
124 0 : fScaleCharmonia(1),
125 0 : fScaleBottomonia(1),
126 0 : fScaleCCbar(1),
127 0 : fScaleBBbar(1),
128 :
129 0 : fJpsiPol(0),
130 0 : fChic1Pol(0),
131 0 : fChic2Pol(0),
132 0 : fPsiPPol(0),
133 0 : fUpsPol(0),
134 0 : fUpsPPol(0),
135 0 : fUpsPPPol(0),
136 0 : fPolFrame(0),
137 :
138 0 : fCMSEnergyTeV(0),
139 0 : fCMSEnergyTeVArray(),
140 0 : fSigmaReaction(0),
141 0 : fSigmaReactionArray(),
142 0 : fSigmaJPsi(0),
143 0 : fSigmaJPsiArray(),
144 0 : fSigmaChic1(0),
145 0 : fSigmaChic1Array(),
146 0 : fSigmaChic2(0),
147 0 : fSigmaChic2Array(),
148 0 : fSigmaPsiP(0),
149 0 : fSigmaPsiPArray(),
150 0 : fSigmaUpsilon(0),
151 0 : fSigmaUpsilonArray(),
152 0 : fSigmaUpsilonP(0),
153 0 : fSigmaUpsilonPArray(),
154 0 : fSigmaUpsilonPP(0),
155 0 : fSigmaUpsilonPPArray(),
156 0 : fSigmaCCbar(0),
157 0 : fSigmaCCbarArray(),
158 0 : fSigmaBBbar(0),
159 0 : fSigmaBBbarArray(),
160 :
161 0 : fSigmaSilent(kFALSE)
162 0 : {
163 : // Constructor
164 :
165 : // x-sections for pp @ 7 TeV:
166 : // -charmonia: 4pi integral of fit function for inclusive J/psi dsigma/dy LHC data
167 : // gives 60 mub; so sigma_prompt = 54 mub, while Ref = R.Vogt_arXiv:1003.3497 (Table 2)
168 : // gives 35 mub. Below we use sigma_direct from the Ref scaled by the factor 54/35.
169 : // -bottomonia: 4pi integral of fit function for inclusive Upsilon1S dsigma/dy LHC data
170 : // gives 0.56 mub, sigmas for 2S & 3S obtained using LHCb data for ratios 2S/1S & 3S/1S
171 : // -ccbar & bbbar: NLO pQCD computations - http://www-alice.gsi.de/ana/MNR/results.html
172 0 : fCMSEnergyTeVArray[0] = 7.00;
173 0 : fSigmaReactionArray[0] = 0.070;
174 0 : fSigmaJPsiArray[0] = 33.6e-6;
175 0 : fSigmaChic1Array[0] = 32.6e-6;
176 0 : fSigmaChic2Array[0] = 53.8e-6;
177 0 : fSigmaPsiPArray[0] = 7.6e-6;
178 0 : fSigmaUpsilonArray[0] = 0.56e-6;
179 0 : fSigmaUpsilonPArray[0] = 0.18e-6;
180 0 : fSigmaUpsilonPPArray[0] = 0.08e-6;
181 0 : fSigmaCCbarArray[0] = 6.91e-3;
182 0 : fSigmaBBbarArray[0] = 0.232e-3;
183 :
184 : //x-sections for pp @ 10 TeV: charmonia and bottomonia from 14 TeV numbers
185 : // scaled down according to ccbar and bbbar cross-sections
186 0 : fCMSEnergyTeVArray[1] = 10.00;
187 0 : fSigmaReactionArray[1] = 0.070;
188 0 : fSigmaJPsiArray[1] = 26.06e-6;
189 0 : fSigmaChic1Array[1] = 25.18e-6;
190 0 : fSigmaChic2Array[1] = 41.58e-6;
191 0 : fSigmaPsiPArray[1] = 5.88e-6;
192 0 : fSigmaUpsilonArray[1] = 0.658e-6;
193 0 : fSigmaUpsilonPArray[1] = 0.218e-6;
194 0 : fSigmaUpsilonPPArray[1] = 0.122e-6;
195 0 : fSigmaCCbarArray[1] = 8.9e-3;
196 0 : fSigmaBBbarArray[1] = 0.33e-3;
197 :
198 : //x-sections for pp @ 14 TeV: charmonia from hep-ph/0311048 Tab.9, page 19,
199 : // bottomonium from hep-ph/0311048 Tab.9, page 19 taken into account that
200 : // feed-down from chib is included
201 0 : fCMSEnergyTeVArray[2] = 14.00;
202 0 : fSigmaReactionArray[2] = 0.070;
203 0 : fSigmaJPsiArray[2] = 32.9e-6;
204 0 : fSigmaChic1Array[2] = 31.8e-6;
205 0 : fSigmaChic2Array[2] = 52.5e-6;
206 0 : fSigmaPsiPArray[2] = 7.43e-6;
207 0 : fSigmaUpsilonArray[2] = 0.989e-6;
208 0 : fSigmaUpsilonPArray[2] = 0.502e-6;
209 0 : fSigmaUpsilonPPArray[2] = 0.228e-6;
210 0 : fSigmaCCbarArray[2] = 11.2e-3;
211 0 : fSigmaBBbarArray[2] = 0.445e-3;
212 :
213 : // x-sections for Min. Bias p-Pb & Pb-p @ 2.76 TeV: charmonia and bottomonia
214 : // from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
215 : // and with Glauber scaling
216 0 : fCMSEnergyTeVArray[3] = 2.00; // for 2.76 TeV
217 0 : fSigmaReactionArray[3] = 2.10;
218 0 : fSigmaJPsiArray[3] = 3.54e-3; // 208*0.507*33.6e-6
219 0 : fSigmaChic1Array[3] = 3.44e-3;
220 0 : fSigmaChic2Array[3] = 5.66e-3;
221 0 : fSigmaPsiPArray[3] = 0.80e-3;
222 0 : fSigmaUpsilonArray[3] = 0.045e-3; // 208*0.384*0.56e-6
223 0 : fSigmaUpsilonPArray[3] = 0.014e-3;
224 0 : fSigmaUpsilonPPArray[3] = 0.006e-3;
225 0 : fSigmaCCbarArray[3] = 0.73; // 208*3.50e-3
226 0 : fSigmaBBbarArray[3] = 0.019; // 208*0.089e-3
227 :
228 0 : fCMSEnergyTeVArray[4] = -fCMSEnergyTeVArray[3];
229 0 : fSigmaReactionArray[4] = fSigmaReactionArray[3];
230 0 : fSigmaJPsiArray[4] = fSigmaJPsiArray[3];
231 0 : fSigmaChic1Array[4] = fSigmaChic1Array[3];
232 0 : fSigmaChic2Array[4] = fSigmaChic2Array[3];
233 0 : fSigmaPsiPArray[4] = fSigmaPsiPArray[3];
234 0 : fSigmaUpsilonArray[4] = fSigmaUpsilonArray[3];
235 0 : fSigmaUpsilonPArray[4] = fSigmaUpsilonPArray[3];
236 0 : fSigmaUpsilonPPArray[4] = fSigmaUpsilonPPArray[3];
237 0 : fSigmaCCbarArray[4] = fSigmaCCbarArray[3];
238 0 : fSigmaBBbarArray[4] = fSigmaBBbarArray[3];
239 :
240 : // x-sections for Min. Bias p-Pb & Pb-p @ 4.4 TeV: charmonia and bottomonia
241 : // from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
242 : // and with Glauber scaling
243 0 : fCMSEnergyTeVArray[5] = 4.00; // for 4.4 TeV
244 0 : fSigmaReactionArray[5] = 2.10;
245 0 : fSigmaJPsiArray[5] = 5.00e-3; // 208*0.715*33.6e-6
246 0 : fSigmaChic1Array[5] = 4.86e-3;
247 0 : fSigmaChic2Array[5] = 7.99e-3;
248 0 : fSigmaPsiPArray[5] = 1.12e-3;
249 0 : fSigmaUpsilonArray[5] = 0.074e-3; // 208*0.629*0.56e-6
250 0 : fSigmaUpsilonPArray[5] = 0.023e-3;
251 0 : fSigmaUpsilonPPArray[5] = 0.010e-3;
252 0 : fSigmaCCbarArray[5] = 1.03; // 208*4.94e-3
253 0 : fSigmaBBbarArray[5] = 0.030; // 208*0.146e-3
254 :
255 0 : fCMSEnergyTeVArray[6] = -fCMSEnergyTeVArray[5];
256 0 : fSigmaReactionArray[6] = fSigmaReactionArray[5];
257 0 : fSigmaJPsiArray[6] = fSigmaJPsiArray[5];
258 0 : fSigmaChic1Array[6] = fSigmaChic1Array[5];
259 0 : fSigmaChic2Array[6] = fSigmaChic2Array[5];
260 0 : fSigmaPsiPArray[6] = fSigmaPsiPArray[5];
261 0 : fSigmaUpsilonArray[6] = fSigmaUpsilonArray[5];
262 0 : fSigmaUpsilonPArray[6] = fSigmaUpsilonPArray[5];
263 0 : fSigmaUpsilonPPArray[6] = fSigmaUpsilonPPArray[5];
264 0 : fSigmaCCbarArray[6] = fSigmaCCbarArray[5];
265 0 : fSigmaBBbarArray[6] = fSigmaBBbarArray[5];
266 :
267 : // x-sections for Min. Bias p-Pb & Pb-p @ 5.03 TeV: charmonia and bottomonia
268 : // from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
269 : // and with Glauber scaling
270 0 : fCMSEnergyTeVArray[7] = 5.00; // for 5.03 TeV
271 0 : fSigmaReactionArray[7] = 2.10;
272 0 : fSigmaJPsiArray[7] = 5.50e-3; // 208*0.787*33.6e-6
273 0 : fSigmaChic1Array[7] = 5.35e-3;
274 0 : fSigmaChic2Array[7] = 8.79e-3;
275 0 : fSigmaPsiPArray[7] = 1.23e-3;
276 0 : fSigmaUpsilonArray[7] = 0.083e-3; // 208*0.716*0.56e-6
277 0 : fSigmaUpsilonPArray[7] = 0.026e-3;
278 0 : fSigmaUpsilonPPArray[7] = 0.011e-3;
279 0 : fSigmaCCbarArray[7] = 1.13; // 208*5.44e-3
280 0 : fSigmaBBbarArray[7] = 0.035; // 208*0.166e-3
281 :
282 0 : fCMSEnergyTeVArray[8] = -fCMSEnergyTeVArray[7];
283 0 : fSigmaReactionArray[8] = fSigmaReactionArray[7];
284 0 : fSigmaJPsiArray[8] = fSigmaJPsiArray[7];
285 0 : fSigmaChic1Array[8] = fSigmaChic1Array[7];
286 0 : fSigmaChic2Array[8] = fSigmaChic2Array[7];
287 0 : fSigmaPsiPArray[8] = fSigmaPsiPArray[7];
288 0 : fSigmaUpsilonArray[8] = fSigmaUpsilonArray[7];
289 0 : fSigmaUpsilonPArray[8] = fSigmaUpsilonPArray[7];
290 0 : fSigmaUpsilonPPArray[8] = fSigmaUpsilonPPArray[7];
291 0 : fSigmaCCbarArray[8] = fSigmaCCbarArray[7];
292 0 : fSigmaBBbarArray[8] = fSigmaBBbarArray[7];
293 :
294 : // x-sections for Min. Bias p-Pb & Pb-p @ 8.8 TeV: charmonia and bottomonia
295 : // from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
296 : // and with Glauber scaling
297 0 : fCMSEnergyTeVArray[9] = 9.00; // for 8.8 TeV
298 0 : fSigmaReactionArray[9] = 2.10;
299 0 : fSigmaJPsiArray[9] = 8.19e-3; // 208*1.172*33.6e-6
300 0 : fSigmaChic1Array[9] = 7.95e-3;
301 0 : fSigmaChic2Array[9] = 13.1e-3;
302 0 : fSigmaPsiPArray[9] = 1.85e-3;
303 0 : fSigmaUpsilonArray[9] = 0.146e-3; // 208*1.25*0.56e-6
304 0 : fSigmaUpsilonPArray[9] = 0.047e-3;
305 0 : fSigmaUpsilonPPArray[9] = 0.021e-3;
306 0 : fSigmaCCbarArray[9] = 1.68; // 208*8.1e-3
307 0 : fSigmaBBbarArray[9] = 0.061; // 208*0.29e-3
308 :
309 0 : fCMSEnergyTeVArray[10] = -fCMSEnergyTeVArray[9];
310 0 : fSigmaReactionArray[10] = fSigmaReactionArray[9];
311 0 : fSigmaJPsiArray[10] = fSigmaJPsiArray[9];
312 0 : fSigmaChic1Array[10] = fSigmaChic1Array[9];
313 0 : fSigmaChic2Array[10] = fSigmaChic2Array[9];
314 0 : fSigmaPsiPArray[10] = fSigmaPsiPArray[9];
315 0 : fSigmaUpsilonArray[10] = fSigmaUpsilonArray[9];
316 0 : fSigmaUpsilonPArray[10] = fSigmaUpsilonPArray[9];
317 0 : fSigmaUpsilonPPArray[10] = fSigmaUpsilonPPArray[9];
318 0 : fSigmaCCbarArray[10] = fSigmaCCbarArray[9];
319 0 : fSigmaBBbarArray[10] = fSigmaBBbarArray[9];
320 :
321 : // x-sections for Min. Bias Pb-Pb @ 2.76 TeV: charmonia and bottomonia
322 : // from 7 TeV numbers scaled according to pQCD ccbar and bbbar x-sections
323 : // and with Glauber scaling
324 0 : fCMSEnergyTeVArray[11] = 3.00; // for 2.76 TeV
325 0 : fSigmaReactionArray[11] = 7.65;
326 0 : fSigmaJPsiArray[11] = 0.737; // 208*208*0.507*33.6e-6
327 0 : fSigmaChic1Array[11] = 0.715;
328 0 : fSigmaChic2Array[11] = 1.179;
329 0 : fSigmaPsiPArray[11] = 0.166;
330 0 : fSigmaUpsilonArray[11] = 0.0093; // 208*208*0.384*0.56e-6
331 0 : fSigmaUpsilonPArray[11] = 0.0030;
332 0 : fSigmaUpsilonPPArray[11] = 0.0013;
333 0 : fSigmaCCbarArray[11] = 151.; // 208*208*3.50e-3
334 0 : fSigmaBBbarArray[11] = 3.8; // 208*208*0.089e-3
335 :
336 0 : }
337 :
338 : //_________________________________________________________________________
339 0 : AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
340 0 : {
341 : // Destructor
342 :
343 0 : }
344 :
345 : //_________________________________________________________________________
346 : void AliGenMUONCocktailpp::SetCMSEnergy(CMSEnergyCode cmsEnergy)
347 : {
348 : // setter for CMSEnergy and corresponding cross-sections
349 0 : fCMSEnergyTeV = fCMSEnergyTeVArray[cmsEnergy];
350 0 : fSigmaReaction = fSigmaReactionArray[cmsEnergy];
351 0 : fSigmaJPsi = fSigmaJPsiArray[cmsEnergy];
352 0 : fSigmaChic1 = fSigmaChic1Array[cmsEnergy];
353 0 : fSigmaChic2 = fSigmaChic2Array[cmsEnergy];
354 0 : fSigmaPsiP = fSigmaPsiPArray[cmsEnergy];
355 0 : fSigmaUpsilon = fSigmaUpsilonArray[cmsEnergy];
356 0 : fSigmaUpsilonP = fSigmaUpsilonPArray[cmsEnergy];
357 0 : fSigmaUpsilonPP = fSigmaUpsilonPPArray[cmsEnergy];
358 0 : fSigmaCCbar = fSigmaCCbarArray[cmsEnergy];
359 0 : fSigmaBBbar = fSigmaBBbarArray[cmsEnergy];
360 0 : }
361 :
362 : //_________________________________________________________________________
363 : void AliGenMUONCocktailpp::SetResPolarization(Double_t JpsiPol, Double_t PsiPPol, Double_t UpsPol,
364 : Double_t UpsPPol, Double_t UpsPPPol, char *PolFrame){
365 : // setter for resonances polarization
366 0 : if (strcmp(PolFrame,"kColSop")==0){
367 0 : fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0;
368 0 : fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0;
369 0 : fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0;
370 0 : fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0;
371 0 : fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
372 0 : fPolFrame = 0;
373 0 : } else if (strcmp(PolFrame,"kHelicity")==0){
374 0 : fJpsiPol = (JpsiPol>=-1 && JpsiPol<=1) ? JpsiPol : 0;
375 0 : fPsiPPol = (PsiPPol>=-1 && PsiPPol<=1) ? PsiPPol : 0;
376 0 : fUpsPol = (UpsPol>=-1 && UpsPol<=1) ? UpsPol : 0;
377 0 : fUpsPPol = (UpsPPol>=-1 && UpsPPol<=1) ? UpsPPol : 0;
378 0 : fUpsPPPol = (UpsPPPol>=-1 && UpsPPPol<=1) ? UpsPPPol : 0;
379 0 : fPolFrame = 1;
380 :
381 0 : } else {
382 0 : AliInfo(Form("The polarization frame is not valid"));
383 0 : AliInfo(Form("No polarization will be set"));
384 0 : fJpsiPol=0.;
385 0 : fPsiPPol=0.;
386 0 : fUpsPol=0.;
387 0 : fUpsPPol=0.;
388 0 : fUpsPPPol=0.;
389 : }
390 0 : }
391 :
392 : //_________________________________________________________________________
393 : void AliGenMUONCocktailpp::CreateCocktail()
394 : {
395 : // create and add resonances and open HF to the coctail
396 0 : Int_t cmsEnergy = Int_t(fCMSEnergyTeV);
397 :
398 : // For temporary use of p-Pb & Pb-p shadowing at 5.03 TeV for lower energies
399 0 : if (cmsEnergy == 2 || cmsEnergy == 4) cmsEnergy = 5;
400 0 : if (cmsEnergy == -2 || cmsEnergy == -4) cmsEnergy = -5;
401 :
402 : // These limits are only used for renormalization of quarkonia cross section
403 : // Pythia events are generated in 4pi
404 0 : Double_t ptMin = fPtMin;
405 0 : Double_t ptMax = fPtMax;
406 0 : Double_t yMin = fYMin;;
407 0 : Double_t yMax = fYMax;;
408 0 : Double_t phiMin = fPhiMin*180./TMath::Pi();
409 0 : Double_t phiMax = fPhiMax*180./TMath::Pi();
410 0 : 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));
411 :
412 : // Cross sections in barns
413 :
414 0 : Double_t sigmajpsi = fSigmaJPsi;
415 0 : Double_t sigmachic1 = fSigmaChic1;
416 0 : Double_t sigmachic2 = fSigmaChic2;
417 0 : Double_t sigmapsiP = fSigmaPsiP;
418 0 : Double_t sigmaupsilon = fSigmaUpsilon;
419 0 : Double_t sigmaupsilonP = fSigmaUpsilonP;
420 0 : Double_t sigmaupsilonPP = fSigmaUpsilonPP;
421 0 : Double_t sigmaccbar = fSigmaCCbar;
422 0 : Double_t sigmabbbar = fSigmaBBbar;
423 :
424 : // Cross sections corrected with the BR in mu+mu-
425 : // (only in case of use of AliDecayerPolarized)
426 :
427 0 : if(TMath::Abs(fJpsiPol) > 1.e-30) {sigmajpsi = fSigmaJPsi*0.0593;}
428 0 : if(TMath::Abs(fChic1Pol) > 1.e-30) {sigmachic1 = fSigmaChic1*0.;} // tb consistent
429 0 : if(TMath::Abs(fChic2Pol) > 1.e-30) {sigmachic2 = fSigmaChic2*0.;} // tb consistent
430 0 : if(TMath::Abs(fPsiPPol) > 1.e-30) {sigmapsiP = fSigmaPsiP*0.0075;}
431 0 : if(TMath::Abs(fUpsPol) > 1.e-30) {sigmaupsilon = fSigmaUpsilon*0.0248;}
432 0 : if(TMath::Abs(fUpsPPol) > 1.e-30) {sigmaupsilonP = fSigmaUpsilonP*0.0193;}
433 0 : if(TMath::Abs(fUpsPPPol) > 1.e-30) {sigmaupsilonPP = fSigmaUpsilonPP*0.0218;}
434 :
435 : // Cross sections scaled to manage the statistics
436 :
437 0 : sigmajpsi *= fScaleJPsi*fScaleCharmonia;
438 0 : sigmachic1 *= fScaleCharmonia;
439 0 : sigmachic2 *= fScaleCharmonia;
440 0 : sigmapsiP *= fScaleCharmonia;
441 0 : sigmaupsilon *= fScaleBottomonia;
442 0 : sigmaupsilonP *= fScaleBottomonia;
443 0 : sigmaupsilonPP *= fScaleBottomonia;
444 0 : sigmaccbar *= fScaleCCbar;
445 0 : sigmabbbar *= fScaleBBbar;
446 :
447 0 : AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
448 :
449 : // Create and add resonances to the generator
450 : AliGenParam * genjpsi=0;
451 : AliGenParam * genchic1=0;
452 : AliGenParam * genchic2=0;
453 : AliGenParam * genpsiP=0;
454 : AliGenParam * genupsilon=0;
455 : AliGenParam * genupsilonP=0;
456 : AliGenParam * genupsilonPP=0;
457 :
458 0 : Char_t nameJpsi[10];
459 0 : Char_t nameChic1[10];
460 0 : Char_t nameChic2[10];
461 0 : Char_t namePsiP[10];
462 0 : Char_t nameUps[10];
463 0 : Char_t nameUpsP[10];
464 0 : Char_t nameUpsPP[10];
465 :
466 0 : snprintf(nameJpsi,10, "Jpsi");
467 0 : snprintf(nameChic1,10, "Chic1");
468 0 : snprintf(nameChic2,10, "Chic2");
469 0 : snprintf(namePsiP,10, "PsiP");
470 0 : snprintf(nameUps,10, "Ups");
471 0 : snprintf(nameUpsP,10, "UpsP");
472 0 : snprintf(nameUpsPP,10, "UpsPP");
473 :
474 0 : Char_t tname[40] = "";
475 0 : if(cmsEnergy == 10) {snprintf(tname, 40, "CDF pp 10");
476 0 : } else if (cmsEnergy == 14){snprintf(tname, 40, "CDF pp");
477 0 : } else if (cmsEnergy == 7) {snprintf(tname, 40, "pp 7");
478 : // } else if (cmsEnergy == 2) {snprintf(tname, 40, "pp 2.76");
479 0 : } else if (cmsEnergy == 5) {snprintf(tname, 40, "pPb 5.03");
480 0 : if (fCentralityBin > 0) snprintf(tname, 40, "pPb 5.03c%d",fCentralityBin);
481 0 : } else if (cmsEnergy == -5){snprintf(tname, 40, "Pbp 5.03");
482 0 : if (fCentralityBin > 0) snprintf(tname, 40, "Pbp 5.03c%d",fCentralityBin);
483 0 : } else if (cmsEnergy == 9) {snprintf(tname, 40, "pPb 8.8");
484 0 : if (fCentralityBin > 0) snprintf(tname, 40, "pPb 8.8c%d",fCentralityBin);
485 0 : } else if (cmsEnergy == -9){snprintf(tname, 40, "Pbp 8.8");
486 0 : if (fCentralityBin > 0) snprintf(tname, 40, "Pbp 8.8c%d",fCentralityBin);
487 0 : } else if (cmsEnergy == 3) {snprintf(tname, 40, "PbPb 2.76");
488 0 : if (fCentralityBin > 0) snprintf(tname, 40, "PbPb 2.76c%d",fCentralityBin);
489 : } else {
490 0 : AliError("Initialisation failed, wrong cmsEnergy");
491 0 : return;
492 : }
493 0 : genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, tname, "Jpsi");
494 0 : genchic1 = new AliGenParam(1, AliGenMUONlib::kChic1, tname, "Chic1");
495 0 : genchic2 = new AliGenParam(1, AliGenMUONlib::kChic2, tname, "Chic2");
496 0 : genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, tname, "PsiP");
497 0 : genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, tname, "Upsilon");
498 0 : genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, tname, "UpsilonP");
499 0 : genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, tname, "UpsilonPP");
500 :
501 : // Hard process yield per pA or AA collision for i-th centrality bin is R*r[i]*shad[i]
502 : // where R is the ratio of hard and geometrical x-sections, r[i] is the ratio of these
503 : // x-section fractions for given centrality and shad[i] is the shadowing factor (in 4pi).
504 : // The latter is assumed to be the same for HF-hadrons & quarkonia of the same flavour.
505 : Int_t i = 0;
506 0 : Double_t chard[20] = {0}; // charm & beauty shadowing factors are different
507 0 : Double_t bhard[20] = {0};
508 0 : chard[0] = 1; // 1st element for pp and min. bias (MB) collisions
509 0 : bhard[0] = 1;
510 :
511 : // 4 centrality bins for p-Pb & Pb-p at 5.03 TeV: 0-20-40-60-100 %
512 0 : if (cmsEnergy == 5 || cmsEnergy == -5) {
513 : const Int_t n5 = 5; // 1st element for MB collisions
514 0 : Double_t r5[n5] = {1, 1.936, 1.473, 0.914, 0.333}; // ratio of hard-over-geo fractions
515 0 : Double_t cshad5[n5] = {0.806, 0.742, 0.796, 0.870, 0.955};// EPS09-LO shadowing factors
516 0 : Double_t bshad5[n5] = {0.917, 0.889, 0.913, 0.944, 0.981};
517 0 : for(i=0; i<n5; i++) {
518 0 : chard[i] = cshad5[i]*r5[i];
519 0 : bhard[i] = bshad5[i]*r5[i];
520 : }
521 0 : }
522 :
523 : // 4 centrality bins for p-Pb & Pb-p at 8.8 TeV: 0-20-40-60-100 %
524 0 : if (cmsEnergy == 9 || cmsEnergy == -9) {
525 : const Int_t n9 = 5; // 1st element for MB collisions
526 0 : Double_t r9[n9] = {1, 1.936, 1.473, 0.914, 0.333}; // ratio of hard-over-geo fractions
527 0 : Double_t cshad9[n9] = {0.785, 0.715, 0.775, 0.856, 0.951};// EKS98 shadowing factors
528 0 : Double_t bshad9[n9] = {0.889, 0.853, 0.884, 0.926, 0.975};
529 0 : for(i=0; i<n9; i++) {
530 0 : chard[i] = cshad9[i]*r9[i];
531 0 : bhard[i] = bshad9[i]*r9[i];
532 : }
533 0 : }
534 :
535 : // 11 centrality bins for Pb-Pb at 2.76 TeV: 0-5-10-20-30-40-50-60-70-80-90-100 %
536 0 : if (cmsEnergy == 3) {
537 : const Int_t n3 = 12; // 1st element for MB collisions
538 0 : Double_t r3[n3] = {1, 4.661, 3.647, 2.551, 1.544, 0.887, 0.474,
539 : 0.235, 0.106, 0.044, 0.017, 0.007}; // ratio of hard-over-geo fractions
540 0 : Double_t cshad3[n3] = {0.662, 0.622, 0.631, 0.650, 0.681, 0.718,
541 : 0.760, 0.805, 0.849, 0.888, 0.918, 0.944};// EKS98 shadowing factors
542 0 : Double_t bshad3[n3] = {0.874, 0.856, 0.861, 0.869, 0.883, 0.898,
543 : 0.915, 0.932, 0.948, 0.962, 0.972, 0.981};
544 0 : for(i=0; i<n3; i++) {
545 0 : chard[i] = cshad3[i]*r3[i];
546 0 : bhard[i] = bshad3[i]*r3[i];
547 : }
548 0 : }
549 :
550 0 : AddReso2Generator(nameJpsi,genjpsi,chard[fCentralityBin]*sigmajpsi,fJpsiPol);
551 0 : AddReso2Generator(nameChic1,genchic1,chard[fCentralityBin]*sigmachic1,fChic1Pol);
552 0 : AddReso2Generator(nameChic2,genchic2,chard[fCentralityBin]*sigmachic2,fChic2Pol);
553 0 : AddReso2Generator(namePsiP,genpsiP,chard[fCentralityBin]*sigmapsiP,fPsiPPol);
554 :
555 0 : AddReso2Generator(nameUps,genupsilon,bhard[fCentralityBin]*sigmaupsilon,fUpsPol);
556 0 : AddReso2Generator(nameUpsP,genupsilonP,bhard[fCentralityBin]*sigmaupsilonP,fUpsPPol);
557 0 : AddReso2Generator(nameUpsPP,genupsilonPP,bhard[fCentralityBin]*sigmaupsilonPP,fUpsPPPol);
558 :
559 : //------------------------------------------------------------------
560 : // Generator of charm
561 0 : AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4, cmsEnergy);
562 0 : gencharm->SetMomentumRange(0,9999);
563 0 : gencharm->SetForceDecay(kAll);
564 0 : Double_t ratioccbar = chard[fCentralityBin]*sigmaccbar/fSigmaReaction;
565 0 : if (!TVirtualMC::GetMC()) gencharm->SetDecayer(fDecayer);
566 0 : gencharm->Init();
567 0 : if (!fSigmaSilent) {
568 0 : AliInfo(Form("c-cbar prod. cross-section in pp %5.3g b",sigmaccbar));
569 0 : AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
570 0 : }
571 0 : AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
572 : //------------------------------------------------------------------
573 : // Generator of beauty
574 0 : AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5, cmsEnergy);
575 0 : genbeauty->SetMomentumRange(0,9999);
576 0 : genbeauty->SetForceDecay(kAll);
577 0 : Double_t ratiobbbar = bhard[fCentralityBin]*sigmabbbar/fSigmaReaction;
578 0 : if (!TVirtualMC::GetMC()) genbeauty->SetDecayer(fDecayer);
579 0 : genbeauty->Init();
580 0 : if (!fSigmaSilent) {
581 0 : AliInfo(Form("b-bbar prod. cross-section in pp %5.3g b",sigmabbbar));
582 0 : AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
583 0 : }
584 0 : AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
585 :
586 : //-------------------------------------------------------------------
587 : // Pythia generator
588 : //
589 : // This has to go into the Config.C
590 : //
591 : // AliGenPythia *pythia = new AliGenPythia(1);
592 : // pythia->SetProcess(kPyMbMSEL1);
593 : // pythia->SetStrucFunc(kCTEQ5L);
594 : // pythia->SetEnergyCMS(14000.);
595 : // AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
596 : // Decay_t dt = gener->GetDecayModePythia();
597 : // pythia->SetForceDecay(dt);
598 : // pythia->SetPtRange(0.,100.);
599 : // pythia->SetYRange(-8.,8.);
600 : // pythia->SetPhiRange(0.,360.);
601 : // pythia->SetPtHard(2.76,-1.0);
602 : // pythia->SwitchHFOff();
603 : // pythia->Init();
604 : // AddGenerator(pythia,"Pythia",1);
605 :
606 0 : }
607 :
608 : //-------------------------------------------------------------------
609 : void AliGenMUONCocktailpp::AddReso2Generator(Char_t* nameReso,
610 : AliGenParam* const genReso,
611 : Double_t sigmaReso,
612 : Double_t polReso)
613 : {
614 : // add resonances to the cocktail
615 0 : Double_t phiMin = fPhiMin*180./TMath::Pi();
616 0 : Double_t phiMax = fPhiMax*180./TMath::Pi();
617 :
618 : // first step: generation in 4pi
619 0 : genReso->SetPtRange(0.,100.);
620 0 : genReso->SetYRange(-8.,8.);
621 0 : genReso->SetPhiRange(0.,360.);
622 0 : genReso->SetForceDecay(fDecayModeResonance);
623 0 : if (!TVirtualMC::GetMC()) genReso->SetDecayer(fDecayer);
624 0 : genReso->Init(); // generation in 4pi
625 : // Ratios with respect to the reaction cross-section in the
626 : // kinematics limit of the MUONCocktail
627 0 : Double_t ratioReso = sigmaReso / fSigmaReaction * genReso->GetRelativeArea(fPtMin,fPtMax,fYMin,fYMax,phiMin,phiMax);
628 0 : if (!fSigmaSilent) {
629 0 : AliInfo(Form("%s prod. cross-section in pp %5.3g b",nameReso,sigmaReso));
630 0 : AliInfo(Form("%s prod. probability per collision in acceptance %5.3g",nameReso,ratioReso));
631 0 : }
632 : // second step: generation in selected kinematical range
633 0 : genReso->SetPtRange(fPtMin, fPtMax);
634 0 : genReso->SetYRange(fYMin, fYMax);
635 0 : genReso->SetPhiRange(phiMin, phiMax);
636 0 : genReso->Init(); // generation in selected kinematical range
637 :
638 : TVirtualMCDecayer *decReso = 0;
639 0 : if(TMath::Abs(polReso) > 1.e-30){
640 0 : AliInfo(Form("******Setting polarized decayer for %s''",nameReso));
641 0 : if(fPolFrame==0) {
642 0 : decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kColSop,AliDecayerPolarized::kMuon);
643 0 : AliInfo(Form("******Reference frame: %s, alpha: %f","Collins-Soper",polReso));
644 0 : }
645 0 : if(fPolFrame==1) {
646 0 : decReso = new AliDecayerPolarized(polReso,AliDecayerPolarized::kHelicity,AliDecayerPolarized::kMuon);
647 0 : AliInfo(Form("******Reference frame: %s, alpha: %f","Helicity",polReso));
648 0 : }
649 0 : if (decReso) {
650 0 : decReso->SetForceDecay(kAll);
651 0 : decReso->Init();
652 0 : genReso->SetDecayer(decReso);
653 0 : }
654 : }
655 :
656 0 : AddGenerator(genReso,nameReso,ratioReso); // Adding Generator
657 0 : }
658 :
659 : //-------------------------------------------------------------------
660 : void AliGenMUONCocktailpp::Init()
661 : {
662 : // Initialisation
663 0 : TIter next(fEntries);
664 : AliGenCocktailEntry *entry;
665 0 : if (fStack) {
666 0 : while((entry = (AliGenCocktailEntry*)next())) {
667 0 : entry->Generator()->SetStack(fStack);
668 : }
669 : }
670 0 : }
671 :
672 : //_________________________________________________________________________
673 : void AliGenMUONCocktailpp::Generate()
674 : {
675 : // Generate event
676 0 : TIter next(fEntries);
677 : AliGenCocktailEntry *entry = 0;
678 : AliGenerator* gen = 0;
679 :
680 0 : if (fHeader) delete fHeader;
681 0 : fHeader = new AliGenCocktailEventHeader("MUON Cocktail Header");
682 :
683 0 : const TObjArray *partArray = gAlice->GetMCApp()->Particles();
684 :
685 : // Generate the vertex position used by all generators
686 0 : if(fVertexSmear == kPerEvent) Vertex();
687 :
688 : // Loop on primordialTrigger:
689 : // minimum muon multiplicity above a pt cut in a theta acceptance region
690 :
691 : Bool_t primordialTrigger = kFALSE;
692 0 : while(!primordialTrigger) {
693 : //Reseting stack (if fMuonMultiplicity > 0 : S. Grigoryan, June 2012)
694 0 : AliRunLoader * runloader = AliRunLoader::Instance();
695 : // if (runloader)
696 0 : if (runloader && fMuonMultiplicity > 0)
697 0 : if (runloader->Stack())
698 0 : runloader->Stack()->Clean();
699 : // Loop over generators and generate events
700 : Int_t igen = 0;
701 : Int_t npart = 0;
702 : const char* genName = 0;
703 0 : while((entry = (AliGenCocktailEntry*)next())) {
704 0 : gen = entry->Generator();
705 0 : genName = entry->GetName();
706 0 : gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
707 0 : gen->SetTime(fTime);
708 :
709 0 : npart = (strcmp(genName,"Pythia") == 0) ? 1 :
710 0 : gRandom->Poisson(entry->Rate());
711 :
712 0 : if (npart > 0) {
713 0 : igen++;
714 0 : if (igen == 1) entry->SetFirst(0);
715 0 : else entry->SetFirst((partArray->GetEntriesFast())+1);
716 :
717 0 : gen->SetNumberParticles(npart);
718 0 : gen->Generate();
719 0 : entry->SetLast(partArray->GetEntriesFast());
720 0 : }
721 : }
722 0 : next.Reset();
723 :
724 :
725 : // Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
726 : // in the muon spectrometer acceptance
727 : Int_t iPart;
728 0 : fNGenerated++;
729 0 : Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
730 0 : for(iPart=0; iPart<maxPart; iPart++){
731 :
732 0 : TParticle *part = gAlice->GetMCApp()->Particle(iPart);
733 0 : if ( TMath::Abs(part->GetPdgCode()) == 13 ){
734 0 : if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
735 0 : (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
736 0 : (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
737 0 : (part->Pt()>fMuonPtCut) &&
738 0 : (part->P()>fMuonPCut)) {
739 0 : numberOfMuons++;
740 0 : }
741 : }
742 : }
743 0 : if (numberOfMuons >= fMuonMultiplicity) {
744 : primordialTrigger = kTRUE;
745 0 : fHeader->SetNProduced(maxPart);
746 : }
747 :
748 : }
749 0 : fNSucceded++;
750 :
751 0 : TArrayF eventVertex;
752 0 : eventVertex.Set(3);
753 0 : for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
754 :
755 0 : fHeader->SetPrimaryVertex(eventVertex);
756 0 : fHeader->SetInteractionTime(fTime);
757 :
758 0 : gAlice->SetGenEventHeader(fHeader);
759 :
760 : // AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
761 0 : AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
762 0 : }
763 :
764 :
|
