--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id$ */
+
+// Implementation of AliDecayer using Pythia8
+// Author: andreas.morsch@cern.ch
+#include <TMath.h>
+#include <TPDGCode.h>
+#include <TPythia8.h>
+#include "AliDecayerPythia8.h"
+#include "ParticleData.h"
+
+ClassImp(AliDecayerPythia8)
+
+AliDecayerPythia8::AliDecayerPythia8():
+ TPythia8Decayer(),
+ fDecay(kAll),
+ fHeavyFlavour(kTRUE)
+{
+ // Constructor
+}
+
+void AliDecayerPythia8::ForceDecay()
+{
+ //
+// Force a particle decay mode
+// Switch heavy flavour production off if requested
+ if (!fHeavyFlavour) SwitchOffHeavyFlavour();
+//
+ Decay_t decay = fDecay;
+ TPythia8::Instance()->ReadString("HadronLevel:Decay = on");
+
+ if (decay == kNoDecayHeavy) return;
+
+//
+// select mode
+ Int_t products[2];
+ Int_t mult[2];
+ Int_t products1[3];
+ Int_t mult1[3];
+
+ switch (decay)
+ {
+ case kHardMuons:
+ products1[0] = 13;
+ products1[1] = 443;
+ products1[2] = 100443;
+ mult1[0] = 1;
+ mult1[1] = 1;
+ mult1[2] = 1;
+ ForceParticleDecay( 511, products1, mult1, 3);
+ ForceParticleDecay( 521, products1, mult1, 3);
+ ForceParticleDecay( 531, products1, mult1, 3);
+ ForceParticleDecay( 5122, products1, mult1, 3);
+ ForceParticleDecay( 5132, products1, mult1, 3);
+ ForceParticleDecay( 5232, products1, mult1, 3);
+ ForceParticleDecay( 5332, products1, mult1, 3);
+ ForceParticleDecay( 100443, 443, 1); // Psi' -> J/Psi X
+ ForceParticleDecay( 443, 13, 2); // J/Psi -> mu+ mu-
+ ForceParticleDecay( 411,13,1); // D+/-
+ ForceParticleDecay( 421,13,1); // D0
+ ForceParticleDecay( 431,13,1); // D_s
+ ForceParticleDecay( 4122,13,1); // Lambda_c
+ ForceParticleDecay( 4132,13,1); // Xsi_c
+ ForceParticleDecay( 4232,13,1); // Sigma_c
+ ForceParticleDecay( 4332,13,1); // Omega_c
+ break;
+ case kChiToJpsiGammaToMuonMuon:
+ products[0] = 443;
+ products[1] = 22;
+ mult[0] = 1;
+ mult[1] = 1;
+ ForceParticleDecay( 20443, products, mult, 2); // Chi_1c -> J/Psi Gamma
+ ForceParticleDecay( 445, products, mult, 2); // Chi_2c -> J/Psi Gamma
+ ForceParticleDecay( 443, 13, 2); // J/Psi -> mu+ mu-
+ break;
+ case kChiToJpsiGammaToElectronElectron:
+ products[0] = 443;
+ products[1] = 22;
+ mult[0] = 1;
+ mult[1] = 1;
+ ForceParticleDecay( 20443, products, mult, 2); // Chi_1c -> J/Psi Gamma
+ ForceParticleDecay( 445, products, mult, 2); // Chi_2c -> J/Psi Gamma
+ ForceParticleDecay( 443, 11, 2); // J/Psi -> e+ e-
+ break;
+
+ case kBSemiMuonic:
+ ForceParticleDecay( 511,13,1); // B0
+ ForceParticleDecay( 521,13,1); // B+/-
+ ForceParticleDecay( 531,13,1); // B_s
+ ForceParticleDecay( 5122,13,1); // Lambda_b
+ ForceParticleDecay( 5132,13,1); // Xsi_b
+ ForceParticleDecay( 5232,13,1); // Sigma_b
+ ForceParticleDecay( 5332,13,1); // Omega_b
+ break;
+ case kSemiMuonic:
+ ForceParticleDecay( 411,13,1); // D+/-
+ ForceParticleDecay( 421,13,1); // D0
+ ForceParticleDecay( 431,13,1); // D_s
+ ForceParticleDecay( 4122,13,1); // Lambda_c
+ ForceParticleDecay( 4132,13,1); // Xsi_c
+ ForceParticleDecay( 4232,13,1); // Sigma_c
+ ForceParticleDecay( 4332,13,1); // Omega_c
+ ForceParticleDecay( 511,13,1); // B0
+ ForceParticleDecay( 521,13,1); // B+/-
+ ForceParticleDecay( 531,13,1); // B_s
+ ForceParticleDecay( 5122,13,1); // Lambda_b
+ ForceParticleDecay( 5132,13,1); // Xsi_b
+ ForceParticleDecay( 5232,13,1); // Sigma_b
+ ForceParticleDecay( 5332,13,1); // Omega_b
+ break;
+ case kDiMuon:
+ ForceParticleDecay( 113,13,2); // rho
+ ForceParticleDecay( 221,13,2); // eta
+ ForceParticleDecay( 223,13,2); // omega
+ ForceParticleDecay( 333,13,2); // phi
+ ForceParticleDecay( 443,13,2); // J/Psi
+ ForceParticleDecay(100443,13,2);// Psi'
+ ForceParticleDecay( 553,13,2); // Upsilon
+ ForceParticleDecay(100553,13,2);// Upsilon'
+ ForceParticleDecay(200553,13,2);// Upsilon''
+ break;
+ case kBSemiElectronic:
+ ForceParticleDecay( 511,11,1); // B0
+ ForceParticleDecay( 521,11,1); // B+/-
+ ForceParticleDecay( 531,11,1); // B_s
+ ForceParticleDecay( 5122,11,1); // Lambda_b
+ ForceParticleDecay( 5132,11,1); // Xsi_b
+ ForceParticleDecay( 5232,11,1); // Sigma_b
+ ForceParticleDecay( 5332,11,1); // Omega_b
+ break;
+ case kSemiElectronic:
+ ForceParticleDecay( 411,11,1); // D+/-
+ ForceParticleDecay( 421,11,1); // D0
+ ForceParticleDecay( 431,11,1); // D_s
+ ForceParticleDecay( 4122,11,1); // Lambda_c
+ ForceParticleDecay( 4132,11,1); // Xsi_c
+ ForceParticleDecay( 4232,11,1); // Sigma_c
+ ForceParticleDecay( 4332,11,1); // Omega_c
+ ForceParticleDecay( 511,11,1); // B0
+ ForceParticleDecay( 521,11,1); // B+/-
+ ForceParticleDecay( 531,11,1); // B_s
+ ForceParticleDecay( 5122,11,1); // Lambda_b
+ ForceParticleDecay( 5132,11,1); // Xsi_b
+ ForceParticleDecay( 5232,11,1); // Sigma_b
+ ForceParticleDecay( 5332,11,1); // Omega_b
+ break;
+ case kDiElectron:
+ ForceParticleDecay( 113,11,2); // rho
+ ForceParticleDecay( 333,11,2); // phi
+ ForceParticleDecay( 221,11,2); // eta
+ ForceParticleDecay( 223,11,2); // omega
+ ForceParticleDecay( 443,11,2); // J/Psi
+ ForceParticleDecay(100443,11,2);// Psi'
+ ForceParticleDecay( 553,11,2); // Upsilon
+ ForceParticleDecay(100553,11,2);// Upsilon'
+ ForceParticleDecay(200553,11,2);// Upsilon''
+ break;
+ case kBJpsiDiMuon:
+
+ products[0] = 443;
+ products[1] = 100443;
+ mult[0] = 1;
+ mult[1] = 1;
+
+ ForceParticleDecay( 511, products, mult, 2); // B0 -> J/Psi (Psi') X
+ ForceParticleDecay( 521, products, mult, 2); // B+/- -> J/Psi (Psi') X
+ ForceParticleDecay( 531, products, mult, 2); // B_s -> J/Psi (Psi') X
+ ForceParticleDecay( 5122, products, mult, 2); // Lambda_b -> J/Psi (Psi') X
+ ForceParticleDecay( 100443, 443, 1); // Psi' -> J/Psi X
+ ForceParticleDecay( 443,13,2); // J/Psi -> mu+ mu-
+ break;
+ case kBPsiPrimeDiMuon:
+ ForceParticleDecay( 511,100443,1); // B0
+ ForceParticleDecay( 521,100443,1); // B+/-
+ ForceParticleDecay( 531,100443,1); // B_s
+ ForceParticleDecay( 5122,100443,1); // Lambda_b
+ ForceParticleDecay(100443,13,2); // Psi'
+ break;
+ case kBJpsiDiElectron:
+ ForceParticleDecay( 511,443,1); // B0
+ ForceParticleDecay( 521,443,1); // B+/-
+ ForceParticleDecay( 531,443,1); // B_s
+ ForceParticleDecay( 5122,443,1); // Lambda_b
+ ForceParticleDecay( 443,11,2); // J/Psi
+ break;
+ case kBJpsi:
+ ForceParticleDecay( 511,443,1); // B0
+ ForceParticleDecay( 521,443,1); // B+/-
+ ForceParticleDecay( 531,443,1); // B_s
+ ForceParticleDecay( 5122,443,1); // Lambda_b
+ break;
+ case kBPsiPrimeDiElectron:
+ ForceParticleDecay( 511,100443,1); // B0
+ ForceParticleDecay( 521,100443,1); // B+/-
+ ForceParticleDecay( 531,100443,1); // B_s
+ ForceParticleDecay( 5122,100443,1); // Lambda_b
+ ForceParticleDecay(100443,11,2); // Psi'
+ break;
+ case kPiToMu:
+ ForceParticleDecay(211,13,1); // pi->mu
+ break;
+ case kKaToMu:
+ ForceParticleDecay(321,13,1); // K->mu
+ break;
+ case kAllMuonic:
+ ForceParticleDecay(211,13,1); // pi->mu
+ ForceParticleDecay(321,13,1); // K->mu
+ break;
+ case kWToMuon:
+ ForceParticleDecay( 24, 13,1); // W -> mu
+ break;
+ case kWToCharm:
+ ForceParticleDecay( 24, 4,1); // W -> c
+ break;
+ case kWToCharmToMuon:
+ ForceParticleDecay( 24, 4,1); // W -> c
+ ForceParticleDecay( 411,13,1); // D+/- -> mu
+ ForceParticleDecay( 421,13,1); // D0 -> mu
+ ForceParticleDecay( 431,13,1); // D_s -> mu
+ ForceParticleDecay( 4122,13,1); // Lambda_c
+ ForceParticleDecay( 4132,13,1); // Xsi_c
+ ForceParticleDecay( 4232,13,1); // Sigma_c
+ ForceParticleDecay( 4332,13,1); // Omega_c
+ break;
+ case kZDiMuon:
+ ForceParticleDecay( 23, 13,2); // Z -> mu+ mu-
+ break;
+ case kZDiElectron:
+ ForceParticleDecay( 23, 11,2); // Z -> e+ e-
+ break;
+ case kHadronicD:
+ ForceHadronicD(1);
+ break;
+ case kHadronicDWithout4Bodies:
+ ForceHadronicD(0);
+ break;
+ case kPhiKK:
+ ForceParticleDecay(333,321,2); // Phi->K+K-
+ break;
+ case kOmega:
+// ForceOmega();
+ case kAll:
+ break;
+ case kNoDecay:
+ TPythia8::Instance()->ReadString("HadronLevel:Decay = off");
+ break;
+ case kNoDecayHeavy:
+ case kNeutralPion:
+ break;
+ }
+}
+
+Float_t AliDecayerPythia8::GetPartialBranchingRatio(Int_t ipart)
+{
+ // Get the partial branching ration for the forced decay channels
+
+ Pythia8::Pythia* thePythia = TPythia8::Instance()->Pythia8();
+ Pythia8::ParticleDataTable table = thePythia->particleData;
+ Pythia8::ParticleDataEntry* pd = table.particleDataPtr(ipart);
+ Pythia8::DecayTable decays = pd->decay;
+
+
+ Int_t nc = decays.size();
+ Float_t br = 0.;
+//
+// Loop over decay channels
+ for (Int_t ic = 0; ic < nc; ic++) {
+ Pythia8::DecayChannel& decCh = decays[ic];
+ for (Int_t i = 1; i <= decCh.multiplicity(); i++) {
+ br += decCh.bRatio();
+ }
+ }
+ return (br);
+}
+
+
+Int_t AliDecayerPythia8::CountProducts(Pythia8::DecayChannel& channel , Int_t particle)
+{
+// Count decay products of a given type
+ Int_t np = 0;
+ for (Int_t i = 1; i <= channel.multiplicity(); i++) {
+ if (TMath::Abs(channel.product(i)) == particle) np++;
+ }
+ return np;
+}
+
+
+
+void AliDecayerPythia8::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
+{
+//
+// Force decay of particle into products with multiplicity mult
+
+ Pythia8::Pythia* thePythia = TPythia8::Instance()->Pythia8();
+ Pythia8::ParticleDataTable table = thePythia->particleData;
+ Pythia8::ParticleDataEntry* pd = table.particleDataPtr(particle);
+ pd->setMayDecay(true);
+ Pythia8::DecayTable decays = pd->decay;
+
+
+ Int_t nc = decays.size();
+//
+// Loop over decay channels
+ for (Int_t ic = 0; ic < nc; ic++) {
+ Pythia8::DecayChannel& decCh = decays[ic];
+ if (CountProducts(decCh, product) >= mult) {
+ decCh.onMode(1);
+ } else {
+ decCh.onMode(0);
+ }
+ }
+}
+
+void AliDecayerPythia8::ForceParticleDecay(Int_t particle, Int_t* products, Int_t* mult, Int_t npart)
+{
+//
+// Force decay of particle into products with multiplicity mult
+
+ Pythia8::Pythia* thePythia = TPythia8::Instance()->Pythia8();
+ Pythia8::ParticleDataTable table = thePythia->particleData;
+ Pythia8::ParticleDataEntry* pd = table.particleDataPtr(particle);
+ pd->setMayDecay(true);
+ Pythia8::DecayTable decays = pd->decay;
+
+ Int_t nc = decays.size();
+//
+// Loop over decay channels
+ for (Int_t ic = 0; ic < nc; ic++) {
+ Int_t nprod = 0;
+ Pythia8::DecayChannel& decCh = decays[ic];
+
+ for (Int_t i = 0; i < npart; i++) {
+ nprod += (CountProducts(decCh, products[i]) >= mult[i]);
+ }
+
+ if (nprod) {
+ decCh.onMode(1);
+ } else {
+ decCh.onMode(0);
+ }
+ }
+}
+
+
+Float_t AliDecayerPythia8::GetLifetime(Int_t kf)
+{
+ // Return lifetime of particle
+ Pythia8::Pythia* thePythia = TPythia8::Instance()->Pythia8();
+ Pythia8::ParticleDataTable table = thePythia->particleData;
+ Float_t tau = table.tau0(kf);
+ return ( tau);
+}
+
+void AliDecayerPythia8::SwitchOffHeavyFlavour()
+{
+ // Switch off heavy flavour production
+ //
+// Maximum number of quark flavours used in pdf
+ TPythia8::Instance()->ReadString("PDFinProcess:nQuarkIn = 3");
+// Maximum number of flavors that can be used in showers
+ TPythia8::Instance()->ReadString("SpaceShower:nQuarkIn = 3");
+ TPythia8::Instance()->ReadString("TimeShower:nGammaToQuark = 3");
+ TPythia8::Instance()->ReadString("TimeShower:nGluonToQuark = 3");
+}
+
+
+void AliDecayerPythia8::ForceHadronicD(Int_t optUse4Bodies)
+{
+//
+// Force golden D decay modes
+//
+ const Int_t kNHadrons = 5;
+ Int_t hadron[kNHadrons] = {411, 421, 431, 4112, 4122};
+
+ // for D+ -> K0* (-> K- pi+) pi+
+ Int_t iKstar0 = 313;
+ Int_t iKstarbar0 = -313;
+ Int_t products[2] = {kKPlus, kPiMinus}, mult[2] = {1, 1};
+ ForceParticleDecay(iKstar0, products, mult, 2);
+ // for Ds -> Phi pi+
+ Int_t iPhi = 333;
+ ForceParticleDecay(iPhi, kKPlus, 2); // Phi->K+K-
+ // for D0 -> rho0 pi+ k-
+ Int_t iRho0=113;
+ ForceParticleDecay(iRho0, kPiPlus, 2); // Rho0->pi+pi-
+ // for Lambda_c -> Delta++ K-
+ Int_t iDeltaPP = 2224;
+ Int_t productsD[2] = {kProton, kPiPlus}, multD[2] = {1, 1};
+ ForceParticleDecay(iDeltaPP, productsD, multD, 2);
+
+
+ Int_t decayP1[kNHadrons][4] =
+ {
+ {kKMinus, kPiPlus, kPiPlus, 0},
+ {kKMinus, kPiPlus, 0 , 0},
+ {kKPlus , iKstarbar0, 0 , 0},
+ {-1 , -1 , -1 , -1},
+ {kProton, iKstarbar0, 0 , 0}
+ };
+ Int_t decayP2[kNHadrons][4] =
+ {
+ {iKstarbar0, kPiPlus, 0 , 0},
+ {kKMinus , kPiPlus, kPiPlus, kPiMinus},
+ {iPhi , kPiPlus, 0 , 0},
+ {-1 , -1 , -1 , -1},
+ {iDeltaPP , kKMinus, 0 , 0}
+ };
+ Int_t decayP3[kNHadrons][4] =
+ {
+ {-1 , -1 , -1 , -1},
+ {kKMinus , kPiPlus, iRho0 , 0 },
+ {-1 , -1 , -1 , -1},
+ {-1 , -1 , -1 , -1},
+ {kProton , kKMinus, kPiPlus , 0}
+ };
+ if(optUse4Bodies==0){
+ for(Int_t iDau=0;iDau<4;iDau++){
+ decayP2[1][iDau]=-1;
+ decayP3[1][iDau]=-1;
+ }
+ }
+
+
+ Pythia8::Pythia* thePythia = TPythia8::Instance()->Pythia8();
+ Pythia8::ParticleDataTable table = thePythia->particleData;
+
+ for (Int_t ihadron = 0; ihadron < kNHadrons; ihadron++)
+ {
+ Pythia8::ParticleDataEntry* pd = table.particleDataPtr(hadron[ihadron]);
+ pd->setMayDecay(true);
+ Pythia8::DecayTable decays = pd->decay;
+
+
+ for (Int_t ic = 0; ic < decays.size(); ic++) {
+ Pythia8::DecayChannel& decCh = decays[ic];
+ if ((
+ decCh.product(0) == decayP1[ihadron][0] &&
+ decCh.product(1) == decayP1[ihadron][1] &&
+ decCh.product(2) == decayP1[ihadron][2] &&
+ decCh.product(3) == decayP1[ihadron][3] &&
+ decCh.product(4) == 0
+ )
+ || (
+ decCh.product(0) == decayP2[ihadron][0] &&
+ decCh.product(1) == decayP2[ihadron][1] &&
+ decCh.product(2) == decayP2[ihadron][2] &&
+ decCh.product(3) == decayP2[ihadron][3] &&
+ decCh.product(4) == 0
+ )
+ || (
+ decCh.product(0) == decayP3[ihadron][0] &&
+ decCh.product(1) == decayP3[ihadron][1] &&
+ decCh.product(2) == decayP3[ihadron][2] &&
+ decCh.product(3) == decayP3[ihadron][3] &&
+ decCh.product(4) == 0
+ ))
+ {
+ decCh.onMode(1);
+ } else {
+ decCh.onMode(0);
+ } // selected channel ?
+ } // decay channels
+ } // hadrons
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff