AliDecayerPythia::AliDecayerPythia():
fPythia(AliPythia::Instance()),
- fDecay(kAll)
+ fDecay(kAll),
+ fHeavyFlavour(kTRUE)
{
// Default Constructor
for (Int_t i=0; i< 501; i++) fBraPart[i]= 1.;
AliDecayerPythia::AliDecayerPythia(const AliDecayerPythia &decayer):
AliDecayer(decayer),
fPythia(0),
- fDecay(kAll)
+ fDecay(kAll),
+ fHeavyFlavour(kTRUE)
{
// Copy Constructor
decayer.Copy(*this);
}
}
//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
-
+
if (fDecay != kNeutralPion) {
- fPythia->SetMDCY(fPythia->Pycomp(111) ,1,0);
+ fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
} else {
- fPythia->SetMDCY(fPythia->Pycomp(111) ,1,1);
+ fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
}
-
-
+
+
fPythia->SetMDCY(fPythia->Pycomp(310) ,1,0);
fPythia->SetMDCY(fPythia->Pycomp(3122),1,0);
fPythia->SetMDCY(fPythia->Pycomp(3112),1,0);
fPythia->SetMDCY(fPythia->Pycomp(3312),1,0);
fPythia->SetMDCY(fPythia->Pycomp(3322),1,0);
fPythia->SetMDCY(fPythia->Pycomp(3334),1,0);
+
// .. Force decay channels
ForceDecay();
}
Lu1Ent(0, idpart, energy, theta, phi);
fPythia->GetPrimaries();
-// fPythia->Pylist(1);
-
}
Int_t AliDecayerPythia::ImportParticles(TClonesArray *particles)
void AliDecayerPythia::ForceDecay()
{
// Force a particle decay mode
+// Switch heavy flavour production off if requested
+ if (!fHeavyFlavour) SwitchOffHeavyFlavour();
+//
Decay_t decay=fDecay;
fPythia->SetMSTJ(21,2);
if (decay == kNoDecayHeavy) return;
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( 443, 13, 2); // J/Psi -> mu+ mu-
ForceParticleDecay( 411,13,1); // D+/-
ForceParticleDecay( 421,13,1); // D0
ForceParticleDecay( 431,13,1); // D_s
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+/-
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 kZDiMuon:
ForceParticleDecay( 23, 13,2); // Z -> mu+ mu-
break;
+ case kZDiElectron:
+ ForceParticleDecay( 23, 11,2); // Z -> e+ e-
+ break;
case kHadronicD:
- ForceHadronicD();
+ ForceHadronicD(1);
+ break;
+ case kHadronicDWithout4Bodies:
+ ForceHadronicD(0);
break;
case kPhiKK:
ForceParticleDecay(333,321,2); // Phi->K+K-
case kNoDecayHeavy:
case kNeutralPion:
break;
-
}
}
+void AliDecayerPythia::SwitchOffHeavyFlavour()
+{
+ // Switch off heavy flavour production
+ //
+ // Maximum number of quark flavours used in pdf
+ fPythia->SetMSTP(58, 3);
+ // Maximum number of flavors that can be used in showers
+ fPythia->SetMSTJ(45, 3);
+ // Switch off g->QQbar splitting in decay table
+ for (Int_t i = 156; i <= 160; i++) fPythia->SetMDME(i, 1, 0);
+}
+
void AliDecayerPythia::Lu1Ent(Int_t flag, Int_t idpart,
Double_t mom, Double_t theta, Double_t phi)
{
}
-void AliDecayerPythia::ForceHadronicD()
+void AliDecayerPythia::ForceHadronicD(Int_t optUse4Bodies)
{
//
// Force golden D decay modes
//
- const Int_t kNHadrons = 4;
- Int_t channel;
- Int_t hadron[kNHadrons] = {411, 421, 431, 4112};
- // 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-
- Int_t decayP1[kNHadrons][3] =
-
- {
- {kKMinus, kPiPlus, kPiPlus},
- {kKMinus, kPiPlus, 0 },
- {kKPlus , iKstarbar0, 0 },
- {-1 , -1 , -1 }
- };
- Int_t decayP2[kNHadrons][3] =
- {
- {iKstarbar0, kPiPlus, 0 },
- {-1 , -1 , -1 },
- {iPhi , kPiPlus, 0 },
- {-1 , -1 , -1 }
- };
-
-
- for (Int_t ihadron = 0; ihadron < kNHadrons; ihadron++)
+ const Int_t kNHadrons = 5;
+ Int_t channel;
+ 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;
+ }
+ }
+
+ for (Int_t ihadron = 0; ihadron < kNHadrons; ihadron++)
{
- Int_t kc = fPythia->Pycomp(hadron[ihadron]);
- fPythia->SetMDCY(kc,1,1);
- Int_t ifirst = fPythia->GetMDCY(kc,2);
- Int_t ilast = ifirst + fPythia->GetMDCY(kc,3)-1;
-
- for (channel = ifirst; channel <= ilast; channel++) {
- if ((
- fPythia->GetKFDP(channel,1) == decayP1[ihadron][0] &&
- fPythia->GetKFDP(channel,2) == decayP1[ihadron][1] &&
- fPythia->GetKFDP(channel,3) == decayP1[ihadron][2] &&
- fPythia->GetKFDP(channel,4) == 0
- ) || (
- fPythia->GetKFDP(channel,1) == decayP2[ihadron][0] &&
- fPythia->GetKFDP(channel,2) == decayP2[ihadron][1] &&
- fPythia->GetKFDP(channel,3) == decayP2[ihadron][2] &&
- fPythia->GetKFDP(channel,4) == 0
- ))
-
- {
- fPythia->SetMDME(channel,1,1);
- } else {
- fPythia->SetMDME(channel,1,0);
- fBraPart[kc] -= fPythia->GetBRAT(channel);
- } // selected channel ?
- } // decay channels
+ Int_t kc = fPythia->Pycomp(hadron[ihadron]);
+ fPythia->SetMDCY(kc,1,1);
+ Int_t ifirst = fPythia->GetMDCY(kc,2);
+ Int_t ilast = ifirst + fPythia->GetMDCY(kc,3)-1;
+
+ for (channel = ifirst; channel <= ilast; channel++) {
+ if ((
+ fPythia->GetKFDP(channel,1) == decayP1[ihadron][0] &&
+ fPythia->GetKFDP(channel,2) == decayP1[ihadron][1] &&
+ fPythia->GetKFDP(channel,3) == decayP1[ihadron][2] &&
+ fPythia->GetKFDP(channel,4) == decayP1[ihadron][3] &&
+ fPythia->GetKFDP(channel,5) == 0
+ ) || (
+ fPythia->GetKFDP(channel,1) == decayP2[ihadron][0] &&
+ fPythia->GetKFDP(channel,2) == decayP2[ihadron][1] &&
+ fPythia->GetKFDP(channel,3) == decayP2[ihadron][2] &&
+ fPythia->GetKFDP(channel,4) == decayP2[ihadron][3] &&
+ fPythia->GetKFDP(channel,5) == 0
+ ) || (
+ fPythia->GetKFDP(channel,1) == decayP3[ihadron][0] &&
+ fPythia->GetKFDP(channel,2) == decayP3[ihadron][1] &&
+ fPythia->GetKFDP(channel,3) == decayP3[ihadron][2] &&
+ fPythia->GetKFDP(channel,4) == decayP3[ihadron][3] &&
+ fPythia->GetKFDP(channel,5) == 0
+ ))
+
+ {
+ fPythia->SetMDME(channel,1,1);
+ } else {
+ fPythia->SetMDME(channel,1,0);
+ fBraPart[kc] -= fPythia->GetBRAT(channel);
+ } // selected channel ?
+ } // decay channels
} // hadrons
}
+
void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
{
//