PYTHIA6/AliPythia.cxx

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