/*
$Log$
+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.
#include "AliDecayerPythia.h"
#include "AliPythia.h"
+#include <TPDGCode.h>
#include <TLorentzVector.h>
+#include <TClonesArray.h>
ClassImp(AliDecayerPythia)
#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;
+
AliDecayerPythia::AliDecayerPythia()
{
// Default Constructor
fPythia=AliPythia::Instance();
- for (Int_t i=0; i< 501; i++) fBraPart[i]=1;
+ for (Int_t i=0; i< 501; i++) fBraPart[i]= 1.;
}
void AliDecayerPythia::Init()
{
// Initialisation
//
+ if (!fgInit) {
+ fgInit = kTRUE;
+ fPythia->SetDecayTable();
+ }
+
// Switch on heavy flavor decays
- Int_t kc, i, j;
- Int_t heavy[8] = {411, 421, 431, 4122, 511, 521, 531, 5122};
- for (j=0; j < 8; j++) {
+ Int_t kc, i, 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++) {
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);
+ }
}
}
+
ForceDecay();
}
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;
-
-//
-// Make clean
-// AllowAllDecays();
+ fPythia->SetMSTJ(21,2);
+ if (decay == kNoDecayHeavy) return;
+
//
// select mode
switch (decay)
{
- case semimuonic:
+ 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( 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;
- case dimuon:
-// ForceParticleDecay( 41,13,2); // phi
+ 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 semielectronic:
-
+ 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 dielectron:
-
- ForceParticleDecay( 41,11,2); // phi
+ 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 b_jpsi_dimuon:
+ 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
break;
- case b_psip_dimuon:
+ 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(30443,13,2); // Psi'
+ ForceParticleDecay(100443,13,2); // Psi'
break;
- case b_jpsi_dielectron:
+ 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 b_psip_dielectron:
+ 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(30443,11,2); // Psi'
+ ForceParticleDecay(100443,11,2); // Psi'
break;
- case pitomu:
+ case kPiToMu:
ForceParticleDecay(211,13,1); // pi->mu
break;
- case katomu:
+ case kKaToMu:
ForceParticleDecay(321,13,1); // K->mu
break;
- case hadronicD:
+ case kHadronicD:
ForceHadronicD();
- case all:
break;
- case nodecay:
+ 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)
}
}
-
-void AliDecayerPythia::AllowAllDecays()
-{
-// Reset decay flags
- Int_t i;
- for (i=1; i<= 2000; i++) {
- fPythia->SetMDME(i,1,1);
- }
-//
- for (i=0; i<501; i++){
- fBraPart[i]=1;
- }
-}
-
-
void AliDecayerPythia::DefineParticles()
{
//
// 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)
{
// Stream an object of class AliDecayerPythia.
R__b.WriteArray(fBraPart, 501);
}
}
-
+#endif
void AliDecayerPythia::Copy(AliDecayerPythia &decayer) const
{
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff