Use new pythia version.
[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"
8d2cd130 20
21ClassImp(AliPythia)
22
23#ifndef WIN32
24# define pyclus pyclus_
25# define pycell pycell_
452af8c7 26# define pyshow pyshow_
27# define pyrobo pyrobo_
8d2cd130 28# define type_of_call
29#else
30# define pyclus PYCLUS
31# define pycell PYCELL
452af8c7 32# define pyrobo PYROBO
8d2cd130 33# define type_of_call _stdcall
34#endif
35
36extern "C" void type_of_call pyclus(Int_t & );
37extern "C" void type_of_call pycell(Int_t & );
452af8c7 38extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
39extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
8d2cd130 40
41//_____________________________________________________________________________
42
43AliPythia* AliPythia::fgAliPythia=NULL;
44
45AliPythia::AliPythia()
46{
47// Default Constructor
48//
49// Set random number
7cdba479 50 if (!AliPythiaRndm::GetPythiaRandom())
51 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 52
53}
54
55void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
56{
57// Initialise the process to generate
7cdba479 58 if (!AliPythiaRndm::GetPythiaRandom())
59 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 60
61 fProcess = process;
62 fEcms = energy;
63 fStrucFunc = strucfunc;
64// don't decay p0
65 SetMDCY(Pycomp(111),1,0);
66// select structure function
67 SetMSTP(52,2);
68 SetMSTP(51,strucfunc);
69//
70// Pythia initialisation for selected processes//
71//
72// Make MSEL clean
73//
74 for (Int_t i=1; i<= 200; i++) {
75 SetMSUB(i,0);
76 }
77// select charm production
78 switch (process)
79 {
80 case kPyCharm:
81 SetMSEL(4);
82//
83// heavy quark masses
84
85 SetPMAS(4,1,1.2);
86 SetMSTU(16,2);
87//
88// primordial pT
89 SetMSTP(91,1);
90 SetPARP(91,1.);
91 SetPARP(93,5.);
92//
93 break;
94 case kPyBeauty:
95 SetMSEL(5);
96 SetPMAS(5,1,4.75);
97 SetMSTU(16,2);
98 break;
99 case kPyJpsi:
100 SetMSEL(0);
101// gg->J/Psi g
102 SetMSUB(86,1);
103 break;
104 case kPyJpsiChi:
105 SetMSEL(0);
106// gg->J/Psi g
107 SetMSUB(86,1);
108// gg-> chi_0c g
109 SetMSUB(87,1);
110// gg-> chi_1c g
111 SetMSUB(88,1);
112// gg-> chi_2c g
113 SetMSUB(89,1);
114 break;
115 case kPyCharmUnforced:
116 SetMSEL(0);
117// gq->qg
118 SetMSUB(28,1);
119// gg->qq
120 SetMSUB(53,1);
121// gg->gg
122 SetMSUB(68,1);
123 break;
124 case kPyBeautyUnforced:
125 SetMSEL(0);
126// gq->qg
127 SetMSUB(28,1);
128// gg->qq
129 SetMSUB(53,1);
130// gg->gg
131 SetMSUB(68,1);
132 break;
133 case kPyMb:
134// Minimum Bias pp-Collisions
135//
136//
137// select Pythia min. bias model
138 SetMSEL(0);
139 SetMSUB(92,1); // single diffraction AB-->XB
140 SetMSUB(93,1); // single diffraction AB-->AX
141 SetMSUB(94,1); // double diffraction
142 SetMSUB(95,1); // low pt production
143 SetMSTP(81,1); // multiple interactions switched on
144 SetMSTP(82,3); // model with varying impact param. & a single Gaussian
145 SetPARP(82,3.47); // set value pT_0 for turn-off of the cross section of
146 // multiple interaction at a reference energy = 14000 GeV
147 SetPARP(89,14000.); // reference energy for the above parameter
148 SetPARP(90,0.174); // set exponent for energy dependence of pT_0
149 case kPyMbNonDiffr:
150// Minimum Bias pp-Collisions
151//
152//
153// select Pythia min. bias model
154 SetMSEL(0);
155 SetMSUB(95,1); // low pt production
156 SetMSTP(81,1); // multiple interactions switched on
157 SetMSTP(82,3); // model with varying impact param. & a single Gaussian
158 SetPARP(82,3.47); // set value pT_0 for turn-off of the cross section of
159 // multiple interaction at a reference energy = 14000 GeV
160 SetPARP(89,14000.); // reference energy for the above parameter
161 SetPARP(90,0.174); // set exponent for energy dependence of pT_0
162
163 break;
164 case kPyJets:
165//
166// QCD Jets
167//
168 SetMSEL(1);
169 break;
170 case kPyDirectGamma:
171 SetMSEL(10);
172 break;
adf4d898 173 case kPyCharmPbPbMNR:
174 case kPyD0PbPbMNR:
8d2cd130 175 // Tuning of Pythia parameters aimed to get a resonable agreement
176 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
177 // c-cbar single inclusive and double differential distributions.
178 // This parameter settings are meant to work with Pb-Pb collisions
adf4d898 179 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
8d2cd130 180 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
181 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
182
183 // All QCD processes
184 SetMSEL(1);
185
186 // No multiple interactions
187 SetMSTP(81,0);
188 SetPARP(81,0.0);
189 SetPARP(82,0.0);
190
191 // Initial/final parton shower on (Pythia default)
192 SetMSTP(61,1);
193 SetMSTP(71,1);
194
195 // 2nd order alpha_s
196 SetMSTP(2,2);
197
198 // QCD scales
199 SetMSTP(32,2);
200 SetPARP(34,1.0);
201
adf4d898 202 // Intrinsic <kT>
8d2cd130 203 SetMSTP(91,1);
204 SetPARP(91,1.304);
205 SetPARP(93,6.52);
206
207 // Set c-quark mass
208 SetPMAS(4,1,1.2);
209
210 break;
adf4d898 211 case kPyCharmpPbMNR:
212 case kPyD0pPbMNR:
213 // Tuning of Pythia parameters aimed to get a resonable agreement
214 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
215 // c-cbar single inclusive and double differential distributions.
216 // This parameter settings are meant to work with p-Pb collisions
217 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
218 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
219 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
220
221 // All QCD processes
222 SetMSEL(1);
223
224 // No multiple interactions
225 SetMSTP(81,0);
226 SetPARP(81,0.0);
227 SetPARP(82,0.0);
228
229 // Initial/final parton shower on (Pythia default)
230 SetMSTP(61,1);
231 SetMSTP(71,1);
232
233 // 2nd order alpha_s
234 SetMSTP(2,2);
235
236 // QCD scales
237 SetMSTP(32,2);
238 SetPARP(34,1.0);
239
240 // Intrinsic <kT>
241 SetMSTP(91,1);
242 SetPARP(91,1.16);
243 SetPARP(93,5.8);
244
245 // Set c-quark mass
246 SetPMAS(4,1,1.2);
247
248 break;
249 case kPyCharmppMNR:
250 case kPyD0ppMNR:
251 // Tuning of Pythia parameters aimed to get a resonable agreement
252 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
253 // c-cbar single inclusive and double differential distributions.
254 // This parameter settings are meant to work with pp collisions
255 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
256 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
257 // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
258
259 // All QCD processes
260 SetMSEL(1);
261
262 // No multiple interactions
263 SetMSTP(81,0);
264 SetPARP(81,0.0);
265 SetPARP(82,0.0);
266
267 // Initial/final parton shower on (Pythia default)
268 SetMSTP(61,1);
269 SetMSTP(71,1);
270
271 // 2nd order alpha_s
272 SetMSTP(2,2);
273
274 // QCD scales
275 SetMSTP(32,2);
276 SetPARP(34,1.0);
277
278 // Intrinsic <kT^2>
279 SetMSTP(91,1);
280 SetPARP(91,1.);
281 SetPARP(93,5.);
282
283 // Set c-quark mass
284 SetPMAS(4,1,1.2);
285
286 break;
287 case kPyBeautyPbPbMNR:
8d2cd130 288 // Tuning of Pythia parameters aimed to get a resonable agreement
289 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
290 // b-bbar single inclusive and double differential distributions.
291 // This parameter settings are meant to work with Pb-Pb 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.75GeV. Example in ConfigBeautyPPR.C.
295
296 // All QCD processes
297 SetMSEL(1);
298
299 // No multiple interactions
300 SetMSTP(81,0);
301 SetPARP(81,0.0);
302 SetPARP(82,0.0);
303
304 // Initial/final parton shower on (Pythia default)
305 SetMSTP(61,1);
306 SetMSTP(71,1);
307
308 // 2nd order alpha_s
309 SetMSTP(2,2);
310
311 // QCD scales
312 SetMSTP(32,2);
313 SetPARP(34,1.0);
314 SetPARP(67,1.0);
315 SetPARP(71,1.0);
316
adf4d898 317 // Intrinsic <kT>
8d2cd130 318 SetMSTP(91,1);
319 SetPARP(91,2.035);
320 SetPARP(93,10.17);
321
322 // Set b-quark mass
323 SetPMAS(5,1,4.75);
324
325 break;
adf4d898 326 case kPyBeautypPbMNR:
327 // Tuning of Pythia parameters aimed to get a resonable agreement
328 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
329 // b-bbar single inclusive and double differential distributions.
330 // This parameter settings are meant to work with p-Pb collisions
331 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
332 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
333 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
334
335 // All QCD processes
336 SetMSEL(1);
337
338 // No multiple interactions
339 SetMSTP(81,0);
340 SetPARP(81,0.0);
341 SetPARP(82,0.0);
342
343 // Initial/final parton shower on (Pythia default)
344 SetMSTP(61,1);
345 SetMSTP(71,1);
346
347 // 2nd order alpha_s
348 SetMSTP(2,2);
349
350 // QCD scales
351 SetMSTP(32,2);
352 SetPARP(34,1.0);
353 SetPARP(67,1.0);
354 SetPARP(71,1.0);
355
356 // Intrinsic <kT>
357 SetMSTP(91,1);
358 SetPARP(91,1.60);
359 SetPARP(93,8.00);
360
361 // Set b-quark mass
362 SetPMAS(5,1,4.75);
363
364 break;
365 case kPyBeautyppMNR:
366 // Tuning of Pythia parameters aimed to get a resonable agreement
367 // between with the NLO calculation by Mangano, Nason, Ridolfi for the
368 // b-bbar single inclusive and double differential distributions.
369 // This parameter settings are meant to work with pp collisions
370 // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
371 // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
372 // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
373
374 // All QCD processes
375 SetMSEL(1);
376
377 // No multiple interactions
378 SetMSTP(81,0);
379 SetPARP(81,0.0);
380 SetPARP(82,0.0);
381
382 // Initial/final parton shower on (Pythia default)
383 SetMSTP(61,1);
384 SetMSTP(71,1);
385
386 // 2nd order alpha_s
387 SetMSTP(2,2);
388
389 // QCD scales
390 SetMSTP(32,2);
391 SetPARP(34,1.0);
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);
402
403 break;
8d2cd130 404 }
405//
406// Initialize PYTHIA
407 SetMSTP(41,1); // all resonance decays switched on
408
409 Initialize("CMS","p","p",fEcms);
410
411}
412
413Int_t AliPythia::CheckedLuComp(Int_t kf)
414{
415// Check Lund particle code (for debugging)
416 Int_t kc=Pycomp(kf);
417 printf("\n Lucomp kf,kc %d %d",kf,kc);
418 return kc;
419}
420
421void AliPythia::SetNuclei(Int_t a1, Int_t a2)
422{
423// Treat protons as inside nuclei with mass numbers a1 and a2
424// The MSTP array in the PYPARS common block is used to enable and
425// select the nuclear structure functions.
426// MSTP(52) : (D=1) choice of proton and nuclear structure-function library
427// =1: internal PYTHIA acording to MSTP(51)
428// =2: PDFLIB proton s.f., with MSTP(51) = 1000xNGROUP+NSET
429// If the following mass number both not equal zero, nuclear corrections of the stf are used.
430// MSTP(192) : Mass number of nucleus side 1
431// MSTP(193) : Mass number of nucleus side 2
432 SetMSTP(52,2);
433 SetMSTP(192, a1);
434 SetMSTP(193, a2);
435}
436
437
438AliPythia* AliPythia::Instance()
439{
440// Set random number generator
441 if (fgAliPythia) {
442 return fgAliPythia;
443 } else {
444 fgAliPythia = new AliPythia();
445 return fgAliPythia;
446 }
447}
448
449void AliPythia::PrintParticles()
450{
451// Print list of particl properties
452 Int_t np = 0;
c31f1d37 453 char* name = new char[16];
8d2cd130 454 for (Int_t kf=0; kf<1000000; kf++) {
455 for (Int_t c = 1; c > -2; c-=2) {
8d2cd130 456 Int_t kc = Pycomp(c*kf);
457 if (kc) {
458 Float_t mass = GetPMAS(kc,1);
459 Float_t width = GetPMAS(kc,2);
460 Float_t tau = GetPMAS(kc,4);
c31f1d37 461
8d2cd130 462 Pyname(kf,name);
463
464 np++;
465
466 printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
467 c*kf, name, mass, width, tau);
468 }
469 }
470 }
471 printf("\n Number of particles %d \n \n", np);
472}
473
474void AliPythia::ResetDecayTable()
475{
476// Set default values for pythia decay switches
477 Int_t i;
478 for (i = 1; i < 501; i++) SetMDCY(i,1,fDefMDCY[i]);
479 for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
480}
481
482void AliPythia::SetDecayTable()
483{
484// Set default values for pythia decay switches
485//
486 Int_t i;
487 for (i = 1; i < 501; i++) fDefMDCY[i] = GetMDCY(i,1);
488 for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
489}
490
491void AliPythia::Pyclus(Int_t& njet)
492{
493// Call Pythia clustering algorithm
494//
495 pyclus(njet);
496}
497
498void AliPythia::Pycell(Int_t& njet)
499{
500// Call Pythia jet reconstruction algorithm
501//
502 pycell(njet);
503}
504
452af8c7 505void AliPythia::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
506{
507// Call Pythia jet reconstruction algorithm
508//
509 Int_t numpart = fPyjets->N;
510 for (Int_t i = 0; i < numpart; i++)
511 {
512 if (fPyjets->K[2][i] == 7) ip1 = i+1;
513 if (fPyjets->K[2][i] == 8) ip2 = i+1;
514 }
515
516
517 qmax = 2. * GetVINT(51);
518 printf("Pyshow %d %d %f", ip1, ip2, qmax);
519
520 pyshow(ip1, ip2, qmax);
521}
522
523void AliPythia::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
524{
525 pyrobo(imi, ima, the, phi, bex, bey, bez);
526}
527
528
529
530void AliPythia::Quench()
531{
532//
533//
534// Simple Jet Quenching routine:
535// =============================
536// The jet formed by all final state partons radiated by the parton created
537// in the hard collisions is quenched by a factor z using:
538// (E + p_z)new = (1-z) (E + p_z)old
539//
540// The lost momentum is first balanced by one gluon with virtuality > 0.
541// Subsequently the gluon splits to yield two gluons with E = p.
542//
543 Float_t p0[2][5];
544 Float_t p1[2][5];
545 Float_t p2[2][5];
546 Int_t klast[2] = {-1, -1};
547 Int_t kglu[2];
548 for (Int_t i = 0; i < 4; i++)
549 {
550 p0[0][i] = 0.;
551 p0[1][i] = 0.;
552 p1[0][i] = 0.;
553 p1[1][i] = 0.;
554 p2[0][i] = 0.;
555 p2[1][i] = 0.;
556 }
557
558 Int_t numpart = fPyjets->N;
559
560 for (Int_t i = 0; i < numpart; i++)
561 {
562 Int_t imo = fPyjets->K[2][i];
563 Int_t kst = fPyjets->K[0][i];
564 Int_t pdg = fPyjets->K[1][i];
565
566// Quarks and gluons only
567 if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
568
569// Particles from hard scattering only
570
571
572 Float_t px = fPyjets->P[0][i];
573 Float_t py = fPyjets->P[1][i];
574 Float_t pz = fPyjets->P[2][i];
575 Float_t e = fPyjets->P[3][i];
576 Float_t m = fPyjets->P[4][i];
577 Float_t pt = TMath::Sqrt(px * px + py * py);
578// Skip comment lines
579 if (kst != 1 && kst != 2) continue;
580
581 Float_t mt = TMath::Sqrt(px * px + py * py + m * m);
582
583 //
584 // Some cuts to be in a save kinematic region
585 //
586 if (imo != 7 && imo != 8) continue;
587 Int_t index = imo - 7;
588 klast[index] = i;
589
590 p0[index][0] += px;
591 p0[index][1] += py;
592 p0[index][2] += pz;
593 p0[index][3] += e;
594
595 //
596 // Fix z
597 //
598
599 Float_t z = 0.2;
600 Float_t eppzOld = e + pz;
601 Float_t empzOld = e - pz;
602
603
604 //
605 // Kinematics of the original parton
606 //
607
608 Float_t eppzNew = (1. - z) * eppzOld;
609 Float_t empzNew = empzOld - mt * mt * z / eppzOld;
610 Float_t eNew0 = 0.5 * (eppzNew + empzNew);
611 Float_t pzNew0 = 0.5 * (eppzNew - empzNew);
612 //
613 // Skip if pt too small
614 //
615 if (m * m > eppzNew * empzNew) continue;
616 Float_t ptNew = TMath::Sqrt(eppzNew * empzNew - m * m);
617 Float_t pxNew0 = ptNew / pt * px;
618 Float_t pyNew0 = ptNew / pt * py;
619
620 p1[index][0] += pxNew0;
621 p1[index][1] += pyNew0;
622 p1[index][2] += pzNew0;
623 p1[index][3] += eNew0;
624 //
625 // Update event record
626 //
627 fPyjets->P[0][i] = pxNew0;
628 fPyjets->P[1][i] = pyNew0;
629 fPyjets->P[2][i] = pzNew0;
630 fPyjets->P[3][i] = eNew0;
631
632 }
633
634 //
635 // Gluons
636 //
637
638 for (Int_t k = 0; k < 2; k++)
639 {
640 for (Int_t j = 0; j < 4; j++)
641 {
642 p2[k][j] = p0[k][j] - p1[k][j];
643 }
644 p2[k][4] = p2[k][3] * p2[k][3] - p2[k][0] * p2[k][0] - p2[k][1] * p2[k][1] - p2[k][2] * p2[k][2];
645
646 if (p2[k][4] > 0.)
647 {
648
649 //
650 // Bring gluon back to mass shell via momentum scaling
651 // (momentum will not be conserved, but energy)
652 //
653 // not used anymore
654/*
655 Float_t psq = p2[k][0] * p2[k][0] + p2[k][1] * p2[k][1] + p2[k][2] * p2[k][2];
656 Float_t fact = TMath::Sqrt(1. + p2[k][4] / psq);
657 p2[k][0] *= fact;
658 p2[k][1] *= fact;
659 p2[k][2] *= fact;
660 p2[k][3] = TMath::Sqrt(psq) * fact;
661 p2[k][4] = 0.;
662*/
663 }
664 }
665
666 if (p2[0][4] > 0.) {
667 p2[0][4] = TMath::Sqrt(p2[0][4]);
668 } else {
669 printf("Warning negative mass squared ! \n");
670 }
671
672 if (p2[1][4] > 0.) {
673 p2[1][4] = TMath::Sqrt(p2[1][4]);
674 } else {
675 printf("Warning negative mass squared ! \n");
676 }
677
678 //
679 // Add the gluons
680 //
681
682
683 for (Int_t i = 0; i < 2; i++) {
684 Int_t ish, jmin, jmax, iGlu, iNew;
685 Int_t in = klast[i];
686 ish = 0;
687
688 if (in == -1) continue;
689 if (i == 1 && klast[1] > klast[0]) in += ish;
690
691 jmin = in - 1;
692 ish = 1;
693
694 if (p2[i][4] > 0) ish = 2;
695
696 iGlu = in;
697 iNew = in + ish;
698 jmax = numpart + ish - 1;
699
700 if (fPyjets->K[0][in-1] == 1 || fPyjets->K[0][in-1] == 21 || fPyjets->K[0][in-1] == 11) {
701 jmin = in;
702 iGlu = in + 1;
703 iNew = in;
704 }
705
706 kglu[i] = iGlu;
707
708 for (Int_t j = jmax; j > jmin; j--)
709 {
710 for (Int_t k = 0; k < 5; k++) {
711 fPyjets->K[k][j] = fPyjets->K[k][j-ish];
712 fPyjets->P[k][j] = fPyjets->P[k][j-ish];
713 fPyjets->V[k][j] = fPyjets->V[k][j-ish];
714 }
715 } // end shifting
716 numpart += ish;
717 (fPyjets->N) += ish;
718
719 if (ish == 1) {
720 fPyjets->P[0][iGlu] = p2[i][0];
721 fPyjets->P[1][iGlu] = p2[i][1];
722 fPyjets->P[2][iGlu] = p2[i][2];
723 fPyjets->P[3][iGlu] = p2[i][3];
724 fPyjets->P[4][iGlu] = p2[i][4];
725
726 fPyjets->K[0][iGlu] = 2;
727 fPyjets->K[1][iGlu] = 21;
728 fPyjets->K[2][iGlu] = fPyjets->K[2][iNew];
729 fPyjets->K[3][iGlu] = -1;
730 fPyjets->K[4][iGlu] = -1;
731 } else {
732 //
733 // Split gluon in rest frame.
734 //
735 Double_t bx = p2[i][0] / p2[i][3];
736 Double_t by = p2[i][1] / p2[i][3];
737 Double_t bz = p2[i][2] / p2[i][3];
738
739 Float_t pst = p2[i][4] / 2.;
740
741 Float_t cost = 2. * gRandom->Rndm() - 1.;
742 Float_t sint = TMath::Sqrt(1. - cost * cost);
743 Float_t phi = 2. * TMath::Pi() * gRandom->Rndm();
744
745 Float_t pz1 = pst * cost;
746 Float_t pz2 = -pst * cost;
747 Float_t pt1 = pst * sint;
748 Float_t pt2 = -pst * sint;
749 Float_t px1 = pt1 * TMath::Cos(phi);
750 Float_t py1 = pt1 * TMath::Sin(phi);
751 Float_t px2 = pt2 * TMath::Cos(phi);
752 Float_t py2 = pt2 * TMath::Sin(phi);
753
754 fPyjets->P[0][iGlu] = px1;
755 fPyjets->P[1][iGlu] = py1;
756 fPyjets->P[2][iGlu] = pz1;
757 fPyjets->P[3][iGlu] = pst;
758 fPyjets->P[4][iGlu] = 0.;
759
760 fPyjets->K[0][iGlu] = 2;
761 fPyjets->K[1][iGlu] = 21;
762 fPyjets->K[2][iGlu] = fPyjets->K[2][iNew];
763 fPyjets->K[3][iGlu] = -1;
764 fPyjets->K[4][iGlu] = -1;
765
766 fPyjets->P[0][iGlu+1] = px2;
767 fPyjets->P[1][iGlu+1] = py2;
768 fPyjets->P[2][iGlu+1] = pz2;
769 fPyjets->P[3][iGlu+1] = pst;
770 fPyjets->P[4][iGlu+1] = 0.;
771
772 fPyjets->K[0][iGlu+1] = 2;
773 fPyjets->K[1][iGlu+1] = 21;
774 fPyjets->K[2][iGlu+1] = fPyjets->K[2][iNew];
775 fPyjets->K[3][iGlu+1] = -1;
776 fPyjets->K[4][iGlu+1] = -1;
777 SetMSTU(1,0);
778 SetMSTU(2,0);
779
780 //
781 // Boost back
782 //
783 Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
784
785 }
786 } // end adding gluons
787} // end quench
788
789