/*
$Log$
+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.
#include "AliDecayerPythia.h"
#include "AliPythia.h"
+#include <TPDGCode.h>
#include <TLorentzVector.h>
#include <TClonesArray.h>
#ifndef WIN32
# define py1ent py1ent_
+# define opendecaytable opendecaytable_
# define type_of_call
#else
# define lu1ent PY1ENT
+# define opendecaytable OPENDECAYTABLE
# define type_of_call _stdcall
#endif
extern "C" void type_of_call
py1ent(Int_t&, Int_t&, Double_t&, Double_t&, Double_t&);
+extern "C" void type_of_call
+ opendecaytable(Int_t&);
+
Bool_t AliDecayerPythia::fgInit = kFALSE;
// Switch on heavy flavor decays
Int_t kc, i, j;
- Int_t heavy[8] = {411, 421, 431, 4122, 511, 521, 531, 5122};
+ Int_t heavy[14] = {411, 421, 431, 4122, 4132, 4232, 4332, 511, 521, 531, 5122, 5132, 5232, 5332};
fPythia->ResetDecayTable();
- for (j=0; j < 8; j++) {
+ for (j=0; j < 14; j++) {
kc=fPythia->Pycomp(heavy[j]);
- fPythia->SetMDCY(kc,1,1);
- for (i=fPythia->GetMDCY(kc,2);
- i<fPythia->GetMDCY(kc,2)+fPythia->GetMDCY(kc,3);
- i++) {
+ if (fDecay == kNoDecayHeavy) {
+ fPythia->SetMDCY(kc,1,0);
+ } else {
+ fPythia->SetMDCY(kc,1,1);
+ for (i=fPythia->GetMDCY(kc,2);
+ i<fPythia->GetMDCY(kc,2)+fPythia->GetMDCY(kc,3);
+ i++) {
fPythia->SetMDME(i,1,1);
+ }
}
}
Lu1Ent(0, idpart, energy, theta, phi);
fPythia->GetPrimaries();
+// fPythia->Pylist(1);
+
}
Int_t AliDecayerPythia::ImportParticles(TClonesArray *particles)
{
// Force a particle decay mode
Decay_t decay=fDecay;
-
+ fPythia->SetMSTJ(21,2);
+ if (decay == kNoDecayHeavy) return;
+
//
// select mode
ForceParticleDecay( 411,13,1); // D+/-
ForceParticleDecay( 421,13,1); // D0
ForceParticleDecay( 431,13,1); // D_s
- ForceParticleDecay( 4122,13,1); // Lambda_c
+ 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
- break;
+ ForceParticleDecay( 5132,13,1); // Xsi_b
+ ForceParticleDecay( 5232,13,1); // Sigma_b
+ ForceParticleDecay( 5332,13,1); // Omega_b
+ break;
case kDiMuon:
+ ForceParticleDecay( 221,13,2); // eta
+ ForceParticleDecay( 223,13,2); // omega
ForceParticleDecay( 333,13,2); // phi
ForceParticleDecay( 443,13,2); // J/Psi
- ForceParticleDecay(20443,13,2); // Psi'
+ ForceParticleDecay(100443,13,2); // Psi'
ForceParticleDecay( 553,13,2); // Upsilon
- ForceParticleDecay(20553,13,2); // Upsilon'
- ForceParticleDecay(30553,13,2); // Upsilon''
+ ForceParticleDecay(100553,13,2); // Upsilon'
+ ForceParticleDecay(200553,13,2); // Upsilon''
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( 333,11,2); // phi
+ ForceParticleDecay( 221,11,2); // eta
+ ForceParticleDecay( 223,11,2); // omega
ForceParticleDecay( 443,11,2); // J/Psi
- ForceParticleDecay(30443,11,2); // Psi'
+ ForceParticleDecay(100443,11,2); // Psi'
ForceParticleDecay( 553,11,2); // Upsilon
- ForceParticleDecay(30553,11,2); // Upsilon'
+ ForceParticleDecay(100553,11,2); // Upsilon'
+ ForceParticleDecay(200553,11,2); // Upsilon''
break;
case kBJpsiDiMuon:
ForceParticleDecay( 511,443,1); // B0
ForceParticleDecay( 521,30443,1); // B+/-
ForceParticleDecay( 531,30443,1); // B_s
ForceParticleDecay( 5122,30443,1); // Lambda_b
- ForceParticleDecay(30443,13,2); // Psi'
+ ForceParticleDecay(100443,13,2); // Psi'
break;
case kBJpsiDiElectron:
ForceParticleDecay( 511,443,1); // B0
ForceParticleDecay( 521,30443,1); // B+/-
ForceParticleDecay( 531,30443,1); // B_s
ForceParticleDecay( 5122,30443,1); // Lambda_b
- ForceParticleDecay(30443,11,2); // Psi'
+ ForceParticleDecay(100443,11,2); // Psi'
break;
case kPiToMu:
ForceParticleDecay(211,13,1); // pi->mu
break;
case kHadronicD:
ForceHadronicD();
+ break;
+ case kPhiKK:
+ ForceParticleDecay(333,321,2); // Phi->K+K-
+ break;
+ case kOmega:
+ ForceOmega();
case kAll:
break;
case kNoDecay:
+ fPythia->SetMSTJ(21,0);
+ break;
+ case kNoDecayHeavy:
break;
}
}
void AliDecayerPythia::ForceHadronicD()
{
+//
// Force golden D decay modes
//
-
+ const Int_t kNHadrons = 4;
Int_t channel;
+ Int_t hadron[kNHadrons] = {411, 421, 431, 4112};
+ Int_t decayP[kNHadrons][3] =
+ {
+ {kKMinus, kPiPlus, kPiPlus},
+ {kKMinus, kPiPlus, 0 },
+ {-1 , -1 , -1 },
+ {-1 , -1 , -1 }
+ };
-// D+ -> K- pi+ pi+
- Int_t kc=fPythia->Pycomp(411);
- 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 (channel==837) {
- fPythia->SetMDME(channel,1,1);
- } else {
- fPythia->SetMDME(channel,1,0);
- fBraPart[kc]-=fPythia->GetBRAT(channel);
- }
- }
-
-// D0 -> K- pi+
- kc=fPythia->Pycomp(421);
- fPythia->SetMDCY(kc,1,1);
- ifirst=fPythia->GetMDCY(kc,2);
- ilast=ifirst+fPythia->GetMDCY(kc,3)-1;
- for (channel=ifirst; channel<=ilast;channel++) {
- if (channel==881) {
- fPythia->SetMDME(channel,1,1);
- } else {
- fPythia->SetMDME(channel,1,0);
- fBraPart[kc]-=fPythia->GetBRAT(channel);
- }
- }
-
-// no D_s decays
- kc=fPythia->Pycomp(431);
- fPythia->SetMDCY(kc,1,1);
- ifirst=fPythia->GetMDCY(kc,2);
- ilast=ifirst+fPythia->GetMDCY(kc,3)-1;
- for (channel=ifirst; channel<=ilast;channel++) {
- fPythia->SetMDME(channel,1,0);
- fBraPart[kc]-=fPythia->GetBRAT(channel);
- }
-// no Lambda_c decays
- kc=fPythia->Pycomp(4122);
- fPythia->SetMDCY(kc,1,1);
- ifirst=fPythia->GetMDCY(kc,2);
- ilast=ifirst+fPythia->GetMDCY(kc,3)-1;
- for (channel=ifirst; channel<=ilast;channel++) {
- fPythia->SetMDME(channel,1,0);
- fBraPart[kc]-=fPythia->GetBRAT(channel);
- }
+ 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) == decayP[ihadron][0] &&
+ fPythia->GetKFDP(channel,2) == decayP[ihadron][1] &&
+ fPythia->GetKFDP(channel,3) == decayP[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
+ } // hadrons
}
void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
// gMC->Gspart(113,"Phi",3,mass,0,tlife);
}
+void AliDecayerPythia::ForceOmega()
+{
+ // Force Omega -> Lambda K- Decay
+ Int_t kc=fPythia->Pycomp(3334);
+ fPythia->SetMDCY(kc,1,1);
+ Int_t ifirst = fPythia->GetMDCY(kc,2);
+ Int_t ilast = ifirst + fPythia->GetMDCY(kc,3)-1;
+ for (Int_t channel = ifirst; channel <= ilast; channel++) {
+ if (
+ fPythia->GetKFDP(channel,1) == kLambda0 &&
+ fPythia->GetKFDP(channel,2) == kKMinus &&
+ fPythia->GetKFDP(channel,3) == 0
+ )
+ {
+ fPythia->SetMDME(channel,1,1);
+ } else {
+ fPythia->SetMDME(channel,1,0);
+ } // selected channel ?
+ } // decay channels
+}
+
+
+
Float_t AliDecayerPythia::GetPartialBranchingRatio(Int_t kf)
{
// Get branching ratio
return fBraPart[kc];
}
+Float_t AliDecayerPythia::GetLifetime(Int_t kf)
+{
+// Get branching ratio
+ Int_t kc=fPythia->Pycomp(TMath::Abs(kf));
+ return fPythia->GetPMAS(kc,4)*3.3333e-12;
+}
+
+void AliDecayerPythia::WriteDecayTable()
+{
+//
+// Write the decay table
+ fPythia->Pyupda(1,15);
+}
+
+void AliDecayerPythia::ReadDecayTable()
+{
+//
+// Read the decay table
+ printf("Reading Deacy table\n \n ");
+
+ Int_t lun = 15;
+ opendecaytable(lun);
+ fPythia->Pyupda(2,lun);
+}
+
#ifdef never
void AliDecayerPythia::Streamer(TBuffer &R__b)
{