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