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