From 948d10f49b699507dde9d22dbd3fd866ad20426a Mon Sep 17 00:00:00 2001 From: morsch Date: Thu, 11 Dec 2008 13:29:42 +0000 Subject: [PATCH] Decayer with alice decay options using Pythia8 --- PYTHIA8/AliDecayerPythia8.cxx | 479 ++++++++++++++++++++++++++++++++++ PYTHIA8/AliDecayerPythia8.h | 43 +++ PYTHIA8/AliPythia8LinkDef.h | 1 + PYTHIA8/libAliPythia8.pkg | 2 +- 4 files changed, 524 insertions(+), 1 deletion(-) create mode 100644 PYTHIA8/AliDecayerPythia8.cxx create mode 100644 PYTHIA8/AliDecayerPythia8.h diff --git a/PYTHIA8/AliDecayerPythia8.cxx b/PYTHIA8/AliDecayerPythia8.cxx new file mode 100644 index 00000000000..c3eff9764b1 --- /dev/null +++ b/PYTHIA8/AliDecayerPythia8.cxx @@ -0,0 +1,479 @@ +/************************************************************************** + * 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 +#include +#include +#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 +} + diff --git a/PYTHIA8/AliDecayerPythia8.h b/PYTHIA8/AliDecayerPythia8.h new file mode 100644 index 00000000000..54c97082831 --- /dev/null +++ b/PYTHIA8/AliDecayerPythia8.h @@ -0,0 +1,43 @@ +#ifndef ALIDECAYERPYTHIA8_H +#define ALIDECAYERPYTHIA8_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/* $Id$*/ + +// Implementation of TVirtualMCDecayer using Pythia8 +// Author: andreas.morsch@cern.ch + +#include +#include +#include "AliDecayer.h" + +class AliDecayerPythia8 : public TPythia8Decayer { + public: + AliDecayerPythia8(); + virtual ~AliDecayerPythia8(){;} + virtual void SetForceDecay(Decay_t decay) {fDecay=decay;} + virtual void SetForceDecay(Int_t decay) {SetForceDecay((Decay_t) decay);} + virtual void ForceDecay(); + virtual Float_t GetPartialBranchingRatio(Int_t ipart); + virtual void HeavyFlavourOff() {fHeavyFlavour = kFALSE;} + virtual Float_t GetLifetime(Int_t kf); + private: + Int_t CountProducts(Pythia8::DecayChannel& decCh, Int_t particle); + void ForceParticleDecay(Int_t particle, Int_t product, Int_t mult); + void ForceParticleDecay(Int_t particle, Int_t* products, Int_t* mult, Int_t npart); + void SwitchOffHeavyFlavour(); + void ForceHadronicD(Int_t optUser4Bodies = 1); + protected: + Decay_t fDecay; // Forced decay mode + Bool_t fHeavyFlavour; //! Flag for heavy flavors + ClassDef(AliDecayerPythia8, 1) // Particle Decayer using Pythia8 +}; +#endif + + + + + + + diff --git a/PYTHIA8/AliPythia8LinkDef.h b/PYTHIA8/AliPythia8LinkDef.h index 46111e81fd0..089c96a7491 100644 --- a/PYTHIA8/AliPythia8LinkDef.h +++ b/PYTHIA8/AliPythia8LinkDef.h @@ -3,4 +3,5 @@ #pragma link off all classes; #pragma link off all functions; #pragma link C++ class AliPythia8+; +#pragma link C++ class AliDecayerPythia8+; #endif diff --git a/PYTHIA8/libAliPythia8.pkg b/PYTHIA8/libAliPythia8.pkg index 1072fbd8d65..ab6c42b54b7 100644 --- a/PYTHIA8/libAliPythia8.pkg +++ b/PYTHIA8/libAliPythia8.pkg @@ -1,4 +1,4 @@ -SRCS= AliPythia8.cxx +SRCS= AliPythia8.cxx AliDecayerPythia8.cxx HDRS:= $(SRCS:.cxx=.h) DHDR:= AliPythia8LinkDef.h EINCLUDE:= $(PYTHIA8)/include PYTHIA6 -- 2.43.0