Option kPyCharmPbMNR added. Produce charm pairs in agreement with MNR
[u/mrichter/AliRoot.git] / EVGEN / 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.19  2002/02/20 08:52:20  morsch
19 Correct documentation of SetNuclei method.
20
21 Revision 1.18  2002/02/07 10:43:06  morsch
22 Tuned pp-min.bias settings (M.Monteno, R.Ugoccioni and N.Carrer)
23
24 Revision 1.17  2001/12/19 15:40:43  morsch
25 For kPyJets enforce simple jet topology, i.e no initial or final state
26 gluon radiation and no primordial pT.
27
28 Revision 1.16  2001/10/12 11:13:59  morsch
29 Missing break statements added (thanks to Nicola Carrer)
30
31 Revision 1.15  2001/03/27 10:54:50  morsch
32 Add ResetDecayTable() and SsetDecayTable() methods.
33
34 Revision 1.14  2001/03/09 13:03:40  morsch
35 Process_t and Struc_Func_t moved to AliPythia.h
36
37 Revision 1.13  2000/12/18 08:55:35  morsch
38 Make AliPythia dependent generartors work with new scheme of random number generation
39
40 Revision 1.12  2000/11/30 07:12:50  alibrary
41 Introducing new Rndm and QA classes
42
43 Revision 1.11  2000/10/20 06:30:06  fca
44 Use version 0 to avoid streamer generation
45
46 Revision 1.10  2000/10/06 14:18:44  morsch
47 Upper cut of prim. pT distribution set to 5. GeV
48
49 Revision 1.9  2000/09/18 10:41:35  morsch
50 Add possibility to use nuclear structure functions from PDF library V8.
51
52 Revision 1.8  2000/09/06 14:26:24  morsch
53 Decayer functionality of AliPythia has been moved to AliDecayerPythia.
54 Class is now a singleton.
55
56 Revision 1.7  2000/06/09 20:34:50  morsch
57 All coding rule violations except RS3 corrected
58
59 Revision 1.6  1999/11/09 07:38:48  fca
60 Changes for compatibility with version 2.23 of ROOT
61
62 Revision 1.5  1999/11/03 17:43:20  fca
63 New version from G.Martinez & A.Morsch
64
65 Revision 1.4  1999/09/29 09:24:14  fca
66 Introduction of the Copyright and cvs Log
67
68 */
69
70
71 #include "AliPythia.h"
72
73 ClassImp(AliPythia)
74
75 //_____________________________________________________________________________
76
77 AliPythia* AliPythia::fgAliPythia=NULL;
78
79 AliPythia::AliPythia()
80 {
81 // Default Constructor
82 //
83 //  Set random number
84     if (!sRandom) sRandom=fRandom;
85
86 }
87
88 void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
89 {
90 // Initialise the process to generate 
91     fProcess = process;
92     fEcms = energy;
93     fStrucFunc = strucfunc;
94 //  don't decay p0
95     SetMDCY(Pycomp(111),1,0);
96 //  select structure function 
97     SetMSTP(52,2);
98     SetMSTP(51,strucfunc);
99 //
100 // Pythia initialisation for selected processes//
101 //
102 // Make MSEL clean
103 //
104     for (Int_t i=1; i<= 200; i++) {
105         SetMSUB(i,0);
106     }
107 //  select charm production
108     switch (process) 
109     {
110     case kPyCharm:
111         SetMSEL(4);
112 //
113 //  heavy quark masses
114
115         SetPMAS(4,1,1.2);
116         SetMSTU(16,2);
117 //
118 //    primordial pT
119         SetMSTP(91,1);
120         SetPARP(91,1.);
121         SetPARP(93,5.);
122 //
123         break;
124     case kPyBeauty:
125         SetMSEL(5);
126         SetPMAS(5,1,4.75);
127         SetMSTU(16,2);
128         break;
129     case kPyJpsi:
130         SetMSEL(0);
131 // gg->J/Psi g
132         SetMSUB(86,1);
133         break;
134     case kPyJpsiChi:
135         SetMSEL(0);
136 // gg->J/Psi g
137         SetMSUB(86,1);
138 // gg-> chi_0c g
139         SetMSUB(87,1);
140 // gg-> chi_1c g
141         SetMSUB(88,1);
142 // gg-> chi_2c g
143         SetMSUB(89,1);  
144         break;
145     case kPyCharmUnforced:
146         SetMSEL(0);
147 // gq->qg   
148         SetMSUB(28,1);
149 // gg->qq
150         SetMSUB(53,1);
151 // gg->gg
152         SetMSUB(68,1);
153         break;
154     case kPyBeautyUnforced:
155         SetMSEL(0);
156 // gq->qg   
157         SetMSUB(28,1);
158 // gg->qq
159         SetMSUB(53,1);
160 // gg->gg
161         SetMSUB(68,1);
162         break;
163     case kPyMb:
164 // Minimum Bias pp-Collisions
165 //
166 //   
167 //      select Pythia min. bias model
168         SetMSEL(0);
169         SetMSUB(92,1);      // single diffraction AB-->XB
170         SetMSUB(93,1);      // single diffraction AB-->AX
171         SetMSUB(94,1);      // double diffraction
172         SetMSUB(95,1);      // low pt production
173         SetMSTP(81,1);      // multiple interactions switched on
174         SetMSTP(82,3);      // model with varying impact param. & a single Gaussian
175         SetPARP(82,3.47);   // set value pT_0  for turn-off of the cross section of                  
176                             // multiple interaction at a reference energy = 14000 GeV
177         SetPARP(89,14000.); // reference energy for the above parameter
178         SetPARP(90,0.174);  // set exponent for energy dependence of pT_0 
179         break;
180     case kPyJets:
181         SetMSEL(1);
182 // no initial state radiation   
183         SetMSTP(61,0);
184 // no final state radiation
185         SetMSTP(71,0);
186 // no primordial pT
187         SetMSTP(91,0);
188 //      SetMSTP(111,0); 
189         SetMSTU(16,1);  
190         SetMSTJ(1,1);
191         
192         break;
193     case kPyDirectGamma:
194         SetMSEL(10);
195         break;
196     case kPyCharmPbMNR:
197       // Tuning of Pythia parameters aimed to get a resonable agreement
198       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
199       // c-cbar single inclusive and double differential distributions.
200       // This parameter settings are meant to work with Pb-Pb collisions
201       // (AliGenPythia::SetNuclei) and with kCTEQ_4L PDFs.
202       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
203       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
204
205       // All QCD processes
206       SetMSEL(1);
207
208       // No multiple interactions
209       SetMSTP(81,0);
210       SetPARP(81,0.0);
211       SetPARP(82,0.0);
212
213       // Initial/final parton shower on (Pythia default)
214       SetMSTP(61,1);
215       SetMSTP(71,1);
216
217       // 2nd order alpha_s
218       SetMSTP(2,2);
219
220       // QCD scales
221       SetMSTP(32,2);
222       SetPARP(34,1.0);
223
224       // Intrinsic <kT^2>
225       SetMSTP(91,1);
226       SetPARP(91,1.304);
227       SetPARP(93,6.52);
228
229       // Set c-quark mass
230       SetPMAS(4,1,1.2);
231
232       break;
233     }
234 //
235 //  Initialize PYTHIA
236     SetMSTP(41,1);   // all resonance decays switched on
237
238     Initialize("CMS","p","p",fEcms);
239
240 }
241
242 Int_t AliPythia::CheckedLuComp(Int_t kf)
243 {
244 // Check Lund particle code (for debugging)
245     Int_t kc=Pycomp(kf);
246     printf("\n Lucomp kf,kc %d %d",kf,kc);
247     return kc;
248 }
249
250 void AliPythia::SetNuclei(Int_t a1, Int_t a2)
251 {
252 // Treat protons as inside nuclei with mass numbers a1 and a2  
253 //    The MSTP array in the PYPARS common block is used to enable and 
254 //    select the nuclear structure functions. 
255 //    MSTP(52)  : (D=1) choice of proton and nuclear structure-function library
256 //            =1: internal PYTHIA acording to MSTP(51) 
257 //            =2: PDFLIB proton  s.f., with MSTP(51)  = 1000xNGROUP+NSET
258 //    If the following mass number both not equal zero, nuclear corrections of the stf are used.
259 //    MSTP(192) : Mass number of nucleus side 1
260 //    MSTP(193) : Mass number of nucleus side 2
261     SetMSTP(52,2);
262     SetMSTP(192, a1);
263     SetMSTP(193, a2);  
264 }
265         
266
267 AliPythia* AliPythia::Instance()
268
269 // Set random number generator 
270     if (fgAliPythia) {
271         return fgAliPythia;
272     } else {
273         fgAliPythia = new AliPythia();
274         return fgAliPythia;
275     }
276 }
277
278 void AliPythia::PrintParticles()
279
280 // Print list of particl properties
281     Int_t np = 0;
282     
283     for (Int_t kf=0; kf<1000000; kf++) {
284         for (Int_t c = 1;  c > -2; c-=2) {
285             
286             Int_t kc = Pycomp(c*kf);
287             if (kc) {
288                 Float_t mass  = GetPMAS(kc,1);
289                 Float_t width = GetPMAS(kc,2);  
290                 Float_t tau   = GetPMAS(kc,4);
291                 
292                 char*   name = new char[8];
293                 Pyname(kf,name);
294         
295                 np++;
296                 
297                 printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e", 
298                        c*kf, name, mass, width, tau);
299             }
300         }
301     }
302     printf("\n Number of particles %d \n \n", np);
303 }
304
305 void  AliPythia::ResetDecayTable()
306 {
307 //  Set default values for pythia decay switches
308     Int_t i;
309     for (i = 1; i <  501; i++) SetMDCY(i,1,fDefMDCY[i]);
310     for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
311 }
312
313 void  AliPythia::SetDecayTable()
314 {
315 //  Set default values for pythia decay switches
316 //
317     Int_t i;
318     for (i = 1; i <  501; i++) fDefMDCY[i] = GetMDCY(i,1);
319     for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
320 }
321
322
323 #ifndef WIN32
324 #define pyr    pyr_
325 #define pyrset pyrset_
326 #define pyrget pyrget_
327 #else
328 #define pyr    PYR
329 #define pyrset PYRSET
330 #define pyrget PYRGET
331 #endif
332
333 extern "C" {
334   Double_t pyr(Int_t*) {return sRandom->Rndm();}
335   void pyrset(Int_t*,Int_t*) {}
336   void pyrget(Int_t*,Int_t*) {}
337 }
338
339
340
341