Changes related to the initialization of random numbers generators. Now one can use...
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia.cxx
CommitLineData
8d2cd130 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
7cdba479 16/* $Id$ */
8d2cd130 17
18#include "AliPythia.h"
7cdba479 19#include "AliPythiaRndm.h"
8d2cd130 20
21ClassImp(AliPythia)
22
23#ifndef WIN32
24# define pyclus pyclus_
25# define pycell pycell_
26# define type_of_call
27#else
28# define pyclus PYCLUS
29# define pycell PYCELL
30# define type_of_call _stdcall
31#endif
32
33extern "C" void type_of_call pyclus(Int_t & );
34extern "C" void type_of_call pycell(Int_t & );
35
36//_____________________________________________________________________________
37
38AliPythia* AliPythia::fgAliPythia=NULL;
39
40AliPythia::AliPythia()
41{
42// Default Constructor
43//
44// Set random number
7cdba479 45 if (!AliPythiaRndm::GetPythiaRandom())
46 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 47
48}
49
50void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
51{
52// Initialise the process to generate
7cdba479 53 if (!AliPythiaRndm::GetPythiaRandom())
54 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 55
56 fProcess = process;
57 fEcms = energy;
58 fStrucFunc = strucfunc;
59// don't decay p0
60 SetMDCY(Pycomp(111),1,0);
61// select structure function
62 SetMSTP(52,2);
63 SetMSTP(51,strucfunc);
64//
65// Pythia initialisation for selected processes//
66//
67// Make MSEL clean
68//
69 for (Int_t i=1; i<= 200; i++) {
70 SetMSUB(i,0);
71 }
72// select charm production
73 switch (process)
74 {
75 case kPyCharm:
76 SetMSEL(4);
77//
78// heavy quark masses
79
80 SetPMAS(4,1,1.2);
81 SetMSTU(16,2);
82//
83// primordial pT
84 SetMSTP(91,1);
85 SetPARP(91,1.);
86 SetPARP(93,5.);
87//
88 break;
89 case kPyBeauty:
90 SetMSEL(5);
91 SetPMAS(5,1,4.75);
92 SetMSTU(16,2);
93 break;
94 case kPyJpsi:
95 SetMSEL(0);
96// gg->J/Psi g
97 SetMSUB(86,1);
98 break;
99 case kPyJpsiChi:
100 SetMSEL(0);
101// gg->J/Psi g
102 SetMSUB(86,1);
103// gg-> chi_0c g
104 SetMSUB(87,1);
105// gg-> chi_1c g
106 SetMSUB(88,1);
107// gg-> chi_2c g
108 SetMSUB(89,1);
109 break;
110 case kPyCharmUnforced:
111 SetMSEL(0);
112// gq->qg
113 SetMSUB(28,1);
114// gg->qq
115 SetMSUB(53,1);
116// gg->gg
117 SetMSUB(68,1);
118 break;
119 case kPyBeautyUnforced:
120 SetMSEL(0);
121// gq->qg
122 SetMSUB(28,1);
123// gg->qq
124 SetMSUB(53,1);
125// gg->gg
126 SetMSUB(68,1);
127 break;
128 case kPyMb:
129// Minimum Bias pp-Collisions
130//
131//
132// select Pythia min. bias model
133 SetMSEL(0);
134 SetMSUB(92,1); // single diffraction AB-->XB
135 SetMSUB(93,1); // single diffraction AB-->AX
136 SetMSUB(94,1); // double diffraction
137 SetMSUB(95,1); // low pt production
138 SetMSTP(81,1); // multiple interactions switched on
139 SetMSTP(82,3); // model with varying impact param. & a single Gaussian
140 SetPARP(82,3.47); // set value pT_0 for turn-off of the cross section of
141 // multiple interaction at a reference energy = 14000 GeV
142 SetPARP(89,14000.); // reference energy for the above parameter
143 SetPARP(90,0.174); // set exponent for energy dependence of pT_0
144 case kPyMbNonDiffr:
145// Minimum Bias pp-Collisions
146//
147//
148// select Pythia min. bias model
149 SetMSEL(0);
150 SetMSUB(95,1); // low pt production
151 SetMSTP(81,1); // multiple interactions switched on
152 SetMSTP(82,3); // model with varying impact param. & a single Gaussian
153 SetPARP(82,3.47); // set value pT_0 for turn-off of the cross section of
154 // multiple interaction at a reference energy = 14000 GeV
155 SetPARP(89,14000.); // reference energy for the above parameter
156 SetPARP(90,0.174); // set exponent for energy dependence of pT_0
157
158 break;
159 case kPyJets:
160//
161// QCD Jets
162//
163 SetMSEL(1);
164 break;
165 case kPyDirectGamma:
166 SetMSEL(10);
167 break;
adf4d898 168 case kPyCharmPbPbMNR:
169 case kPyD0PbPbMNR:
8d2cd130 170 // Tuning of Pythia parameters aimed to get a resonable agreement
171 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
172 // c-cbar single inclusive and double differential distributions.
173 // This parameter settings are meant to work with Pb-Pb collisions
adf4d898 174 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
8d2cd130 175 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
176 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
177
178 // All QCD processes
179 SetMSEL(1);
180
181 // No multiple interactions
182 SetMSTP(81,0);
183 SetPARP(81,0.0);
184 SetPARP(82,0.0);
185
186 // Initial/final parton shower on (Pythia default)
187 SetMSTP(61,1);
188 SetMSTP(71,1);
189
190 // 2nd order alpha_s
191 SetMSTP(2,2);
192
193 // QCD scales
194 SetMSTP(32,2);
195 SetPARP(34,1.0);
196
adf4d898 197 // Intrinsic <kT>
8d2cd130 198 SetMSTP(91,1);
199 SetPARP(91,1.304);
200 SetPARP(93,6.52);
201
202 // Set c-quark mass
203 SetPMAS(4,1,1.2);
204
205 break;
adf4d898 206 case kPyCharmpPbMNR:
207 case kPyD0pPbMNR:
208 // Tuning of Pythia parameters aimed to get a resonable agreement
209 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
210 // c-cbar single inclusive and double differential distributions.
211 // This parameter settings are meant to work with p-Pb collisions
212 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
213 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
214 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
215
216 // All QCD processes
217 SetMSEL(1);
218
219 // No multiple interactions
220 SetMSTP(81,0);
221 SetPARP(81,0.0);
222 SetPARP(82,0.0);
223
224 // Initial/final parton shower on (Pythia default)
225 SetMSTP(61,1);
226 SetMSTP(71,1);
227
228 // 2nd order alpha_s
229 SetMSTP(2,2);
230
231 // QCD scales
232 SetMSTP(32,2);
233 SetPARP(34,1.0);
234
235 // Intrinsic <kT>
236 SetMSTP(91,1);
237 SetPARP(91,1.16);
238 SetPARP(93,5.8);
239
240 // Set c-quark mass
241 SetPMAS(4,1,1.2);
242
243 break;
244 case kPyCharmppMNR:
245 case kPyD0ppMNR:
246 // Tuning of Pythia parameters aimed to get a resonable agreement
247 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
248 // c-cbar single inclusive and double differential distributions.
249 // This parameter settings are meant to work with pp collisions
250 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
251 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
252 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
253
254 // All QCD processes
255 SetMSEL(1);
256
257 // No multiple interactions
258 SetMSTP(81,0);
259 SetPARP(81,0.0);
260 SetPARP(82,0.0);
261
262 // Initial/final parton shower on (Pythia default)
263 SetMSTP(61,1);
264 SetMSTP(71,1);
265
266 // 2nd order alpha_s
267 SetMSTP(2,2);
268
269 // QCD scales
270 SetMSTP(32,2);
271 SetPARP(34,1.0);
272
273 // Intrinsic <kT^2>
274 SetMSTP(91,1);
275 SetPARP(91,1.);
276 SetPARP(93,5.);
277
278 // Set c-quark mass
279 SetPMAS(4,1,1.2);
280
281 break;
282 case kPyBeautyPbPbMNR:
8d2cd130 283 // Tuning of Pythia parameters aimed to get a resonable agreement
284 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
285 // b-bbar single inclusive and double differential distributions.
286 // This parameter settings are meant to work with Pb-Pb collisions
287 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
288 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
289 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
290
291 // All QCD processes
292 SetMSEL(1);
293
294 // No multiple interactions
295 SetMSTP(81,0);
296 SetPARP(81,0.0);
297 SetPARP(82,0.0);
298
299 // Initial/final parton shower on (Pythia default)
300 SetMSTP(61,1);
301 SetMSTP(71,1);
302
303 // 2nd order alpha_s
304 SetMSTP(2,2);
305
306 // QCD scales
307 SetMSTP(32,2);
308 SetPARP(34,1.0);
309 SetPARP(67,1.0);
310 SetPARP(71,1.0);
311
adf4d898 312 // Intrinsic <kT>
8d2cd130 313 SetMSTP(91,1);
314 SetPARP(91,2.035);
315 SetPARP(93,10.17);
316
317 // Set b-quark mass
318 SetPMAS(5,1,4.75);
319
320 break;
adf4d898 321 case kPyBeautypPbMNR:
322 // Tuning of Pythia parameters aimed to get a resonable agreement
323 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
324 // b-bbar single inclusive and double differential distributions.
325 // This parameter settings are meant to work with p-Pb collisions
326 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
327 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
328 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
329
330 // All QCD processes
331 SetMSEL(1);
332
333 // No multiple interactions
334 SetMSTP(81,0);
335 SetPARP(81,0.0);
336 SetPARP(82,0.0);
337
338 // Initial/final parton shower on (Pythia default)
339 SetMSTP(61,1);
340 SetMSTP(71,1);
341
342 // 2nd order alpha_s
343 SetMSTP(2,2);
344
345 // QCD scales
346 SetMSTP(32,2);
347 SetPARP(34,1.0);
348 SetPARP(67,1.0);
349 SetPARP(71,1.0);
350
351 // Intrinsic <kT>
352 SetMSTP(91,1);
353 SetPARP(91,1.60);
354 SetPARP(93,8.00);
355
356 // Set b-quark mass
357 SetPMAS(5,1,4.75);
358
359 break;
360 case kPyBeautyppMNR:
361 // Tuning of Pythia parameters aimed to get a resonable agreement
362 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
363 // b-bbar single inclusive and double differential distributions.
364 // This parameter settings are meant to work with pp collisions
365 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
366 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
367 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
368
369 // All QCD processes
370 SetMSEL(1);
371
372 // No multiple interactions
373 SetMSTP(81,0);
374 SetPARP(81,0.0);
375 SetPARP(82,0.0);
376
377 // Initial/final parton shower on (Pythia default)
378 SetMSTP(61,1);
379 SetMSTP(71,1);
380
381 // 2nd order alpha_s
382 SetMSTP(2,2);
383
384 // QCD scales
385 SetMSTP(32,2);
386 SetPARP(34,1.0);
387 SetPARP(67,1.0);
388 SetPARP(71,1.0);
389
390 // Intrinsic <kT>
391 SetMSTP(91,1);
392 SetPARP(91,1.);
393 SetPARP(93,5.);
394
395 // Set b-quark mass
396 SetPMAS(5,1,4.75);
397
398 break;
8d2cd130 399 }
400//
401// Initialize PYTHIA
402 SetMSTP(41,1); // all resonance decays switched on
403
404 Initialize("CMS","p","p",fEcms);
405
406}
407
408Int_t AliPythia::CheckedLuComp(Int_t kf)
409{
410// Check Lund particle code (for debugging)
411 Int_t kc=Pycomp(kf);
412 printf("\n Lucomp kf,kc %d %d",kf,kc);
413 return kc;
414}
415
416void AliPythia::SetNuclei(Int_t a1, Int_t a2)
417{
418// Treat protons as inside nuclei with mass numbers a1 and a2
419// The MSTP array in the PYPARS common block is used to enable and
420// select the nuclear structure functions.
421// MSTP(52) : (D=1) choice of proton and nuclear structure-function library
422// =1: internal PYTHIA acording to MSTP(51)
423// =2: PDFLIB proton s.f., with MSTP(51) = 1000xNGROUP+NSET
424// If the following mass number both not equal zero, nuclear corrections of the stf are used.
425// MSTP(192) : Mass number of nucleus side 1
426// MSTP(193) : Mass number of nucleus side 2
427 SetMSTP(52,2);
428 SetMSTP(192, a1);
429 SetMSTP(193, a2);
430}
431
432
433AliPythia* AliPythia::Instance()
434{
435// Set random number generator
436 if (fgAliPythia) {
437 return fgAliPythia;
438 } else {
439 fgAliPythia = new AliPythia();
440 return fgAliPythia;
441 }
442}
443
444void AliPythia::PrintParticles()
445{
446// Print list of particl properties
447 Int_t np = 0;
448
449 for (Int_t kf=0; kf<1000000; kf++) {
450 for (Int_t c = 1; c > -2; c-=2) {
451
452 Int_t kc = Pycomp(c*kf);
453 if (kc) {
454 Float_t mass = GetPMAS(kc,1);
455 Float_t width = GetPMAS(kc,2);
456 Float_t tau = GetPMAS(kc,4);
457
458 char* name = new char[8];
459 Pyname(kf,name);
460
461 np++;
462
463 printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
464 c*kf, name, mass, width, tau);
465 }
466 }
467 }
468 printf("\n Number of particles %d \n \n", np);
469}
470
471void AliPythia::ResetDecayTable()
472{
473// Set default values for pythia decay switches
474 Int_t i;
475 for (i = 1; i < 501; i++) SetMDCY(i,1,fDefMDCY[i]);
476 for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
477}
478
479void AliPythia::SetDecayTable()
480{
481// Set default values for pythia decay switches
482//
483 Int_t i;
484 for (i = 1; i < 501; i++) fDefMDCY[i] = GetMDCY(i,1);
485 for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
486}
487
488void AliPythia::Pyclus(Int_t& njet)
489{
490// Call Pythia clustering algorithm
491//
492 pyclus(njet);
493}
494
495void AliPythia::Pycell(Int_t& njet)
496{
497// Call Pythia jet reconstruction algorithm
498//
499 pycell(njet);
500}
501
502
503