Classed imported from EVGEN.
[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 /*
17 $Log$
18 Revision 1.28  2002/12/09 08:22:56  morsch
19 UA1 jet finder (Pycell) for software triggering added.
20
21 Revision 1.27  2002/11/15 00:39:37  morsch
22 - Correct initialisation of sRandom.
23 - QCD Jets with initial and final state gluon radiation is default
24 - pt kick for jets default
25 - Interface to Pyclus added.
26
27 Revision 1.26  2002/11/14 00:37:32  morsch
28 Warning message for kPyJets added.
29
30 Revision 1.25  2002/10/14 14:55:35  hristov
31 Merging the VirtualMC branch to the main development branch (HEAD)
32
33 Revision 1.20.6.1  2002/06/10 14:57:41  hristov
34 Merged with v3-08-02
35
36 Revision 1.24  2002/05/22 13:22:53  morsch
37 Process kPyMbNonDiffr added.
38
39 Revision 1.23  2002/05/06 07:17:29  morsch
40 Pyr gives random number r in interval 0 < r < 1.
41
42 Revision 1.22  2002/04/26 10:28:48  morsch
43 Option kPyBeautyPbMNR added (N. Carrer).
44
45 Revision 1.21  2002/03/25 14:46:16  morsch
46 Case  kPyD0PbMNR added (N. Carrer).
47
48 Revision 1.20  2002/03/03 13:48:50  morsch
49 Option  kPyCharmPbMNR added. Produce charm pairs in agreement with MNR
50 NLO calculations (Nicola Carrer).
51
52 Revision 1.19  2002/02/20 08:52:20  morsch
53 Correct documentation of SetNuclei method.
54
55 Revision 1.18  2002/02/07 10:43:06  morsch
56 Tuned pp-min.bias settings (M.Monteno, R.Ugoccioni and N.Carrer)
57
58 Revision 1.17  2001/12/19 15:40:43  morsch
59 For kPyJets enforce simple jet topology, i.e no initial or final state
60 gluon radiation and no primordial pT.
61
62 Revision 1.16  2001/10/12 11:13:59  morsch
63 Missing break statements added (thanks to Nicola Carrer)
64
65 Revision 1.15  2001/03/27 10:54:50  morsch
66 Add ResetDecayTable() and SsetDecayTable() methods.
67
68 Revision 1.14  2001/03/09 13:03:40  morsch
69 Process_t and Struc_Func_t moved to AliPythia.h
70
71 Revision 1.13  2000/12/18 08:55:35  morsch
72 Make AliPythia dependent generartors work with new scheme of random number generation
73
74 Revision 1.12  2000/11/30 07:12:50  alibrary
75 Introducing new Rndm and QA classes
76
77 Revision 1.11  2000/10/20 06:30:06  fca
78 Use version 0 to avoid streamer generation
79
80 Revision 1.10  2000/10/06 14:18:44  morsch
81 Upper cut of prim. pT distribution set to 5. GeV
82
83 Revision 1.9  2000/09/18 10:41:35  morsch
84 Add possibility to use nuclear structure functions from PDF library V8.
85
86 Revision 1.8  2000/09/06 14:26:24  morsch
87 Decayer functionality of AliPythia has been moved to AliDecayerPythia.
88 Class is now a singleton.
89
90 Revision 1.7  2000/06/09 20:34:50  morsch
91 All coding rule violations except RS3 corrected
92
93 Revision 1.6  1999/11/09 07:38:48  fca
94 Changes for compatibility with version 2.23 of ROOT
95
96 Revision 1.5  1999/11/03 17:43:20  fca
97 New version from G.Martinez & A.Morsch
98
99 Revision 1.4  1999/09/29 09:24:14  fca
100 Introduction of the Copyright and cvs Log
101
102 */
103
104
105 #include "AliPythia.h"
106
107 ClassImp(AliPythia)
108
109 #ifndef WIN32
110 # define pyclus pyclus_
111 # define pycell pycell_
112 # define type_of_call
113 #else
114 # define pyclus PYCLUS
115 # define pycell PYCELL
116 # define type_of_call _stdcall
117 #endif
118
119 extern "C" void type_of_call pyclus(Int_t & );
120 extern "C" void type_of_call pycell(Int_t & );
121
122 //_____________________________________________________________________________
123
124 AliPythia* AliPythia::fgAliPythia=NULL;
125
126 AliPythia::AliPythia()
127 {
128 // Default Constructor
129 //
130 //  Set random number
131     if (!sRandom) sRandom=fRandom;
132
133 }
134
135 void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
136 {
137 // Initialise the process to generate 
138     if (!sRandom) sRandom = gRandom;
139     
140     fProcess = process;
141     fEcms = energy;
142     fStrucFunc = strucfunc;
143 //  don't decay p0
144     SetMDCY(Pycomp(111),1,0);
145 //  select structure function 
146     SetMSTP(52,2);
147     SetMSTP(51,strucfunc);
148 //
149 // Pythia initialisation for selected processes//
150 //
151 // Make MSEL clean
152 //
153     for (Int_t i=1; i<= 200; i++) {
154         SetMSUB(i,0);
155     }
156 //  select charm production
157     switch (process) 
158     {
159     case kPyCharm:
160         SetMSEL(4);
161 //
162 //  heavy quark masses
163
164         SetPMAS(4,1,1.2);
165         SetMSTU(16,2);
166 //
167 //    primordial pT
168         SetMSTP(91,1);
169         SetPARP(91,1.);
170         SetPARP(93,5.);
171 //
172         break;
173     case kPyBeauty:
174         SetMSEL(5);
175         SetPMAS(5,1,4.75);
176         SetMSTU(16,2);
177         break;
178     case kPyJpsi:
179         SetMSEL(0);
180 // gg->J/Psi g
181         SetMSUB(86,1);
182         break;
183     case kPyJpsiChi:
184         SetMSEL(0);
185 // gg->J/Psi g
186         SetMSUB(86,1);
187 // gg-> chi_0c g
188         SetMSUB(87,1);
189 // gg-> chi_1c g
190         SetMSUB(88,1);
191 // gg-> chi_2c g
192         SetMSUB(89,1);  
193         break;
194     case kPyCharmUnforced:
195         SetMSEL(0);
196 // gq->qg   
197         SetMSUB(28,1);
198 // gg->qq
199         SetMSUB(53,1);
200 // gg->gg
201         SetMSUB(68,1);
202         break;
203     case kPyBeautyUnforced:
204         SetMSEL(0);
205 // gq->qg   
206         SetMSUB(28,1);
207 // gg->qq
208         SetMSUB(53,1);
209 // gg->gg
210         SetMSUB(68,1);
211         break;
212     case kPyMb:
213 // Minimum Bias pp-Collisions
214 //
215 //   
216 //      select Pythia min. bias model
217         SetMSEL(0);
218         SetMSUB(92,1);      // single diffraction AB-->XB
219         SetMSUB(93,1);      // single diffraction AB-->AX
220         SetMSUB(94,1);      // double diffraction
221         SetMSUB(95,1);      // low pt production
222         SetMSTP(81,1);      // multiple interactions switched on
223         SetMSTP(82,3);      // model with varying impact param. & a single Gaussian
224         SetPARP(82,3.47);   // set value pT_0  for turn-off of the cross section of                  
225                             // multiple interaction at a reference energy = 14000 GeV
226         SetPARP(89,14000.); // reference energy for the above parameter
227         SetPARP(90,0.174);  // set exponent for energy dependence of pT_0
228     case kPyMbNonDiffr:
229 // Minimum Bias pp-Collisions
230 //
231 //   
232 //      select Pythia min. bias model
233         SetMSEL(0);
234         SetMSUB(95,1);      // low pt production
235         SetMSTP(81,1);      // multiple interactions switched on
236         SetMSTP(82,3);      // model with varying impact param. & a single Gaussian
237         SetPARP(82,3.47);   // set value pT_0  for turn-off of the cross section of                  
238                             // multiple interaction at a reference energy = 14000 GeV
239         SetPARP(89,14000.); // reference energy for the above parameter
240         SetPARP(90,0.174);  // set exponent for energy dependence of pT_0
241  
242         break;
243     case kPyJets:
244 //
245 //  QCD Jets
246 //
247         SetMSEL(1);
248         break;
249     case kPyDirectGamma:
250         SetMSEL(10);
251         break;
252     case kPyCharmPbMNR:
253     case kPyD0PbMNR:
254       // Tuning of Pythia parameters aimed to get a resonable agreement
255       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
256       // c-cbar single inclusive and double differential distributions.
257       // This parameter settings are meant to work with Pb-Pb collisions
258       // (AliGenPythia::SetNuclei) and with kCTEQ_4L PDFs.
259       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
260       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
261
262       // All QCD processes
263       SetMSEL(1);
264
265       // No multiple interactions
266       SetMSTP(81,0);
267       SetPARP(81,0.0);
268       SetPARP(82,0.0);
269
270       // Initial/final parton shower on (Pythia default)
271       SetMSTP(61,1);
272       SetMSTP(71,1);
273
274       // 2nd order alpha_s
275       SetMSTP(2,2);
276
277       // QCD scales
278       SetMSTP(32,2);
279       SetPARP(34,1.0);
280
281       // Intrinsic <kT^2>
282       SetMSTP(91,1);
283       SetPARP(91,1.304);
284       SetPARP(93,6.52);
285
286       // Set c-quark mass
287       SetPMAS(4,1,1.2);
288
289       break;
290     case kPyBeautyPbMNR:
291       // Tuning of Pythia parameters aimed to get a resonable agreement
292       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
293       // b-bbar single inclusive and double differential distributions.
294       // This parameter settings are meant to work with Pb-Pb collisions
295       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
296       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
297       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
298
299       // All QCD processes
300       SetMSEL(1);
301
302       // No multiple interactions
303       SetMSTP(81,0);
304       SetPARP(81,0.0);
305       SetPARP(82,0.0);
306
307       // Initial/final parton shower on (Pythia default)
308       SetMSTP(61,1);
309       SetMSTP(71,1);
310
311       // 2nd order alpha_s
312       SetMSTP(2,2);
313
314       // QCD scales
315       SetMSTP(32,2);
316       SetPARP(34,1.0);
317       SetPARP(67,1.0);
318       SetPARP(71,1.0);
319
320       // Intrinsic <kT^2>
321       SetMSTP(91,1);
322       SetPARP(91,2.035);
323       SetPARP(93,10.17);
324
325       // Set b-quark mass
326       SetPMAS(5,1,4.75);
327
328       break;
329     }
330 //
331 //  Initialize PYTHIA
332     SetMSTP(41,1);   // all resonance decays switched on
333
334     Initialize("CMS","p","p",fEcms);
335
336 }
337
338 Int_t AliPythia::CheckedLuComp(Int_t kf)
339 {
340 // Check Lund particle code (for debugging)
341     Int_t kc=Pycomp(kf);
342     printf("\n Lucomp kf,kc %d %d",kf,kc);
343     return kc;
344 }
345
346 void AliPythia::SetNuclei(Int_t a1, Int_t a2)
347 {
348 // Treat protons as inside nuclei with mass numbers a1 and a2  
349 //    The MSTP array in the PYPARS common block is used to enable and 
350 //    select the nuclear structure functions. 
351 //    MSTP(52)  : (D=1) choice of proton and nuclear structure-function library
352 //            =1: internal PYTHIA acording to MSTP(51) 
353 //            =2: PDFLIB proton  s.f., with MSTP(51)  = 1000xNGROUP+NSET
354 //    If the following mass number both not equal zero, nuclear corrections of the stf are used.
355 //    MSTP(192) : Mass number of nucleus side 1
356 //    MSTP(193) : Mass number of nucleus side 2
357     SetMSTP(52,2);
358     SetMSTP(192, a1);
359     SetMSTP(193, a2);  
360 }
361         
362
363 AliPythia* AliPythia::Instance()
364
365 // Set random number generator 
366     if (fgAliPythia) {
367         return fgAliPythia;
368     } else {
369         fgAliPythia = new AliPythia();
370         return fgAliPythia;
371     }
372 }
373
374 void AliPythia::PrintParticles()
375
376 // Print list of particl properties
377     Int_t np = 0;
378     
379     for (Int_t kf=0; kf<1000000; kf++) {
380         for (Int_t c = 1;  c > -2; c-=2) {
381             
382             Int_t kc = Pycomp(c*kf);
383             if (kc) {
384                 Float_t mass  = GetPMAS(kc,1);
385                 Float_t width = GetPMAS(kc,2);  
386                 Float_t tau   = GetPMAS(kc,4);
387                 
388                 char*   name = new char[8];
389                 Pyname(kf,name);
390         
391                 np++;
392                 
393                 printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e", 
394                        c*kf, name, mass, width, tau);
395             }
396         }
397     }
398     printf("\n Number of particles %d \n \n", np);
399 }
400
401 void  AliPythia::ResetDecayTable()
402 {
403 //  Set default values for pythia decay switches
404     Int_t i;
405     for (i = 1; i <  501; i++) SetMDCY(i,1,fDefMDCY[i]);
406     for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
407 }
408
409 void  AliPythia::SetDecayTable()
410 {
411 //  Set default values for pythia decay switches
412 //
413     Int_t i;
414     for (i = 1; i <  501; i++) fDefMDCY[i] = GetMDCY(i,1);
415     for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
416 }
417
418 void  AliPythia::Pyclus(Int_t& njet)
419 {
420 //  Call Pythia clustering algorithm
421 //
422     pyclus(njet);
423 }
424
425 void  AliPythia::Pycell(Int_t& njet)
426 {
427 //  Call Pythia jet reconstruction algorithm
428 //
429     pycell(njet);
430 }
431
432
433
434 #ifndef WIN32
435 #define pyr    pyr_
436 #define pyrset pyrset_
437 #define pyrget pyrget_
438 #define pyclus pyclus_
439 #define pycell pycell_
440 #else
441 #define pyr    PYR
442 #define pyrset PYRSET
443 #define pyrget PYRGET
444 #define pyclus PYCLUS
445 #define pycell PYCELL
446 #endif
447
448 extern "C" {
449   Double_t pyr(Int_t*) 
450 {
451       Float_t r;
452       do r=sRandom->Rndm(); while(0 >= r || r >= 1);
453       return r;
454 }
455   void pyrset(Int_t*,Int_t*) {}
456   void pyrget(Int_t*,Int_t*) {}
457 }
458
459
460
461
462