Change W50511 common block in pythia-6.4.25
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia6.cxx
CommitLineData
39d810c8 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
16/* $Id: AliPythia.cxx,v 1.40 2007/10/09 08:43:24 morsch Exp $ */
17
18#include "AliPythia6.h"
19#include "AliStack.h"
20#include "AliPythiaRndm.h"
21#include "AliFastGlauber.h"
22#include "AliQuenchingWeights.h"
23
24#include "TVector3.h"
25#include "TParticle.h"
26#include "PyquenCommon.h"
27
28ClassImp(AliPythia6)
29
30#ifndef WIN32
31# define pyclus pyclus_
32# define pycell pycell_
33# define pyshow pyshow_
f3cbba5f 34# define pyshowq pyshowq_
39d810c8 35# define pyrobo pyrobo_
36# define pyquen pyquen_
37# define pyevnw pyevnw_
8b9bb87a 38# define pyjoin pyjoin_
f3cbba5f 39# define qpygin0 qpygin0_
1cdfb96f 40# define setpowwght setpowwght_
39d810c8 41# define type_of_call
42#else
43# define pyclus PYCLUS
44# define pycell PYCELL
f3cbba5f 45# define pyshow PYSHOW
46# define pyshowq PYSHOWQ
39d810c8 47# define pyrobo PYROBO
48# define pyquen PYQUEN
49# define pyevnw PYEVNW
8b9bb87a 50# define pyjoin PYJOIN
f3cbba5f 51# define qpygin0 QPYGIN0
1cdfb96f 52# define setpowwght SETPOWWGHT
39d810c8 53# define type_of_call _stdcall
54#endif
55
8b9bb87a 56extern "C" void type_of_call pyjoin(Int_t &, Int_t * );
39d810c8 57extern "C" void type_of_call pyclus(Int_t & );
58extern "C" void type_of_call pycell(Int_t & );
59extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
f3cbba5f 60extern "C" void type_of_call pyshowq(Int_t &, Int_t &, Double_t &);
61extern "C" void type_of_call qpygin0();
39d810c8 62extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
63extern "C" void type_of_call pyquen(Double_t &, Int_t &, Double_t &);
64extern "C" void type_of_call pyevnw();
1cdfb96f 65extern "C" void type_of_call setpowwght(Double_t &);
39d810c8 66
67
68//_____________________________________________________________________________
69
70AliPythia6* AliPythia6::fgAliPythia=NULL;
71
72AliPythia6::AliPythia6():
73 TPythia6(),
74 AliPythiaBase(),
75 fProcess(kPyMb),
76 fEcms(0.),
77 fStrucFunc(kCTEQ5L),
e28ccdaf 78 fProjectile("p"),
79 fTarget("p"),
39d810c8 80 fXJet(0.),
81 fYJet(0.),
82 fNGmax(30),
83 fZmax(0.97),
84 fGlauber(0),
85 fQuenchingWeights(0)
86{
87// Default Constructor
88//
89// Set random number
a23e397a 90 Int_t i;
91 for (i = 0; i < 501; i++) fDefMDCY[i] = 0;
92 for (i = 0; i < 2001; i++) fDefMDME[i] = 0;
93 for (i = 0; i < 4; i++) fZQuench[i] = 0;
94
39d810c8 95 if (!AliPythiaRndm::GetPythiaRandom())
96 AliPythiaRndm::SetPythiaRandom(GetRandom());
97 fGlauber = 0;
98 fQuenchingWeights = 0;
99}
100
101AliPythia6::AliPythia6(const AliPythia6& pythia):
102 TPythia6(),
103 AliPythiaBase(),
104 fProcess(kPyMb),
105 fEcms(0.),
106 fStrucFunc(kCTEQ5L),
e28ccdaf 107 fProjectile("p"),
108 fTarget("p"),
39d810c8 109 fXJet(0.),
110 fYJet(0.),
111 fNGmax(30),
112 fZmax(0.97),
113 fGlauber(0),
114 fQuenchingWeights(0)
115{
116 // Copy Constructor
a23e397a 117 Int_t i;
118 for (i = 0; i < 501; i++) fDefMDCY[i] = 0;
119 for (i = 0; i < 2001; i++) fDefMDME[i] = 0;
120 for (i = 0; i < 4; i++) fZQuench[i] = 0;
39d810c8 121 pythia.Copy(*this);
122}
123
75d4f39e 124void AliPythia6::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc, Int_t /*tune*/)
39d810c8 125{
126// Initialise the process to generate
127 if (!AliPythiaRndm::GetPythiaRandom())
128 AliPythiaRndm::SetPythiaRandom(GetRandom());
129
130 fProcess = process;
131 fEcms = energy;
132 fStrucFunc = strucfunc;
133//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
134 SetMDCY(Pycomp(111) ,1,0);
135 SetMDCY(Pycomp(310) ,1,0);
136 SetMDCY(Pycomp(3122),1,0);
137 SetMDCY(Pycomp(3112),1,0);
138 SetMDCY(Pycomp(3212),1,0);
139 SetMDCY(Pycomp(3222),1,0);
140 SetMDCY(Pycomp(3312),1,0);
141 SetMDCY(Pycomp(3322),1,0);
142 SetMDCY(Pycomp(3334),1,0);
143 // Select structure function
144 SetMSTP(52,2);
145 SetMSTP(51,AliStructFuncType::PDFsetIndex(strucfunc));
146 // Particles produced in string fragmentation point directly to either of the two endpoints
147 // of the string (depending in the side they were generated from).
148 SetMSTU(16,2);
149
150//
151// Pythia initialisation for selected processes//
152//
153// Make MSEL clean
154//
155 for (Int_t i=1; i<= 200; i++) {
156 SetMSUB(i,0);
157 }
158// select charm production
159 switch (process)
160 {
161 case kPyOldUEQ2ordered: //Old underlying events with Q2 ordered QCD processes
162// Multiple interactions on.
163 SetMSTP(81,1);
164// Double Gaussian matter distribution.
165 SetMSTP(82,4);
166 SetPARP(83,0.5);
167 SetPARP(84,0.4);
168// pT0.
169 SetPARP(82,2.0);
170// Reference energy for pT0 and energy rescaling pace.
171 SetPARP(89,1800);
172 SetPARP(90,0.25);
173// String drawing almost completely minimizes string length.
174 SetPARP(85,0.9);
175 SetPARP(86,0.95);
176// ISR and FSR activity.
177 SetPARP(67,4);
178 SetPARP(71,4);
179// Lambda_FSR scale.
180 SetPARJ(81,0.29);
181 break;
182 case kPyOldUEQ2ordered2:
183// Old underlying events with Q2 ordered QCD processes
184// Multiple interactions on.
185 SetMSTP(81,1);
186// Double Gaussian matter distribution.
187 SetMSTP(82,4);
188 SetPARP(83,0.5);
189 SetPARP(84,0.4);
190// pT0.
191 SetPARP(82,2.0);
192// Reference energy for pT0 and energy rescaling pace.
193 SetPARP(89,1800);
194 SetPARP(90,0.16); // here is the difference with kPyOldUEQ2ordered
195// String drawing almost completely minimizes string length.
196 SetPARP(85,0.9);
197 SetPARP(86,0.95);
198// ISR and FSR activity.
199 SetPARP(67,4);
200 SetPARP(71,4);
201// Lambda_FSR scale.
202 SetPARJ(81,0.29);
203 break;
204 case kPyOldPopcorn:
205// Old production mechanism: Old Popcorn
206 SetMSEL(1);
207 SetMSTJ(12,3);
208// (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon
209 SetMSTP(88,2);
210// (D=1)see can be used to form baryons (BARYON JUNCTION)
211 SetMSTJ(1,1);
212 AtlasTuning();
213 break;
214 case kPyCharm:
215 SetMSEL(4);
216// heavy quark masses
217
218 SetPMAS(4,1,1.2);
219//
220// primordial pT
221 SetMSTP(91,1);
222 SetPARP(91,1.);
223 SetPARP(93,5.);
224//
225 break;
226 case kPyBeauty:
227 SetMSEL(5);
228 SetPMAS(5,1,4.75);
229 break;
230 case kPyJpsi:
231 SetMSEL(0);
232// gg->J/Psi g
233 SetMSUB(86,1);
234 break;
235 case kPyJpsiChi:
236 SetMSEL(0);
237// gg->J/Psi g
238 SetMSUB(86,1);
239// gg-> chi_0c g
240 SetMSUB(87,1);
241// gg-> chi_1c g
242 SetMSUB(88,1);
243// gg-> chi_2c g
244 SetMSUB(89,1);
245 break;
246 case kPyCharmUnforced:
247 SetMSEL(0);
248// gq->qg
249 SetMSUB(28,1);
250// gg->qq
251 SetMSUB(53,1);
252// gg->gg
253 SetMSUB(68,1);
254 break;
255 case kPyBeautyUnforced:
256 SetMSEL(0);
257// gq->qg
258 SetMSUB(28,1);
259// gg->qq
260 SetMSUB(53,1);
261// gg->gg
262 SetMSUB(68,1);
263 break;
264 case kPyMb:
265// Minimum Bias pp-Collisions
266//
267//
268// select Pythia min. bias model
269 SetMSEL(0);
270 SetMSUB(92,1); // single diffraction AB-->XB
271 SetMSUB(93,1); // single diffraction AB-->AX
272 SetMSUB(94,1); // double diffraction
273 SetMSUB(95,1); // low pt production
274
275 AtlasTuning();
276 break;
67b7bb0e 277 case kPyMbAtlasTuneMC09:
278// Minimum Bias pp-Collisions
279//
280//
281// select Pythia min. bias model
282 SetMSEL(0);
283 SetMSUB(92,1); // single diffraction AB-->XB
284 SetMSUB(93,1); // single diffraction AB-->AX
285 SetMSUB(94,1); // double diffraction
286 SetMSUB(95,1); // low pt production
287
b4d68b2a 288 AtlasTuningMC09();
67b7bb0e 289 break;
04081a91 290
291 case kPyMbWithDirectPhoton:
292// Minimum Bias pp-Collisions with direct photon processes added
293//
294//
295// select Pythia min. bias model
296 SetMSEL(0);
297 SetMSUB(92,1); // single diffraction AB-->XB
298 SetMSUB(93,1); // single diffraction AB-->AX
299 SetMSUB(94,1); // double diffraction
300 SetMSUB(95,1); // low pt production
301
302 SetMSUB(14,1); //
303 SetMSUB(18,1); //
304 SetMSUB(29,1); //
305 SetMSUB(114,1); //
306 SetMSUB(115,1); //
307
308
309 AtlasTuning();
310 break;
311
39d810c8 312 case kPyMbDefault:
313// Minimum Bias pp-Collisions
314//
315//
316// select Pythia min. bias model
317 SetMSEL(0);
318 SetMSUB(92,1); // single diffraction AB-->XB
319 SetMSUB(93,1); // single diffraction AB-->AX
320 SetMSUB(94,1); // double diffraction
321 SetMSUB(95,1); // low pt production
322
323 break;
324 case kPyLhwgMb:
325// Les Houches Working Group 05 Minimum Bias pp-Collisions: hep-ph/0604120
326// -> Pythia 6.3 or above is needed
327//
328 SetMSEL(0);
329 SetMSUB(92,1); // single diffraction AB-->XB
330 SetMSUB(93,1); // single diffraction AB-->AX
331 SetMSUB(94,1); // double diffraction
332 SetMSUB(95,1); // low pt production
333 SetMSTP(51,AliStructFuncType::PDFsetIndex(kCTEQ6ll));
334 SetMSTP(52,2);
335 SetMSTP(68,1);
336 SetMSTP(70,2);
337 SetMSTP(81,1); // Multiple Interactions ON
338 SetMSTP(82,4); // Double Gaussian Model
339 SetMSTP(88,1);
340
341 SetPARP(82,2.3); // [GeV] PT_min at Ref. energy
342 SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
343 SetPARP(84,0.5); // Core radius
344 SetPARP(85,0.9); // Regulates gluon prod. mechanism
345 SetPARP(90,0.2); // 2*epsilon (exponent in power law)
346
347 break;
348 case kPyMbNonDiffr:
349// Minimum Bias pp-Collisions
350//
351//
352// select Pythia min. bias model
353 SetMSEL(0);
354 SetMSUB(95,1); // low pt production
355
356 AtlasTuning();
357 break;
358 case kPyMbMSEL1:
359 ConfigHeavyFlavor();
360// Intrinsic <kT^2>
361 SetMSTP(91,1);// Width (1=gaussian) primordial kT dist. inside hadrons
362 SetPARP(91,1.); // <kT^2> = PARP(91,1.)^2
363 SetPARP(93,5.); // Upper cut-off
364// Set Q-quark mass
365 SetPMAS(4,1,1.2); // Charm quark mass
366 SetPMAS(5,1,4.78); // Beauty quark mass
367 SetPARP(71,4.); // Defaut value
368// Atlas Tuning
369 AtlasTuning();
370 break;
371 case kPyJets:
372//
373// QCD Jets
374//
375 SetMSEL(1);
376 // Pythia Tune A (CDF)
377 //
378 SetPARP(67,2.5); // Regulates Initial State Radiation (value from best fit to D0 dijet analysis)
379 SetMSTP(82,4); // Double Gaussian Model
380 SetPARP(82,2.0); // [GeV] PT_min at Ref. energy
381 SetPARP(84,0.4); // Core radius
382 SetPARP(85,0.90) ; // Regulates gluon prod. mechanism
383 SetPARP(86,0.95); // Regulates gluon prod. mechanism
384 SetPARP(89,1800.); // [GeV] Ref. energy
385 SetPARP(90,0.25); // 2*epsilon (exponent in power law)
386 break;
387 case kPyDirectGamma:
388 SetMSEL(10);
389 break;
390 case kPyCharmPbPbMNR:
391 case kPyD0PbPbMNR:
392 case kPyDPlusPbPbMNR:
393 case kPyDPlusStrangePbPbMNR:
394 // Tuning of Pythia parameters aimed to get a resonable agreement
395 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
396 // c-cbar single inclusive and double differential distributions.
397 // This parameter settings are meant to work with Pb-Pb collisions
398 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
399 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
400 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
401 ConfigHeavyFlavor();
402 // Intrinsic <kT>
403 SetMSTP(91,1);
404 SetPARP(91,1.304);
405 SetPARP(93,6.52);
406 // Set c-quark mass
407 SetPMAS(4,1,1.2);
408 break;
409 case kPyCharmpPbMNR:
410 case kPyD0pPbMNR:
411 case kPyDPluspPbMNR:
412 case kPyDPlusStrangepPbMNR:
413 // Tuning of Pythia parameters aimed to get a resonable agreement
414 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
415 // c-cbar single inclusive and double differential distributions.
416 // This parameter settings are meant to work with p-Pb collisions
417 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
418 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
419 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
420 ConfigHeavyFlavor();
421 // Intrinsic <kT>
422 SetMSTP(91,1);
423 SetPARP(91,1.16);
424 SetPARP(93,5.8);
425
426 // Set c-quark mass
427 SetPMAS(4,1,1.2);
428 break;
429 case kPyCharmppMNR:
430 case kPyD0ppMNR:
431 case kPyDPlusppMNR:
432 case kPyDPlusStrangeppMNR:
68504d42 433 case kPyLambdacppMNR:
39d810c8 434 // Tuning of Pythia parameters aimed to get a resonable agreement
435 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
436 // c-cbar single inclusive and double differential distributions.
437 // This parameter settings are meant to work with pp collisions
438 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
439 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
440 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
441 ConfigHeavyFlavor();
442 // Intrinsic <kT^2>
443 SetMSTP(91,1);
444 SetPARP(91,1.);
445 SetPARP(93,5.);
446
447 // Set c-quark mass
448 SetPMAS(4,1,1.2);
449 break;
450 case kPyCharmppMNRwmi:
451 // Tuning of Pythia parameters aimed to get a resonable agreement
452 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
453 // c-cbar single inclusive and double differential distributions.
454 // This parameter settings are meant to work with pp collisions
455 // and with kCTEQ5L PDFs.
456 // Added multiple interactions according to ATLAS tune settings.
457 // To get a "reasonable" agreement with MNR results, events have to be
458 // generated with the minimum ptHard (AliGenPythia::SetPtHard)
459 // set to 2.76 GeV.
460 // To get a "perfect" agreement with MNR results, events have to be
461 // generated in four ptHard bins with the following relative
462 // normalizations:
463 // 2.76-3 GeV: 25%
464 // 3-4 GeV: 40%
465 // 4-8 GeV: 29%
466 // >8 GeV: 6%
467 ConfigHeavyFlavor();
468 // Intrinsic <kT^2>
469 SetMSTP(91,1);
470 SetPARP(91,1.);
471 SetPARP(93,5.);
472
473 // Set c-quark mass
474 SetPMAS(4,1,1.2);
475 AtlasTuning();
476 break;
477 case kPyBeautyPbPbMNR:
478 // Tuning of Pythia parameters aimed to get a resonable agreement
479 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
480 // b-bbar single inclusive and double differential distributions.
481 // This parameter settings are meant to work with Pb-Pb collisions
482 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
483 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
484 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
485 ConfigHeavyFlavor();
486 // QCD scales
487 SetPARP(67,1.0);
488 SetPARP(71,1.0);
489 // Intrinsic <kT>
490 SetMSTP(91,1);
491 SetPARP(91,2.035);
492 SetPARP(93,10.17);
493 // Set b-quark mass
494 SetPMAS(5,1,4.75);
495 break;
496 case kPyBeautypPbMNR:
497 // Tuning of Pythia parameters aimed to get a resonable agreement
498 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
499 // b-bbar single inclusive and double differential distributions.
500 // This parameter settings are meant to work with p-Pb collisions
501 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
502 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
503 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
504 ConfigHeavyFlavor();
505 // QCD scales
506 SetPARP(67,1.0);
507 SetPARP(71,1.0);
508 // Intrinsic <kT>
509 SetMSTP(91,1);
510 SetPARP(91,1.60);
511 SetPARP(93,8.00);
512 // Set b-quark mass
513 SetPMAS(5,1,4.75);
514 break;
515 case kPyBeautyppMNR:
516 // Tuning of Pythia parameters aimed to get a resonable agreement
517 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
518 // b-bbar single inclusive and double differential distributions.
519 // This parameter settings are meant to work with pp collisions
520 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
521 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
522 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
523 ConfigHeavyFlavor();
524 // QCD scales
525 SetPARP(67,1.0);
526 SetPARP(71,1.0);
527
528 // Intrinsic <kT>
529 SetMSTP(91,1);
530 SetPARP(91,1.);
531 SetPARP(93,5.);
532
533 // Set b-quark mass
534 SetPMAS(5,1,4.75);
535 break;
70574ff8 536 case kPyBeautyJets:
39d810c8 537 case kPyBeautyppMNRwmi:
538 // Tuning of Pythia parameters aimed to get a resonable agreement
539 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
540 // b-bbar single inclusive and double differential distributions.
541 // This parameter settings are meant to work with pp collisions
542 // and with kCTEQ5L PDFs.
543 // Added multiple interactions according to ATLAS tune settings.
544 // To get a "reasonable" agreement with MNR results, events have to be
545 // generated with the minimum ptHard (AliGenPythia::SetPtHard)
546 // set to 2.76 GeV.
547 // To get a "perfect" agreement with MNR results, events have to be
548 // generated in four ptHard bins with the following relative
549 // normalizations:
550 // 2.76-4 GeV: 5%
551 // 4-6 GeV: 31%
552 // 6-8 GeV: 28%
553 // >8 GeV: 36%
554 ConfigHeavyFlavor();
555 // QCD scales
556 SetPARP(67,1.0);
557 SetPARP(71,1.0);
558
559 // Intrinsic <kT>
560 SetMSTP(91,1);
561 SetPARP(91,1.);
562 SetPARP(93,5.);
563
564 // Set b-quark mass
565 SetPMAS(5,1,4.75);
566
567 AtlasTuning();
568 break;
569 case kPyW:
570
571 //Inclusive production of W+/-
572 SetMSEL(0);
573 //f fbar -> W+
574 SetMSUB(2,1);
575 // //f fbar -> g W+
576 // SetMSUB(16,1);
577 // //f fbar -> gamma W+
578 // SetMSUB(20,1);
579 // //f g -> f W+
580 // SetMSUB(31,1);
581 // //f gamma -> f W+
582 // SetMSUB(36,1);
583
584 // Initial/final parton shower on (Pythia default)
585 // With parton showers on we are generating "W inclusive process"
586 SetMSTP(61,1); //Initial QCD & QED showers on
587 SetMSTP(71,1); //Final QCD & QED showers on
588
589 break;
590
591 case kPyZ:
592
593 //Inclusive production of Z
594 SetMSEL(0);
595 //f fbar -> Z/gamma
596 SetMSUB(1,1);
597
598 // // f fbar -> g Z/gamma
599 // SetMSUB(15,1);
600 // // f fbar -> gamma Z/gamma
601 // SetMSUB(19,1);
602 // // f g -> f Z/gamma
603 // SetMSUB(30,1);
604 // // f gamma -> f Z/gamma
605 // SetMSUB(35,1);
606
607 //only Z included, not gamma
608 SetMSTP(43,2);
609
610 // Initial/final parton shower on (Pythia default)
611 // With parton showers on we are generating "Z inclusive process"
612 SetMSTP(61,1); //Initial QCD & QED showers on
613 SetMSTP(71,1); //Final QCD & QED showers on
9a8774a1 614 break;
df607629 615 case kPyZgamma:
616 //Inclusive production of Z
617 SetMSEL(0);
618 //f fbar -> Z/gamma
619 SetMSUB(1,1);
620 // Initial/final parton shower on (Pythia default)
621 // With parton showers on we are generating "Z inclusive process"
622 SetMSTP(61,1); //Initial QCD & QED showers on
623 SetMSTP(71,1); //Final QCD & QED showers on
624 break;
9a8774a1 625 case kPyMBRSingleDiffraction:
626 case kPyMBRDoubleDiffraction:
627 case kPyMBRCentralDiffraction:
39d810c8 628 break;
64da86aa 629 case kPyJetsPWHG:
630 // N.B.
631 // ====
632 // For the case of jet production the following parameter setting
633 // limits the transverse momentum of secondary scatterings, due
634 // to multiple parton interactions, to be less than that of the
635 // primary interaction (see POWHEG Dijet paper arXiv:1012.3380
636 // [hep-ph] sec. 4.1 and also the PYTHIA Manual).
637 SetMSTP(86,1);
638 // maximum number of errors before pythia aborts (def=10)
639 SetMSTU(22,10);
640 // number of warnings printed on the shell
641 SetMSTU(26,20);
642 break;
a6d31e74 643 case kPyCharmPWHG:
644 case kPyBeautyPWHG:
df607629 645 case kPyWPWHG:
a6d31e74 646 // number of warnings printed on the shell
647 SetMSTU(26,20);
648
39d810c8 649 }
650//
651// Initialize PYTHIA
652 SetMSTP(41,1); // all resonance decays switched on
e136a92a 653 if (process == kPyJetsPWHG || process == kPyCharmPWHG || process == kPyBeautyPWHG || process == kPyWPWHG) {
64da86aa 654 Initialize("USER","","",0.);
655 } else {
656 Initialize("CMS",fProjectile,fTarget,fEcms);
657 }
39d810c8 658}
659
660Int_t AliPythia6::CheckedLuComp(Int_t kf)
661{
662// Check Lund particle code (for debugging)
663 Int_t kc=Pycomp(kf);
664 return kc;
665}
666
667void AliPythia6::SetNuclei(Int_t a1, Int_t a2)
668{
669// Treat protons as inside nuclei with mass numbers a1 and a2
670// The MSTP array in the PYPARS common block is used to enable and
671// select the nuclear structure functions.
672// MSTP(52) : (D=1) choice of proton and nuclear structure-function library
673// =1: internal PYTHIA acording to MSTP(51)
674// =2: PDFLIB proton s.f., with MSTP(51) = 1000xNGROUP+NSET
675// If the following mass number both not equal zero, nuclear corrections of the stf are used.
676// MSTP(192) : Mass number of nucleus side 1
677// MSTP(193) : Mass number of nucleus side 2
678 SetMSTP(52,2);
679 SetMSTP(192, a1);
680 SetMSTP(193, a2);
681}
682
683
684AliPythia6* AliPythia6::Instance()
685{
686// Set random number generator
687 if (fgAliPythia) {
688 return fgAliPythia;
689 } else {
690 fgAliPythia = new AliPythia6();
691 return fgAliPythia;
692 }
693}
694
695void AliPythia6::PrintParticles()
696{
697// Print list of particl properties
698 Int_t np = 0;
699 char* name = new char[16];
700 for (Int_t kf=0; kf<1000000; kf++) {
701 for (Int_t c = 1; c > -2; c-=2) {
702 Int_t kc = Pycomp(c*kf);
703 if (kc) {
704 Float_t mass = GetPMAS(kc,1);
705 Float_t width = GetPMAS(kc,2);
706 Float_t tau = GetPMAS(kc,4);
707
708 Pyname(kf,name);
709
710 np++;
711
712 printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
713 c*kf, name, mass, width, tau);
714 }
715 }
716 }
717 printf("\n Number of particles %d \n \n", np);
718}
719
720void AliPythia6::ResetDecayTable()
721{
722// Set default values for pythia decay switches
723 Int_t i;
724 for (i = 1; i < 501; i++) SetMDCY(i,1,fDefMDCY[i]);
725 for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
726}
727
728void AliPythia6::SetDecayTable()
729{
730// Set default values for pythia decay switches
731//
732 Int_t i;
733 for (i = 1; i < 501; i++) fDefMDCY[i] = GetMDCY(i,1);
734 for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
735}
736
8b9bb87a 737void AliPythia6::Pyjoin(Int_t& npart, Int_t *ipart)
738{
739// Call Pythia join alogorithm to set up a string between
740// npart partons, given by indices in array ipart[npart]
741//
742 pyjoin(npart, ipart);
743}
744
f3cbba5f 745void AliPythia6::Pyshowq(Int_t ip1, Int_t ip2, Double_t qmax)
746{
747// Call qPythia showering
748//
749 pyshowq(ip1, ip2, qmax);
750}
751
752void AliPythia6::Qpygin0()
753{
754 //position of the hard scattering in the nuclear overlapping area.
755 //just for qpythia.
756 qpygin0();
757}
758
39d810c8 759void AliPythia6::Pyclus(Int_t& njet)
760{
761// Call Pythia clustering algorithm
762//
763 pyclus(njet);
764}
765
766void AliPythia6::Pycell(Int_t& njet)
767{
768// Call Pythia jet reconstruction algorithm
769//
770 pycell(njet);
771}
772
773void AliPythia6::GetJet(Int_t i, Float_t& px, Float_t& py, Float_t& pz, Float_t& e)
774{
775 // Get jet number i
776 Int_t n = GetN();
777 px = GetPyjets()->P[0][n+i];
778 py = GetPyjets()->P[1][n+i];
779 pz = GetPyjets()->P[2][n+i];
780 e = GetPyjets()->P[3][n+i];
781}
782
783void AliPythia6::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
784{
785// Call Pythia showering
786//
787 pyshow(ip1, ip2, qmax);
788}
789
790void AliPythia6::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
791{
792 pyrobo(imi, ima, the, phi, bex, bey, bez);
793}
794
795
796
797void AliPythia6::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECMethod, Float_t zmax, Int_t ngmax)
798{
799// Initializes
800// (1) The quenching model using quenching weights according to C. Salgado and U. Wiedemann
801// (2) The nuclear geometry using the Glauber Model
802//
803
18b7a4a1 804 fGlauber = AliFastGlauber::Instance();
39d810c8 805 fGlauber->Init(2);
806 fGlauber->SetCentralityClass(cMin, cMax);
807
808 fQuenchingWeights = new AliQuenchingWeights();
809 fQuenchingWeights->InitMult();
810 fQuenchingWeights->SetK(k);
811 fQuenchingWeights->SetECMethod(AliQuenchingWeights::kECMethod(iECMethod));
812 fNGmax = ngmax;
813 fZmax = zmax;
814
815}
816
817
818void AliPythia6::Quench()
819{
820//
821//
822// Simple Jet Quenching routine:
823// =============================
824// The jet formed by all final state partons radiated by the parton created
825// in the hard collisions is quenched by a factor (1-z) using light cone variables in
826// the initial parton reference frame:
827// (E + p_z)new = (1-z) (E + p_z)old
828//
829//
830//
831//
832// The lost momentum is first balanced by one gluon with virtuality > 0.
833// Subsequently the gluon splits to yield two gluons with E = p.
834//
835//
836//
837 static Float_t eMean = 0.;
838 static Int_t icall = 0;
839
840 Double_t p0[4][5];
841 Double_t p1[4][5];
842 Double_t p2[4][5];
843 Int_t klast[4] = {-1, -1, -1, -1};
844
845 Int_t numpart = fPyjets->N;
846 Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0., phi = 0.;
847 Double_t pxq[4], pyq[4], pzq[4], eq[4], yq[4], mq[4], pq[4], phiq[4], thetaq[4], ptq[4];
848 Bool_t quenched[4];
a23e397a 849 Double_t wjtKick[4] = {0., 0., 0., 0.};
39d810c8 850 Int_t nGluon[4];
851 Int_t qPdg[4];
852 Int_t imo, kst, pdg;
853
854//
855// Sore information about Primary partons
856//
857// j =
858// 0, 1 partons from hard scattering
859// 2, 3 partons from initial state radiation
860//
861 for (Int_t i = 2; i <= 7; i++) {
862 Int_t j = 0;
863 // Skip gluons that participate in hard scattering
864 if (i == 4 || i == 5) continue;
865 // Gluons from hard Scattering
866 if (i == 6 || i == 7) {
867 j = i - 6;
868 pxq[j] = fPyjets->P[0][i];
869 pyq[j] = fPyjets->P[1][i];
870 pzq[j] = fPyjets->P[2][i];
871 eq[j] = fPyjets->P[3][i];
872 mq[j] = fPyjets->P[4][i];
873 } else {
874 // Gluons from initial state radiation
875 //
876 // Obtain 4-momentum vector from difference between original parton and parton after gluon
877 // radiation. Energy is calculated independently because initial state radition does not
878 // conserve strictly momentum and energy for each partonic system independently.
879 //
880 // Not very clean. Should be improved !
881 //
882 //
883 j = i;
884 pxq[j] = fPyjets->P[0][i] - fPyjets->P[0][i+2];
885 pyq[j] = fPyjets->P[1][i] - fPyjets->P[1][i+2];
886 pzq[j] = fPyjets->P[2][i] - fPyjets->P[2][i+2];
887 mq[j] = fPyjets->P[4][i];
888 eq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j] + mq[j] * mq[j]);
889 }
890//
891// Calculate some kinematic variables
892//
893 yq[j] = 0.5 * TMath::Log((eq[j] + pzq[j] + 1.e-14) / (eq[j] - pzq[j] + 1.e-14));
894 pq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j]);
895 phiq[j] = TMath::Pi()+TMath::ATan2(-pyq[j], -pxq[j]);
896 ptq[j] = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j]);
897 thetaq[j] = TMath::ATan2(ptq[j], pzq[j]);
898 qPdg[j] = fPyjets->K[1][i];
899 }
900
901 Double_t int0[4];
902 Double_t int1[4];
903
904 fGlauber->GetI0I1ForPythiaAndXY(4, phiq, int0, int1, fXJet, fYJet, 15.);
905
906 for (Int_t j = 0; j < 4; j++) {
907 //
908 // Quench only central jets and with E > 10.
909 //
910
911
912 Int_t itype = (qPdg[j] == 21) ? 2 : 1;
5c843db7 913 // Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]);
914 Double_t eloss = fQuenchingWeights->GetELossRandomK(itype, int0[j], int1[j], eq[j]);
39d810c8 915
916 if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
917 fZQuench[j] = 0.;
918 } else {
919 if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) {
920 icall ++;
921 eMean += eloss;
922 }
923 //
924 // Extra pt
925 Double_t l = fQuenchingWeights->CalcLk(int0[j], int1[j]);
926 wjtKick[j] = TMath::Sqrt(l * fQuenchingWeights->CalcQk(int0[j], int1[j]));
927 //
928 // Fractional energy loss
929 fZQuench[j] = eloss / eq[j];
930 //
931 // Avoid complete loss
932 //
1044c4d8 933 if (fZQuench[j] > fZmax) fZQuench[j] = fZmax;
39d810c8 934 //
935 // Some debug printing
936
937
938// 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",
939// j, itype, eq[j], phiq[j], l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]);
940
941// fZQuench[j] = 0.8;
942// while (fZQuench[j] >= 0.95) fZQuench[j] = gRandom->Exp(0.2);
943 }
944
945 quenched[j] = (fZQuench[j] > 0.01);
946 } // primary partons
947
948
949
950 Double_t pNew[1000][4];
951 Int_t kNew[1000];
952 Int_t icount = 0;
953 Double_t zquench[4];
954
955//
956// System Loop
957 for (Int_t isys = 0; isys < 4; isys++) {
958// Skip to next system if not quenched.
959 if (!quenched[isys]) continue;
960
961 nGluon[isys] = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
962 if (nGluon[isys] > fNGmax) nGluon[isys] = fNGmax;
963 zquench[isys] = 1. - TMath::Power(1. - fZQuench[isys], 1./Double_t(nGluon[isys]));
964 wjtKick[isys] = wjtKick[isys] / TMath::Sqrt(Double_t(nGluon[isys]));
965
966
967
968 Int_t igMin = -1;
969 Int_t igMax = -1;
970 Double_t pg[4] = {0., 0., 0., 0.};
971
972//
973// Loop on radiation events
974
975 for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) {
976 while (1) {
977 icount = 0;
978 for (Int_t k = 0; k < 4; k++)
979 {
980 p0[isys][k] = 0.;
981 p1[isys][k] = 0.;
982 p2[isys][k] = 0.;
983 }
984// Loop over partons
985 for (Int_t i = 0; i < numpart; i++)
986 {
987 imo = fPyjets->K[2][i];
988 kst = fPyjets->K[0][i];
989 pdg = fPyjets->K[1][i];
990
991
992
993// Quarks and gluons only
994 if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
995// Particles from hard scattering only
996
997 if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
998 Int_t imom = imo % 1000;
999 if ((isys == 0 || isys == 1) && ((imom != (isys + 7)))) continue;
1000 if ((isys == 2 || isys == 3) && ((imom != (isys + 1)))) continue;
1001
1002
1003// Skip comment lines
1004 if (kst != 1 && kst != 2) continue;
1005//
1006// Parton kinematic
1007 px = fPyjets->P[0][i];
1008 py = fPyjets->P[1][i];
1009 pz = fPyjets->P[2][i];
1010 e = fPyjets->P[3][i];
1011 m = fPyjets->P[4][i];
1012 pt = TMath::Sqrt(px * px + py * py);
1013 p = TMath::Sqrt(px * px + py * py + pz * pz);
1014 phi = TMath::Pi() + TMath::ATan2(-py, -px);
1015 theta = TMath::ATan2(pt, pz);
1016
1017//
1018// Save 4-momentum sum for balancing
1019 Int_t index = isys;
1020
1021 p0[index][0] += px;
1022 p0[index][1] += py;
1023 p0[index][2] += pz;
1024 p0[index][3] += e;
1025
1026 klast[index] = i;
1027
1028//
1029// Fractional energy loss
1030 Double_t z = zquench[index];
1031
1032
1033// Don't fully quench radiated gluons
1034//
1035 if (imo > 1000) {
1036// This small factor makes sure that the gluons are not too close in phase space to avoid recombination
1037//
1038
1039 z = 0.02;
1040 }
1041// printf("z: %d %f\n", imo, z);
1042
1043
1044//
1045
1046 //
1047 //
1048 // Transform into frame in which initial parton is along z-axis
1049 //
1050 TVector3 v(px, py, pz);
1051 v.RotateZ(-phiq[index]); v.RotateY(-thetaq[index]);
1052 Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pl = v.Z();
1053
1054 Double_t jt = TMath::Sqrt(pxs * pxs + pys * pys);
1055 Double_t mt2 = jt * jt + m * m;
1056 Double_t zmax = 1.;
1057 //
1058 // Kinematic limit on z
1059 //
1060 if (m > 0.) zmax = 1. - m / TMath::Sqrt(m * m + jt * jt);
1061 //
1062 // Change light-cone kinematics rel. to initial parton
1063 //
1064 Double_t eppzOld = e + pl;
1065 Double_t empzOld = e - pl;
1066
1067 Double_t eppzNew = (1. - z) * eppzOld;
1068 Double_t empzNew = empzOld - mt2 * z / eppzOld;
1069 Double_t eNew = 0.5 * (eppzNew + empzNew);
1070 Double_t plNew = 0.5 * (eppzNew - empzNew);
1071
1072 Double_t jtNew;
1073 //
1074 // if mt very small (or sometimes even < 0 for numerical reasons) set it to 0
1075 Double_t mt2New = eppzNew * empzNew;
1076 if (mt2New < 1.e-8) mt2New = 0.;
1077 if (z < zmax) {
1078 if (m * m > mt2New) {
1079 //
1080 // This should not happen
1081 //
1082 Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
1083 jtNew = 0;
1084 } else {
1085 jtNew = TMath::Sqrt(mt2New - m * m);
1086 }
1087 } else {
1088 // If pT is to small (probably a leading massive particle) we scale only the energy
1089 // This can cause negative masses of the radiated gluon
1090 // Let's hope for the best ...
1091 jtNew = jt;
1092 eNew = TMath::Sqrt(plNew * plNew + mt2);
1093
1094 }
1095 //
1096 // Calculate new px, py
1097 //
1044c4d8 1098 Double_t pxNew = 0;
1099 Double_t pyNew = 0;
1100
1101 if (jt > 0.) {
1102 pxNew = jtNew / jt * pxs;
1103 pyNew = jtNew / jt * pys;
1104 }
39d810c8 1105
1106// Double_t dpx = pxs - pxNew;
1107// Double_t dpy = pys - pyNew;
1108// Double_t dpz = pl - plNew;
1109// Double_t de = e - eNew;
1110// Double_t dmass2 = de * de - dpx * dpx - dpy * dpy - dpz * dpz;
1111// printf("New mass (1) %e %e %e %e %e %e %e \n", dmass2, jt, jtNew, pl, plNew, e, eNew);
1112// printf("New mass (2) %e %e \n", pxNew, pyNew);
1113 //
1114 // Rotate back
1115 //
1116 TVector3 w(pxNew, pyNew, plNew);
1117 w.RotateY(thetaq[index]); w.RotateZ(phiq[index]);
1118 pxNew = w.X(); pyNew = w.Y(); plNew = w.Z();
1119
1120 p1[index][0] += pxNew;
1121 p1[index][1] += pyNew;
1122 p1[index][2] += plNew;
1123 p1[index][3] += eNew;
1124 //
1125 // Updated 4-momentum vectors
1126 //
1127 pNew[icount][0] = pxNew;
1128 pNew[icount][1] = pyNew;
1129 pNew[icount][2] = plNew;
1130 pNew[icount][3] = eNew;
1131 kNew[icount] = i;
1132 icount++;
1133 } // parton loop
1134 //
1135 // Check if there was phase-space for quenching
1136 //
1137
1138 if (icount == 0) quenched[isys] = kFALSE;
1139 if (!quenched[isys]) break;
1140
1141 for (Int_t j = 0; j < 4; j++)
1142 {
1143 p2[isys][j] = p0[isys][j] - p1[isys][j];
1144 }
1145 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];
1146 if (p2[isys][4] > 0.) {
1147 p2[isys][4] = TMath::Sqrt(p2[isys][4]);
1148 break;
1149 } else {
1150 printf("Warning negative mass squared in system %d %f ! \n", isys, zquench[isys]);
1151 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]);
1152 if (p2[isys][4] < -0.01) {
1153 printf("Negative mass squared !\n");
1154 // Here we have to put the gluon back to mass shell
1155 // This will lead to a small energy imbalance
1156 p2[isys][4] = 0.;
1157 p2[isys][3] = TMath::Sqrt(p2[isys][0] * p2[isys][0] + p2[isys][1] * p2[isys][1] + p2[isys][2] * p2[isys][2]);
1158 break;
1159 } else {
1160 p2[isys][4] = 0.;
1161 break;
1162 }
1163 }
1164 /*
1165 zHeavy *= 0.98;
1166 printf("zHeavy lowered to %f\n", zHeavy);
1167 if (zHeavy < 0.01) {
1168 printf("No success ! \n");
1169 icount = 0;
1170 quenched[isys] = kFALSE;
1171 break;
1172 }
1173 */
1174 } // iteration on z (while)
1175
1176// Update event record
1177 for (Int_t k = 0; k < icount; k++) {
1178// 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] );
1179 fPyjets->P[0][kNew[k]] = pNew[k][0];
1180 fPyjets->P[1][kNew[k]] = pNew[k][1];
1181 fPyjets->P[2][kNew[k]] = pNew[k][2];
1182 fPyjets->P[3][kNew[k]] = pNew[k][3];
1183 }
1184 //
1185 // Add the gluons
1186 //
1187 Int_t ish = 0;
1188 Int_t iGlu;
1189 if (!quenched[isys]) continue;
1190//
1191// Last parton from shower i
1192 Int_t in = klast[isys];
1193//
1194// Continue if no parton in shower i selected
1195 if (in == -1) continue;
1196//
1197// If this is the second initial parton and it is behind the first move pointer by previous ish
1198 if (isys == 1 && klast[1] > klast[0]) in += ish;
1199//
1200// Starting index
1201
1202// jmin = in - 1;
1203// How many additional gluons will be generated
1204 ish = 1;
1205 if (p2[isys][4] > 0.05) ish = 2;
1206//
1207// Position of gluons
1208 iGlu = numpart;
1209 if (iglu == 0) igMin = iGlu;
1210 igMax = iGlu;
1211 numpart += ish;
1212 (fPyjets->N) += ish;
1213
1214 if (ish == 1) {
1215 fPyjets->P[0][iGlu] = p2[isys][0];
1216 fPyjets->P[1][iGlu] = p2[isys][1];
1217 fPyjets->P[2][iGlu] = p2[isys][2];
1218 fPyjets->P[3][iGlu] = p2[isys][3];
1219 fPyjets->P[4][iGlu] = p2[isys][4];
1220
1221 fPyjets->K[0][iGlu] = 1;
1222 if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu] = 1;
1223 fPyjets->K[1][iGlu] = 21;
1224 fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
1225 fPyjets->K[3][iGlu] = -1;
1226 fPyjets->K[4][iGlu] = -1;
1227
1228 pg[0] += p2[isys][0];
1229 pg[1] += p2[isys][1];
1230 pg[2] += p2[isys][2];
1231 pg[3] += p2[isys][3];
1232 } else {
1233 //
1234 // Split gluon in rest frame.
1235 //
1236 Double_t bx = p2[isys][0] / p2[isys][3];
1237 Double_t by = p2[isys][1] / p2[isys][3];
1238 Double_t bz = p2[isys][2] / p2[isys][3];
1239 Double_t pst = p2[isys][4] / 2.;
1240 //
1241 // Isotropic decay ????
1242 Double_t cost = 2. * gRandom->Rndm() - 1.;
60e55aee 1243 Double_t sint = TMath::Sqrt((1.-cost)*(1.+cost));
2ab330c9 1244 Double_t phis = 2. * TMath::Pi() * gRandom->Rndm();
39d810c8 1245
1246 Double_t pz1 = pst * cost;
1247 Double_t pz2 = -pst * cost;
1248 Double_t pt1 = pst * sint;
1249 Double_t pt2 = -pst * sint;
2ab330c9 1250 Double_t px1 = pt1 * TMath::Cos(phis);
1251 Double_t py1 = pt1 * TMath::Sin(phis);
1252 Double_t px2 = pt2 * TMath::Cos(phis);
1253 Double_t py2 = pt2 * TMath::Sin(phis);
39d810c8 1254
1255 fPyjets->P[0][iGlu] = px1;
1256 fPyjets->P[1][iGlu] = py1;
1257 fPyjets->P[2][iGlu] = pz1;
1258 fPyjets->P[3][iGlu] = pst;
1259 fPyjets->P[4][iGlu] = 0.;
1260
1261 fPyjets->K[0][iGlu] = 1 ;
1262 fPyjets->K[1][iGlu] = 21;
1263 fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
1264 fPyjets->K[3][iGlu] = -1;
1265 fPyjets->K[4][iGlu] = -1;
1266
1267 fPyjets->P[0][iGlu+1] = px2;
1268 fPyjets->P[1][iGlu+1] = py2;
1269 fPyjets->P[2][iGlu+1] = pz2;
1270 fPyjets->P[3][iGlu+1] = pst;
1271 fPyjets->P[4][iGlu+1] = 0.;
1272
1273 fPyjets->K[0][iGlu+1] = 1;
1274 if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu+1] = 1;
1275 fPyjets->K[1][iGlu+1] = 21;
1276 fPyjets->K[2][iGlu+1] = fPyjets->K[2][in] + 1000;
1277 fPyjets->K[3][iGlu+1] = -1;
1278 fPyjets->K[4][iGlu+1] = -1;
1279 SetMSTU(1,0);
1280 SetMSTU(2,0);
1281 //
1282 // Boost back
1283 //
1284 Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
1285 }
1286/*
1287 for (Int_t ig = iGlu; ig < iGlu+ish; ig++) {
1288 Double_t px, py, pz;
1289 px = fPyjets->P[0][ig];
1290 py = fPyjets->P[1][ig];
1291 pz = fPyjets->P[2][ig];
1292 TVector3 v(px, py, pz);
1293 v.RotateZ(-phiq[isys]);
1294 v.RotateY(-thetaq[isys]);
1295 Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pzs = v.Z();
1296 Double_t r = AliPythiaRndm::GetPythiaRandom()->Rndm();
1297 Double_t jtKick = 0.3 * TMath::Sqrt(-TMath::Log(r));
1298 if (ish == 2) jtKick = wjtKick[i] * TMath::Sqrt(-TMath::Log(r)) / TMath::Sqrt(2.);
1299 Double_t phiKick = 2. * TMath::Pi() * AliPythiaRndm::GetPythiaRandom()->Rndm();
1300 pxs += jtKick * TMath::Cos(phiKick);
1301 pys += jtKick * TMath::Sin(phiKick);
1302 TVector3 w(pxs, pys, pzs);
1303 w.RotateY(thetaq[isys]);
1304 w.RotateZ(phiq[isys]);
1305 fPyjets->P[0][ig] = w.X();
1306 fPyjets->P[1][ig] = w.Y();
1307 fPyjets->P[2][ig] = w.Z();
1308 fPyjets->P[2][ig] = w.Mag();
1309 }
1310*/
1311 } // kGluon
1312
1313
1314 // Check energy conservation
1315 Double_t pxs = 0.;
1316 Double_t pys = 0.;
1317 Double_t pzs = 0.;
1318 Double_t es = 14000.;
1319
1320 for (Int_t i = 0; i < numpart; i++)
1321 {
1322 kst = fPyjets->K[0][i];
1323 if (kst != 1 && kst != 2) continue;
1324 pxs += fPyjets->P[0][i];
1325 pys += fPyjets->P[1][i];
1326 pzs += fPyjets->P[2][i];
1327 es -= fPyjets->P[3][i];
1328 }
1329 if (TMath::Abs(pxs) > 1.e-2 ||
1330 TMath::Abs(pys) > 1.e-2 ||
1331 TMath::Abs(pzs) > 1.e-1) {
1332 printf("%e %e %e %e\n", pxs, pys, pzs, es);
1333// Fatal("Quench()", "4-Momentum non-conservation");
1334 }
1335
1336 } // end quenching loop (systems)
1337// Clean-up
1338 for (Int_t i = 0; i < numpart; i++)
1339 {
1340 imo = fPyjets->K[2][i];
1341 if (imo > 1000) {
1342 fPyjets->K[2][i] = fPyjets->K[2][i] % 1000;
1343 }
1344 }
1345// this->Pylist(1);
1346} // end quench
1347
1348
1349void AliPythia6::Pyquen(Double_t a, Int_t ibf, Double_t b)
1350{
1351 // Igor Lokthine's quenching routine
1352 // http://lokhtin.web.cern.ch/lokhtin/pyquen/pyquen.txt
1353
1354 pyquen(a, ibf, b);
1355}
1356
1357void AliPythia6::SetPyquenParameters(Double_t t0, Double_t tau0, Int_t nf, Int_t iengl, Int_t iangl)
1358{
1359 // Set the parameters for the PYQUEN package.
1360 // See comments in PyquenCommon.h
1361
1362
1363 PYQPAR.t0 = t0;
1364 PYQPAR.tau0 = tau0;
1365 PYQPAR.nf = nf;
1366 PYQPAR.iengl = iengl;
1367 PYQPAR.iangl = iangl;
1368}
1369
1370void AliPythia6::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1371{
1372//
1373// Load event into Pythia Common Block
1374//
1375
1376 Int_t npart = stack -> GetNprimary();
1377 Int_t n0 = 0;
1378
1379 if (!flag) {
1380 GetPyjets()->N = npart;
1381 } else {
1382 n0 = GetPyjets()->N;
1383 GetPyjets()->N = n0 + npart;
1384 }
1385
1386
1387 for (Int_t part = 0; part < npart; part++) {
1388 TParticle *mPart = stack->Particle(part);
1389
1390 Int_t kf = mPart->GetPdgCode();
1391 Int_t ks = mPart->GetStatusCode();
1392 Int_t idf = mPart->GetFirstDaughter();
1393 Int_t idl = mPart->GetLastDaughter();
1394
1395 if (reHadr) {
1396 if (ks == 11 || ks == 12) {
1397 ks -= 10;
1398 idf = -1;
1399 idl = -1;
1400 }
1401 }
1402
1403 Float_t px = mPart->Px();
1404 Float_t py = mPart->Py();
1405 Float_t pz = mPart->Pz();
1406 Float_t e = mPart->Energy();
1407 Float_t m = mPart->GetCalcMass();
1408
1409
1410 (GetPyjets())->P[0][part+n0] = px;
1411 (GetPyjets())->P[1][part+n0] = py;
1412 (GetPyjets())->P[2][part+n0] = pz;
1413 (GetPyjets())->P[3][part+n0] = e;
1414 (GetPyjets())->P[4][part+n0] = m;
1415
1416 (GetPyjets())->K[1][part+n0] = kf;
1417 (GetPyjets())->K[0][part+n0] = ks;
1418 (GetPyjets())->K[3][part+n0] = idf + 1;
1419 (GetPyjets())->K[4][part+n0] = idl + 1;
1420 (GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1421 }
1422}
1423
1424
1425void AliPythia6::Pyevnw()
1426{
1427 // New multiple interaction scenario
1428 pyevnw();
1429}
1430
1431void AliPythia6::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])
1432{
1433 // Return event specific quenching parameters
1434 xp = fXJet;
1435 yp = fYJet;
1436 for (Int_t i = 0; i < 4; i++) z[i] = fZQuench[i];
1437
1438}
1439
1440void AliPythia6::ConfigHeavyFlavor()
1441{
1442 //
1443 // Default configuration for Heavy Flavor production
1444 //
1445 // All QCD processes
1446 //
1447 SetMSEL(1);
1448
1449 // No multiple interactions
1450 SetMSTP(81,0);
1451 SetPARP(81, 0.);
1452 SetPARP(82, 0.);
1453 // Initial/final parton shower on (Pythia default)
1454 SetMSTP(61,1);
1455 SetMSTP(71,1);
1456
1457 // 2nd order alpha_s
1458 SetMSTP(2,2);
1459
1460 // QCD scales
1461 SetMSTP(32,2);
1462 SetPARP(34,1.0);
1463}
1464
1465void AliPythia6::AtlasTuning()
1466{
1467 //
1468 // Configuration for the ATLAS tuning
1469 SetMSTP(51,AliStructFuncType::PDFsetIndex(kCTEQ5L));
1470 SetMSTP(81,1); // Multiple Interactions ON
1471 SetMSTP(82,4); // Double Gaussian Model
1472 SetPARP(81,1.9); // Min. pt for multiple interactions (default in 6.2-14)
1473 SetPARP(82,1.8); // [GeV] PT_min at Ref. energy
1474 SetPARP(89,1000.); // [GeV] Ref. energy
1475 SetPARP(90,0.16); // 2*epsilon (exponent in power law)
1476 SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
1477 SetPARP(84,0.5); // Core radius
1478 SetPARP(85,0.33); // Regulates gluon prod. mechanism
1479 SetPARP(86,0.66); // Regulates gluon prod. mechanism
1480 SetPARP(67,1); // Regulates Initial State Radiation
1481}
1482
b4d68b2a 1483void AliPythia6::AtlasTuningMC09()
67b7bb0e 1484{
1485 //
1486 // Configuration for the ATLAS tuning
1487 printf("ATLAS New TUNE MC09\n");
1488 SetMSTP(81,21); // treatment for MI, ISR, FSR and beam remnants: MI on, new model
1489 SetMSTP(82, 4); // Double Gaussian Model
1490 SetMSTP(52, 2); // External PDF
1491 SetMSTP(51, 20650); // MRST LO*
1492
1493
1494 SetMSTP(70, 0); // (was 2: def manual 1, def code 0) virtuality scale for ISR
1495 SetMSTP(72, 1); // (was 0: def 1) maximum scale for FSR
1496 SetMSTP(88, 1); // (was 0: def 1) strategy for qq junction to di-quark or baryon in beam remnant
1497 SetMSTP(90, 0); // (was 1: def 0) strategy of compensate the primordial kT
1498
1499 SetPARP(78, 0.3); // the amount of color reconnection in the final state
1500 SetPARP(80, 0.1); // probability of color partons kicked out from beam remnant
1501 SetPARP(82, 2.3); // [GeV] PT_min at Ref. energy
1502 SetPARP(83, 0.8); // Core density in proton matter distribution (def.value)
1503 SetPARP(84, 0.7); // Core radius
1504 SetPARP(90, 0.25); // 2*epsilon (exponent in power law)
1505 SetPARJ(81, 0.29); // (was 0.14: def 0.29) Labmda value in running alpha_s for parton showers
1506
1507 SetMSTP(95, 6);
1508 SetPARJ(41, 0.3); // a and b parameters of the symmm. Lund FF
1509 SetPARJ(42, 0.58);
1510 SetPARJ(46, 0.75); // mod. of the Lund FF for heavy end-point quarks
1511 SetPARP(89,1800.); // [GeV] Ref. energy
1512}
1513
1cdfb96f 1514void AliPythia6::SetWeightPower(Double_t pow)
1515{
1516 setpowwght(pow);
1517 SetMSTP(142, 1); // Tell Pythia to use pyevwt to calculate event wghts
1518}
1519
39d810c8 1520void AliPythia6::SetPtHardRange(Float_t ptmin, Float_t ptmax)
1521{
1522 // Set the pt hard range
1523 SetCKIN(3, ptmin);
1524 SetCKIN(4, ptmax);
1525}
1526
1527void AliPythia6::SetYHardRange(Float_t ymin, Float_t ymax)
1528{
1529 // Set the y hard range
1530 SetCKIN(7, ymin);
1531 SetCKIN(8, ymax);
1532}
1533
1534
1535void AliPythia6::SetFragmentation(Int_t flag)
1536{
1537 // Switch fragmentation on/off
1538 SetMSTP(111, flag);
1539}
1540
1541void AliPythia6::SetInitialAndFinalStateRadiation(Int_t flag1, Int_t flag2)
1542{
1543// initial state radiation
1544 SetMSTP(61, flag1);
1545// final state radiation
1546 SetMSTP(71, flag2);
1547}
1548
1549void AliPythia6::SetIntrinsicKt(Float_t kt)
1550{
1551 // Set the inreinsic kt
1552 if (kt > 0.) {
1553 SetMSTP(91,1);
1554 SetPARP(91,kt);
1555 SetPARP(93, 4. * kt);
1556 } else {
1557 SetMSTP(91,0);
1558 }
1559}
1560
1561void AliPythia6::SwitchHFOff()
1562{
1563 // Switch off heavy flavor
1564 // Maximum number of quark flavours used in pdf
1565 SetMSTP(58, 3);
1566 // Maximum number of flavors that can be used in showers
1567 SetMSTJ(45, 3);
1568}
1569
1570void AliPythia6::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
1571 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
1572{
1573// Set pycell parameters
1574 SetPARU(51, etamax);
1575 SetMSTU(51, neta);
1576 SetMSTU(52, nphi);
1577 SetPARU(58, thresh);
1578 SetPARU(52, etseed);
1579 SetPARU(53, minet);
1580 SetPARU(54, r);
1581 SetMSTU(54, 2);
1582}
1583
1584void AliPythia6::ModifiedSplitting()
1585{
1586 // Modified splitting probability as a model for quenching
1587 SetPARJ(200, 0.8);
1588 SetMSTJ(41, 1); // QCD radiation only
1589 SetMSTJ(42, 2); // angular ordering
1590 SetMSTJ(44, 2); // option to run alpha_s
1591 SetMSTJ(47, 0); // No correction back to hard scattering element
1592 SetMSTJ(50, 0); // No coherence in first branching
1593 SetPARJ(82, 1.); // Cut off for parton showers
1594}
1595
1596void AliPythia6::SwitchHadronisationOff()
1597{
1598 // Switch off hadronisarion
1599 SetMSTJ(1, 0);
1600}
1601
1602void AliPythia6::SwitchHadronisationOn()
1603{
1604 // Switch on hadronisarion
1605 SetMSTJ(1, 1);
1606}
1607
1608
1609void AliPythia6::GetXandQ(Float_t& x1, Float_t& x2, Float_t& q)
1610{
1611 // Get x1, x2 and Q for this event
1612
1613 q = GetVINT(51);
1614 x1 = GetVINT(41);
1615 x2 = GetVINT(42);
1616}
1617
1618Float_t AliPythia6::GetXSection()
1619{
1620 // Get the total cross-section
1621 return (GetPARI(1));
1622}
1623
1624Float_t AliPythia6::GetPtHard()
1625{
1626 // Get the pT hard for this event
1627 return GetVINT(47);
1628}
1629
1630Int_t AliPythia6::ProcessCode()
1631{
1632 // Get the subprocess code
1633 return GetMSTI(1);
1634}
1635
1636void AliPythia6::PrintStatistics()
1637{
1638 // End of run statistics
1639 Pystat(1);
1640}
1641
1642void AliPythia6::EventListing()
1643{
1644 // End of run statistics
1645 Pylist(2);
1646}
1647
1648AliPythia6& AliPythia6::operator=(const AliPythia6& rhs)
1649{
1650// Assignment operator
1651 rhs.Copy(*this);
1652 return *this;
1653}
1654
1655 void AliPythia6::Copy(TObject&) const
1656{
1657 //
1658 // Copy
1659 //
1660 Fatal("Copy","Not implemented!\n");
1661}