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