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