0d6a1c7e15f002f861437f1756dda443f9d96484
[u/mrichter/AliRoot.git] / PYTHIA8 / AliPythia8.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          *
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 /* $Id$ */
18 #include <TString.h>
19 #include <TVector3.h>
20 #include <TMath.h>
21
22 #include "AliPythia8.h"
23 #include "AliLog.h"
24 #include "AliStack.h"
25 #include "AliPythiaRndm.h"
26
27
28 ClassImp(AliPythia8)
29     
30 // Particles produced in string fragmentation point directly to either of the two endpoints
31 // of the string (depending in the side they were generated from).
32 //    SetMSTU(16,2); // ????
33 //  String drawing almost completely minimizes string length.
34 //  Probability that an additional interaction gives two gluons
35 //  ... with color connection to nearest neighbours
36 //    SetPARP(85,0.9);
37 //  ... as closed gluon loop
38 //    SetPARP(86,0.95);
39 // Lambda_FSR scale.
40 //      SetPARJ(81, 0.29);
41 // Baryon production model
42 //      SetMSTJ(12,3); 
43 // String fragmentation
44 //      SetMSTJ(1,1);  
45 // sea quarks can be used for baryon formation
46 //      SetMSTP(88,2); 
47 // choice of max. virtuality for ISR
48 //      SetMSTP(68,1);
49 // regularisation scheme of ISR 
50 //      SetMSTP(70,2);
51 // all resonance decays switched on
52 //    SetMSTP(41,1);   
53 AliPythia8* AliPythia8::fgAliPythia8=NULL;
54
55 AliPythia8::AliPythia8():
56     AliTPythia8(),
57     AliPythiaBase(),
58     fProcess(kPyMb),
59     fEcms(0.),
60     fStrucFunc(kCTEQ5L),
61     fEtSeed(0.),
62     fMinEtJet(0.),
63     fRJet(0.),
64     fYScale(0.),
65     fPtScale(0.),
66     fNJetMin(0),
67     fNJetMax(0)
68 {
69 // Default Constructor
70 //
71 //  Set random number
72     if (!AliPythiaRndm::GetPythiaRandom()) 
73         AliPythiaRndm::SetPythiaRandom(GetRandom());
74 }
75
76 AliPythia8::AliPythia8(const AliPythia8& pythia):
77     AliTPythia8(), 
78     AliPythiaBase(),
79     fProcess(kPyMb),
80     fEcms(0.),
81     fStrucFunc(kCTEQ5L),
82     fEtSeed(0.),
83     fMinEtJet(0.),
84     fRJet(0.),
85     fYScale(0.),
86     fPtScale(0.),
87     fNJetMin(0),
88     fNJetMax(0)
89 {
90     // Copy Constructor
91     pythia.Copy(*this);
92 }
93
94 void AliPythia8::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
95 {
96 // Initialise the process to generate 
97     if (!AliPythiaRndm::GetPythiaRandom()) 
98       AliPythiaRndm::SetPythiaRandom(GetRandom());
99     
100     fProcess   = process;
101     fEcms      = energy;
102     fStrucFunc = strucfunc;
103 //...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
104     ReadString("111:mayDecay  = off");
105     ReadString("310:mayDecay  = off");
106     ReadString("3122:mayDecay = off");
107     ReadString("3112:mayDecay = off");
108     ReadString("3212:mayDecay = off");
109     ReadString("3222:mayDecay = off");
110     ReadString("3312:mayDecay = off");
111     ReadString("3322:mayDecay = off");
112     ReadString("3334:mayDecay = off");
113     // Select structure function 
114
115     ReadString("PDF:useLHAPDF = on");
116     ReadString(Form("PDF:LHAPDFset = %s", AliStructFuncType::PDFsetName(fStrucFunc).Data()));
117
118     // Particles produced in string fragmentation point directly to either of the two endpoints
119     // of the string (depending in the side they were generated from).
120
121 //    SetMSTU(16,2); // ????
122
123 //
124 // Pythia initialisation for selected processes//
125 //
126     switch (process) 
127     {
128     case kPyOldUEQ2ordered:  //Old underlying events with Q2 ordered QCD processes
129 //        Multiple interactions on.
130         ReadString("PartonLevel:MI = on");
131 // Double Gaussian matter distribution.
132         ReadString("MultipleInteractions:bProfile = 2");
133         ReadString("MultipleInteractions:coreFraction = 0.5");
134         ReadString("MultipleInteractions:coreRadius = 0.4");
135 //  pT0.
136         ReadString("MultipleInteractions:pTmin = 2.0");
137 //  Reference energy for pT0 and energy rescaling pace.
138         ReadString("MultipleInteractions:ecmRef = 1800.");
139         ReadString("MultipleInteractions:ecmPow = 0.25");
140 //  String drawing almost completely minimizes string length.
141 //      SetPARP(85,0.9);
142 //      SetPARP(86,0.95);
143 // ISR and FSR activity.
144 // Q^2 scale of the hard scattering
145         ReadString("SigmaProcess:factorMultFac = 4.");
146 // Lambda_FSR scale.
147 //      SetPARJ(81, 0.29);
148         break;
149     case kPyOldUEQ2ordered2:   
150 // Old underlying events with Q2 ordered QCD processes
151 // Multiple interactions on.
152         ReadString("PartonLevel:MI = on");
153 // Double Gaussian matter distribution.
154         ReadString("MultipleInteractions:bProfile = 2");
155         ReadString("MultipleInteractions:coreFraction = 0.5");
156         ReadString("MultipleInteractions:coreRadius = 0.4");
157 // pT0.
158         ReadString("MultipleInteractions:pTmin = 2.0");
159 // Reference energy for pT0 and energy rescaling pace.
160         ReadString("MultipleInteractions:ecmRef = 1800.");
161         ReadString("MultipleInteractions:ecmPow = 0.16");
162 // String drawing almost completely minimizes string length.
163 //      SetPARP(85,0.9);
164 //      SetPARP(86,0.95);
165 // ISR and FSR activity.
166         ReadString("SigmaProcess:factorMultFac = 4.");
167 // Lambda_FSR scale.
168 //      SetPARJ(81,0.29);       
169         break;
170     case kPyOldPopcorn:  
171 // Old production mechanism: Old Popcorn
172         ReadString("HardQCD:all = on");
173 //      SetMSTJ(12,3); 
174 // (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon
175 //      SetMSTP(88,2); 
176 // (D=1)see can be used to form  baryons (BARYON JUNCTION)
177 //      SetMSTJ(1,1);  
178         AtlasTuning();
179         break;
180     case kPyCharm:
181         ReadString("HardQCD:gg2ccbar = on");
182         ReadString("HardQCD:qqbar2ccbar = on");
183 //  heavy quark masses
184         ReadString("ParticleData:mcRun = 1.2");
185 //
186 //    primordial pT
187         ReadString("BeamRemnants:primordialKT = on");
188         ReadString("BeamRemnants:primordialKTsoft = 0.");
189         ReadString("BeamRemnants:primordialKThard = 1.");
190         ReadString("BeamRemnants:halfScaleForKT = 0.");
191         ReadString("BeamRemnants:halfMassForKT = 0.");
192         break;
193     case kPyBeauty:
194         ReadString("HardQCD:gg2bbbar = on");
195         ReadString("HardQCD:qqbar2bbbar = on");
196         ReadString("ParticleData:mbRun = 4.75");
197         break;
198     case kPyJpsi:
199 // gg->J/Psi g
200         ReadString("Charmonium:gg2QQbar[3S1(1)]g = on");
201         break;
202     case kPyJpsiChi:
203         ReadString("Charmonium:all = on");
204         break;
205     case kPyCharmUnforced:
206 // gq->qg   
207         ReadString("HardQCD:gq2qg = on");
208 // gg->qq
209         ReadString("HardQCD:gg2qq = on");
210 // gg->gg
211         ReadString("HardQCD:gg2gg = on");
212         break;
213     case kPyBeautyUnforced:
214 // gq->qg   
215         ReadString("HardQCD:gq2qg = on");
216 // gg->qq
217         ReadString("HardQCD:gg2qq = on");
218 // gg->gg
219         ReadString("HardQCD:gg2gg = on");
220         break;
221     case kPyMb:
222 // Minimum Bias pp-Collisions
223 //
224 //   
225 //      select Pythia min. bias model
226 // single diffraction AB-->XB
227         ReadString("SoftQCD:minBias = on");
228         ReadString("SoftQCD:singleDiffractive = on");
229         ReadString("SoftQCD:doubleDiffractive = on");
230         AtlasTuning();
231         break;
232     case kPyMbDefault:
233 // Minimum Bias pp-Collisions
234 //
235 //   
236 //      select Pythia min. bias model
237         ReadString("SoftQCD:minBias = on");
238         ReadString("SoftQCD:singleDiffractive = on");
239         ReadString("SoftQCD:doubleDiffractive = on");
240         break;
241     case kPyLhwgMb:
242 // Les Houches Working Group 05 Minimum Bias pp-Collisions: hep-ph/0604120
243 //  -> Pythia 6.3 or above is needed
244 //   
245         ReadString("SoftQCD:minBias = on");
246         ReadString("SoftQCD:singleDiffractive = on");
247         ReadString("SoftQCD:doubleDiffractive = on");
248         ReadString(Form("PDF:LHAPDFset = %s", AliStructFuncType::PDFsetName(kCTEQ6ll).Data()));
249
250 //      SetMSTP(68,1);
251 //      SetMSTP(70,2);
252 //      ReadString("PartonLevel:MI = on");
253 // Double Gaussian matter distribution.
254         ReadString("MultipleInteractions:bProfile = 2");
255         ReadString("MultipleInteractions:coreFraction = 0.5");
256         ReadString("MultipleInteractions:coreRadius = 0.5");
257         ReadString("MultipleInteractions:expPow = 0.16");
258         ReadString("MultipleInteractions:pTmin = 2.3");
259 //      SetMSTP(88,1);
260 //      SetPARP(85,0.9);           // Regulates gluon prod. mechanism
261         break;
262     case kPyMbNonDiffr:
263 // Minimum Bias pp-Collisions
264 //
265 //   
266 //      select Pythia min. bias model
267         ReadString("SoftQCD:minBias = on");
268         AtlasTuning();
269         break;
270     case kPyMbMSEL1:
271         ConfigHeavyFlavor();
272 // Intrinsic <kT^2>
273         ReadString("BeamRemnants:primordialKT = on");
274         ReadString("BeamRemnants:primordialKTsoft = 0.");
275         ReadString("BeamRemnants:primordialKThard = 1.");
276         ReadString("BeamRemnants:halfScaleForKT = 0.");
277         ReadString("BeamRemnants:halfMassForKT = 0.");
278 // Set Q-quark mass
279         ReadString("ParticleData:mcRun = 1.20");
280         ReadString("ParticleData:mbRun = 4.78");
281 // Atlas Tuning
282         AtlasTuning();
283         break;
284     case kPyJets:
285 //
286 //  QCD Jets
287 //
288         ReadString("HardQCD:all = on");
289 //
290 //  Pythia Tune A (CDF)
291 //
292         ReadString("PartonLevel:MI = on");
293         ReadString("MultipleInteractions:pTmin = 2.0");
294         ReadString("MultipleInteractions:pT0Ref = 2.8");
295         ReadString("MultipleInteractions:ecmRef = 1800.");
296         ReadString("MultipleInteractions:expPow = 0.25");
297         ReadString("MultipleInteractions:bProfile = 2");
298         ReadString("MultipleInteractions:coreFraction = 0.16");
299         ReadString("MultipleInteractions:coreRadius = 0.4");
300         ReadString("SigmaProcess:factorMultFac = 2.5");
301 //      SetPARP(85,0.90) ;         // Regulates gluon prod. mechanism
302 //      SetPARP(86,0.95);          // Regulates gluon prod. mechanism
303        break;
304     case kPyDirectGamma:
305         ReadString("PromptPhoton:all = on");
306         break;
307     case kPyCharmPbPbMNR:
308     case kPyD0PbPbMNR:
309     case kPyDPlusPbPbMNR:
310     case kPyDPlusStrangePbPbMNR:
311       // Tuning of Pythia parameters aimed to get a resonable agreement
312       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
313       // c-cbar single inclusive and double differential distributions.
314       // This parameter settings are meant to work with Pb-Pb collisions
315       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
316       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
317       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
318         ConfigHeavyFlavor();
319       // Intrinsic <kT>
320         ReadString("BeamRemnants:primordialKT = on");
321         ReadString("BeamRemnants:primordialKTsoft = 0.");
322         ReadString("BeamRemnants:primordialKThard = 1.304");
323         ReadString("BeamRemnants:halfScaleForKT = 0.");
324         ReadString("BeamRemnants:halfMassForKT = 0.");
325       // Set c-quark mass
326         ReadString("ParticleData:mcRun = 1.20");
327         break;
328     case kPyCharmpPbMNR:
329     case kPyD0pPbMNR:
330     case kPyDPluspPbMNR:
331     case kPyDPlusStrangepPbMNR:
332       // Tuning of Pythia parameters aimed to get a resonable agreement
333       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
334       // c-cbar single inclusive and double differential distributions.
335       // This parameter settings are meant to work with p-Pb collisions
336       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
337       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
338       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
339         ConfigHeavyFlavor();
340       // Intrinsic <kT>
341         ReadString("BeamRemnants:primordialKT = on");
342         ReadString("BeamRemnants:primordialKTsoft = 0.");
343         ReadString("BeamRemnants:primordialKThard = 1.16");
344         ReadString("BeamRemnants:halfScaleForKT = 0.");
345         ReadString("BeamRemnants:halfMassForKT = 0.");
346       // Set c-quark mass
347         ReadString("ParticleData:mcRun = 1.20");
348       break;
349     case kPyCharmppMNR:
350     case kPyD0ppMNR:
351     case kPyDPlusppMNR:
352     case kPyDPlusStrangeppMNR:
353     case kPyLambdacppMNR:
354       // Tuning of Pythia parameters aimed to get a resonable agreement
355       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
356       // c-cbar single inclusive and double differential distributions.
357       // This parameter settings are meant to work with pp collisions
358       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
359       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
360       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
361         ConfigHeavyFlavor();
362       // Intrinsic <kT^2>
363         ReadString("BeamRemnants:primordialKT = on");
364         ReadString("BeamRemnants:primordialKTsoft = 0.");
365         ReadString("BeamRemnants:primordialKThard = 1.");
366         ReadString("BeamRemnants:halfScaleForKT = 0.");
367         ReadString("BeamRemnants:halfMassForKT = 0.");
368       // Set c-quark mass
369         ReadString("ParticleData:mcRun = 1.20");
370       break;
371     case kPyCharmppMNRwmi:
372       // Tuning of Pythia parameters aimed to get a resonable agreement
373       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
374       // c-cbar single inclusive and double differential distributions.
375       // This parameter settings are meant to work with pp collisions
376       // and with kCTEQ5L PDFs.
377       // Added multiple interactions according to ATLAS tune settings.
378       // To get a "reasonable" agreement with MNR results, events have to be 
379       // generated with the minimum ptHard (AliGenPythia::SetPtHard)
380       // set to 2.76 GeV.
381       // To get a "perfect" agreement with MNR results, events have to be 
382       // generated in four ptHard bins with the following relative 
383       // normalizations:
384       // 2.76-3 GeV: 25%
385       //    3-4 GeV: 40%
386       //    4-8 GeV: 29%
387       //     >8 GeV:  6%
388         ConfigHeavyFlavor();
389       // Intrinsic <kT^2>
390         ReadString("BeamRemnants:primordialKT = on");
391         ReadString("BeamRemnants:primordialKTsoft = 0.");
392         ReadString("BeamRemnants:primordialKThard = 1.");
393         ReadString("BeamRemnants:halfScaleForKT = 0.");
394         ReadString("BeamRemnants:halfMassForKT = 0.");
395       // Set c-quark mass
396         ReadString("ParticleData:mcRun = 1.20");
397         AtlasTuning();
398         break;
399     case kPyBeautyPbPbMNR:
400       // Tuning of Pythia parameters aimed to get a resonable agreement
401       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
402       // b-bbar single inclusive and double differential distributions.
403       // This parameter settings are meant to work with Pb-Pb collisions
404       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
405       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
406       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
407         ConfigHeavyFlavor();
408       // QCD scales
409         ReadString("SigmaProcess:factorMultFac = 1.");
410       // Intrinsic <kT>
411         ReadString("BeamRemnants:primordialKT = on");
412         ReadString("BeamRemnants:primordialKTsoft = 0.");
413         ReadString("BeamRemnants:primordialKThard = 2.035");
414         ReadString("BeamRemnants:halfScaleForKT = 0.");
415         ReadString("BeamRemnants:halfMassForKT = 0.");
416       // Set b-quark mass
417         ReadString("ParticleData:mbRun = 4.75");
418       break;
419     case kPyBeautypPbMNR:
420       // Tuning of Pythia parameters aimed to get a resonable agreement
421       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
422       // b-bbar single inclusive and double differential distributions.
423       // This parameter settings are meant to work with p-Pb collisions
424       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
425       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
426       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
427         ConfigHeavyFlavor();
428       // QCD scales
429         ReadString("SigmaProcess:factorMultFac = 1.");
430       // Intrinsic <kT>
431         ReadString("BeamRemnants:primordialKT = on");
432         ReadString("BeamRemnants:primordialKTsoft = 0.");
433         ReadString("BeamRemnants:primordialKThard = 1.6");
434         ReadString("BeamRemnants:halfScaleForKT = 0.");
435         ReadString("BeamRemnants:halfMassForKT = 0.");
436       // Set b-quark mass
437         ReadString("ParticleData:mbRun = 4.75");
438       break;
439     case kPyBeautyppMNR:
440       // Tuning of Pythia parameters aimed to get a resonable agreement
441       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
442       // b-bbar single inclusive and double differential distributions.
443       // This parameter settings are meant to work with pp collisions
444       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
445       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
446       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
447         ConfigHeavyFlavor();
448         // QCD scales
449         ReadString("SigmaProcess:factorMultFac = 1.");
450         // Intrinsic <kT>
451         ReadString("BeamRemnants:primordialKT = on");
452         ReadString("BeamRemnants:primordialKTsoft = 0.");
453         ReadString("BeamRemnants:primordialKThard = 1.0");
454         ReadString("BeamRemnants:halfScaleForKT = 0.");
455         ReadString("BeamRemnants:halfMassForKT = 0.");
456         // Set b-quark mass
457         ReadString("ParticleData:mbRun = 4.75");
458       break;
459      case kPyBeautyppMNRwmi:
460       // Tuning of Pythia parameters aimed to get a resonable agreement
461       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
462       // b-bbar single inclusive and double differential distributions.
463       // This parameter settings are meant to work with pp collisions
464       // and with kCTEQ5L PDFs.
465       // Added multiple interactions according to ATLAS tune settings.
466       // To get a "reasonable" agreement with MNR results, events have to be 
467       // generated with the minimum ptHard (AliGenPythia::SetPtHard)
468       // set to 2.76 GeV.
469       // To get a "perfect" agreement with MNR results, events have to be 
470       // generated in four ptHard bins with the following relative 
471       // normalizations:
472       // 2.76-4 GeV:  5% 
473       //    4-6 GeV: 31%
474       //    6-8 GeV: 28%
475       //     >8 GeV: 36%
476          ConfigHeavyFlavor();
477          // QCD scales
478          ReadString("SigmaProcess:factorMultFac = 1.");
479          // Intrinsic <kT>
480         ReadString("BeamRemnants:primordialKT = on");
481         ReadString("BeamRemnants:primordialKTsoft = 0.");
482         ReadString("BeamRemnants:primordialKThard = 1.0");
483         ReadString("BeamRemnants:halfScaleForKT = 0.");
484         ReadString("BeamRemnants:halfMassForKT = 0.");
485         // Set b-quark mass
486         ReadString("ParticleData:mbRun = 4.75");
487         AtlasTuning();
488         break; 
489     case kPyW:
490         //Inclusive production of W+/-
491         //f fbar -> W+ 
492         ReadString("WeakSingleBoson:ffbar2W = on");
493         // Initial/final parton shower on (Pythia default)
494         // With parton showers on we are generating "W inclusive process"
495         ReadString("PartonLevel:ISR = on");
496         ReadString("PartonLevel:FSR = on");
497         break;  
498     case kPyZ:
499         //Inclusive production of Z
500         //f fbar -> Z/gamma
501         ReadString("WeakSingleBoson:ffbar2gmZ = on");
502         //only Z included, not gamma 
503         ReadString("WeakZ0:gmZmode = 2");
504         // Initial/final parton shower on (Pythia default)
505         // With parton showers on we are generating "Z inclusive process"
506         ReadString("PartonLevel:ISR = on");
507         ReadString("PartonLevel:FSR = on");
508         break;  
509     case kPyMbWithDirectPhoton:
510     case kPyBeautyJets:
511     case kPyMbAtlasTuneMC09: 
512       break;
513     }
514 //
515 //  Initialize PYTHIA
516 //    SetMSTP(41,1);   // all resonance decays switched on
517     Initialize(2212, 2212, fEcms);
518 }
519
520 void AliPythia8::SetNuclei(Int_t /*a1*/, Int_t /*a2*/)
521 {
522 // Treat protons as inside nuclei with mass numbers a1 and a2  
523 //    The MSTP array in the PYPARS common block is used to enable and 
524 //    select the nuclear structure functions. 
525 //    MSTP(52)  : (D=1) choice of proton and nuclear structure-function library
526 //            =1: internal PYTHIA acording to MSTP(51) 
527 //            =2: PDFLIB proton  s.f., with MSTP(51)  = 1000xNGROUP+NSET
528 //    If the following mass number both not equal zero, nuclear corrections of the stf are used.
529 //    MSTP(192) : Mass number of nucleus side 1
530 //    MSTP(193) : Mass number of nucleus side 2
531 //    SetMSTP(52,2);
532 //    SetMSTP(192, a1);
533 //    SetMSTP(193, a2);  
534 }
535         
536
537 AliPythia8* AliPythia8::Instance()
538
539 // Set random number generator 
540     if (fgAliPythia8) {
541         return fgAliPythia8;
542     } else {
543         fgAliPythia8 = new AliPythia8();
544         return fgAliPythia8;
545     }
546 }
547
548 void AliPythia8::PrintParticles()
549
550 // Print list of particl properties
551     ReadString("Main:showAllParticleData");
552 }
553
554 void  AliPythia8::ResetDecayTable()
555 {
556 //  Set default values for pythia decay switches
557 //    Int_t i;
558 //    for (i = 1; i <  501; i++) SetMDCY(i,1,fDefMDCY[i]);
559 //    for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
560 }
561
562 void  AliPythia8::SetDecayTable()
563 {
564 //  Set default values for pythia decay switches
565 //
566 //    Int_t i;
567 //    for (i = 1; i <  501; i++) fDefMDCY[i] = GetMDCY(i,1);
568 //    for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
569 }
570
571 void  AliPythia8::Pyclus(Int_t& njet)
572 {
573 //  Call Pythia clustering algorithm
574 //
575     Bool_t ok = fClusterJet.analyze(Pythia8()->event, fYScale, fPtScale, fNJetMin, fNJetMax);
576     njet = 0;
577     if (ok) njet = fClusterJet.size();
578 }
579
580 void  AliPythia8::Pycell(Int_t& njet)
581 {
582 //  Call Pythia jet reconstruction algorithm
583 //
584     Bool_t ok = fCellJet.analyze(Pythia8()->event, fMinEtJet, fRJet, fEtSeed);
585     njet = 0;
586     if (ok) njet = fCellJet.size();
587 }
588
589 void AliPythia8::GetJet(Int_t i, Float_t& px, Float_t& py, Float_t& pz, Float_t& e)
590 {
591     // Get jet number i
592     Float_t et = fCellJet.eT(i);
593     px = et * TMath::Cos(fCellJet.phiWeighted(i));
594     py = et * TMath::Sin(fCellJet.phiWeighted(i));
595     pz = et * TMath::SinH(fCellJet.etaWeighted(i));
596     e  = et * TMath::CosH(fCellJet.etaWeighted(i));
597 }
598
599 void AliPythia8::GenerateEvent()
600 {
601     // Generate one event
602     AliTPythia8::GenerateEvent();
603 }
604
605 void AliPythia8::GenerateMIEvent()
606 {
607     // New multiple interaction scenario
608     AliWarning("Not implemented. No event will be generated");
609 }
610
611 void AliPythia8::PrintStatistics()
612 {
613     // End of run statistics
614     AliTPythia8::PrintStatistics();
615 }
616
617 void AliPythia8::EventListing()
618 {
619     // End of run statistics
620     AliTPythia8::EventListing();
621 }
622
623 Int_t AliPythia8::ProcessCode()
624 {
625     // Returns the subprocess code for the current event
626     return Pythia8()->info.code();
627 }
628
629 void AliPythia8::ConfigHeavyFlavor()
630 {
631     //
632     // Default configuration for Heavy Flavor production
633     //
634     // All QCD processes
635     //
636     ReadString("HardQCD:all = on");
637
638     // No multiple interactions
639     ReadString("PartonLevel:MI = off");
640     ReadString("MultipleInteractions:pTmin = 0.0");
641     ReadString("MultipleInteractions:pT0Ref = 0.0");
642
643     // Initial/final parton shower on (Pythia default)
644     ReadString("PartonLevel:ISR = on");
645     ReadString("PartonLevel:FSR = on");
646     
647     // 2nd order alpha_s
648     ReadString("SigmaProcess:alphaSorder = 2");
649
650     // QCD scales 
651     ReadString("SigmaProcess:renormScale2 = 2");
652     ReadString("SigmaProcess:renormMultFac = 1.");
653 }
654
655 void AliPythia8::AtlasTuning()
656 {
657     //
658     // Configuration for the ATLAS tuning
659     ReadString(Form("PDF:LHAPDFset = %s", AliStructFuncType::PDFsetName(kCTEQ5L).Data()));
660     ReadString("PartonLevel:MI = on");
661     ReadString("MultipleInteractions:pTmin = 1.9");
662     ReadString("MultipleInteractions:pT0Ref = 1.8");
663     ReadString("MultipleInteractions:ecmRef = 1000.");
664     ReadString("MultipleInteractions:expPow = 0.16");
665     ReadString("MultipleInteractions:bProfile = 2");
666     ReadString("MultipleInteractions:coreFraction = 0.16");
667     ReadString("MultipleInteractions:coreRadius = 0.5");
668 //      SetPARP(85,0.33);          // Regulates gluon prod. mechanism
669 //      SetPARP(86,0.66);          // Regulates gluon prod. mechanism
670     ReadString("SigmaProcess:factorMultFac = 1.");
671 }
672
673 void AliPythia8::SetPtHardRange(Float_t ptmin, Float_t ptmax)
674 {
675     // Set the pt hard range
676     ReadString(Form("PhaseSpace:pTHatMin = %13.3f", ptmin));
677     ReadString(Form("PhaseSpace:pTHatMax = %13.3f", ptmax));
678 }
679
680 void AliPythia8::SetYHardRange(Float_t /*ymin*/, Float_t /*ymax*/)
681 {
682     // Set the y hard range
683     printf("YHardRange not implemented in Pythia8 !!!\n");
684     
685 }
686
687
688 void AliPythia8::SetFragmentation(Int_t flag)
689 {
690     // Switch fragmentation on/off
691     if (flag) {
692         ReadString("HadronLevel:Hadronize = on");
693     } else {
694         ReadString("HadronLevel:Hadronize = off");
695     }
696 }
697
698 void AliPythia8::SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2)
699 {
700 //  initial state radiation    
701     if (flag1) {
702         ReadString("PartonLevel:ISR = on");
703     } else {
704         ReadString("PartonLevel:ISR = off");
705     }
706 // final state radiation    
707     if (flag2) {
708         ReadString("PartonLevel:FSR = on");
709     } else {
710         ReadString("PartonLevel:FSR = off");
711     }
712 }
713
714 void AliPythia8::SetIntrinsicKt(Float_t kt)
715 {
716         ReadString("BeamRemnants:primordialKT = on");
717         ReadString("BeamRemnants:primordialKTsoft = 0.");
718         ReadString(Form("BeamRemnants:primordialKThard = %13.3f", kt));
719         ReadString("BeamRemnants:halfScaleForKT = 0.");
720         ReadString("BeamRemnants:halfMassForKT = 0.");
721 }
722
723 void AliPythia8::SwitchHFOff()
724 {
725     // Switch off heavy flavor
726     // Maximum number of quark flavours used in pdf 
727     ReadString("PDFinProcess:nQuarkIn = 3");
728     // Maximum number of flavors that can be used in showers
729     ReadString("TimeShower:nGluonToQuark = 3");
730     ReadString("SpaceShower:nQuarkIn = 3");
731     
732
733 }
734
735 void AliPythia8::SetPycellParameters(Float_t etaMax, Int_t nEta, Int_t nPhi,
736                                        Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
737 {
738 // Set pycell parameters
739     fCellJet  = Pythia8::CellJet( etaMax, nEta, nPhi, 2, 0, 0., 0., thresh);
740     fEtSeed   = etseed;
741     fMinEtJet = minet;
742     fRJet     = r;
743 }
744
745 void AliPythia8::ModifiedSplitting()
746 {
747 //
748 // We have to see how to implement this in Pythia8 !!!
749 //
750     // Modified splitting probability as a model for quenching
751 //    SetPARJ(200, 0.8);
752 //    SetMSTJ(41, 1);  // QCD radiation only
753 //    SetMSTJ(42, 2);  // angular ordering
754 //    SetMSTJ(44, 2);  // option to run alpha_s
755 //    SetMSTJ(47, 0);  // No correction back to hard scattering element
756 //    SetMSTJ(50, 0);  // No coherence in first branching
757 //    SetPARJ(82, 1.); // Cut off for parton showers
758 }
759
760     
761 void AliPythia8::InitQuenching(Float_t /*cMin*/, Float_t /*cMax*/, Float_t /*k*/, Int_t /*iECMethod*/, Float_t /*zmax*/, Int_t /*ngmax*/)
762 {
763     //
764     //
765     AliWarning("Not implemented !");
766 }
767
768 void AliPythia8::SwitchHadronisationOff()
769 {
770     // Switch off hadronisation
771     ReadString("HadronLevel:Hadronize = off");
772 }
773
774 void AliPythia8::SwitchHadronisationOn()
775 {
776     // Switch on hadronisarion
777     ReadString("HadronLevel:Hadronize = on");
778 }
779
780
781 void AliPythia8::GetXandQ(Float_t& x1, Float_t& x2, Float_t& q)
782 {
783     // Get x1, x2 and Q for this event
784     
785     q  = Pythia8()->info.QFac();
786     x1 = Pythia8()->info.x1();
787     x2 = Pythia8()->info.x2();
788     
789 }
790
791 Float_t AliPythia8::GetXSection()
792 {
793     // Get the total cross-section
794     return Pythia8()->info.sigmaGen();
795 }
796
797 Float_t AliPythia8::GetPtHard()
798 {
799     // Get the pT hard for this event
800     return Pythia8()->info.pTHat();
801 }
802
803
804
805
806 AliPythia8& AliPythia8::operator=(const  AliPythia8& rhs)
807 {
808 // Assignment operator
809     rhs.Copy(*this);
810     return *this;
811 }
812
813  void AliPythia8::Copy(TObject&) const
814 {
815     //
816     // Copy 
817     //
818     Fatal("Copy","Not implemented!\n");
819 }
820
821 //
822 // To be implemented
823 //
824 void AliPythia8::SetNumberOfParticles(Int_t /*i*/)
825 {
826     AliWarning("Not implemented");
827 }
828
829 void AliPythia8::EditEventList(Int_t /*i*/)
830 {
831     AliWarning("Not implemented");
832 }
833
834 void AliPythia8::Pyquen(Double_t /*a*/, Int_t /*b*/, Double_t /*c*/)
835 {
836     AliWarning("Cannot be used with Pythia8");
837 }
838
839 void AliPythia8::HadronizeEvent()
840 {
841     // Needs access to HadronLevel ?
842     AliWarning("Not yet implemented");
843 }
844
845 void AliPythia8::GetQuenchingParameters(Double_t& /*xp*/, Double_t& /*yp*/, Double_t* /*z[4]*/)
846 {
847     AliWarning("Not yet implemented");
848 }
849
850 void AliPythia8::LoadEvent(AliStack* /*stack*/, Int_t /*flag*/, Int_t /*reHadr*/)
851 {
852     AliWarning("Not yet implemented");
853 }