Changes related to the initialization of random numbers generators. Now one can use...
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia.cxx
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 #include "AliPythia.h"
19 #include "AliPythiaRndm.h"
20
21 ClassImp(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
33 extern "C" void type_of_call pyclus(Int_t & );
34 extern "C" void type_of_call pycell(Int_t & );
35
36 //_____________________________________________________________________________
37
38 AliPythia* AliPythia::fgAliPythia=NULL;
39
40 AliPythia::AliPythia()
41 {
42 // Default Constructor
43 //
44 //  Set random number
45     if (!AliPythiaRndm::GetPythiaRandom()) 
46       AliPythiaRndm::SetPythiaRandom(GetRandom());
47
48 }
49
50 void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
51 {
52 // Initialise the process to generate 
53     if (!AliPythiaRndm::GetPythiaRandom()) 
54       AliPythiaRndm::SetPythiaRandom(GetRandom());
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;
168     case kPyCharmPbPbMNR:
169     case kPyD0PbPbMNR:
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
174       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
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
197       // Intrinsic <kT>
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;
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:
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
312       // Intrinsic <kT>
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;
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;
399     }
400 //
401 //  Initialize PYTHIA
402     SetMSTP(41,1);   // all resonance decays switched on
403
404     Initialize("CMS","p","p",fEcms);
405
406 }
407
408 Int_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
416 void 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
433 AliPythia* 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
444 void 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
471 void  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
479 void  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
488 void  AliPythia::Pyclus(Int_t& njet)
489 {
490 //  Call Pythia clustering algorithm
491 //
492     pyclus(njet);
493 }
494
495 void  AliPythia::Pycell(Int_t& njet)
496 {
497 //  Call Pythia jet reconstruction algorithm
498 //
499     pycell(njet);
500 }
501
502
503