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