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