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