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 #include "../FASTSIM/AliFastGlauber.h"
  21 #include "../FASTSIM/AliQuenchingWeights.h"
  22 #include "TVector3.h"
  23
  24 ClassImp(AliPythia)
  25
  26 #ifndef WIN32
  27 # define pyclus pyclus_
  28 # define pycell pycell_
  29 # define pyshow pyshow_
  30 # define pyrobo pyrobo_
  31 # define pyquen pyquen_
  32 # define pyevnw pyevnw_
  33 # define type_of_call
  34 #else
  35 # define pyclus PYCLUS
  36 # define pycell PYCELL
  37 # define pyrobo PYROBO
  38 # define pyquen PYQUEN
  39 # define pyevnw PYEVNW
  40 # define type_of_call _stdcall
  41 #endif
  42
  43 extern "C" void type_of_call pyclus(Int_t & );
  44 extern "C" void type_of_call pycell(Int_t & );
  45 extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
  46 extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
  47 extern "C" void type_of_call pyquen(Double_t &, Int_t &, Double_t &);
  48 extern "C" void type_of_call pyevnw();
  49
  50 //_____________________________________________________________________________
  51
  52 AliPythia* AliPythia::fgAliPythia=NULL;
  53
  54 AliPythia::AliPythia()
  55 {
  56 // Default Constructor
  57 //
  58 //  Set random number
  59     if (!AliPythiaRndm::GetPythiaRandom())
  60       AliPythiaRndm::SetPythiaRandom(GetRandom());
  61     fGlauber          = 0;
  62     fQuenchingWeights = 0;
  63 }
  64
  65 void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
  66 {
  67 // Initialise the process to generate
  68     if (!AliPythiaRndm::GetPythiaRandom())
  69       AliPythiaRndm::SetPythiaRandom(GetRandom());
  70
  71     fProcess = process;
  72     fEcms = energy;
  73     fStrucFunc = strucfunc;
  74 //...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
  75     SetMDCY(Pycomp(111) ,1,0);
  76     SetMDCY(Pycomp(310) ,1,0);
  77     SetMDCY(Pycomp(3122),1,0);
  78     SetMDCY(Pycomp(3112),1,0);
  79     SetMDCY(Pycomp(3212),1,0);
  80     SetMDCY(Pycomp(3222),1,0);
  81     SetMDCY(Pycomp(3312),1,0);
  82     SetMDCY(Pycomp(3322),1,0);
  83     SetMDCY(Pycomp(3334),1,0);
  84     //  select structure function
  85     SetMSTP(52,2);
  86     SetMSTP(51,strucfunc);
  87 //
  88 // Pythia initialisation for selected processes//
  89 //
  90 // Make MSEL clean
  91 //
  92     for (Int_t i=1; i<= 200; i++) {
  93         SetMSUB(i,0);
  94     }
  95 //  select charm production
  96     switch (process)
  97     {
  98     case kPyOldUEQ2ordered:  //Old underlying events with Q2 ordered QCD processes
  99 //        Multiple interactions on.
 100         SetMSTP(81,1);
 101 // Double Gaussian matter distribution.
 102         SetMSTP(82,4);
 103         SetPARP(83,0.5);
 104         SetPARP(84,0.4);
 105 //  pT0.
 106         SetPARP(82,2.0);
 107 //  Reference energy for pT0 and energy rescaling pace.
 108         SetPARP(89,1800);
 109         SetPARP(90,0.25);
 110 //  String drawing almost completely minimizes string length.
 111         SetPARP(85,0.9);
 112         SetPARP(86,0.95);
 113 // ISR and FSR activity.
 114         SetPARP(67,4);
 115         SetPARP(71,4);
 116 // Lambda_FSR scale.
 117         SetPARJ(81,0.29);
 118         break;
 119     case kPyOldUEQ2ordered2:
 120 // Old underlying events with Q2 ordered QCD processes
 121 // Multiple interactions on.
 122         SetMSTP(81,1);
 123 // Double Gaussian matter distribution.
 124         SetMSTP(82,4);
 125         SetPARP(83,0.5);
 126         SetPARP(84,0.4);
 127 // pT0.
 128         SetPARP(82,2.0);
 129 // Reference energy for pT0 and energy rescaling pace.
 130         SetPARP(89,1800);
 131         SetPARP(90,0.16);  // here is the difference with  kPyOldUEQ2ordered
 132 // String drawing almost completely minimizes string length.
 133         SetPARP(85,0.9);
 134         SetPARP(86,0.95);
 135 // ISR and FSR activity.
 136         SetPARP(67,4);
 137         SetPARP(71,4);
 138 // Lambda_FSR scale.
 139         SetPARJ(81,0.29);
 140         break;
 141     case kPyOldPopcorn:
 142 // Old production mechanism: Old Popcorn
 143         SetMSEL(1);
 144         SetMSTJ(12,3);
 145 // (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon
 146         SetMSTP(88,2);
 147 // (D=1)see can be used to form  baryons (BARYON JUNCTION)
 148         SetMSTJ(1,1);
 149         SetMSTP(51,kCTEQ5L);// CTEQ 5L        ! CTEQ5L pdf
 150         SetMSTP(81,1);      // Multiple Interactions ON
 151         SetMSTP(82,4);      // Double Gaussian Model
 152         SetPARP(82,1.8);    // [GeV]    PT_min at Ref. energy
 153         SetPARP(89,1000.);  // [GeV]   Ref. energy
 154         SetPARP(90,0.16);   // 2*epsilon (exponent in power law)
 155         SetPARP(83,0.5);    // Core density in proton matter dist. (def.value)
 156         SetPARP(84,0.5);    // Core radius
 157         SetPARP(85,0.33);   // Regulates gluon prod. mechanism
 158         SetPARP(86,0.66);   // Regulates gluon prod. mechanism
 159         SetPARP(67,1);      // Regulate gluon prod. mechanism
 160         break;
 161     case kPyCharm:
 162         SetMSEL(4);
 163 //  heavy quark masses
 164
 165         SetPMAS(4,1,1.2);
 166         SetMSTU(16,2);
 167 //
 168 //    primordial pT
 169         SetMSTP(91,1);
 170         SetPARP(91,1.);
 171         SetPARP(93,5.);
 172 //
 173         break;
 174     case kPyBeauty:
 175         SetMSEL(5);
 176         SetPMAS(5,1,4.75);
 177         SetMSTU(16,2);
 178         break;
 179     case kPyJpsi:
 180         SetMSEL(0);
 181 // gg->J/Psi g
 182         SetMSUB(86,1);
 183         break;
 184     case kPyJpsiChi:
 185         SetMSEL(0);
 186 // gg->J/Psi g
 187         SetMSUB(86,1);
 188 // gg-> chi_0c g
 189         SetMSUB(87,1);
 190 // gg-> chi_1c g
 191         SetMSUB(88,1);
 192 // gg-> chi_2c g
 193         SetMSUB(89,1);
 194         break;
 195     case kPyCharmUnforced:
 196         SetMSEL(0);
 197 // gq->qg
 198         SetMSUB(28,1);
 199 // gg->qq
 200         SetMSUB(53,1);
 201 // gg->gg
 202         SetMSUB(68,1);
 203         break;
 204     case kPyBeautyUnforced:
 205         SetMSEL(0);
 206 // gq->qg
 207         SetMSUB(28,1);
 208 // gg->qq
 209         SetMSUB(53,1);
 210 // gg->gg
 211         SetMSUB(68,1);
 212         break;
 213     case kPyMb:
 214 // Minimum Bias pp-Collisions
 215 //
 216 //
 217 //      select Pythia min. bias model
 218         SetMSEL(0);
 219         SetMSUB(92,1);             // single diffraction AB-->XB
 220         SetMSUB(93,1);             // single diffraction AB-->AX
 221         SetMSUB(94,1);             // double diffraction
 222         SetMSUB(95,1);             // low pt production
 223
 224 //
 225 // ATLAS Tuning
 226 //
 227
 228         SetMSTP(51, kCTEQ5L);      // CTEQ5L pdf
 229         SetMSTP(81,1);             // Multiple Interactions ON
 230         SetMSTP(82,4);             // Double Gaussian Model
 231
 232         SetPARP(82,1.8);           // [GeV]    PT_min at Ref. energy
 233         SetPARP(89,1000.);         // [GeV]   Ref. energy
 234         SetPARP(90,0.16);          // 2*epsilon (exponent in power law)
 235         SetPARP(83,0.5);           // Core density in proton matter distribution (def.value)
 236         SetPARP(84,0.5);           // Core radius
 237         SetPARP(85,0.33);          // Regulates gluon prod. mechanism
 238         SetPARP(86,0.66);          // Regulates gluon prod. mechanism
 239         SetPARP(67,1);             // Regulates Initial State Radiation
 240         break;
 241     case kPyMbNonDiffr:
 242 // Minimum Bias pp-Collisions
 243 //
 244 //
 245 //      select Pythia min. bias model
 246         SetMSEL(0);
 247         SetMSUB(95,1);             // low pt production
 248
 249 //
 250 // ATLAS Tuning
 251 //
 252
 253         SetMSTP(51,kCTEQ5L);       // CTEQ5L pdf
 254         SetMSTP(81,1);             // Multiple Interactions ON
 255         SetMSTP(82,4);             // Double Gaussian Model
 256
 257         SetPARP(82,1.8);           // [GeV]    PT_min at Ref. energy
 258         SetPARP(89,1000.);         // [GeV]   Ref. energy
 259         SetPARP(90,0.16);          // 2*epsilon (exponent in power law)
 260         SetPARP(83,0.5);           // Core density in proton matter distribution (def.value)
 261         SetPARP(84,0.5);           // Core radius
 262         SetPARP(85,0.33);          // Regulates gluon prod. mechanism
 263         SetPARP(86,0.66);          // Regulates gluon prod. mechanism
 264         SetPARP(67,1);             // Regulates Initial State Radiation
 265         break;
 266     case kPyJets:
 267 //
 268 //  QCD Jets
 269 //
 270         SetMSEL(1);
 271  // Pythia Tune A (CDF)
 272  //
 273        SetPARP(67,4.);            // Regulates Initial State Radiation
 274        SetMSTP(82,4);             // Double Gaussian Model
 275        SetPARP(82,2.0);           // [GeV]    PT_min at Ref. energy
 276        SetPARP(84,0.4);           // Core radius
 277        SetPARP(85,0.90) ;         // Regulates gluon prod. mechanism
 278        SetPARP(86,0.95);          // Regulates gluon prod. mechanism
 279        SetPARP(89,1800.);         // [GeV]   Ref. energy
 280        SetPARP(90,0.25);          // 2*epsilon (exponent in power law)
 281        break;
 282     case kPyDirectGamma:
 283         SetMSEL(10);
 284         break;
 285     case kPyCharmPbPbMNR:
 286     case kPyD0PbPbMNR:
 287       // Tuning of Pythia parameters aimed to get a resonable agreement
 288       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 289       // c-cbar single inclusive and double differential distributions.
 290       // This parameter settings are meant to work with Pb-Pb collisions
 291       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 292       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 293       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
 294
 295       // All QCD processes
 296       SetMSEL(1);
 297
 298       // No multiple interactions
 299       SetMSTP(81,0);
 300       SetPARP(81,0.0);
 301       SetPARP(82,0.0);
 302
 303       // Initial/final parton shower on (Pythia default)
 304       SetMSTP(61,1);
 305       SetMSTP(71,1);
 306
 307       // 2nd order alpha_s
 308       SetMSTP(2,2);
 309
 310       // QCD scales
 311       SetMSTP(32,2);
 312       SetPARP(34,1.0);
 313
 314       // Intrinsic <kT>
 315       SetMSTP(91,1);
 316       SetPARP(91,1.304);
 317       SetPARP(93,6.52);
 318
 319       // Set c-quark mass
 320       SetPMAS(4,1,1.2);
 321
 322       break;
 323     case kPyDPlusPbPbMNR:
 324       // Tuning of Pythia parameters aimed to get a resonable agreement
 325       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 326       // c-cbar single inclusive and double differential distributions.
 327       // This parameter settings are meant to work with Pb-Pb collisions
 328       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 329       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 330       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
 331
 332       // All QCD processes
 333       SetMSEL(1);
 334
 335       // No multiple interactions
 336       SetMSTP(81,0);
 337       SetPARP(81,0.0);
 338       SetPARP(82,0.0);
 339
 340       // Initial/final parton shower on (Pythia default)
 341       SetMSTP(61,1);
 342       SetMSTP(71,1);
 343
 344       // 2nd order alpha_s
 345       SetMSTP(2,2);
 346
 347       // QCD scales
 348       SetMSTP(32,2);
 349       SetPARP(34,1.0);
 350
 351       // Intrinsic <kT>
 352       SetMSTP(91,1);
 353       SetPARP(91,1.304);
 354       SetPARP(93,6.52);
 355
 356       // Set c-quark mass
 357       SetPMAS(4,1,1.2);
 358
 359       break;
 360     case kPyCharmpPbMNR:
 361     case kPyD0pPbMNR:
 362       // Tuning of Pythia parameters aimed to get a resonable agreement
 363       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 364       // c-cbar single inclusive and double differential distributions.
 365       // This parameter settings are meant to work with p-Pb collisions
 366       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 367       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 368       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
 369
 370       // All QCD processes
 371       SetMSEL(1);
 372
 373       // No multiple interactions
 374       SetMSTP(81,0);
 375       SetPARP(81,0.0);
 376       SetPARP(82,0.0);
 377
 378       // Initial/final parton shower on (Pythia default)
 379       SetMSTP(61,1);
 380       SetMSTP(71,1);
 381
 382       // 2nd order alpha_s
 383       SetMSTP(2,2);
 384
 385       // QCD scales
 386       SetMSTP(32,2);
 387       SetPARP(34,1.0);
 388
 389       // Intrinsic <kT>
 390       SetMSTP(91,1);
 391       SetPARP(91,1.16);
 392       SetPARP(93,5.8);
 393
 394       // Set c-quark mass
 395       SetPMAS(4,1,1.2);
 396
 397       break;
 398     case kPyDPluspPbMNR:
 399       // Tuning of Pythia parameters aimed to get a resonable agreement
 400       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 401       // c-cbar single inclusive and double differential distributions.
 402       // This parameter settings are meant to work with p-Pb collisions
 403       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 404       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 405       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
 406
 407       // All QCD processes
 408       SetMSEL(1);
 409
 410       // No multiple interactions
 411       SetMSTP(81,0);
 412       SetPARP(81,0.0);
 413       SetPARP(82,0.0);
 414
 415       // Initial/final parton shower on (Pythia default)
 416       SetMSTP(61,1);
 417       SetMSTP(71,1);
 418
 419       // 2nd order alpha_s
 420       SetMSTP(2,2);
 421
 422       // QCD scales
 423       SetMSTP(32,2);
 424       SetPARP(34,1.0);
 425
 426       // Intrinsic <kT>
 427       SetMSTP(91,1);
 428       SetPARP(91,1.16);
 429       SetPARP(93,5.8);
 430
 431       // Set c-quark mass
 432       SetPMAS(4,1,1.2);
 433
 434       break;
 435     case kPyCharmppMNR:
 436     case kPyD0ppMNR:
 437       // Tuning of Pythia parameters aimed to get a resonable agreement
 438       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 439       // c-cbar single inclusive and double differential distributions.
 440       // This parameter settings are meant to work with pp collisions
 441       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 442       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 443       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
 444
 445       // All QCD processes
 446       SetMSEL(1);
 447
 448       // No multiple interactions
 449       SetMSTP(81,0);
 450       SetPARP(81,0.0);
 451       SetPARP(82,0.0);
 452
 453       // Initial/final parton shower on (Pythia default)
 454       SetMSTP(61,1);
 455       SetMSTP(71,1);
 456
 457       // 2nd order alpha_s
 458       SetMSTP(2,2);
 459
 460       // QCD scales
 461       SetMSTP(32,2);
 462       SetPARP(34,1.0);
 463
 464       // Intrinsic <kT^2>
 465       SetMSTP(91,1);
 466       SetPARP(91,1.);
 467       SetPARP(93,5.);
 468
 469       // Set c-quark mass
 470       SetPMAS(4,1,1.2);
 471
 472       break;
 473     case kPyDPlusppMNR:
 474       // Tuning of Pythia parameters aimed to get a resonable agreement
 475       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 476       // c-cbar single inclusive and double differential distributions.
 477       // This parameter settings are meant to work with pp collisions
 478       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 479       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 480       // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
 481
 482       // All QCD processes
 483       SetMSEL(1);
 484
 485       // No multiple interactions
 486       SetMSTP(81,0);
 487       SetPARP(81,0.0);
 488       SetPARP(82,0.0);
 489
 490       // Initial/final parton shower on (Pythia default)
 491       SetMSTP(61,1);
 492       SetMSTP(71,1);
 493
 494       // 2nd order alpha_s
 495       SetMSTP(2,2);
 496
 497       // QCD scales
 498       SetMSTP(32,2);
 499       SetPARP(34,1.0);
 500
 501       // Intrinsic <kT^2>
 502       SetMSTP(91,1);
 503       SetPARP(91,1.);
 504       SetPARP(93,5.);
 505
 506       // Set c-quark mass
 507       SetPMAS(4,1,1.2);
 508
 509       break;
 510     case kPyBeautyPbPbMNR:
 511       // Tuning of Pythia parameters aimed to get a resonable agreement
 512       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 513       // b-bbar single inclusive and double differential distributions.
 514       // This parameter settings are meant to work with Pb-Pb collisions
 515       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 516       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 517       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
 518
 519       // All QCD processes
 520       SetMSEL(1);
 521
 522       // No multiple interactions
 523       SetMSTP(81,0);
 524       SetPARP(81,0.0);
 525       SetPARP(82,0.0);
 526
 527       // Initial/final parton shower on (Pythia default)
 528       SetMSTP(61,1);
 529       SetMSTP(71,1);
 530
 531       // 2nd order alpha_s
 532       SetMSTP(2,2);
 533
 534       // QCD scales
 535       SetMSTP(32,2);
 536       SetPARP(34,1.0);
 537       SetPARP(67,1.0);
 538       SetPARP(71,1.0);
 539
 540       // Intrinsic <kT>
 541       SetMSTP(91,1);
 542       SetPARP(91,2.035);
 543       SetPARP(93,10.17);
 544
 545       // Set b-quark mass
 546       SetPMAS(5,1,4.75);
 547
 548       break;
 549     case kPyBeautypPbMNR:
 550       // Tuning of Pythia parameters aimed to get a resonable agreement
 551       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 552       // b-bbar single inclusive and double differential distributions.
 553       // This parameter settings are meant to work with p-Pb collisions
 554       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 555       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 556       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
 557
 558       // All QCD processes
 559       SetMSEL(1);
 560
 561       // No multiple interactions
 562       SetMSTP(81,0);
 563       SetPARP(81,0.0);
 564       SetPARP(82,0.0);
 565
 566       // Initial/final parton shower on (Pythia default)
 567       SetMSTP(61,1);
 568       SetMSTP(71,1);
 569
 570       // 2nd order alpha_s
 571       SetMSTP(2,2);
 572
 573       // QCD scales
 574       SetMSTP(32,2);
 575       SetPARP(34,1.0);
 576       SetPARP(67,1.0);
 577       SetPARP(71,1.0);
 578
 579       // Intrinsic <kT>
 580       SetMSTP(91,1);
 581       SetPARP(91,1.60);
 582       SetPARP(93,8.00);
 583
 584       // Set b-quark mass
 585       SetPMAS(5,1,4.75);
 586
 587       break;
 588     case kPyBeautyppMNR:
 589       // Tuning of Pythia parameters aimed to get a resonable agreement
 590       // between with the NLO calculation by Mangano, Nason, Ridolfi for the
 591       // b-bbar single inclusive and double differential distributions.
 592       // This parameter settings are meant to work with pp collisions
 593       // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
 594       // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
 595       // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
 596
 597       // All QCD processes
 598       SetMSEL(1);
 599
 600       // No multiple interactions
 601       SetMSTP(81,0);
 602       SetPARP(81,0.0);
 603       SetPARP(82,0.0);
 604
 605       // Initial/final parton shower on (Pythia default)
 606       SetMSTP(61,1);
 607       SetMSTP(71,1);
 608
 609       // 2nd order alpha_s
 610       SetMSTP(2,2);
 611
 612       // QCD scales
 613       SetMSTP(32,2);
 614       SetPARP(34,1.0);
 615       SetPARP(67,1.0);
 616       SetPARP(71,1.0);
 617
 618       // Intrinsic <kT>
 619       SetMSTP(91,1);
 620       SetPARP(91,1.);
 621       SetPARP(93,5.);
 622
 623       // Set b-quark mass
 624       SetPMAS(5,1,4.75);
 625
 626       break;
 627
 628     case kPyW:
 629
 630       //Inclusive production of W+/-
 631       SetMSEL(0);
 632       //f fbar -> W+
 633       SetMSUB(2,1);
 634       //        //f fbar -> g W+
 635       //        SetMSUB(16,1);
 636       //        //f fbar -> gamma W+
 637       //        SetMSUB(20,1);
 638       //        //f g -> f W+
 639       //        SetMSUB(31,1);
 640       //        //f gamma -> f W+
 641       //        SetMSUB(36,1);
 642
 643       // Initial/final parton shower on (Pythia default)
 644       // With parton showers on we are generating "W inclusive process"
 645       SetMSTP(61,1); //Initial QCD & QED showers on
 646       SetMSTP(71,1); //Final QCD & QED showers on
 647
 648       break;
 649
 650     case kPyZ:
 651
 652       //Inclusive production of Z
 653       SetMSEL(0);
 654       //f fbar -> Z/gamma
 655       SetMSUB(1,1);
 656
 657       //       // f fbar -> g Z/gamma
 658       //       SetMSUB(15,1);
 659       //       // f fbar -> gamma Z/gamma
 660       //       SetMSUB(19,1);
 661       //       // f g -> f Z/gamma
 662       //       SetMSUB(30,1);
 663       //       // f gamma -> f Z/gamma
 664       //       SetMSUB(35,1);
 665
 666       //only Z included, not gamma
 667       SetMSTP(43,2);
 668
 669       // Initial/final parton shower on (Pythia default)
 670       // With parton showers on we are generating "Z inclusive process"
 671       SetMSTP(61,1); //Initial QCD & QED showers on
 672       SetMSTP(71,1); //Final QCD & QED showers on
 673
 674       break;
 675
 676     }
 677 //
 678 //  Initialize PYTHIA
 679     SetMSTP(41,1);   // all resonance decays switched on
 680
 681     Initialize("CMS","p","p",fEcms);
 682
 683 }
 684
 685 Int_t AliPythia::CheckedLuComp(Int_t kf)
 686 {
 687 // Check Lund particle code (for debugging)
 688     Int_t kc=Pycomp(kf);
 689     printf("\n Lucomp kf,kc %d %d",kf,kc);
 690     return kc;
 691 }
 692
 693 void AliPythia::SetNuclei(Int_t a1, Int_t a2)
 694 {
 695 // Treat protons as inside nuclei with mass numbers a1 and a2
 696 //    The MSTP array in the PYPARS common block is used to enable and
 697 //    select the nuclear structure functions.
 698 //    MSTP(52)  : (D=1) choice of proton and nuclear structure-function library
 699 //            =1: internal PYTHIA acording to MSTP(51)
 700 //            =2: PDFLIB proton  s.f., with MSTP(51)  = 1000xNGROUP+NSET
 701 //    If the following mass number both not equal zero, nuclear corrections of the stf are used.
 702 //    MSTP(192) : Mass number of nucleus side 1
 703 //    MSTP(193) : Mass number of nucleus side 2
 704     SetMSTP(52,2);
 705     SetMSTP(192, a1);
 706     SetMSTP(193, a2);
 707 }
 708
 709
 710 AliPythia* AliPythia::Instance()
 711 {
 712 // Set random number generator
 713     if (fgAliPythia) {
 714         return fgAliPythia;
 715     } else {
 716         fgAliPythia = new AliPythia();
 717         return fgAliPythia;
 718     }
 719 }
 720
 721 void AliPythia::PrintParticles()
 722 {
 723 // Print list of particl properties
 724     Int_t np = 0;
 725     char*   name = new char[16];
 726     for (Int_t kf=0; kf<1000000; kf++) {
 727         for (Int_t c = 1;  c > -2; c-=2) {
 728             Int_t kc = Pycomp(c*kf);
 729             if (kc) {
 730                 Float_t mass  = GetPMAS(kc,1);
 731                 Float_t width = GetPMAS(kc,2);
 732                 Float_t tau   = GetPMAS(kc,4);
 733
 734                 Pyname(kf,name);
 735
 736                 np++;
 737
 738                 printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
 739                        c*kf, name, mass, width, tau);
 740             }
 741         }
 742     }
 743     printf("\n Number of particles %d \n \n", np);
 744 }
 745
 746 void  AliPythia::ResetDecayTable()
 747 {
 748 //  Set default values for pythia decay switches
 749     Int_t i;
 750     for (i = 1; i <  501; i++) SetMDCY(i,1,fDefMDCY[i]);
 751     for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
 752 }
 753
 754 void  AliPythia::SetDecayTable()
 755 {
 756 //  Set default values for pythia decay switches
 757 //
 758     Int_t i;
 759     for (i = 1; i <  501; i++) fDefMDCY[i] = GetMDCY(i,1);
 760     for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
 761 }
 762
 763 void  AliPythia::Pyclus(Int_t& njet)
 764 {
 765 //  Call Pythia clustering algorithm
 766 //
 767     pyclus(njet);
 768 }
 769
 770 void  AliPythia::Pycell(Int_t& njet)
 771 {
 772 //  Call Pythia jet reconstruction algorithm
 773 //
 774     pycell(njet);
 775 }
 776
 777 void  AliPythia::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
 778 {
 779 //  Call Pythia jet reconstruction algorithm
 780 //
 781     pyshow(ip1, ip2, qmax);
 782 }
 783
 784 void AliPythia::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
 785 {
 786     pyrobo(imi, ima, the, phi, bex, bey, bez);
 787 }
 788
 789
 790
 791 void AliPythia::InitQuenching(Float_t cMin, Float_t cMax, Float_t k, Int_t iECMethod)
 792 {
 793 // Initializes
 794 // (1) The quenching model using quenching weights according to C. Salgado and U. Wiedemann
 795 // (2) The nuclear geometry using the Glauber Model
 796 //
 797
 798
 799     fGlauber = new AliFastGlauber();
 800     fGlauber->Init(2);
 801     fGlauber->SetCentralityClass(cMin, cMax);
 802
 803     fQuenchingWeights = new AliQuenchingWeights();
 804     fQuenchingWeights->InitMult();
 805     fQuenchingWeights->SetK(k);
 806     fQuenchingWeights->SetECMethod(AliQuenchingWeights::kECMethod(iECMethod));
 807 }
 808
 809
 810 void  AliPythia::Quench()
 811 {
 812 //
 813 //
 814 //  Simple Jet Quenching routine:
 815 //  =============================
 816 //  The jet formed by all final state partons radiated by the parton created
 817 //  in the hard collisions is quenched by a factor (1-z) using light cone variables in
 818 //  the initial parton reference frame:
 819 //  (E + p_z)new = (1-z) (E + p_z)old
 820 //
 821 //
 822 //
 823 //
 824 //   The lost momentum is first balanced by one gluon with virtuality > 0.
 825 //   Subsequently the gluon splits to yield two gluons with E = p.
 826 //
 827 //
 828 //
 829     static Float_t eMean = 0.;
 830     static Int_t   icall = 0;
 831
 832     Double_t p0[4][5];
 833     Double_t p1[4][5];
 834     Double_t p2[4][5];
 835     Int_t   klast[4] = {-1, -1, -1, -1};
 836
 837     Int_t numpart   = fPyjets->N;
 838     Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0., phi = 0.;
 839     Double_t pxq[4], pyq[4], pzq[4], eq[4], yq[4], mq[4], pq[4], phiq[4], thetaq[4], ptq[4];
 840     Bool_t  quenched[4];
 841     Double_t wjtKick[4];
 842     Int_t nGluon[4];
 843     Int_t qPdg[4];
 844     Int_t   imo, kst, pdg;
 845
 846 //
 847 //  Sore information about Primary partons
 848 //
 849 //  j =
 850 //  0, 1 partons from hard scattering
 851 //  2, 3 partons from initial state radiation
 852 //
 853     for (Int_t i = 2; i <= 7; i++) {
 854         Int_t j = 0;
 855         // Skip gluons that participate in hard scattering
 856         if (i == 4 || i == 5) continue;
 857         // Gluons from hard Scattering
 858         if (i == 6 || i == 7) {
 859             j = i - 6;
 860             pxq[j]    = fPyjets->P[0][i];
 861             pyq[j]    = fPyjets->P[1][i];
 862             pzq[j]    = fPyjets->P[2][i];
 863             eq[j]     = fPyjets->P[3][i];
 864             mq[j]     = fPyjets->P[4][i];
 865         } else {
 866             // Gluons from initial state radiation
 867             //
 868             // Obtain 4-momentum vector from difference between original parton and parton after gluon
 869             // radiation. Energy is calculated independently because initial state radition does not
 870             // conserve strictly momentum and energy for each partonic system independently.
 871             //
 872             // Not very clean. Should be improved !
 873             //
 874             //
 875             j = i;
 876             pxq[j]    = fPyjets->P[0][i] - fPyjets->P[0][i+2];
 877             pyq[j]    = fPyjets->P[1][i] - fPyjets->P[1][i+2];
 878             pzq[j]    = fPyjets->P[2][i] - fPyjets->P[2][i+2];
 879             mq[j]     = fPyjets->P[4][i];
 880             eq[j]     = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j] + mq[j] * mq[j]);
 881         }
 882 //
 883 //  Calculate some kinematic variables
 884 //
 885         yq[j]     = 0.5 * TMath::Log((eq[j] + pzq[j] + 1.e-14) / (eq[j] - pzq[j] + 1.e-14));
 886         pq[j]     = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j]);
 887         phiq[j]   = TMath::Pi()+TMath::ATan2(-pyq[j], -pxq[j]);
 888         ptq[j]    = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j]);
 889         thetaq[j] = TMath::ATan2(ptq[j], pzq[j]);
 890         qPdg[j]   =  fPyjets->K[1][i];
 891     }
 892
 893     Double_t int0[4];
 894     Double_t int1[4];
 895
 896     fGlauber->GetI0I1ForPythiaAndXY(4, phiq, int0, int1, fXJet, fYJet, 15.);
 897
 898     for (Int_t j = 0; j < 4; j++) {
 899         //
 900         // Quench only central jets and with E > 10.
 901         //
 902
 903
 904         Int_t itype = (qPdg[j] == 21) ? 2 : 1;
 905         Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]);
 906
 907         if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
 908             fZQuench[j] = 0.;
 909         } else {
 910             if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) {
 911                 icall ++;
 912                 eMean += eloss;
 913             }
 914             //
 915             // Extra pt
 916             Double_t l =   fQuenchingWeights->CalcLk(int0[j], int1[j]);
 917             wjtKick[j] = TMath::Sqrt(l *  fQuenchingWeights->CalcQk(int0[j], int1[j]));
 918             //
 919             // Fractional energy loss
 920             fZQuench[j] = eloss / eq[j];
 921             //
 922             // Avoid complete loss
 923             //
 924             if (fZQuench[j] == 1.) fZQuench[j] = 0.95;
 925             //
 926             // Some debug printing
 927
 928
 929 //          printf("Initial parton # %3d, Type %3d Energy %10.3f Phi %10.3f Length %10.3f Loss %10.3f Kick %10.3f Mean: %10.3f %10.3f\n",
 930 //                 j, itype, eq[j], phiq[j], l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]);
 931
 932 //          fZQuench[j] = 0.8;
 933 //          while (fZQuench[j] >= 0.95)  fZQuench[j] = gRandom->Exp(0.2);
 934         }
 935
 936         quenched[j] = (fZQuench[j] > 0.01);
 937     } // primary partons
 938
 939
 940
 941     Double_t pNew[1000][4];
 942     Int_t    kNew[1000];
 943     Int_t icount = 0;
 944     Double_t zquench[4];
 945
 946 //
 947 //  System Loop
 948     for (Int_t isys = 0; isys < 4; isys++) {
 949 //      Skip to next system if not quenched.
 950         if (!quenched[isys]) continue;
 951
 952         nGluon[isys]   = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
 953         if (nGluon[isys] > 6) nGluon[isys] = 6;
 954         zquench[isys] = 1. - TMath::Power(1. - fZQuench[isys], 1./Double_t(nGluon[isys]));
 955         wjtKick[isys]  = wjtKick[isys] / TMath::Sqrt(Double_t(nGluon[isys]));
 956
 957
 958
 959         Int_t igMin = -1;
 960         Int_t igMax = -1;
 961         Double_t pg[4] = {0., 0., 0., 0.};
 962
 963 //
 964 // Loop on radiation events
 965
 966         for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) {
 967             while (1) {
 968                 icount = 0;
 969                 for (Int_t k = 0; k < 4; k++)
 970                 {
 971                     p0[isys][k] = 0.;
 972                     p1[isys][k] = 0.;
 973                     p2[isys][k] = 0.;
 974                 }
 975 //      Loop over partons
 976                 for (Int_t i = 0; i < numpart; i++)
 977                 {
 978                     imo =  fPyjets->K[2][i];
 979                     kst =  fPyjets->K[0][i];
 980                     pdg =  fPyjets->K[1][i];
 981
 982
 983
 984 //      Quarks and gluons only
 985                     if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
 986 //      Particles from hard scattering only
 987
 988                     if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
 989                     Int_t imom = imo % 1000;
 990                     if ((isys == 0 || isys == 1) && ((imom != (isys + 7)))) continue;
 991                     if ((isys == 2 || isys == 3) && ((imom != (isys + 1)))) continue;
 992
 993
 994 //      Skip comment lines
 995                     if (kst != 1 && kst != 2) continue;
 996 //
 997 //      Parton kinematic
 998                     px    = fPyjets->P[0][i];
 999                     py    = fPyjets->P[1][i];
1000                     pz    = fPyjets->P[2][i];
1001                     e     = fPyjets->P[3][i];
1002                     m     = fPyjets->P[4][i];
1003                     pt    = TMath::Sqrt(px * px + py * py);
1004                     p     = TMath::Sqrt(px * px + py * py + pz * pz);
1005                     phi   = TMath::Pi() + TMath::ATan2(-py, -px);
1006                     theta = TMath::ATan2(pt, pz);
1007
1008 //
1009 //      Save 4-momentum sum for balancing
1010                     Int_t index = isys;
1011
1012                     p0[index][0] += px;
1013                     p0[index][1] += py;
1014                     p0[index][2] += pz;
1015                     p0[index][3] += e;
1016
1017                     klast[index] = i;
1018
1019 //
1020 //      Fractional energy loss
1021                     Double_t z = zquench[index];
1022
1023
1024 //      Don't fully quench radiated gluons
1025 //
1026                     if (imo > 1000) {
1027 //      This small factor makes sure that the gluons are not too close in phase space to avoid recombination
1028 //
1029
1030                         z = 0.02;
1031                     }
1032 //                  printf("z: %d %f\n", imo, z);
1033
1034
1035 //
1036
1037                     //
1038                     //
1039                     //      Transform into frame in which initial parton is along z-axis
1040                     //
1041                     TVector3 v(px, py, pz);
1042                     v.RotateZ(-phiq[index]);  v.RotateY(-thetaq[index]);
1043                     Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pl  = v.Z();
1044
1045                     Double_t jt  = TMath::Sqrt(pxs * pxs + pys * pys);
1046                     Double_t mt2 = jt * jt + m * m;
1047                     Double_t zmax = 1.;
1048                     //
1049                     // Kinematic limit on z
1050                     //
1051                     if (m > 0.) zmax = 1. - m / TMath::Sqrt(m * m + jt * jt);
1052                     //
1053                     // Change light-cone kinematics rel. to initial parton
1054                     //
1055                     Double_t eppzOld = e + pl;
1056                     Double_t empzOld = e - pl;
1057
1058                     Double_t eppzNew = (1. - z) * eppzOld;
1059                     Double_t empzNew = empzOld - mt2 * z / eppzOld;
1060                     Double_t eNew    = 0.5 * (eppzNew + empzNew);
1061                     Double_t plNew   = 0.5 * (eppzNew - empzNew);
1062
1063                     Double_t jtNew;
1064                     //
1065                     // if mt very small (or sometimes even < 0 for numerical reasons) set it to 0
1066                     Double_t mt2New = eppzNew * empzNew;
1067                     if (mt2New < 1.e-8) mt2New = 0.;
1068                     if (z < zmax) {
1069                         if (m * m > mt2New) {
1070                             //
1071                             // This should not happen
1072                             //
1073                             Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
1074                             jtNew = 0;
1075                         } else {
1076                             jtNew    = TMath::Sqrt(mt2New - m * m);
1077                         }
1078                     } else {
1079                         // If pT is to small (probably a leading massive particle) we scale only the energy
1080                         // This can cause negative masses of the radiated gluon
1081                         // Let's hope for the best ...
1082                         jtNew = jt;
1083                         eNew  = TMath::Sqrt(plNew * plNew + mt2);
1084
1085                     }
1086                     //
1087                     //     Calculate new px, py
1088                     //
1089                     Double_t pxNew   = jtNew / jt * pxs;
1090                     Double_t pyNew   = jtNew / jt * pys;
1091
1092 //                  Double_t dpx = pxs - pxNew;
1093 //                  Double_t dpy = pys - pyNew;
1094 //                  Double_t dpz = pl  - plNew;
1095 //                  Double_t de  = e   - eNew;
1096 //                  Double_t dmass2 = de * de  - dpx * dpx - dpy * dpy - dpz * dpz;
1097 //                  printf("New mass (1) %e %e %e %e %e %e %e \n", dmass2, jt, jtNew, pl, plNew, e, eNew);
1098 //                  printf("New mass (2) %e %e \n", pxNew, pyNew);
1099                     //
1100                     //      Rotate back
1101                     //
1102                     TVector3 w(pxNew, pyNew, plNew);
1103                     w.RotateY(thetaq[index]); w.RotateZ(phiq[index]);
1104                     pxNew = w.X(); pyNew = w.Y(); plNew = w.Z();
1105
1106                     p1[index][0] += pxNew;
1107                     p1[index][1] += pyNew;
1108                     p1[index][2] += plNew;
1109                     p1[index][3] += eNew;
1110                     //
1111                     // Updated 4-momentum vectors
1112                     //
1113                     pNew[icount][0]  = pxNew;
1114                     pNew[icount][1]  = pyNew;
1115                     pNew[icount][2]  = plNew;
1116                     pNew[icount][3]  = eNew;
1117                     kNew[icount]     = i;
1118                     icount++;
1119                 } // parton loop
1120                 //
1121                 // Check if there was phase-space for quenching
1122                 //
1123
1124                 if (icount == 0) quenched[isys] = kFALSE;
1125                 if (!quenched[isys]) break;
1126
1127                 for (Int_t j = 0; j < 4; j++)
1128                 {
1129                     p2[isys][j] = p0[isys][j] - p1[isys][j];
1130                 }
1131                 p2[isys][4] = p2[isys][3] * p2[isys][3] - p2[isys][0] * p2[isys][0] - p2[isys][1] * p2[isys][1] - p2[isys][2] * p2[isys][2];
1132                 if (p2[isys][4] > 0.) {
1133                     p2[isys][4] = TMath::Sqrt(p2[isys][4]);
1134                     break;
1135                 } else {
1136                     printf("Warning negative mass squared in system %d %f ! \n", isys, zquench[isys]);
1137                     printf("4-Momentum: %10.3e %10.3e %10.3e %10.3e %10.3e \n", p2[isys][0], p2[isys][1], p2[isys][2], p2[isys][3], p2[isys][4]);
1138                     if (p2[isys][4] < -0.01) {
1139                         printf("Negative mass squared !\n");
1140                         // Here we have to put the gluon back to mass shell
1141                         // This will lead to a small energy imbalance
1142                         p2[isys][4]  = 0.;
1143                         p2[isys][3]  = TMath::Sqrt(p2[isys][0] * p2[isys][0] + p2[isys][1] * p2[isys][1] + p2[isys][2] * p2[isys][2]);
1144                         break;
1145                     } else {
1146                         p2[isys][4] = 0.;
1147                         break;
1148                     }
1149                 }
1150                 /*
1151                 zHeavy *= 0.98;
1152                 printf("zHeavy lowered to %f\n", zHeavy);
1153                 if (zHeavy < 0.01) {
1154                     printf("No success ! \n");
1155                     icount = 0;
1156                     quenched[isys] = kFALSE;
1157                     break;
1158                 }
1159                 */
1160             } // iteration on z (while)
1161
1162 //          Update  event record
1163             for (Int_t k = 0; k < icount; k++) {
1164 //              printf("%6d %6d %10.3e %10.3e %10.3e %10.3e\n", k, kNew[k], pNew[k][0],pNew[k][1], pNew[k][2], pNew[k][3] );
1165                 fPyjets->P[0][kNew[k]] = pNew[k][0];
1166                 fPyjets->P[1][kNew[k]] = pNew[k][1];
1167                 fPyjets->P[2][kNew[k]] = pNew[k][2];
1168                 fPyjets->P[3][kNew[k]] = pNew[k][3];
1169             }
1170             //
1171             // Add the gluons
1172             //
1173             Int_t ish = 0;
1174             Int_t iGlu;
1175             if (!quenched[isys]) continue;
1176 //
1177 //      Last parton from shower i
1178             Int_t in = klast[isys];
1179 //
1180 //      Continue if no parton in shower i selected
1181             if (in == -1) continue;
1182 //
1183 //      If this is the second initial parton and it is behind the first move pointer by previous ish
1184             if (isys == 1 && klast[1] > klast[0]) in += ish;
1185 //
1186 //      Starting index
1187
1188 //          jmin = in - 1;
1189 // How many additional gluons will be generated
1190             ish  = 1;
1191             if (p2[isys][4] > 0.05) ish = 2;
1192 //
1193 //      Position of gluons
1194             iGlu = numpart;
1195             if (iglu == 0) igMin = iGlu;
1196             igMax = iGlu;
1197             numpart += ish;
1198             (fPyjets->N) += ish;
1199
1200             if (ish == 1) {
1201                 fPyjets->P[0][iGlu] = p2[isys][0];
1202                 fPyjets->P[1][iGlu] = p2[isys][1];
1203                 fPyjets->P[2][iGlu] = p2[isys][2];
1204                 fPyjets->P[3][iGlu] = p2[isys][3];
1205                 fPyjets->P[4][iGlu] = p2[isys][4];
1206
1207                 fPyjets->K[0][iGlu] = 1;
1208                 if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu] = 1;
1209                 fPyjets->K[1][iGlu] = 21;
1210                 fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
1211                 fPyjets->K[3][iGlu] = -1;
1212                 fPyjets->K[4][iGlu] = -1;
1213
1214                 pg[0] += p2[isys][0];
1215                 pg[1] += p2[isys][1];
1216                 pg[2] += p2[isys][2];
1217                 pg[3] += p2[isys][3];
1218             } else {
1219                 //
1220                 // Split gluon in rest frame.
1221                 //
1222                 Double_t bx   =  p2[isys][0] / p2[isys][3];
1223                 Double_t by   =  p2[isys][1] / p2[isys][3];
1224                 Double_t bz   =  p2[isys][2] / p2[isys][3];
1225                 Double_t pst  =  p2[isys][4] / 2.;
1226                 //
1227                 // Isotropic decay ????
1228                 Double_t cost = 2. * gRandom->Rndm() - 1.;
1229                 Double_t sint = TMath::Sqrt(1. - cost * cost);
1230                 Double_t phi =  2. * TMath::Pi() * gRandom->Rndm();
1231
1232                 Double_t pz1 =   pst * cost;
1233                 Double_t pz2 =  -pst * cost;
1234                 Double_t pt1 =   pst * sint;
1235                 Double_t pt2 =  -pst * sint;
1236                 Double_t px1 =   pt1 * TMath::Cos(phi);
1237                 Double_t py1 =   pt1 * TMath::Sin(phi);
1238                 Double_t px2 =   pt2 * TMath::Cos(phi);
1239                 Double_t py2 =   pt2 * TMath::Sin(phi);
1240
1241                 fPyjets->P[0][iGlu] = px1;
1242                 fPyjets->P[1][iGlu] = py1;
1243                 fPyjets->P[2][iGlu] = pz1;
1244                 fPyjets->P[3][iGlu] = pst;
1245                 fPyjets->P[4][iGlu] = 0.;
1246
1247                 fPyjets->K[0][iGlu] = 1 ;
1248                 fPyjets->K[1][iGlu] = 21;
1249                 fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
1250                 fPyjets->K[3][iGlu] = -1;
1251                 fPyjets->K[4][iGlu] = -1;
1252
1253                 fPyjets->P[0][iGlu+1] = px2;
1254                 fPyjets->P[1][iGlu+1] = py2;
1255                 fPyjets->P[2][iGlu+1] = pz2;
1256                 fPyjets->P[3][iGlu+1] = pst;
1257                 fPyjets->P[4][iGlu+1] = 0.;
1258
1259                 fPyjets->K[0][iGlu+1] = 1;
1260                 if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu+1] = 1;
1261                 fPyjets->K[1][iGlu+1] = 21;
1262                 fPyjets->K[2][iGlu+1] = fPyjets->K[2][in] + 1000;
1263                 fPyjets->K[3][iGlu+1] = -1;
1264                 fPyjets->K[4][iGlu+1] = -1;
1265                 SetMSTU(1,0);
1266                 SetMSTU(2,0);
1267                 //
1268                 // Boost back
1269                 //
1270                 Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
1271             }
1272 /*
1273             for (Int_t ig = iGlu; ig < iGlu+ish; ig++) {
1274                 Double_t px, py, pz;
1275                 px = fPyjets->P[0][ig];
1276                 py = fPyjets->P[1][ig];
1277                 pz = fPyjets->P[2][ig];
1278                 TVector3 v(px, py, pz);
1279                 v.RotateZ(-phiq[isys]);
1280                 v.RotateY(-thetaq[isys]);
1281                 Double_t pxs     = v.X(); Double_t pys = v.Y(); Double_t pzs  = v.Z();
1282                 Double_t r       = AliPythiaRndm::GetPythiaRandom()->Rndm();
1283                 Double_t jtKick  = 0.3 * TMath::Sqrt(-TMath::Log(r));
1284                 if (ish == 2)   jtKick  = wjtKick[i] * TMath::Sqrt(-TMath::Log(r)) / TMath::Sqrt(2.);
1285                 Double_t phiKick = 2. * TMath::Pi() * AliPythiaRndm::GetPythiaRandom()->Rndm();
1286                 pxs += jtKick * TMath::Cos(phiKick);
1287                 pys += jtKick * TMath::Sin(phiKick);
1288                 TVector3 w(pxs, pys, pzs);
1289                 w.RotateY(thetaq[isys]);
1290                 w.RotateZ(phiq[isys]);
1291                 fPyjets->P[0][ig] = w.X();
1292                 fPyjets->P[1][ig] = w.Y();
1293                 fPyjets->P[2][ig] = w.Z();
1294                 fPyjets->P[2][ig] = w.Mag();
1295             }
1296 */
1297         } // kGluon
1298
1299
1300     // Check energy conservation
1301         Double_t pxs = 0.;
1302         Double_t pys = 0.;
1303         Double_t pzs = 0.;
1304         Double_t es  = 14000.;
1305
1306         for (Int_t i = 0; i < numpart; i++)
1307         {
1308             kst =  fPyjets->K[0][i];
1309             if (kst != 1 && kst != 2) continue;
1310             pxs += fPyjets->P[0][i];
1311             pys += fPyjets->P[1][i];
1312             pzs += fPyjets->P[2][i];
1313             es  -= fPyjets->P[3][i];
1314         }
1315         if (TMath::Abs(pxs) > 1.e-2 ||
1316             TMath::Abs(pys) > 1.e-2 ||
1317             TMath::Abs(pzs) > 1.e-1) {
1318             printf("%e %e %e %e\n", pxs, pys, pzs, es);
1319 //              Fatal("Quench()", "4-Momentum non-conservation");
1320         }
1321
1322     } // end quenching loop (systems)
1323 // Clean-up
1324     for (Int_t i = 0; i < numpart; i++)
1325     {
1326         imo =  fPyjets->K[2][i];
1327         if (imo > 1000) {
1328             fPyjets->K[2][i] = fPyjets->K[2][i] % 1000;
1329         }
1330     }
1331 //      this->Pylist(1);
1332 } // end quench
1333
1334
1335 void AliPythia::Pyquen(Double_t a, Int_t ibf, Double_t b)
1336 {
1337     // Igor Lokthine's quenching routine
1338     pyquen(a, ibf, b);
1339 }
1340
1341 void AliPythia::Pyevnw()
1342 {
1343     // New multiple interaction scenario
1344     pyevnw();
1345 }
1346
1347 void AliPythia::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])
1348 {
1349     // Return event specific quenching parameters
1350     xp = fXJet;
1351     yp = fYJet;
1352     for (Int_t i = 0; i < 4; i++) z[i] = fZQuench[i];
1353
1354 }
1355
1356