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
 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;
 160     case kPyMbNonDiffr:
 161 // Minimum Bias pp-Collisions
 162 //
 163 //
 164 //      select Pythia min. bias model
 165         SetMSEL(0);
 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
 180         break;
 181     case kPyJets:
 182 //
 183 //  QCD Jets
 184 //
 185         SetMSEL(1);
 186         break;
 187     case kPyDirectGamma:
 188         SetMSEL(10);
 189         break;
 190     case kPyCharmPbPbMNR:
 191     case kPyD0PbPbMNR:
 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
 196       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 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
 219       // Intrinsic <kT>
 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;
 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:
 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
 334       // Intrinsic <kT>
 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
 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
 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
 430 Int_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
 438 void 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
 455 AliPythia* 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
 466 void AliPythia::PrintParticles()
 467 {
 468 // Print list of particl properties
 469     Int_t np = 0;
 470     char*   name = new char[16];
 471     for (Int_t kf=0; kf<1000000; kf++) {
 472         for (Int_t c = 1;  c > -2; c-=2) {
 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);
 478
 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
 491 void  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
 499 void  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
 508 void  AliPythia::Pyclus(Int_t& njet)
 509 {
 510 //  Call Pythia clustering algorithm
 511 //
 512     pyclus(njet);
 513 }
 514
 515 void  AliPythia::Pycell(Int_t& njet)
 516 {
 517 //  Call Pythia jet reconstruction algorithm
 518 //
 519     pycell(njet);
 520 }
 521
 522 void  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
 540 void 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
 547 void  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;
 610         p0[index][3] += e;
 611 //
 612 // Fix z
 613 //
 614
 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