Al layer: positioning correction
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia.cxx
CommitLineData
8d2cd130 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
7cdba479 16/* $Id$ */
8d2cd130 17
18#include "AliPythia.h"
7cdba479 19#include "AliPythiaRndm.h"
0f482ae4 20#include "../FASTSIM/AliFastGlauber.h"
21#include "../FASTSIM/AliQuenchingWeights.h"
22#include "TVector3.h"
8d2cd130 23
24ClassImp(AliPythia)
25
26#ifndef WIN32
27# define pyclus pyclus_
28# define pycell pycell_
452af8c7 29# define pyshow pyshow_
30# define pyrobo pyrobo_
992f2843 31# define pyquen pyquen_
16a82508 32# define pyevnw pyevnw_
8d2cd130 33# define type_of_call
34#else
35# define pyclus PYCLUS
36# define pycell PYCELL
452af8c7 37# define pyrobo PYROBO
992f2843 38# define pyquen PYQUEN
16a82508 39# define pyevnw PYEVNW
8d2cd130 40# define type_of_call _stdcall
41#endif
42
43extern "C" void type_of_call pyclus(Int_t & );
44extern "C" void type_of_call pycell(Int_t & );
452af8c7 45extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
46extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
992f2843 47extern "C" void type_of_call pyquen(Double_t &, Int_t &, Double_t &);
0a2cfc0a 48extern "C" void type_of_call pyevnw(){;}
8d2cd130 49
50//_____________________________________________________________________________
51
52AliPythia* AliPythia::fgAliPythia=NULL;
53
54AliPythia::AliPythia()
55{
56// Default Constructor
57//
58// Set random number
7cdba479 59 if (!AliPythiaRndm::GetPythiaRandom())
60 AliPythiaRndm::SetPythiaRandom(GetRandom());
0f482ae4 61 fGlauber = 0;
62 fQuenchingWeights = 0;
8d2cd130 63}
64
65void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
66{
67// Initialise the process to generate
7cdba479 68 if (!AliPythiaRndm::GetPythiaRandom())
69 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 70
71 fProcess = process;
72 fEcms = energy;
73 fStrucFunc = strucfunc;
1d5b1b20 74//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
75 SetMDCY(Pycomp(111) ,1,0);
76 SetMDCY(Pycomp(310) ,1,0);
77 SetMDCY(Pycomp(3122),1,0);
78 SetMDCY(Pycomp(3112),1,0);
79 SetMDCY(Pycomp(3212),1,0);
80 SetMDCY(Pycomp(3222),1,0);
81 SetMDCY(Pycomp(3312),1,0);
82 SetMDCY(Pycomp(3322),1,0);
83 SetMDCY(Pycomp(3334),1,0);
84 // select structure function
8d2cd130 85 SetMSTP(52,2);
86 SetMSTP(51,strucfunc);
87//
88// Pythia initialisation for selected processes//
89//
90// Make MSEL clean
91//
92 for (Int_t i=1; i<= 200; i++) {
93 SetMSUB(i,0);
94 }
95// select charm production
96 switch (process)
97 {
65f2626c 98 case kPyOldUEQ2ordered: //Old underlying events with Q2 ordered QCD processes
99// Multiple interactions on.
100 SetMSTP(81,1);
101// Double Gaussian matter distribution.
102 SetMSTP(82,4);
103 SetPARP(83,0.5);
104 SetPARP(84,0.4);
105// pT0.
106 SetPARP(82,2.0);
107// Reference energy for pT0 and energy rescaling pace.
108 SetPARP(89,1800);
109 SetPARP(90,0.25);
110// String drawing almost completely minimizes string length.
111 SetPARP(85,0.9);
112 SetPARP(86,0.95);
113// ISR and FSR activity.
114 SetPARP(67,4);
115 SetPARP(71,4);
116// Lambda_FSR scale.
117 SetPARJ(81,0.29);
118 break;
119 case kPyOldUEQ2ordered2:
120// Old underlying events with Q2 ordered QCD processes
121// Multiple interactions on.
122 SetMSTP(81,1);
123// Double Gaussian matter distribution.
124 SetMSTP(82,4);
125 SetPARP(83,0.5);
126 SetPARP(84,0.4);
127// pT0.
128 SetPARP(82,2.0);
129// Reference energy for pT0 and energy rescaling pace.
130 SetPARP(89,1800);
131 SetPARP(90,0.16); // here is the difference with kPyOldUEQ2ordered
132// String drawing almost completely minimizes string length.
133 SetPARP(85,0.9);
134 SetPARP(86,0.95);
135// ISR and FSR activity.
136 SetPARP(67,4);
137 SetPARP(71,4);
138// Lambda_FSR scale.
139 SetPARJ(81,0.29);
140 break;
141 case kPyOldPopcorn:
142// Old production mechanism: Old Popcorn
143 SetMSEL(1);
144 SetMSTJ(12,3);
145// (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon
146 SetMSTP(88,2);
147// (D=1)see can be used to form baryons (BARYON JUNCTION)
148 SetMSTJ(1,1);
e0e89f40 149 AtlasTuning();
65f2626c 150 break;
8d2cd130 151 case kPyCharm:
152 SetMSEL(4);
8d2cd130 153// heavy quark masses
154
155 SetPMAS(4,1,1.2);
156 SetMSTU(16,2);
157//
158// primordial pT
159 SetMSTP(91,1);
160 SetPARP(91,1.);
161 SetPARP(93,5.);
162//
163 break;
164 case kPyBeauty:
165 SetMSEL(5);
166 SetPMAS(5,1,4.75);
167 SetMSTU(16,2);
168 break;
169 case kPyJpsi:
170 SetMSEL(0);
171// gg->J/Psi g
172 SetMSUB(86,1);
173 break;
174 case kPyJpsiChi:
175 SetMSEL(0);
176// gg->J/Psi g
177 SetMSUB(86,1);
178// gg-> chi_0c g
179 SetMSUB(87,1);
180// gg-> chi_1c g
181 SetMSUB(88,1);
182// gg-> chi_2c g
183 SetMSUB(89,1);
184 break;
185 case kPyCharmUnforced:
186 SetMSEL(0);
187// gq->qg
188 SetMSUB(28,1);
189// gg->qq
190 SetMSUB(53,1);
191// gg->gg
192 SetMSUB(68,1);
193 break;
194 case kPyBeautyUnforced:
195 SetMSEL(0);
196// gq->qg
197 SetMSUB(28,1);
198// gg->qq
199 SetMSUB(53,1);
200// gg->gg
201 SetMSUB(68,1);
202 break;
203 case kPyMb:
204// Minimum Bias pp-Collisions
205//
206//
207// select Pythia min. bias model
208 SetMSEL(0);
511db649 209 SetMSUB(92,1); // single diffraction AB-->XB
210 SetMSUB(93,1); // single diffraction AB-->AX
211 SetMSUB(94,1); // double diffraction
212 SetMSUB(95,1); // low pt production
213
e0e89f40 214 AtlasTuning();
511db649 215 break;
8d2cd130 216 case kPyMbNonDiffr:
217// Minimum Bias pp-Collisions
218//
219//
220// select Pythia min. bias model
221 SetMSEL(0);
511db649 222 SetMSUB(95,1); // low pt production
0f482ae4 223
e0e89f40 224 AtlasTuning();
8d2cd130 225 break;
226 case kPyJets:
227//
228// QCD Jets
229//
230 SetMSEL(1);
65f2626c 231 // Pythia Tune A (CDF)
232 //
233 SetPARP(67,4.); // Regulates Initial State Radiation
234 SetMSTP(82,4); // Double Gaussian Model
235 SetPARP(82,2.0); // [GeV] PT_min at Ref. energy
236 SetPARP(84,0.4); // Core radius
237 SetPARP(85,0.90) ; // Regulates gluon prod. mechanism
238 SetPARP(86,0.95); // Regulates gluon prod. mechanism
239 SetPARP(89,1800.); // [GeV] Ref. energy
240 SetPARP(90,0.25); // 2*epsilon (exponent in power law)
241 break;
8d2cd130 242 case kPyDirectGamma:
243 SetMSEL(10);
244 break;
adf4d898 245 case kPyCharmPbPbMNR:
246 case kPyD0PbPbMNR:
90d7b703 247 case kPyDPlusPbPbMNR:
e0e89f40 248 case kPyDPlusStrangePbPbMNR:
90d7b703 249 // Tuning of Pythia parameters aimed to get a resonable agreement
250 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
251 // c-cbar single inclusive and double differential distributions.
252 // This parameter settings are meant to work with Pb-Pb collisions
253 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
254 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
255 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
3dc3ec94 256 ConfigHeavyFlavor();
90d7b703 257 // Intrinsic <kT>
258 SetMSTP(91,1);
259 SetPARP(91,1.304);
260 SetPARP(93,6.52);
90d7b703 261 // Set c-quark mass
262 SetPMAS(4,1,1.2);
8d2cd130 263 break;
adf4d898 264 case kPyCharmpPbMNR:
265 case kPyD0pPbMNR:
90d7b703 266 case kPyDPluspPbMNR:
e0e89f40 267 case kPyDPlusStrangepPbMNR:
90d7b703 268 // Tuning of Pythia parameters aimed to get a resonable agreement
269 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
270 // c-cbar single inclusive and double differential distributions.
271 // This parameter settings are meant to work with p-Pb collisions
272 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
273 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
274 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
3dc3ec94 275 ConfigHeavyFlavor();
90d7b703 276 // Intrinsic <kT>
3dc3ec94 277 SetMSTP(91,1);
278 SetPARP(91,1.16);
279 SetPARP(93,5.8);
280
90d7b703 281 // Set c-quark mass
3dc3ec94 282 SetPMAS(4,1,1.2);
adf4d898 283 break;
284 case kPyCharmppMNR:
285 case kPyD0ppMNR:
90d7b703 286 case kPyDPlusppMNR:
e0e89f40 287 case kPyDPlusStrangeppMNR:
90d7b703 288 // Tuning of Pythia parameters aimed to get a resonable agreement
289 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
290 // c-cbar single inclusive and double differential distributions.
291 // This parameter settings are meant to work with pp collisions
292 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
293 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
294 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
3dc3ec94 295 ConfigHeavyFlavor();
90d7b703 296 // Intrinsic <kT^2>
3dc3ec94 297 SetMSTP(91,1);
298 SetPARP(91,1.);
299 SetPARP(93,5.);
300
90d7b703 301 // Set c-quark mass
3dc3ec94 302 SetPMAS(4,1,1.2);
adf4d898 303 break;
e0e89f40 304 case kPyCharmppMNRwmi:
305 // Tuning of Pythia parameters aimed to get a resonable agreement
306 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
307 // c-cbar single inclusive and double differential distributions.
308 // This parameter settings are meant to work with pp collisions
309 // and with kCTEQ5L PDFs.
310 // Added multiple interactions according to ATLAS tune settings.
311 // To get a "reasonable" agreement with MNR results, events have to be
312 // generated with the minimum ptHard (AliGenPythia::SetPtHard)
313 // set to 2.76 GeV.
314 // To get a "perfect" agreement with MNR results, events have to be
315 // generated in four ptHard bins with the following relative
316 // normalizations:
317 // 2.76-3 GeV: 25%
318 // 3-4 GeV: 40%
319 // 4-8 GeV: 29%
320 // >8 GeV: 6%
321 ConfigHeavyFlavor();
322 // Intrinsic <kT^2>
323 SetMSTP(91,1);
324 SetPARP(91,1.);
325 SetPARP(93,5.);
326
327 // Set c-quark mass
328 SetPMAS(4,1,1.2);
329 AtlasTuning();
330 break;
adf4d898 331 case kPyBeautyPbPbMNR:
8d2cd130 332 // Tuning of Pythia parameters aimed to get a resonable agreement
333 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
334 // b-bbar single inclusive and double differential distributions.
335 // This parameter settings are meant to work with Pb-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.75GeV. Example in ConfigBeautyPPR.C.
3dc3ec94 339 ConfigHeavyFlavor();
8d2cd130 340 // QCD scales
3dc3ec94 341 SetPARP(67,1.0);
342 SetPARP(71,1.0);
adf4d898 343 // Intrinsic <kT>
3dc3ec94 344 SetMSTP(91,1);
345 SetPARP(91,2.035);
346 SetPARP(93,10.17);
8d2cd130 347 // Set b-quark mass
3dc3ec94 348 SetPMAS(5,1,4.75);
adf4d898 349 break;
350 case kPyBeautypPbMNR:
351 // Tuning of Pythia parameters aimed to get a resonable agreement
352 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
353 // b-bbar single inclusive and double differential distributions.
354 // This parameter settings are meant to work with p-Pb collisions
355 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
356 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
357 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
3dc3ec94 358 ConfigHeavyFlavor();
adf4d898 359 // QCD scales
3dc3ec94 360 SetPARP(67,1.0);
361 SetPARP(71,1.0);
adf4d898 362 // Intrinsic <kT>
3dc3ec94 363 SetMSTP(91,1);
364 SetPARP(91,1.60);
365 SetPARP(93,8.00);
adf4d898 366 // Set b-quark mass
3dc3ec94 367 SetPMAS(5,1,4.75);
adf4d898 368 break;
369 case kPyBeautyppMNR:
370 // Tuning of Pythia parameters aimed to get a resonable agreement
371 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
372 // b-bbar single inclusive and double differential distributions.
373 // This parameter settings are meant to work with pp collisions
374 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
375 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
376 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
3dc3ec94 377 ConfigHeavyFlavor();
adf4d898 378 // QCD scales
3dc3ec94 379 SetPARP(67,1.0);
380 SetPARP(71,1.0);
381
382 // Intrinsic <kT>
383 SetMSTP(91,1);
384 SetPARP(91,1.);
385 SetPARP(93,5.);
386
387 // Set b-quark mass
388 SetPMAS(5,1,4.75);
8d2cd130 389 break;
e0e89f40 390 case kPyBeautyppMNRwmi:
391 // Tuning of Pythia parameters aimed to get a resonable agreement
392 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
393 // b-bbar single inclusive and double differential distributions.
394 // This parameter settings are meant to work with pp collisions
395 // and with kCTEQ5L PDFs.
396 // Added multiple interactions according to ATLAS tune settings.
397 // To get a "reasonable" agreement with MNR results, events have to be
398 // generated with the minimum ptHard (AliGenPythia::SetPtHard)
399 // set to 2.76 GeV.
400 // To get a "perfect" agreement with MNR results, events have to be
401 // generated in four ptHard bins with the following relative
402 // normalizations:
403 // 2.76-4 GeV: 5%
404 // 4-6 GeV: 31%
405 // 6-8 GeV: 28%
406 // >8 GeV: 36%
407 ConfigHeavyFlavor();
408 // QCD scales
409 SetPARP(67,1.0);
410 SetPARP(71,1.0);
411
412 // Intrinsic <kT>
413 SetMSTP(91,1);
414 SetPARP(91,1.);
415 SetPARP(93,5.);
416
417 // Set b-quark mass
418 SetPMAS(5,1,4.75);
419
420 AtlasTuning();
421 break;
589380c6 422 case kPyW:
423
424 //Inclusive production of W+/-
425 SetMSEL(0);
426 //f fbar -> W+
427 SetMSUB(2,1);
428 // //f fbar -> g W+
429 // SetMSUB(16,1);
430 // //f fbar -> gamma W+
431 // SetMSUB(20,1);
432 // //f g -> f W+
433 // SetMSUB(31,1);
434 // //f gamma -> f W+
435 // SetMSUB(36,1);
436
437 // Initial/final parton shower on (Pythia default)
438 // With parton showers on we are generating "W inclusive process"
439 SetMSTP(61,1); //Initial QCD & QED showers on
440 SetMSTP(71,1); //Final QCD & QED showers on
441
442 break;
0f6ee828 443
444 case kPyZ:
445
446 //Inclusive production of Z
447 SetMSEL(0);
448 //f fbar -> Z/gamma
449 SetMSUB(1,1);
450
451 // // f fbar -> g Z/gamma
452 // SetMSUB(15,1);
453 // // f fbar -> gamma Z/gamma
454 // SetMSUB(19,1);
455 // // f g -> f Z/gamma
456 // SetMSUB(30,1);
457 // // f gamma -> f Z/gamma
458 // SetMSUB(35,1);
459
460 //only Z included, not gamma
461 SetMSTP(43,2);
462
463 // Initial/final parton shower on (Pythia default)
464 // With parton showers on we are generating "Z inclusive process"
465 SetMSTP(61,1); //Initial QCD & QED showers on
466 SetMSTP(71,1); //Final QCD & QED showers on
467
468 break;
469
8d2cd130 470 }
471//
472// Initialize PYTHIA
473 SetMSTP(41,1); // all resonance decays switched on
474
475 Initialize("CMS","p","p",fEcms);
476
477}
478
479Int_t AliPythia::CheckedLuComp(Int_t kf)
480{
481// Check Lund particle code (for debugging)
482 Int_t kc=Pycomp(kf);
483 printf("\n Lucomp kf,kc %d %d",kf,kc);
484 return kc;
485}
486
487void AliPythia::SetNuclei(Int_t a1, Int_t a2)
488{
489// Treat protons as inside nuclei with mass numbers a1 and a2
490// The MSTP array in the PYPARS common block is used to enable and
491// select the nuclear structure functions.
492// MSTP(52) : (D=1) choice of proton and nuclear structure-function library
493// =1: internal PYTHIA acording to MSTP(51)
494// =2: PDFLIB proton s.f., with MSTP(51) = 1000xNGROUP+NSET
495// If the following mass number both not equal zero, nuclear corrections of the stf are used.
496// MSTP(192) : Mass number of nucleus side 1
497// MSTP(193) : Mass number of nucleus side 2
498 SetMSTP(52,2);
499 SetMSTP(192, a1);
500 SetMSTP(193, a2);
501}
502
503
504AliPythia* AliPythia::Instance()
505{
506// Set random number generator
507 if (fgAliPythia) {
508 return fgAliPythia;
509 } else {
510 fgAliPythia = new AliPythia();
511 return fgAliPythia;
512 }
513}
514
515void AliPythia::PrintParticles()
516{
517// Print list of particl properties
518 Int_t np = 0;
c31f1d37 519 char* name = new char[16];
8d2cd130 520 for (Int_t kf=0; kf<1000000; kf++) {
521 for (Int_t c = 1; c > -2; c-=2) {
8d2cd130 522 Int_t kc = Pycomp(c*kf);
523 if (kc) {
524 Float_t mass = GetPMAS(kc,1);
525 Float_t width = GetPMAS(kc,2);
526 Float_t tau = GetPMAS(kc,4);
c31f1d37 527
8d2cd130 528 Pyname(kf,name);
529
530 np++;
531
532 printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
533 c*kf, name, mass, width, tau);
534 }
535 }
536 }
537 printf("\n Number of particles %d \n \n", np);
538}
539
540void AliPythia::ResetDecayTable()
541{
542// Set default values for pythia decay switches
543 Int_t i;
544 for (i = 1; i < 501; i++) SetMDCY(i,1,fDefMDCY[i]);
545 for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
546}
547
548void AliPythia::SetDecayTable()
549{
550// Set default values for pythia decay switches
551//
552 Int_t i;
553 for (i = 1; i < 501; i++) fDefMDCY[i] = GetMDCY(i,1);
554 for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
555}
556
557void AliPythia::Pyclus(Int_t& njet)
558{
559// Call Pythia clustering algorithm
560//
561 pyclus(njet);
562}
563
564void AliPythia::Pycell(Int_t& njet)
565{
566// Call Pythia jet reconstruction algorithm
567//
568 pycell(njet);
569}
570
452af8c7 571void AliPythia::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
572{
573// Call Pythia jet reconstruction algorithm
574//
452af8c7 575 pyshow(ip1, ip2, qmax);
576}
577
578void AliPythia::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
579{
580 pyrobo(imi, ima, the, phi, bex, bey, bez);
581}
582
583
584
86b6ad68 585void AliPythia::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECMethod)
0f482ae4 586{
587// Initializes
588// (1) The quenching model using quenching weights according to C. Salgado and U. Wiedemann
589// (2) The nuclear geometry using the Glauber Model
590//
591
592
593 fGlauber = new AliFastGlauber();
594 fGlauber->Init(2);
595 fGlauber->SetCentralityClass(cMin, cMax);
596
597 fQuenchingWeights = new AliQuenchingWeights();
598 fQuenchingWeights->InitMult();
86b6ad68 599 fQuenchingWeights->SetK(k);
0f482ae4 600 fQuenchingWeights->SetECMethod(AliQuenchingWeights::kECMethod(iECMethod));
0f482ae4 601}
602
603
452af8c7 604void AliPythia::Quench()
605{
606//
607//
608// Simple Jet Quenching routine:
609// =============================
610// The jet formed by all final state partons radiated by the parton created
0f482ae4 611// in the hard collisions is quenched by a factor (1-z) using light cone variables in
612// the initial parton reference frame:
452af8c7 613// (E + p_z)new = (1-z) (E + p_z)old
614//
0f482ae4 615//
616//
617//
452af8c7 618// The lost momentum is first balanced by one gluon with virtuality > 0.
619// Subsequently the gluon splits to yield two gluons with E = p.
620//
0f482ae4 621//
622//
4e383037 623 static Float_t eMean = 0.;
624 static Int_t icall = 0;
0f482ae4 625
c2c598a3 626 Double_t p0[4][5];
627 Double_t p1[4][5];
628 Double_t p2[4][5];
629 Int_t klast[4] = {-1, -1, -1, -1};
452af8c7 630
631 Int_t numpart = fPyjets->N;
86b6ad68 632 Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0., phi = 0.;
c2c598a3 633 Double_t pxq[4], pyq[4], pzq[4], eq[4], yq[4], mq[4], pq[4], phiq[4], thetaq[4], ptq[4];
634 Bool_t quenched[4];
b280c4cc 635 Double_t wjtKick[4];
c2c598a3 636 Int_t nGluon[4];
86b6ad68 637 Int_t qPdg[4];
0f482ae4 638 Int_t imo, kst, pdg;
b280c4cc 639
511db649 640//
c2c598a3 641// Sore information about Primary partons
642//
643// j =
644// 0, 1 partons from hard scattering
645// 2, 3 partons from initial state radiation
646//
647 for (Int_t i = 2; i <= 7; i++) {
648 Int_t j = 0;
649 // Skip gluons that participate in hard scattering
650 if (i == 4 || i == 5) continue;
651 // Gluons from hard Scattering
652 if (i == 6 || i == 7) {
653 j = i - 6;
654 pxq[j] = fPyjets->P[0][i];
655 pyq[j] = fPyjets->P[1][i];
656 pzq[j] = fPyjets->P[2][i];
657 eq[j] = fPyjets->P[3][i];
658 mq[j] = fPyjets->P[4][i];
659 } else {
660 // Gluons from initial state radiation
661 //
662 // Obtain 4-momentum vector from difference between original parton and parton after gluon
663 // radiation. Energy is calculated independently because initial state radition does not
664 // conserve strictly momentum and energy for each partonic system independently.
665 //
666 // Not very clean. Should be improved !
667 //
668 //
669 j = i;
670 pxq[j] = fPyjets->P[0][i] - fPyjets->P[0][i+2];
671 pyq[j] = fPyjets->P[1][i] - fPyjets->P[1][i+2];
672 pzq[j] = fPyjets->P[2][i] - fPyjets->P[2][i+2];
673 mq[j] = fPyjets->P[4][i];
674 eq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j] + mq[j] * mq[j]);
675 }
676//
677// Calculate some kinematic variables
511db649 678//
4e383037 679 yq[j] = 0.5 * TMath::Log((eq[j] + pzq[j] + 1.e-14) / (eq[j] - pzq[j] + 1.e-14));
0f482ae4 680 pq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j]);
681 phiq[j] = TMath::Pi()+TMath::ATan2(-pyq[j], -pxq[j]);
682 ptq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j]);
683 thetaq[j] = TMath::ATan2(ptq[j], pzq[j]);
86b6ad68 684 qPdg[j] = fPyjets->K[1][i];
685 }
686
687 Double_t int0[4];
688 Double_t int1[4];
86b6ad68 689
b280c4cc 690 fGlauber->GetI0I1ForPythiaAndXY(4, phiq, int0, int1, fXJet, fYJet, 15.);
691
86b6ad68 692 for (Int_t j = 0; j < 4; j++) {
c2c598a3 693 //
694 // Quench only central jets and with E > 10.
695 //
86b6ad68 696
697
698 Int_t itype = (qPdg[j] == 21) ? 2 : 1;
699 Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]);
700
c2c598a3 701 if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
b280c4cc 702 fZQuench[j] = 0.;
0f482ae4 703 } else {
c2c598a3 704 if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) {
4e383037 705 icall ++;
706 eMean += eloss;
707 }
0f482ae4 708 //
709 // Extra pt
86b6ad68 710 Double_t l = fQuenchingWeights->CalcLk(int0[j], int1[j]);
711 wjtKick[j] = TMath::Sqrt(l * fQuenchingWeights->CalcQk(int0[j], int1[j]));
0f482ae4 712 //
713 // Fractional energy loss
b280c4cc 714 fZQuench[j] = eloss / eq[j];
0f482ae4 715 //
716 // Avoid complete loss
717 //
b280c4cc 718 if (fZQuench[j] == 1.) fZQuench[j] = 0.95;
0f482ae4 719 //
720 // Some debug printing
86b6ad68 721
722
bf9bb016 723// printf("Initial parton # %3d, Type %3d Energy %10.3f Phi %10.3f Length %10.3f Loss %10.3f Kick %10.3f Mean: %10.3f %10.3f\n",
724// j, itype, eq[j], phiq[j], l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]);
4e383037 725
b280c4cc 726// fZQuench[j] = 0.8;
727// while (fZQuench[j] >= 0.95) fZQuench[j] = gRandom->Exp(0.2);
0f482ae4 728 }
4e383037 729
b280c4cc 730 quenched[j] = (fZQuench[j] > 0.01);
4e383037 731 } // primary partons
c2c598a3 732
b280c4cc 733
734
6e90ad26 735 Double_t pNew[1000][4];
736 Int_t kNew[1000];
737 Int_t icount = 0;
b280c4cc 738 Double_t zquench[4];
739
6e90ad26 740//
4e383037 741// System Loop
c2c598a3 742 for (Int_t isys = 0; isys < 4; isys++) {
6e90ad26 743// Skip to next system if not quenched.
4e383037 744 if (!quenched[isys]) continue;
745
b280c4cc 746 nGluon[isys] = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
4e383037 747 if (nGluon[isys] > 6) nGluon[isys] = 6;
b280c4cc 748 zquench[isys] = 1. - TMath::Power(1. - fZQuench[isys], 1./Double_t(nGluon[isys]));
4e383037 749 wjtKick[isys] = wjtKick[isys] / TMath::Sqrt(Double_t(nGluon[isys]));
0f482ae4 750
4e383037 751
752
753 Int_t igMin = -1;
754 Int_t igMax = -1;
755 Double_t pg[4] = {0., 0., 0., 0.};
756
757//
758// Loop on radiation events
759
760 for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) {
6e90ad26 761 while (1) {
762 icount = 0;
763 for (Int_t k = 0; k < 4; k++)
764 {
765 p0[isys][k] = 0.;
766 p1[isys][k] = 0.;
767 p2[isys][k] = 0.;
768 }
769// Loop over partons
770 for (Int_t i = 0; i < numpart; i++)
771 {
772 imo = fPyjets->K[2][i];
773 kst = fPyjets->K[0][i];
774 pdg = fPyjets->K[1][i];
775
776
777
0f482ae4 778// Quarks and gluons only
6e90ad26 779 if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
0f482ae4 780// Particles from hard scattering only
c2c598a3 781
6e90ad26 782 if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
c2c598a3 783 Int_t imom = imo % 1000;
784 if ((isys == 0 || isys == 1) && ((imom != (isys + 7)))) continue;
785 if ((isys == 2 || isys == 3) && ((imom != (isys + 1)))) continue;
786
6e90ad26 787
0f482ae4 788// Skip comment lines
6e90ad26 789 if (kst != 1 && kst != 2) continue;
0f482ae4 790//
791// Parton kinematic
6e90ad26 792 px = fPyjets->P[0][i];
793 py = fPyjets->P[1][i];
794 pz = fPyjets->P[2][i];
795 e = fPyjets->P[3][i];
796 m = fPyjets->P[4][i];
797 pt = TMath::Sqrt(px * px + py * py);
798 p = TMath::Sqrt(px * px + py * py + pz * pz);
799 phi = TMath::Pi() + TMath::ATan2(-py, -px);
800 theta = TMath::ATan2(pt, pz);
801
0f482ae4 802//
c2c598a3 803// Save 4-momentum sum for balancing
804 Int_t index = isys;
6e90ad26 805
806 p0[index][0] += px;
807 p0[index][1] += py;
808 p0[index][2] += pz;
809 p0[index][3] += e;
6e90ad26 810
811 klast[index] = i;
812
0f482ae4 813//
814// Fractional energy loss
b280c4cc 815 Double_t z = zquench[index];
4e383037 816
c2c598a3 817
4e383037 818// Don't fully quench radiated gluons
819//
820 if (imo > 1000) {
821// This small factor makes sure that the gluons are not too close in phase space to avoid recombination
822//
823
c2c598a3 824 z = 0.02;
4e383037 825 }
c2c598a3 826// printf("z: %d %f\n", imo, z);
827
4e383037 828
829//
6e90ad26 830
831 //
832 //
833 // Transform into frame in which initial parton is along z-axis
834 //
835 TVector3 v(px, py, pz);
836 v.RotateZ(-phiq[index]); v.RotateY(-thetaq[index]);
837 Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pl = v.Z();
838
839 Double_t jt = TMath::Sqrt(pxs * pxs + pys * pys);
840 Double_t mt2 = jt * jt + m * m;
841 Double_t zmax = 1.;
842 //
843 // Kinematic limit on z
844 //
4e383037 845 if (m > 0.) zmax = 1. - m / TMath::Sqrt(m * m + jt * jt);
6e90ad26 846 //
847 // Change light-cone kinematics rel. to initial parton
848 //
849 Double_t eppzOld = e + pl;
850 Double_t empzOld = e - pl;
851
852 Double_t eppzNew = (1. - z) * eppzOld;
853 Double_t empzNew = empzOld - mt2 * z / eppzOld;
854 Double_t eNew = 0.5 * (eppzNew + empzNew);
855 Double_t plNew = 0.5 * (eppzNew - empzNew);
856
857 Double_t jtNew;
858 //
859 // if mt very small (or sometimes even < 0 for numerical reasons) set it to 0
860 Double_t mt2New = eppzNew * empzNew;
861 if (mt2New < 1.e-8) mt2New = 0.;
4e383037 862 if (z < zmax) {
863 if (m * m > mt2New) {
864 //
865 // This should not happen
866 //
867 Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
868 jtNew = 0;
869 } else {
870 jtNew = TMath::Sqrt(mt2New - m * m);
871 }
6e90ad26 872 } else {
4e383037 873 // If pT is to small (probably a leading massive particle) we scale only the energy
874 // This can cause negative masses of the radiated gluon
875 // Let's hope for the best ...
876 jtNew = jt;
877 eNew = TMath::Sqrt(plNew * plNew + mt2);
878
6e90ad26 879 }
6e90ad26 880 //
881 // Calculate new px, py
882 //
883 Double_t pxNew = jtNew / jt * pxs;
884 Double_t pyNew = jtNew / jt * pys;
885
886// Double_t dpx = pxs - pxNew;
887// Double_t dpy = pys - pyNew;
888// Double_t dpz = pl - plNew;
889// Double_t de = e - eNew;
890// Double_t dmass2 = de * de - dpx * dpx - dpy * dpy - dpz * dpz;
891// printf("New mass (1) %e %e %e %e %e %e %e \n", dmass2, jt, jtNew, pl, plNew, e, eNew);
892// printf("New mass (2) %e %e \n", pxNew, pyNew);
893 //
894 // Rotate back
895 //
896 TVector3 w(pxNew, pyNew, plNew);
897 w.RotateY(thetaq[index]); w.RotateZ(phiq[index]);
898 pxNew = w.X(); pyNew = w.Y(); plNew = w.Z();
899
900 p1[index][0] += pxNew;
901 p1[index][1] += pyNew;
902 p1[index][2] += plNew;
903 p1[index][3] += eNew;
904 //
905 // Updated 4-momentum vectors
906 //
907 pNew[icount][0] = pxNew;
908 pNew[icount][1] = pyNew;
909 pNew[icount][2] = plNew;
910 pNew[icount][3] = eNew;
911 kNew[icount] = i;
912 icount++;
913 } // parton loop
0f482ae4 914 //
6e90ad26 915 // Check if there was phase-space for quenching
0f482ae4 916 //
0f482ae4 917
6e90ad26 918 if (icount == 0) quenched[isys] = kFALSE;
919 if (!quenched[isys]) break;
920
921 for (Int_t j = 0; j < 4; j++)
922 {
923 p2[isys][j] = p0[isys][j] - p1[isys][j];
924 }
925 p2[isys][4] = p2[isys][3] * p2[isys][3] - p2[isys][0] * p2[isys][0] - p2[isys][1] * p2[isys][1] - p2[isys][2] * p2[isys][2];
6e90ad26 926 if (p2[isys][4] > 0.) {
927 p2[isys][4] = TMath::Sqrt(p2[isys][4]);
928 break;
929 } else {
b280c4cc 930 printf("Warning negative mass squared in system %d %f ! \n", isys, zquench[isys]);
4e383037 931 printf("4-Momentum: %10.3e %10.3e %10.3e %10.3e %10.3e \n", p2[isys][0], p2[isys][1], p2[isys][2], p2[isys][3], p2[isys][4]);
6e90ad26 932 if (p2[isys][4] < -0.01) {
4e383037 933 printf("Negative mass squared !\n");
934 // Here we have to put the gluon back to mass shell
935 // This will lead to a small energy imbalance
936 p2[isys][4] = 0.;
937 p2[isys][3] = TMath::Sqrt(p2[isys][0] * p2[isys][0] + p2[isys][1] * p2[isys][1] + p2[isys][2] * p2[isys][2]);
938 break;
6e90ad26 939 } else {
940 p2[isys][4] = 0.;
941 break;
942 }
943 }
6e90ad26 944 /*
6e90ad26 945 zHeavy *= 0.98;
946 printf("zHeavy lowered to %f\n", zHeavy);
947 if (zHeavy < 0.01) {
948 printf("No success ! \n");
949 icount = 0;
950 quenched[isys] = kFALSE;
951 break;
952 }
4e383037 953 */
954 } // iteration on z (while)
955
6e90ad26 956// Update event record
957 for (Int_t k = 0; k < icount; k++) {
958// printf("%6d %6d %10.3e %10.3e %10.3e %10.3e\n", k, kNew[k], pNew[k][0],pNew[k][1], pNew[k][2], pNew[k][3] );
959 fPyjets->P[0][kNew[k]] = pNew[k][0];
960 fPyjets->P[1][kNew[k]] = pNew[k][1];
961 fPyjets->P[2][kNew[k]] = pNew[k][2];
962 fPyjets->P[3][kNew[k]] = pNew[k][3];
0f482ae4 963 }
4e383037 964 //
965 // Add the gluons
966 //
967 Int_t ish = 0;
1837e95c 968 Int_t iGlu;
4e383037 969 if (!quenched[isys]) continue;
0f482ae4 970//
971// Last parton from shower i
4e383037 972 Int_t in = klast[isys];
0f482ae4 973//
974// Continue if no parton in shower i selected
975 if (in == -1) continue;
976//
977// If this is the second initial parton and it is behind the first move pointer by previous ish
4e383037 978 if (isys == 1 && klast[1] > klast[0]) in += ish;
0f482ae4 979//
980// Starting index
452af8c7 981
4e383037 982// jmin = in - 1;
0f482ae4 983// How many additional gluons will be generated
984 ish = 1;
4e383037 985 if (p2[isys][4] > 0.05) ish = 2;
0f482ae4 986//
987// Position of gluons
4e383037 988 iGlu = numpart;
989 if (iglu == 0) igMin = iGlu;
990 igMax = iGlu;
0f482ae4 991 numpart += ish;
992 (fPyjets->N) += ish;
4e383037 993
0f482ae4 994 if (ish == 1) {
4e383037 995 fPyjets->P[0][iGlu] = p2[isys][0];
996 fPyjets->P[1][iGlu] = p2[isys][1];
997 fPyjets->P[2][iGlu] = p2[isys][2];
998 fPyjets->P[3][iGlu] = p2[isys][3];
999 fPyjets->P[4][iGlu] = p2[isys][4];
0f482ae4 1000
4e383037 1001 fPyjets->K[0][iGlu] = 1;
1002 if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu] = 1;
0f482ae4 1003 fPyjets->K[1][iGlu] = 21;
4e383037 1004 fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
0f482ae4 1005 fPyjets->K[3][iGlu] = -1;
1006 fPyjets->K[4][iGlu] = -1;
4e383037 1007
1008 pg[0] += p2[isys][0];
1009 pg[1] += p2[isys][1];
1010 pg[2] += p2[isys][2];
1011 pg[3] += p2[isys][3];
0f482ae4 1012 } else {
1013 //
1014 // Split gluon in rest frame.
1015 //
4e383037 1016 Double_t bx = p2[isys][0] / p2[isys][3];
1017 Double_t by = p2[isys][1] / p2[isys][3];
1018 Double_t bz = p2[isys][2] / p2[isys][3];
1019 Double_t pst = p2[isys][4] / 2.;
0f482ae4 1020 //
1021 // Isotropic decay ????
1022 Double_t cost = 2. * gRandom->Rndm() - 1.;
1023 Double_t sint = TMath::Sqrt(1. - cost * cost);
1024 Double_t phi = 2. * TMath::Pi() * gRandom->Rndm();
1025
1026 Double_t pz1 = pst * cost;
1027 Double_t pz2 = -pst * cost;
1028 Double_t pt1 = pst * sint;
1029 Double_t pt2 = -pst * sint;
1030 Double_t px1 = pt1 * TMath::Cos(phi);
1031 Double_t py1 = pt1 * TMath::Sin(phi);
1032 Double_t px2 = pt2 * TMath::Cos(phi);
1033 Double_t py2 = pt2 * TMath::Sin(phi);
1034
1035 fPyjets->P[0][iGlu] = px1;
1036 fPyjets->P[1][iGlu] = py1;
1037 fPyjets->P[2][iGlu] = pz1;
1038 fPyjets->P[3][iGlu] = pst;
1039 fPyjets->P[4][iGlu] = 0.;
1040
4e383037 1041 fPyjets->K[0][iGlu] = 1 ;
0f482ae4 1042 fPyjets->K[1][iGlu] = 21;
4e383037 1043 fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
0f482ae4 1044 fPyjets->K[3][iGlu] = -1;
1045 fPyjets->K[4][iGlu] = -1;
1046
1047 fPyjets->P[0][iGlu+1] = px2;
1048 fPyjets->P[1][iGlu+1] = py2;
1049 fPyjets->P[2][iGlu+1] = pz2;
1050 fPyjets->P[3][iGlu+1] = pst;
1051 fPyjets->P[4][iGlu+1] = 0.;
1052
4e383037 1053 fPyjets->K[0][iGlu+1] = 1;
1054 if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu+1] = 1;
0f482ae4 1055 fPyjets->K[1][iGlu+1] = 21;
4e383037 1056 fPyjets->K[2][iGlu+1] = fPyjets->K[2][in] + 1000;
0f482ae4 1057 fPyjets->K[3][iGlu+1] = -1;
1058 fPyjets->K[4][iGlu+1] = -1;
1059 SetMSTU(1,0);
1060 SetMSTU(2,0);
1061 //
1062 // Boost back
1063 //
1064 Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
1065 }
4e383037 1066/*
1067 for (Int_t ig = iGlu; ig < iGlu+ish; ig++) {
1068 Double_t px, py, pz;
1069 px = fPyjets->P[0][ig];
1070 py = fPyjets->P[1][ig];
1071 pz = fPyjets->P[2][ig];
1072 TVector3 v(px, py, pz);
1073 v.RotateZ(-phiq[isys]);
1074 v.RotateY(-thetaq[isys]);
1075 Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pzs = v.Z();
1076 Double_t r = AliPythiaRndm::GetPythiaRandom()->Rndm();
1077 Double_t jtKick = 0.3 * TMath::Sqrt(-TMath::Log(r));
1078 if (ish == 2) jtKick = wjtKick[i] * TMath::Sqrt(-TMath::Log(r)) / TMath::Sqrt(2.);
1079 Double_t phiKick = 2. * TMath::Pi() * AliPythiaRndm::GetPythiaRandom()->Rndm();
1080 pxs += jtKick * TMath::Cos(phiKick);
1081 pys += jtKick * TMath::Sin(phiKick);
1082 TVector3 w(pxs, pys, pzs);
1083 w.RotateY(thetaq[isys]);
1084 w.RotateZ(phiq[isys]);
1085 fPyjets->P[0][ig] = w.X();
1086 fPyjets->P[1][ig] = w.Y();
1087 fPyjets->P[2][ig] = w.Z();
1088 fPyjets->P[2][ig] = w.Mag();
1089 }
1090*/
1091 } // kGluon
1092
6e90ad26 1093
4e383037 1094 // Check energy conservation
0f482ae4 1095 Double_t pxs = 0.;
1096 Double_t pys = 0.;
1097 Double_t pzs = 0.;
1098 Double_t es = 14000.;
1099
1100 for (Int_t i = 0; i < numpart; i++)
1101 {
1102 kst = fPyjets->K[0][i];
1103 if (kst != 1 && kst != 2) continue;
1104 pxs += fPyjets->P[0][i];
1105 pys += fPyjets->P[1][i];
1106 pzs += fPyjets->P[2][i];
1107 es -= fPyjets->P[3][i];
1108 }
1109 if (TMath::Abs(pxs) > 1.e-2 ||
1110 TMath::Abs(pys) > 1.e-2 ||
1111 TMath::Abs(pzs) > 1.e-1) {
1112 printf("%e %e %e %e\n", pxs, pys, pzs, es);
4e383037 1113// Fatal("Quench()", "4-Momentum non-conservation");
452af8c7 1114 }
4e383037 1115
1116 } // end quenching loop (systems)
6e90ad26 1117// Clean-up
0f482ae4 1118 for (Int_t i = 0; i < numpart; i++)
1119 {
4e383037 1120 imo = fPyjets->K[2][i];
1121 if (imo > 1000) {
1122 fPyjets->K[2][i] = fPyjets->K[2][i] % 1000;
1123 }
0f482ae4 1124 }
4e383037 1125// this->Pylist(1);
0f482ae4 1126} // end quench
90d7b703 1127
992f2843 1128
1129void AliPythia::Pyquen(Double_t a, Int_t ibf, Double_t b)
1130{
1131 // Igor Lokthine's quenching routine
1132 pyquen(a, ibf, b);
1133}
b280c4cc 1134
16a82508 1135void AliPythia::Pyevnw()
1136{
1137 // New multiple interaction scenario
1138 pyevnw();
1139}
1140
b280c4cc 1141void AliPythia::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])
1142{
1143 // Return event specific quenching parameters
1144 xp = fXJet;
1145 yp = fYJet;
1146 for (Int_t i = 0; i < 4; i++) z[i] = fZQuench[i];
1147
1148}
1149
3dc3ec94 1150void AliPythia::ConfigHeavyFlavor()
1151{
1152 //
1153 // Default configuration for Heavy Flavor production
1154 //
1155 // All QCD processes
1156 //
1157 SetMSEL(1);
1158
1159 // No multiple interactions
1160 SetMSTP(81,0);
1161 SetPARP(81,0.0);
1162 SetPARP(82,0.0);
1163
1164 // Initial/final parton shower on (Pythia default)
1165 SetMSTP(61,1);
1166 SetMSTP(71,1);
1167
1168 // 2nd order alpha_s
1169 SetMSTP(2,2);
1170
1171 // QCD scales
1172 SetMSTP(32,2);
1173 SetPARP(34,1.0);
1174}
e0e89f40 1175
1176void AliPythia::AtlasTuning()
1177{
1178 //
1179 // Configuration for the ATLAS tuning
1180 SetMSTP(51, kCTEQ5L); // CTEQ5L pdf
1181 SetMSTP(81,1); // Multiple Interactions ON
1182 SetMSTP(82,4); // Double Gaussian Model
1183 SetPARP(82,1.8); // [GeV] PT_min at Ref. energy
1184 SetPARP(89,1000.); // [GeV] Ref. energy
1185 SetPARP(90,0.16); // 2*epsilon (exponent in power law)
1186 SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
1187 SetPARP(84,0.5); // Core radius
1188 SetPARP(85,0.33); // Regulates gluon prod. mechanism
1189 SetPARP(86,0.66); // Regulates gluon prod. mechanism
1190 SetPARP(67,1); // Regulates Initial State Radiation
1191}