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