* provided "as is" without express or implied warranty. *
**************************************************************************/
-/*
-$Log$
-Revision 1.2 2003/04/25 09:55:23 morsch
-Initialize from decay table in constructor.
-
-Revision 1.1 2003/03/15 15:00:47 morsch
-Classed imported from EVGEN.
-
-Revision 1.18 2003/02/12 09:09:19 morsch
-Pylist removed.
-
-Revision 1.17 2003/01/31 16:54:38 morsch
-Use TPDCCode.h instead of AliPDG.
-
-Revision 1.16 2003/01/31 15:56:42 morsch
-Forcing of decay channels independent of order in decay table.
-
-Revision 1.15 2002/10/14 14:55:35 hristov
-Merging the VirtualMC branch to the main development branch (HEAD)
-
-Revision 1.14 2002/10/11 10:05:18 morsch
-pdg code for psi' corrected.
-
-Revision 1.10.6.3 2002/10/10 16:40:08 hristov
-Updating VirtualMC to v3-09-02
-
-Revision 1.13 2002/09/16 10:40:48 morsch
-Use correct pdg codes for Upsilon(2S) = 100553 and Upsilon(3S) = 200553.
-
-Revision 1.12 2002/06/05 14:05:46 morsch
-Decayer option kPhiKK for forced phi->K+K- decay added.
-
-Revision 1.11 2002/04/26 10:32:59 morsch
-Option kNoDecayHeavy added.
-
-Revision 1.10 2002/02/22 17:28:05 morsch
-ReadDecayTable() and WriteDecayTable() methods added.
-
-Revision 1.9 2001/12/20 10:37:13 morsch
-- Add omega forced decay.
-- Semileptonic decays for some more B and D baryons.
-
-Revision 1.8 2001/07/04 10:28:20 morsch
-Introduce GetLifetime(Int_T kf) method until functionality provided by
-TParticlePDG.
-
-Revision 1.7 2001/04/12 07:23:28 morsch
-Reactivate forcing option for dimuon and dielectron decay channels of phi (333).
-
-Revision 1.6 2001/03/27 10:53:26 morsch
-Save pythia default decay table at first initialization. Reload at each
-following Init() call.
-
-Revision 1.5 2001/03/09 13:04:06 morsch
-Decay_t moved to AliDecayer.h
-
-Revision 1.4 2001/01/30 09:23:11 hristov
-Streamers removed (R.Brun)
-
-Revision 1.3 2000/12/21 16:24:06 morsch
-Coding convention clean-up
-
-Revision 1.2 2000/09/12 13:58:45 morsch
-SetForceDcay(..) sets the member data storing the forced decay information.
-ForceDecay() executes the change of the decay table.
-
-Revision 1.1 2000/09/06 14:23:43 morsch
-Realisation of AliDecayer using Pythia6
-
-*/
+/* $Id$ */
// Implementation of AliDecayer using Pythia
// Method forwarding to the AliPythia instance.
#ifndef WIN32
# define py1ent py1ent_
# define opendecaytable opendecaytable_
+# define closedecaytable closedecaytable_
# define type_of_call
#else
# define lu1ent PY1ENT
# define opendecaytable OPENDECAYTABLE
+# define closedecaytable CLOSEDECAYTABLE
# define type_of_call _stdcall
#endif
extern "C" void type_of_call
opendecaytable(Int_t&);
+extern "C" void type_of_call
+ closedecaytable(Int_t&);
+
Bool_t AliDecayerPythia::fgInit = kFALSE;
-AliDecayerPythia::AliDecayerPythia()
+AliDecayerPythia::AliDecayerPythia():
+ fPythia(AliPythia::Instance()),
+ fDecay(kAll)
{
// Default Constructor
- fPythia=AliPythia::Instance();
for (Int_t i=0; i< 501; i++) fBraPart[i]= 1.;
ReadDecayTable();
}
+AliDecayerPythia::AliDecayerPythia(const AliDecayerPythia &decayer):
+ AliDecayer(decayer),
+ fPythia(0),
+ fDecay(kAll)
+{
+ // Copy Constructor
+ decayer.Copy(*this);
+}
+
void AliDecayerPythia::Init()
{
// Initialisation
}
}
}
-
+//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
+ fPythia->SetMDCY(fPythia->Pycomp(111) ,1,0);
+ 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(3212),1,0);
+ fPythia->SetMDCY(fPythia->Pycomp(3222),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();
}
//
// 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 kSemiMuonic:
ForceParticleDecay( 411,13,1); // D+/-
ForceParticleDecay( 421,13,1); // D0
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(100443,13,2);// Psi'
ForceParticleDecay( 553,13,2); // Upsilon
- ForceParticleDecay(100553,13,2); // Upsilon'
- ForceParticleDecay(200553,13,2); // Upsilon''
+ ForceParticleDecay(100553,13,2);// Upsilon'
+ ForceParticleDecay(200553,13,2);// Upsilon''
break;
case kSemiElectronic:
ForceParticleDecay( 411,11,1); // D+/-
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(100443,11,2);// Psi'
ForceParticleDecay( 553,11,2); // Upsilon
- ForceParticleDecay(100553,11,2); // Upsilon'
- ForceParticleDecay(200553,11,2); // Upsilon''
+ ForceParticleDecay(100553,11,2);// Upsilon'
+ ForceParticleDecay(200553,11,2);// Upsilon''
break;
case kBJpsiDiMuon:
- ForceParticleDecay( 511,443,1); // B0
- ForceParticleDecay( 521,443,1); // B+/-
- ForceParticleDecay( 531,443,1); // B_s
- ForceParticleDecay( 5122,443,1); // Lambda_b
- ForceParticleDecay( 443,13,2); // J/Psi
+
+ 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,30443,1); // B0
- ForceParticleDecay( 521,30443,1); // B+/-
- ForceParticleDecay( 531,30443,1); // B_s
- ForceParticleDecay( 5122,30443,1); // Lambda_b
+ 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( 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,30443,1); // B0
- ForceParticleDecay( 521,30443,1); // B+/-
- ForceParticleDecay( 531,30443,1); // B_s
- ForceParticleDecay( 5122,30443,1); // Lambda_b
- ForceParticleDecay(100443,11,2); // Psi'
+ 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
case kKaToMu:
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 kHadronicD:
ForceHadronicD();
break;
const Int_t kNHadrons = 4;
Int_t channel;
Int_t hadron[kNHadrons] = {411, 421, 431, 4112};
- Int_t decayP[kNHadrons][3] =
+ // 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 },
- {-1 , -1 , -1 },
- {-1 , -1 , -1 }
+ {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 }
};
Int_t ilast = ifirst + fPythia->GetMDCY(kc,3)-1;
for (channel = ifirst; channel <= ilast; channel++) {
- if (
- fPythia->GetKFDP(channel,1) == decayP[ihadron][0] &&
- fPythia->GetKFDP(channel,2) == decayP[ihadron][1] &&
- fPythia->GetKFDP(channel,3) == decayP[ihadron][2] &&
+ 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 {
} // hadrons
}
+
void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
{
//
-// force decay of particle into products with multiplicity mult
+// Force decay of particle into products with multiplicity mult
Int_t kc=fPythia->Pycomp(particle);
fPythia->SetMDCY(kc,1,1);
}
}
+void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t* products, Int_t* mult, Int_t npart)
+{
+//
+// Force decay of particle into products with multiplicity mult
+
+ Int_t kc=fPythia->Pycomp(particle);
+ fPythia->SetMDCY(kc,1,1);
+ Int_t ifirst=fPythia->GetMDCY(kc,2);
+ Int_t ilast=ifirst+fPythia->GetMDCY(kc,3)-1;
+ fBraPart[kc] = 1;
+//
+// Loop over decay channels
+ for (Int_t channel = ifirst; channel <= ilast; channel++) {
+ Int_t nprod = 0;
+ for (Int_t i = 0; i < npart; i++) {
+ nprod += (CountProducts(channel, products[i]) >= mult[i]);
+ }
+ if (nprod) {
+ fPythia->SetMDME(channel,1,1);
+ } else {
+ fPythia->SetMDME(channel,1,0);
+ fBraPart[kc]-=fPythia->GetBRAT(channel);
+ }
+ }
+}
+
void AliDecayerPythia::DefineParticles()
{
//
Int_t lun = 15;
opendecaytable(lun);
fPythia->Pyupda(3,lun);
+ closedecaytable(lun);
+
}
#ifdef never
}
#endif
-void AliDecayerPythia::Copy(AliDecayerPythia &decayer) const
+void AliDecayerPythia::Copy(TObject &) const
{
//
// Copy *this onto AliDecayerPythia -- not implemented
797 0 42 0.001000 e+ nu_e pi0
798 0 42 0.001000 e+ nu_e eta
799 0 42 0.001000 e+ nu_e eta'
+
800 0 42 0.001000 e+ nu_e rho0
801 0 42 0.001000 e+ nu_e omega
802 1 42 0.070000 mu+ nu_mu Kbar0
2553 1 0 0.610139 W+ e-
*/
+