PYTHIA8/AliDecayerPythia8.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 // Implementation of AliDecayer using Pythia8
  19 // Author: andreas.morsch@cern.ch
  20 #include <TMath.h>
  21 #include <TPDGCode.h>
  22 #include <TPythia8.h>
  23 #include "AliDecayerPythia8.h"
  24 #include "ParticleData.h"
  25
  26 ClassImp(AliDecayerPythia8)
  27
  28 AliDecayerPythia8::AliDecayerPythia8():
  29     TPythia8Decayer(),
  30     fDecay(kAll),
  31     fHeavyFlavour(kTRUE)
  32 {
  33     // Constructor
  34 }
  35
  36 void AliDecayerPythia8::ForceDecay()
  37 {
  38   //
  39 // Force a particle decay mode
  40 // Switch heavy flavour production off if requested
  41     if (!fHeavyFlavour) SwitchOffHeavyFlavour();
  42 //
  43     Decay_t decay = fDecay;
  44     TPythia8::Instance()->ReadString("HadronLevel:Decay = on");
  45
  46     if (decay == kNoDecayHeavy) return;
  47
  48 //
  49 // select mode
  50     Int_t products[2];
  51     Int_t mult[2];
  52     Int_t products1[3];
  53     Int_t mult1[3];
  54
  55     switch (decay)
  56     {
  57     case kHardMuons:
  58         products1[0] =     13;
  59         products1[1] =    443;
  60         products1[2] = 100443;
  61         mult1[0] = 1;
  62         mult1[1] = 1;
  63         mult1[2] = 1;
  64         ForceParticleDecay(  511, products1, mult1, 3);
  65         ForceParticleDecay(  521, products1, mult1, 3);
  66         ForceParticleDecay(  531, products1, mult1, 3);
  67         ForceParticleDecay( 5122, products1, mult1, 3);
  68         ForceParticleDecay( 5132, products1, mult1, 3);
  69         ForceParticleDecay( 5232, products1, mult1, 3);
  70         ForceParticleDecay( 5332, products1, mult1, 3);
  71         ForceParticleDecay( 100443, 443, 1);  // Psi'  -> J/Psi X
  72         ForceParticleDecay(    443,  13, 2);  // J/Psi -> mu+ mu-
  73         ForceParticleDecay(  411,13,1); // D+/-
  74         ForceParticleDecay(  421,13,1); // D0
  75         ForceParticleDecay(  431,13,1); // D_s
  76         ForceParticleDecay( 4122,13,1); // Lambda_c
  77         ForceParticleDecay( 4132,13,1); // Xsi_c
  78         ForceParticleDecay( 4232,13,1); // Sigma_c
  79         ForceParticleDecay( 4332,13,1); // Omega_c
  80         break;
  81    case kChiToJpsiGammaToMuonMuon:
  82         products[0] =    443;
  83         products[1] =     22;
  84         mult[0] = 1;
  85         mult[1] = 1;
  86         ForceParticleDecay( 20443, products, mult, 2);   // Chi_1c  -> J/Psi  Gamma
  87         ForceParticleDecay(   445, products, mult, 2);   // Chi_2c  -> J/Psi  Gamma
  88         ForceParticleDecay(   443, 13, 2);               // J/Psi -> mu+ mu-
  89         break;
  90     case kChiToJpsiGammaToElectronElectron:
  91         products[0] =    443;
  92         products[1] =     22;
  93         mult[0] = 1;
  94         mult[1] = 1;
  95         ForceParticleDecay( 20443, products, mult, 2);   // Chi_1c  -> J/Psi  Gamma
  96         ForceParticleDecay(   445, products, mult, 2);   // Chi_2c  -> J/Psi  Gamma
  97         ForceParticleDecay(   443, 11, 2);               // J/Psi -> e+ e-
  98         break;
  99
 100     case kBSemiMuonic:
 101         ForceParticleDecay(  511,13,1); // B0
 102         ForceParticleDecay(  521,13,1); // B+/-
 103         ForceParticleDecay(  531,13,1); // B_s
 104         ForceParticleDecay( 5122,13,1); // Lambda_b
 105         ForceParticleDecay( 5132,13,1); // Xsi_b
 106         ForceParticleDecay( 5232,13,1); // Sigma_b
 107         ForceParticleDecay( 5332,13,1); // Omega_b
 108         break;
 109     case kSemiMuonic:
 110         ForceParticleDecay(  411,13,1); // D+/-
 111         ForceParticleDecay(  421,13,1); // D0
 112         ForceParticleDecay(  431,13,1); // D_s
 113         ForceParticleDecay( 4122,13,1); // Lambda_c
 114         ForceParticleDecay( 4132,13,1); // Xsi_c
 115         ForceParticleDecay( 4232,13,1); // Sigma_c
 116         ForceParticleDecay( 4332,13,1); // Omega_c
 117         ForceParticleDecay(  511,13,1); // B0
 118         ForceParticleDecay(  521,13,1); // B+/-
 119         ForceParticleDecay(  531,13,1); // B_s
 120         ForceParticleDecay( 5122,13,1); // Lambda_b
 121         ForceParticleDecay( 5132,13,1); // Xsi_b
 122         ForceParticleDecay( 5232,13,1); // Sigma_b
 123         ForceParticleDecay( 5332,13,1); // Omega_b
 124         break;
 125     case kDiMuon:
 126         ForceParticleDecay(  113,13,2); // rho
 127         ForceParticleDecay(  221,13,2); // eta
 128         ForceParticleDecay(  223,13,2); // omega
 129         ForceParticleDecay(  333,13,2); // phi
 130         ForceParticleDecay(  443,13,2); // J/Psi
 131         ForceParticleDecay(100443,13,2);// Psi'
 132         ForceParticleDecay(  553,13,2); // Upsilon
 133         ForceParticleDecay(100553,13,2);// Upsilon'
 134         ForceParticleDecay(200553,13,2);// Upsilon''
 135         break;
 136     case kBSemiElectronic:
 137         ForceParticleDecay(  511,11,1); // B0
 138         ForceParticleDecay(  521,11,1); // B+/-
 139         ForceParticleDecay(  531,11,1); // B_s
 140         ForceParticleDecay( 5122,11,1); // Lambda_b
 141         ForceParticleDecay( 5132,11,1); // Xsi_b
 142         ForceParticleDecay( 5232,11,1); // Sigma_b
 143         ForceParticleDecay( 5332,11,1); // Omega_b
 144         break;
 145     case kSemiElectronic:
 146         ForceParticleDecay(  411,11,1); // D+/-
 147         ForceParticleDecay(  421,11,1); // D0
 148         ForceParticleDecay(  431,11,1); // D_s
 149         ForceParticleDecay( 4122,11,1); // Lambda_c
 150         ForceParticleDecay( 4132,11,1); // Xsi_c
 151         ForceParticleDecay( 4232,11,1); // Sigma_c
 152         ForceParticleDecay( 4332,11,1); // Omega_c
 153         ForceParticleDecay(  511,11,1); // B0
 154         ForceParticleDecay(  521,11,1); // B+/-
 155         ForceParticleDecay(  531,11,1); // B_s
 156         ForceParticleDecay( 5122,11,1); // Lambda_b
 157         ForceParticleDecay( 5132,11,1); // Xsi_b
 158         ForceParticleDecay( 5232,11,1); // Sigma_b
 159         ForceParticleDecay( 5332,11,1); // Omega_b
 160         break;
 161     case kDiElectron:
 162         ForceParticleDecay(  113,11,2); // rho
 163         ForceParticleDecay(  333,11,2); // phi
 164         ForceParticleDecay(  221,11,2); // eta
 165         ForceParticleDecay(  223,11,2); // omega
 166         ForceParticleDecay(  443,11,2); // J/Psi
 167         ForceParticleDecay(100443,11,2);// Psi'
 168         ForceParticleDecay(  553,11,2); // Upsilon
 169         ForceParticleDecay(100553,11,2);// Upsilon'
 170         ForceParticleDecay(200553,11,2);// Upsilon''
 171         break;
 172     case kBJpsiDiMuon:
 173
 174         products[0] =    443;
 175         products[1] = 100443;
 176         mult[0] = 1;
 177         mult[1] = 1;
 178
 179         ForceParticleDecay(  511, products, mult, 2); // B0   -> J/Psi (Psi') X
 180         ForceParticleDecay(  521, products, mult, 2); // B+/- -> J/Psi (Psi') X
 181         ForceParticleDecay(  531, products, mult, 2); // B_s  -> J/Psi (Psi') X
 182         ForceParticleDecay( 5122, products, mult, 2); // Lambda_b -> J/Psi (Psi') X
 183         ForceParticleDecay( 100443, 443, 1);          // Psi'  -> J/Psi X
 184         ForceParticleDecay(    443,13,2);             // J/Psi -> mu+ mu-
 185         break;
 186     case kBPsiPrimeDiMuon:
 187         ForceParticleDecay(  511,100443,1); // B0
 188         ForceParticleDecay(  521,100443,1); // B+/-
 189         ForceParticleDecay(  531,100443,1); // B_s
 190         ForceParticleDecay( 5122,100443,1); // Lambda_b
 191         ForceParticleDecay(100443,13,2);    // Psi'
 192         break;
 193     case kBJpsiDiElectron:
 194         ForceParticleDecay(  511,443,1); // B0
 195         ForceParticleDecay(  521,443,1); // B+/-
 196         ForceParticleDecay(  531,443,1); // B_s
 197         ForceParticleDecay( 5122,443,1); // Lambda_b
 198         ForceParticleDecay(  443,11,2);  // J/Psi
 199         break;
 200     case kBJpsi:
 201         ForceParticleDecay(  511,443,1); // B0
 202         ForceParticleDecay(  521,443,1); // B+/-
 203         ForceParticleDecay(  531,443,1); // B_s
 204         ForceParticleDecay( 5122,443,1); // Lambda_b
 205         break;
 206     case kBPsiPrimeDiElectron:
 207         ForceParticleDecay(  511,100443,1); // B0
 208         ForceParticleDecay(  521,100443,1); // B+/-
 209         ForceParticleDecay(  531,100443,1); // B_s
 210         ForceParticleDecay( 5122,100443,1); // Lambda_b
 211         ForceParticleDecay(100443,11,2);   // Psi'
 212         break;
 213     case kPiToMu:
 214         ForceParticleDecay(211,13,1); // pi->mu
 215         break;
 216     case kKaToMu:
 217         ForceParticleDecay(321,13,1); // K->mu
 218         break;
 219     case kAllMuonic:
 220         ForceParticleDecay(211,13,1); // pi->mu
 221         ForceParticleDecay(321,13,1); // K->mu
 222         break;
 223     case kWToMuon:
 224         ForceParticleDecay(  24, 13,1); // W -> mu
 225         break;
 226     case kWToCharm:
 227         ForceParticleDecay(   24, 4,1); // W -> c
 228         break;
 229     case kWToCharmToMuon:
 230         ForceParticleDecay(   24, 4,1); // W -> c
 231         ForceParticleDecay(  411,13,1); // D+/- -> mu
 232         ForceParticleDecay(  421,13,1); // D0  -> mu
 233         ForceParticleDecay(  431,13,1); // D_s  -> mu
 234         ForceParticleDecay( 4122,13,1); // Lambda_c
 235         ForceParticleDecay( 4132,13,1); // Xsi_c
 236         ForceParticleDecay( 4232,13,1); // Sigma_c
 237         ForceParticleDecay( 4332,13,1); // Omega_c
 238         break;
 239     case kZDiMuon:
 240         ForceParticleDecay(  23, 13,2); // Z -> mu+ mu-
 241         break;
 242     case kZDiElectron:
 243         ForceParticleDecay(  23, 11,2); // Z -> e+ e-
 244         break;
 245     case kHadronicD:
 246         ForceHadronicD(1);
 247         break;
 248     case kHadronicDWithout4Bodies:
 249         ForceHadronicD(0);
 250         break;
 251     case kPhiKK:
 252         ForceParticleDecay(333,321,2); // Phi->K+K-
 253         break;
 254     case kOmega:
 255 //      ForceOmega();
 256     case kAll:
 257         break;
 258     case kNoDecay:
 259         TPythia8::Instance()->ReadString("HadronLevel:Decay = off");
 260         break;
 261     case kNoDecayHeavy:
 262     case kNeutralPion:
 263         break;
 264     }
 265 }
 266
 267 Float_t AliDecayerPythia8::GetPartialBranchingRatio(Int_t ipart)
 268 {
 269     // Get the partial branching ration for the forced decay channels
 270
 271     Pythia8::Pythia* thePythia       = TPythia8::Instance()->Pythia8();
 272     Pythia8::ParticleDataTable table = thePythia->particleData;
 273     Pythia8::ParticleDataEntry* pd   = table.particleDataPtr(ipart);
 274     Pythia8::DecayTable  decays      = pd->decay;
 275
 276
 277     Int_t nc = decays.size();
 278     Float_t br = 0.;
 279 //
 280 //  Loop over decay channels
 281     for (Int_t ic = 0; ic < nc; ic++) {
 282         Pythia8::DecayChannel& decCh = decays[ic];
 283         for (Int_t i = 1; i <= decCh.multiplicity(); i++) {
 284             br += decCh.bRatio();
 285         }
 286     }
 287     return (br);
 288 }
 289
 290
 291 Int_t AliDecayerPythia8::CountProducts(Pythia8::DecayChannel& channel , Int_t particle)
 292 {
 293 // Count decay products of a given type
 294     Int_t np = 0;
 295     for (Int_t i = 1; i <= channel.multiplicity(); i++) {
 296         if (TMath::Abs(channel.product(i)) == particle) np++;
 297     }
 298     return np;
 299 }
 300
 301
 302
 303 void AliDecayerPythia8::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
 304 {
 305 //
 306 //  Force decay of particle into products with multiplicity mult
 307
 308     Pythia8::Pythia* thePythia       = TPythia8::Instance()->Pythia8();
 309     Pythia8::ParticleDataTable table = thePythia->particleData;
 310     Pythia8::ParticleDataEntry* pd   = table.particleDataPtr(particle);
 311     pd->setMayDecay(true);
 312     Pythia8::DecayTable  decays      = pd->decay;
 313
 314
 315     Int_t nc = decays.size();
 316 //
 317 //  Loop over decay channels
 318     for (Int_t ic = 0; ic < nc; ic++) {
 319         Pythia8::DecayChannel& decCh = decays[ic];
 320         if (CountProducts(decCh, product) >= mult) {
 321             decCh.onMode(1);
 322         } else {
 323             decCh.onMode(0);
 324         }
 325     }
 326 }
 327
 328 void AliDecayerPythia8::ForceParticleDecay(Int_t particle, Int_t* products, Int_t* mult, Int_t npart)
 329 {
 330 //
 331 //  Force decay of particle into products with multiplicity mult
 332
 333     Pythia8::Pythia* thePythia       = TPythia8::Instance()->Pythia8();
 334     Pythia8::ParticleDataTable table = thePythia->particleData;
 335     Pythia8::ParticleDataEntry* pd   = table.particleDataPtr(particle);
 336     pd->setMayDecay(true);
 337     Pythia8::DecayTable decays       = pd->decay;
 338
 339     Int_t nc = decays.size();
 340 //
 341 //  Loop over decay channels
 342     for (Int_t ic = 0; ic < nc; ic++) {
 343         Int_t nprod = 0;
 344         Pythia8::DecayChannel& decCh = decays[ic];
 345
 346         for (Int_t i = 0; i < npart; i++) {
 347             nprod += (CountProducts(decCh, products[i]) >= mult[i]);
 348         }
 349
 350         if (nprod) {
 351             decCh.onMode(1);
 352         } else {
 353             decCh.onMode(0);
 354         }
 355     }
 356 }
 357
 358
 359 Float_t AliDecayerPythia8::GetLifetime(Int_t kf)
 360 {
 361     // Return lifetime of particle
 362     Pythia8::Pythia* thePythia       = TPythia8::Instance()->Pythia8();
 363     Pythia8::ParticleDataTable table = thePythia->particleData;
 364     Float_t tau = table.tau0(kf);
 365     return ( tau);
 366 }
 367
 368 void  AliDecayerPythia8::SwitchOffHeavyFlavour()
 369 {
 370     // Switch off heavy flavour production
 371     //
 372 // Maximum number of quark flavours used in pdf
 373     TPythia8::Instance()->ReadString("PDFinProcess:nQuarkIn = 3");
 374 // Maximum number of flavors that can be used in showers
 375     TPythia8::Instance()->ReadString("SpaceShower:nQuarkIn = 3");
 376     TPythia8::Instance()->ReadString("TimeShower:nGammaToQuark = 3");
 377     TPythia8::Instance()->ReadString("TimeShower:nGluonToQuark = 3");
 378 }
 379
 380
 381 void AliDecayerPythia8::ForceHadronicD(Int_t optUse4Bodies)
 382 {
 383 //
 384 // Force golden D decay modes
 385 //
 386     const Int_t kNHadrons = 5;
 387     Int_t hadron[kNHadrons] = {411,  421, 431, 4112, 4122};
 388
 389     // for D+ -> K0* (-> K- pi+) pi+
 390     Int_t iKstar0    =  313;
 391     Int_t iKstarbar0 = -313;
 392     Int_t products[2] = {kKPlus, kPiMinus}, mult[2] = {1, 1};
 393     ForceParticleDecay(iKstar0, products, mult, 2);
 394     // for Ds -> Phi pi+
 395     Int_t iPhi = 333;
 396     ForceParticleDecay(iPhi, kKPlus, 2); // Phi->K+K-
 397     // for D0 -> rho0 pi+ k-
 398     Int_t iRho0=113;
 399     ForceParticleDecay(iRho0, kPiPlus, 2); // Rho0->pi+pi-
 400     // for Lambda_c -> Delta++ K-
 401     Int_t iDeltaPP = 2224;
 402     Int_t productsD[2] = {kProton, kPiPlus}, multD[2] = {1, 1};
 403     ForceParticleDecay(iDeltaPP, productsD, multD, 2);
 404
 405
 406     Int_t decayP1[kNHadrons][4] =
 407         {
 408             {kKMinus, kPiPlus,    kPiPlus, 0},
 409             {kKMinus, kPiPlus,    0      , 0},
 410             {kKPlus , iKstarbar0, 0      , 0},
 411             {-1     , -1        , -1     , -1},
 412             {kProton, iKstarbar0, 0      , 0}
 413         };
 414     Int_t decayP2[kNHadrons][4] =
 415         {
 416             {iKstarbar0, kPiPlus, 0      , 0},
 417             {kKMinus   , kPiPlus, kPiPlus, kPiMinus},
 418             {iPhi      , kPiPlus, 0      , 0},
 419             {-1        , -1     , -1     , -1},
 420             {iDeltaPP  , kKMinus, 0      , 0}
 421         };
 422     Int_t decayP3[kNHadrons][4] =
 423         {
 424             {-1        , -1     , -1      , -1},
 425             {kKMinus   , kPiPlus, iRho0   , 0 },
 426             {-1        , -1     , -1      , -1},
 427             {-1        , -1     , -1      , -1},
 428             {kProton   , kKMinus, kPiPlus , 0}
 429         };
 430     if(optUse4Bodies==0){
 431         for(Int_t iDau=0;iDau<4;iDau++){
 432             decayP2[1][iDau]=-1;
 433             decayP3[1][iDau]=-1;
 434         }
 435     }
 436
 437
 438     Pythia8::Pythia* thePythia       = TPythia8::Instance()->Pythia8();
 439     Pythia8::ParticleDataTable table = thePythia->particleData;
 440
 441     for (Int_t ihadron = 0; ihadron < kNHadrons; ihadron++)
 442     {
 443         Pythia8::ParticleDataEntry* pd   = table.particleDataPtr(hadron[ihadron]);
 444         pd->setMayDecay(true);
 445         Pythia8::DecayTable  decays      = pd->decay;
 446
 447
 448         for (Int_t ic = 0; ic < decays.size(); ic++) {
 449             Pythia8::DecayChannel& decCh = decays[ic];
 450             if ((
 451                     decCh.product(0) == decayP1[ihadron][0] &&
 452                     decCh.product(1) == decayP1[ihadron][1] &&
 453                     decCh.product(2) == decayP1[ihadron][2] &&
 454                     decCh.product(3) == decayP1[ihadron][3] &&
 455                     decCh.product(4) == 0
 456                     )
 457                 || (
 458                     decCh.product(0) == decayP2[ihadron][0] &&
 459                     decCh.product(1) == decayP2[ihadron][1] &&
 460                     decCh.product(2) == decayP2[ihadron][2] &&
 461                     decCh.product(3) == decayP2[ihadron][3] &&
 462                     decCh.product(4) == 0
 463                     )
 464                 || (
 465                     decCh.product(0) == decayP3[ihadron][0] &&
 466                     decCh.product(1) == decayP3[ihadron][1] &&
 467                     decCh.product(2) == decayP3[ihadron][2] &&
 468                     decCh.product(3) == decayP3[ihadron][3] &&
 469                     decCh.product(4) == 0
 470                     ))
 471             {
 472                 decCh.onMode(1);
 473             } else {
 474                 decCh.onMode(0);
 475             } // selected channel ?
 476         } // decay channels
 477     } // hadrons
 478 }
 479