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