798e76a9e8230c73aece941d82d2dfb1cbdaaa0d
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia.cxx
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$ */
17
18 #include "AliPythia.h"
19 #include "AliPythiaRndm.h"
20
21 ClassImp(AliPythia)
22
23 #ifndef WIN32
24 # define pyclus pyclus_
25 # define pycell pycell_
26 # define pyshow pyshow_
27 # define pyrobo pyrobo_
28 # define type_of_call
29 #else
30 # define pyclus PYCLUS
31 # define pycell PYCELL
32 # define pyrobo PYROBO
33 # define type_of_call _stdcall
34 #endif
35
36 extern "C" void type_of_call pyclus(Int_t & );
37 extern "C" void type_of_call pycell(Int_t & );
38 extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
39 extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
40
41 //_____________________________________________________________________________
42
43 AliPythia* AliPythia::fgAliPythia=NULL;
44
45 AliPythia::AliPythia()
46 {
47 // Default Constructor
48 //
49 //  Set random number
50     if (!AliPythiaRndm::GetPythiaRandom()) 
51       AliPythiaRndm::SetPythiaRandom(GetRandom());
52
53 }
54
55 void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
56 {
57 // Initialise the process to generate 
58     if (!AliPythiaRndm::GetPythiaRandom()) 
59       AliPythiaRndm::SetPythiaRandom(GetRandom());
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;
173     case kPyCharmPbPbMNR:
174     case kPyD0PbPbMNR:
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
179       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
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
202       // Intrinsic <kT>
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;
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:
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
317       // Intrinsic <kT>
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;
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;
404     }
405 //
406 //  Initialize PYTHIA
407     SetMSTP(41,1);   // all resonance decays switched on
408
409     Initialize("CMS","p","p",fEcms);
410
411 }
412
413 Int_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
421 void 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
438 AliPythia* 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
449 void AliPythia::PrintParticles()
450
451 // Print list of particl properties
452     Int_t np = 0;
453     char*   name = new char[16];    
454     for (Int_t kf=0; kf<1000000; kf++) {
455         for (Int_t c = 1;  c > -2; c-=2) {
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);
461
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
474 void  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
482 void  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
491 void  AliPythia::Pyclus(Int_t& njet)
492 {
493 //  Call Pythia clustering algorithm
494 //
495     pyclus(njet);
496 }
497
498 void  AliPythia::Pycell(Int_t& njet)
499 {
500 //  Call Pythia jet reconstruction algorithm
501 //
502     pycell(njet);
503 }
504
505 void  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
523 void 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
530 void  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