// Author: andreas.morsch@cern.ch
#include <TMath.h>
#include <TPDGCode.h>
-#include <TPythia8.h>
+#include <TLorentzVector.h>
+#include <TClonesArray.h>
+#include <TParticle.h>
+#include "AliTPythia8.h"
#include "AliDecayerPythia8.h"
#include "ParticleData.h"
ClassImp(AliDecayerPythia8)
+Bool_t AliDecayerPythia8::fgInit = kFALSE;
+
AliDecayerPythia8::AliDecayerPythia8():
- TPythia8Decayer(),
- fDecay(kAll),
- fHeavyFlavour(kTRUE)
+ TVirtualMCDecayer(),
+ fPythia8(AliTPythia8::Instance()),
+ fDebug(0),
+ fDecay(kAll),
+ fHeavyFlavour(kTRUE)
{
// Constructor
+ fPythia8->Pythia8()->readString("SoftQCD:elastic = on");
+ fPythia8->Pythia8()->init();
+}
+
+//___________________________________________________________________________
+void AliDecayerPythia8::Decay(Int_t pdg, TLorentzVector* p)
+{
+ // Decay a single particle
+ ClearEvent();
+ AppendParticle(pdg, p);
+ Int_t idPart = fPythia8->Pythia8()->event[0].id();
+ fPythia8->Pythia8()->particleData.mayDecay(idPart,kTRUE);
+ fPythia8->Pythia8()->moreDecays();
+ if (fDebug > 0) fPythia8->EventListing();
+}
+
+//___________________________________________________________________________
+Int_t AliDecayerPythia8::ImportParticles(TClonesArray *particles)
+{
+ //import the decay products into particles array
+ const Float_t kconvT=0.001/2.999792458e8; // mm/c to seconds conversion
+ const Float_t kconvL=1./10; // mm to cm conversion
+ int np = fPythia8->ImportParticles(particles, "All");
+ // pythia assigns decay time in mm/c, convert to seconds
+ for (int ip=1;ip<np;ip++) {
+ TParticle* prod = (TParticle*)particles->At(ip);
+ if (!prod) continue;
+ prod->SetProductionVertex(prod->Vx()*kconvL,prod->Vy()*kconvL,prod->Vz()*kconvL,kconvT*prod->T());
+ }
+ return np;
+}
+
+
+void AliDecayerPythia8::Init()
+{
+// Initialisation
+//
+ if (!fgInit) {
+ fgInit = kTRUE;
+ // fPythia->SetDecayTable();
+ }
+
+// Switch on heavy flavor decays
+
+ Int_t j;
+ Int_t heavy[14] = {411, 421, 431, 4122, 4132, 4232, 4332, 511, 521, 531, 5122, 5132, 5232, 5332};
+// fPythia->ResetDecayTable();
+ for (j=0; j < 14; j++) {
+ if (fDecay == kNoDecayHeavy) {
+ fPythia8->ReadString(Form("%d:onMode = off", heavy[j]));
+ } else {
+ fPythia8->ReadString(Form("%d:onMode = on", heavy[j]));
+ }
+ }
+
+
+//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
+
+ if (fDecay != kNeutralPion) {
+ fPythia8->ReadString("111:onMode = off");
+ } else {
+ fPythia8->ReadString("111:onMode = on");
+ }
+
+ fPythia8->ReadString("310:onMode = off");
+ fPythia8->ReadString("3112:onMode = off");
+ fPythia8->ReadString("3212:onMode = off");
+ fPythia8->ReadString("3222:onMode = off");
+ fPythia8->ReadString("3312:onMode = off");
+ fPythia8->ReadString("3322:onMode = off");
+ fPythia8->ReadString("3334:onMode = off");
+// .. Force decay channels
+ ForceDecay();
}
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");
+ fPythia8->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
+// B0 -> mu X
+ fPythia8->ReadString("511:onMode = off");
+ fPythia8->ReadString("511:onIfAny = 13 443 100443");
+// B+/- -> mu X
+ fPythia8->ReadString("521:onMode = off");
+ fPythia8->ReadString("521:onIfAny = 13 443 100443");
+// Bs -> mu X
+ fPythia8->ReadString("531:onMode = off");
+ fPythia8->ReadString("531:onIfAny = 13 443 100443");
+// Lambda_b -> mu X
+ fPythia8->ReadString("5122:onMode = off");
+ fPythia8->ReadString("5122:onIfAny = 13 443 100443");
+// Sigma_b- -> mu X
+ fPythia8->ReadString("5132:onMode = off");
+ fPythia8->ReadString("5132:onIfAny = 13 443 100443");
+// Sigma_b0 -> mu X
+ fPythia8->ReadString("5232:onMode = off");
+ fPythia8->ReadString("5232:onIfAny = 13 443 100443");
+// Omega_b -> mu X
+ fPythia8->ReadString("5332:onMode = off");
+ fPythia8->ReadString("5332:onIfAny = 13 443 100443");
+// Psi' -> mu+ mu-
+ fPythia8->ReadString("100443:onMode = off");
+ fPythia8->ReadString("100443:onIfAny = 443");
+// Psi -> mu+ mu-
+ fPythia8->ReadString("443:onMode = off");
+ fPythia8->ReadString("443:onIfAll = 13 13");
+// D+/- -> mu X
+ fPythia8->ReadString("411:onMode = off");
+ fPythia8->ReadString("411:onIfAll = 13");
+// D0 -> mu X
+ fPythia8->ReadString("421:onMode = off");
+ fPythia8->ReadString("421:onIfAll = 13");
+// D_s -> mu X
+ fPythia8->ReadString("431:onMode = off");
+ fPythia8->ReadString("431:onIfAll = 13");
+// Lambda_c -> mu X
+ fPythia8->ReadString("4122:onMode = off");
+ fPythia8->ReadString("4122:onIfAll = 13");
+// Sigma_c -> mu X
+ fPythia8->ReadString("4132:onMode = off");
+ fPythia8->ReadString("4132:onIfAll = 13");
+// Sigma_c+ -> mu X
+ fPythia8->ReadString("4232:onMode = off");
+ fPythia8->ReadString("4232:onIfAll = 13");
+// Omega_c -> mu X
+ fPythia8->ReadString("4332:onMode = off");
+ fPythia8->ReadString("4332:onIfAll = 13");
+
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-
+// Chi_1c -> J/Psi Gamma
+ fPythia8->ReadString("20443:onMode = off");
+ fPythia8->ReadString("20443:onIfAll = 443 22");
+// Chi_2c -> J/Psi Gamma
+ fPythia8->ReadString("445:onMode = off");
+ fPythia8->ReadString("445:onIfAll = 443 22");
+// J/Psi -> mu+ mu-
+ fPythia8->ReadString("443:onMode = off");
+ fPythia8->ReadString("443:onIfAll = 13 13");
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-
+// Chi_1c -> J/Psi Gamma
+ fPythia8->ReadString("20443:onMode = off");
+ fPythia8->ReadString("20443:onIfAll = 443 22");
+// Chi_2c -> J/Psi Gamma
+ fPythia8->ReadString("445:onMode = off");
+ fPythia8->ReadString("445:onIfAll = 443 22");
+// J/Psi -> e+ e-
+ fPythia8->ReadString("443:onMode = off");
+ fPythia8->ReadString("443:onIfAll = 11 11");
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
+// B0 -> mu X
+ fPythia8->ReadString("511:onMode = off");
+ fPythia8->ReadString("511:onIfAny = 13");
+// B+/- -> mu X
+ fPythia8->ReadString("521:onMode = off");
+ fPythia8->ReadString("521:onIfAny = 13");
+// B_s -> mu X
+ fPythia8->ReadString("531:onMode = off");
+ fPythia8->ReadString("531:onIfAny = 13");
+// Lambda_b -> mu X
+ fPythia8->ReadString("5122:onMode = off");
+ fPythia8->ReadString("5122:onIfAny = 13");
+// Sigma_b -> mu X
+ fPythia8->ReadString("5132:onMode = off");
+ fPythia8->ReadString("5132:onIfAny = 13");
+// Sigma_b0 -> mu X
+ fPythia8->ReadString("5232:onMode = off");
+ fPythia8->ReadString("5232:onIfAny = 13");
+// Omega_b -> mu X
+ fPythia8->ReadString("5332:onMode = off");
+ fPythia8->ReadString("5332:onIfAny = 13");
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
+// D+- -> mu X
+ fPythia8->ReadString("411:onMode = off");
+ fPythia8->ReadString("411:onIfAll = 13");
+// D0 -> mu X
+ fPythia8->ReadString("421:onMode = off");
+ fPythia8->ReadString("421:onIfAll = 13");
+// D_s -> mu X
+ fPythia8->ReadString("431:onMode = off");
+ fPythia8->ReadString("431:onIfAll = 13");
+// Lambda_c -> mu X
+ fPythia8->ReadString("4122:onMode = off");
+ fPythia8->ReadString("4122:onIfAll = 13");
+// Sigma_c -> mu X
+ fPythia8->ReadString("4132:onMode = off");
+ fPythia8->ReadString("4132:onIfAll = 13");
+// Sigma -> mu X
+ fPythia8->ReadString("4232:onMode = off");
+ fPythia8->ReadString("4232:onIfAll = 13");
+// Omega_c -> mu X
+ fPythia8->ReadString("4332:onMode = off");
+ fPythia8->ReadString("4332:onIfAll = 13");
+// B0 -> mu X
+ fPythia8->ReadString("511:onMode = off");
+ fPythia8->ReadString("511:onIfAny = 13");
+// B+/- -> mu X
+ fPythia8->ReadString("521:onMode = off");
+ fPythia8->ReadString("521:onIfAny = 13");
+// B_s -> mu X
+ fPythia8->ReadString("531:onMode = off");
+ fPythia8->ReadString("531:onIfAny = 13");
+// Lambda_c -> mu X
+ fPythia8->ReadString("5122:onMode = off");
+ fPythia8->ReadString("5122:onIfAny = 13");
+// Sigma_c -> mu X
+ fPythia8->ReadString("5132:onMode = off");
+ fPythia8->ReadString("5132:onIfAny = 13");
+// Sigma_c -> mu X
+ fPythia8->ReadString("5232:onMode = off");
+ fPythia8->ReadString("5232:onIfAny = 13");
+// Omega_c -> mu X
+ fPythia8->ReadString("5332:onMode = off");
+ fPythia8->ReadString("5332:onIfAny = 13");
+
+ break;
+ case kJpsiDiMuon:
+// J/Psi-> mu+ mu-
+ fPythia8->ReadString("443:onMode = off");
+ fPythia8->ReadString("443:onIfAll = 13 13");
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''
+// Rho -> mu+ mu-
+ fPythia8->ReadString("113:onMode = off");
+ fPythia8->ReadString("113:onIfAll = 13 13");
+// Eta-> mu+ mu-
+ fPythia8->ReadString("221:onMode = off");
+ fPythia8->ReadString("221:onIfAll = 13 13");
+// omega-> mu+ mu-
+ fPythia8->ReadString("223:onMode = off");
+ fPythia8->ReadString("223:onIfAll = 13 13");
+// phi-> mu+ mu-
+ fPythia8->ReadString("333:onMode = off");
+ fPythia8->ReadString("333:onIfAll = 13 13");
+// J/Psi-> mu+ mu-
+ fPythia8->ReadString("443:onMode = off");
+ fPythia8->ReadString("443:onIfAll = 13 13");
+// Psi'-> mu+ mu-
+ fPythia8->ReadString("100443:onMode = off");
+ fPythia8->ReadString("100443:onIfAll = 13 13");
+// Ups-> mu+ mu-
+ fPythia8->ReadString("553:onMode = off");
+ fPythia8->ReadString("553:onIfAll = 13 13");
+// Ups'-> mu+ mu-
+ fPythia8->ReadString("100553:onMode = off");
+ fPythia8->ReadString("100553:onIfAll = 13 13");
+// Ups''-> mu+ mu-
+ fPythia8->ReadString("200553:onMode = off");
+ fPythia8->ReadString("200553:onIfAll = 13 13");
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
+// B0 - > e+ e-
+ fPythia8->ReadString("511:onMode = off");
+ fPythia8->ReadString("511:onIfAny = 11");
+// B+- -> e+ e-
+ fPythia8->ReadString("521:onMode = off");
+ fPythia8->ReadString("521:onIfAny = 11");
+// B_s -> e+ e-
+ fPythia8->ReadString("531:onMode = off");
+ fPythia8->ReadString("531:onIfAny = 11");
+// Lambda_b -> e+ e-
+ fPythia8->ReadString("5122:onMode = off");
+ fPythia8->ReadString("5122:onIfAny = 11");
+// Sigma_b -> e+ e-
+ fPythia8->ReadString("5132:onMode = off");
+ fPythia8->ReadString("5132:onIfAny = 11");
+// Sigma_b -> e+ e-
+ fPythia8->ReadString("5232:onMode = off");
+ fPythia8->ReadString("5232:onIfAny = 11");
+// Omega_b ->e+ e-
+ fPythia8->ReadString("5332:onMode = off");
+ fPythia8->ReadString("5332:onIfAny = 11");
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
+// D+/- -> e X
+ fPythia8->ReadString("411:onMode = off");
+ fPythia8->ReadString("411:onIfAll = 11");
+// D0 -> e X
+ fPythia8->ReadString("421:onMode = off");
+ fPythia8->ReadString("421:onIfAll = 11");
+// D_s ->e X
+ fPythia8->ReadString("431:onMode = off");
+ fPythia8->ReadString("431:onIfAll = 11");
+// Lambda_c -> e X
+ fPythia8->ReadString("4122:onMode = off");
+ fPythia8->ReadString("4122:onIfAll = 11");
+// Sigma_c -> e X
+ fPythia8->ReadString("4132:onMode = off");
+ fPythia8->ReadString("4132:onIfAll = 11");
+// Sigma_c -> e X
+ fPythia8->ReadString("4232:onMode = off");
+ fPythia8->ReadString("4232:onIfAll = 11");
+// Omega_c -> e X
+ fPythia8->ReadString("4332:onMode = off");
+ fPythia8->ReadString("4332:onIfAll = 11");
+// B0 -> e X
+ fPythia8->ReadString("511:onMode = off");
+ fPythia8->ReadString("511:onIfAny = 11");
+// B+/- -> e X
+ fPythia8->ReadString("521:onMode = off");
+ fPythia8->ReadString("521:onIfAny = 11");
+// B_s -> e X
+ fPythia8->ReadString("531:onMode = off");
+ fPythia8->ReadString("531:onIfAny = 11");
+// Lambda_b -> e X
+ fPythia8->ReadString("5122:onMode = off");
+ fPythia8->ReadString("5122:onIfAny = 11");
+// Sigma_b -> e X
+ fPythia8->ReadString("5132:onMode = off");
+ fPythia8->ReadString("5132:onIfAny = 11");
+// Sigma_b -> e X
+ fPythia8->ReadString("5232:onMode = off");
+ fPythia8->ReadString("5232:onIfAny = 11");
+// Omega_b -> e X
+ fPythia8->ReadString("5332:onMode = off");
+ fPythia8->ReadString("5332:onIfAny = 11");
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''
+// Rho -> e+e-
+ fPythia8->ReadString("113:onMode = off");
+ fPythia8->ReadString("113:onIfAll = 11 11");
+// Eta -> e+e-
+ fPythia8->ReadString("221:onMode = off");
+ fPythia8->ReadString("221:onIfAll = 11 11");
+// omega -> e+e-
+ fPythia8->ReadString("223:onMode = off");
+ fPythia8->ReadString("223:onIfAll = 11 11");
+// phi -> e+e-
+ fPythia8->ReadString("333:onMode = off");
+ fPythia8->ReadString("333:onIfAll = 11 11");
+// J/Psi -> e+e-
+ fPythia8->ReadString("443:onMode = off");
+ fPythia8->ReadString("443:onIfAll = 11 11");
+// Psi' -> e+e-
+ fPythia8->ReadString("100443:onMode = off");
+ fPythia8->ReadString("100443:onIfAll = 11 11");
+// Ups -> e+e-
+ fPythia8->ReadString("553:onMode = off");
+ fPythia8->ReadString("553:onIfAll = 11 11");
+// Ups' -> e+e-
+ fPythia8->ReadString("100553:onMode = off");
+ fPythia8->ReadString("100553:onIfAll = 11 11");
+// Ups'' -> e+e-
+ fPythia8->ReadString("200553:onMode = off");
+ fPythia8->ReadString("200553:onIfAll = 11 11");
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-
+// B0 -> J/Psi (Psi') X
+ fPythia8->ReadString("511:onMode = off");
+ fPythia8->ReadString("511:onIfAny = 443 100443");
+// B+/- -> J/Psi (Psi') X
+ fPythia8->ReadString("521:onMode = off");
+ fPythia8->ReadString("521:onIfAny = 443 100443");
+// B_s -> J/Psi (Psi') X
+ fPythia8->ReadString("531:onMode = off");
+ fPythia8->ReadString("531:onIfAny = 443 100443");
+// Lambda_b -> J/Psi (Psi') X
+ fPythia8->ReadString("5122:onMode = off");
+ fPythia8->ReadString("5122:onIfAny = 443 100443");
+//
+// J/Psi -> mu+ mu-
+ fPythia8->ReadString("443:onMode = off");
+ fPythia8->ReadString("443:onIfAll = 13 13");
+// Psi' -> mu+ mu-
+ fPythia8->ReadString("100443:onMode = off");
+ fPythia8->ReadString("100443:onIfAll = 13 13");
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'
+// B0 -> Psi' X
+ fPythia8->ReadString("511:onMode = off");
+ fPythia8->ReadString("511:onIfAny = 100443");
+// B+/- -> Psi' X
+ fPythia8->ReadString("521:onMode = off");
+ fPythia8->ReadString("521:onIfAny = 100443");
+// B_s -> Psi' X
+ fPythia8->ReadString("531:onMode = off");
+ fPythia8->ReadString("531:onIfAny = 100443");
+// Lambda_b -> Psi' X
+ fPythia8->ReadString("5122:onMode = off");
+ fPythia8->ReadString("5122:onIfAny = 100443");
+//
+// Psi' -> mu+ mu-
+ fPythia8->ReadString("100443:onMode = off");
+ fPythia8->ReadString("100443:onIfAll = 13 13");
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
+// B0 -> Psi X
+ fPythia8->ReadString("511:onMode = off");
+ fPythia8->ReadString("511:onIfAny = 443");
+// B+/- -> Psi X
+ fPythia8->ReadString("521:onMode = off");
+ fPythia8->ReadString("521:onIfAny = 443");
+// B_s -> Psi X
+ fPythia8->ReadString("531:onMode = off");
+ fPythia8->ReadString("531:onIfAny = 443");
+// Lambda_b -> Psi X
+ fPythia8->ReadString("5122:onMode = off");
+ fPythia8->ReadString("5122:onIfAny = 443");
+//
+// Psi -> mu+ mu-
+ fPythia8->ReadString("443:onMode = off");
+ fPythia8->ReadString("443:onIfAll = 11 11");
+
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
+// B0 -> Psi X
+ fPythia8->ReadString("511:onMode = off");
+ fPythia8->ReadString("511:onIfAny = 443");
+// B+/- -> Psi X
+ fPythia8->ReadString("521:onMode = off");
+ fPythia8->ReadString("521:onIfAny = 443");
+// B_s -> Psi X
+ fPythia8->ReadString("531:onMode = off");
+ fPythia8->ReadString("531:onIfAny = 443");
+// Lambda_b -> Psi X
+ fPythia8->ReadString("5122:onMode = off");
+ fPythia8->ReadString("5122:onIfAny = 443");
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'
+// B0 -> Psi' X
+ fPythia8->ReadString("511:onMode = off");
+ fPythia8->ReadString("511:onIfAny = 100443");
+// B+/- -> Psi' X
+ fPythia8->ReadString("521:onMode = off");
+ fPythia8->ReadString("521:onIfAny = 100443");
+// B_s -> Psi' X
+ fPythia8->ReadString("531:onMode = off");
+ fPythia8->ReadString("531:onIfAny = 100443");
+// Lambda_b -> Psi' X
+ fPythia8->ReadString("5122:onMode = off");
+ fPythia8->ReadString("5122:onIfAny = 100443");
+//
+// Psi' -> mu+ mu-
+ fPythia8->ReadString("100443:onMode = off");
+ fPythia8->ReadString("100443:onIfAll = 11 11");
break;
case kPiToMu:
- ForceParticleDecay(211,13,1); // pi->mu
+// pi -> mu nu
+ fPythia8->ReadString("211:onMode = off");
+ fPythia8->ReadString("211:onIfAny = 13");
break;
case kKaToMu:
- ForceParticleDecay(321,13,1); // K->mu
+// K -> mu nu
+ fPythia8->ReadString("321:onMode = off");
+ fPythia8->ReadString("321:onIfAny = 13");
break;
case kAllMuonic:
- ForceParticleDecay(211,13,1); // pi->mu
- ForceParticleDecay(321,13,1); // K->mu
+// pi/K -> mu
+ fPythia8->ReadString("211:onMode = off");
+ fPythia8->ReadString("211:onIfAny = 13");
+ fPythia8->ReadString("321:onMode = off");
+ fPythia8->ReadString("321:onIfAny = 13");
break;
case kWToMuon:
- ForceParticleDecay( 24, 13,1); // W -> mu
+// W -> mu X
+ fPythia8->ReadString("24:onMode = off");
+ fPythia8->ReadString("24:onIfAny = 13");
break;
case kWToCharm:
- ForceParticleDecay( 24, 4,1); // W -> c
+// W -> c X
+ fPythia8->ReadString("24:onMode = off");
+ fPythia8->ReadString("24:onIfAny = 4");
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
+// W -> c X
+ fPythia8->ReadString("24:onMode = off");
+ fPythia8->ReadString("24:onIfAny = 4");
+// D+- -> mu X
+ fPythia8->ReadString("411:onMode = off");
+ fPythia8->ReadString("411:onIfAll = 13");
+// D0 -> mu X
+ fPythia8->ReadString("421:onMode = off");
+ fPythia8->ReadString("421:onIfAll = 13");
+// D_s -> mu X
+ fPythia8->ReadString("431:onMode = off");
+ fPythia8->ReadString("431:onIfAll = 13");
+// Lambda_c -> mu X
+ fPythia8->ReadString("4122:onMode = off");
+ fPythia8->ReadString("4122:onIfAll = 13");
+// Sigma_c -> mu X
+ fPythia8->ReadString("4132:onMode = off");
+ fPythia8->ReadString("4132:onIfAll = 13");
+// Sigma_c -> mu X
+ fPythia8->ReadString("4232:onMode = off");
+ fPythia8->ReadString("4232:onIfAll = 13");
+// Omega_c -> mu X
+ fPythia8->ReadString("4332:onMode = off");
+ fPythia8->ReadString("4332:onIfAll = 13");
break;
case kZDiMuon:
- ForceParticleDecay( 23, 13,2); // Z -> mu+ mu-
+// Z -> mu+ mu-
+ fPythia8->ReadString("23:onMode = off");
+ fPythia8->ReadString("23:onIfAll = 13 13");
break;
case kZDiElectron:
- ForceParticleDecay( 23, 11,2); // Z -> e+ e-
+// Z -> e+ e-
+ fPythia8->ReadString("23:onMode = off");
+ fPythia8->ReadString("23:onIfAll = 11 11");
break;
case kHadronicD:
ForceHadronicD(1);
ForceHadronicD(0);
break;
case kPhiKK:
- ForceParticleDecay(333,321,2); // Phi->K+K-
+ // Phi-> K+ K-
+ fPythia8->ReadString("333:onMode = off");
+ fPythia8->ReadString("333:onIfAll = 321 321");
break;
case kOmega:
-// ForceOmega();
+ // Omega -> Lambda K
+ fPythia8->ReadString("3334:onMode = off");
+ fPythia8->ReadString("3334:onIfAll = 3122 321 ");
+ break;
+ case kLambda:
+ // Lambda -> p pi-
+ fPythia8->ReadString("3122:onMode = off");
+ fPythia8->ReadString("3122:onIfAll = 2212 211 ");
+ break;
+ case kBeautyUpgrade:
+ fPythia8->ReadString("5122:onMode = off");
+ fPythia8->ReadString("4122:onMode = off");
+ fPythia8->ReadString("5122:onIfAll = 4122");
+ fPythia8->ReadString("4122:onIfAll = 3122");
+ break;
case kAll:
break;
case kNoDecay:
- TPythia8::Instance()->ReadString("HadronLevel:Decay = off");
+ fPythia8->ReadString("HadronLevel:Decay = off");
break;
case kNoDecayHeavy:
+ case kNoDecayBeauty:
case kNeutralPion:
+ case kPsiPrimeJpsiDiElectron:
+ case kElectronEM:
+ case kGammaEM:
+ case kDiElectronEM:
break;
}
}
{
// 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();
+ Pythia8::Pythia* thePythia = fPythia8->Pythia8();
+ Pythia8::ParticleData & table = thePythia->particleData;
+ Pythia8::ParticleDataEntry* pd = table.particleDataEntryPtr(ipart);
+
+ Int_t nc = pd->sizeChannels();
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++) {
+ Pythia8::DecayChannel& decCh = pd->channel(ic);
+ for (Int_t i = 0; i < decCh.multiplicity(); i++) {
br += decCh.bRatio();
}
}
}
-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;
+ Pythia8::Pythia* thePythia = fPythia8->Pythia8();
+ Pythia8::ParticleData& table = thePythia->particleData;
Float_t tau = table.tau0(kf);
- return ( tau);
+ return (tau);
}
void AliDecayerPythia8::SwitchOffHeavyFlavour()
// Switch off heavy flavour production
//
// Maximum number of quark flavours used in pdf
- TPythia8::Instance()->ReadString("PDFinProcess:nQuarkIn = 3");
+ fPythia8->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");
+ fPythia8->ReadString("SpaceShower:nQuarkIn = 3");
+ fPythia8->ReadString("TimeShower:nGammaToQuark = 3");
+ fPythia8->ReadString("TimeShower:nGluonToQuark = 3");
}
//
// 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);
+ // K* -> K pi
+ fPythia8->ReadString("313:onMode = off");
+ fPythia8->ReadString("313:onIfAll = 321 211");
// for Ds -> Phi pi+
- Int_t iPhi = 333;
- ForceParticleDecay(iPhi, kKPlus, 2); // Phi->K+K-
+ fPythia8->ReadString("333:onMode = off");
+ fPythia8->ReadString("333:onIfAll = 321 321");
// for D0 -> rho0 pi+ k-
- Int_t iRho0=113;
- ForceParticleDecay(iRho0, kPiPlus, 2); // Rho0->pi+pi-
+ fPythia8->ReadString("113:onMode = off");
+ fPythia8->ReadString("113:onIfAll = 211 211");
// 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;
- }
- }
+ fPythia8->ReadString("2224:onMode = off");
+ fPythia8->ReadString("2224:onIfAll = 2212 211");
+ // for Lambda_c -> Lambda(1520) K-
+ fPythia8->ReadString("3124:onMode = off");
+ fPythia8->ReadString("3124:onIfAll = 2212 321");
- Pythia8::Pythia* thePythia = TPythia8::Instance()->Pythia8();
- Pythia8::ParticleDataTable table = thePythia->particleData;
+ fPythia8->ReadString("411:onMode = off");
+ fPythia8->ReadString("421:onMode = off");
+ fPythia8->ReadString("431:onMode = off");
+ fPythia8->ReadString("4112:onMode = off");
+ fPythia8->ReadString("4122:onMode = off");
+
+ // D+/- -> K pi pi
+ fPythia8->ReadString("411:onIfMatch = 321 211 211");
+ // D+/- -> K* pi
+ fPythia8->ReadString("411:onIfMatch = 313 211");
+ // D0 -> K pi
+ fPythia8->ReadString("421:onIfMatch = 321 211");
+
+ if (optUse4Bodies) {
+ // D0 -> K pi pi pi
+ fPythia8->ReadString("421:onIfMatch = 321 211 211 211");
+ // D0 -> K pi rho
+ fPythia8->ReadString("421:onIfMatch = 321 211 113");
+ // D0 -> K*0 pi pi
+ fPythia8->ReadString("421:onIfMatch = 313 211 211");
+ }
- 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
+ // D_s -> K K*
+ fPythia8->ReadString("431:onIfMatch = 321 313");
+ // D_s -> Phi pi
+ fPythia8->ReadString("431:onIfMatch = 333 211");
+
+ // Lambda_c -> p K*
+ fPythia8->ReadString("4122:onIfMatch = 2212 313");
+ // Lambda_c -> Delta K
+ fPythia8->ReadString("4122:onIfMatch = 2224 321");
+ // Lambda_c -> Lambda(1520) pi
+ fPythia8->ReadString("4122:onIfMatch = 3124 211");
+ // Lambda_c -> p K pi
+ fPythia8->ReadString("4122:onIfMatch = 2212 321 211");
+ // Lambda_c -> Lambda pi
+ fPythia8->ReadString("4122:onIfMatch = 3122 211");
+
+}
+
+//___________________________________________________________________________
+void AliDecayerPythia8::ReadDecayTable()
+{
+ //to read a decay table (not yet implemented)
+}
+
+
+//___________________________________________________________________________
+void AliDecayerPythia8::AppendParticle(Int_t pdg, TLorentzVector* p)
+{
+ // Append a particle to the stack
+ fPythia8->Pythia8()->event.append(pdg, 11, 0, 0, p->Px(), p->Py(), p->Pz(), p->E(), p->M());
+}
+
+
+//___________________________________________________________________________
+void AliDecayerPythia8::ClearEvent()
+{
+ // Clear the event stack
+ fPythia8->Pythia8()->event.clear();
}