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