#include "AliDecayerPythia.h"
#include "AliPythia.h"
+#include "AliLog.h"
#include <TPDGCode.h>
#include <TLorentzVector.h>
#include <TClonesArray.h>
+#include <TParticle.h>
ClassImp(AliDecayerPythia)
Bool_t AliDecayerPythia::fgInit = kFALSE;
-AliDecayerPythia::AliDecayerPythia()
+AliDecayerPythia::AliDecayerPythia():
+ fPythia(AliPythia::Instance()),
+ fDecay(kAll),
+ fHeavyFlavour(kTRUE),
+ fLongLived(kFALSE),
+ fPatchOmegaDalitz(0),
+ fPi0(1)
{
// 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),
+ fHeavyFlavour(kTRUE),
+ fLongLived(kFALSE),
+ fPatchOmegaDalitz(0),
+ fPi0(1)
+{
+ // Copy Constructor
+ decayer.Copy(*this);
+ for (Int_t i = 0; i < 501; i++) fBraPart[i] = 0.;
+}
+
void AliDecayerPythia::Init()
{
+
// Initialisation
//
if (!fgInit) {
}
}
}
+//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
+
+/*
+ if (fDecay != kNeutralPion) {
+ fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
+ } else {
+ fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
+ }
+*/
+ if (fPi0) fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
+ Int_t isw = 0;
+ if (fLongLived) isw = 1;
+
+ fPythia->SetMDCY(fPythia->Pycomp(310) ,1, isw);
+ fPythia->SetMDCY(fPythia->Pycomp(3122),1, isw);
+ fPythia->SetMDCY(fPythia->Pycomp(3112),1, isw);
+// fPythia->SetMDCY(fPythia->Pycomp(3212),1, isw); // Sigma0 decays elem.
+ fPythia->SetMDCY(fPythia->Pycomp(3222),1, isw);
+ fPythia->SetMDCY(fPythia->Pycomp(3312),1, isw);
+ fPythia->SetMDCY(fPythia->Pycomp(3322),1, isw);
+ fPythia->SetMDCY(fPythia->Pycomp(3334),1, isw);
+
+// .. Force decay channels
ForceDecay();
}
{
// Decay a particle
//
- Float_t energy = p->Energy();
- Float_t theta = p->Theta();
- Float_t phi = p->Phi();
-
+ Float_t energy = p->Energy();
+ Float_t theta = p->Theta();
+ Float_t phi = p->Phi();
+
+ if(!fDecayerExodus) {
Lu1Ent(0, idpart, energy, theta, phi);
- fPythia->GetPrimaries();
-// fPythia->Pylist(1);
-
+ } else {
+
+ // EXODUS decayer
+ if(idpart == 111){
+ fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 0);
+ Lu1Ent(0, idpart, energy, theta, phi);
+ fPythia->PizeroDalitz();
+ fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 1);
+ }
+ else if(idpart == 221){
+ fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 0);
+ Lu1Ent(0, idpart, energy, theta, phi);
+ fPythia->EtaDalitz();
+ fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 1);
+ }
+ else if(idpart == 113){
+ Lu1Ent(0, idpart, energy, theta, phi);
+ fPythia->RhoDirect();
+ }
+ else if(idpart == 223){
+ fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
+ Lu1Ent(0, idpart, energy, theta, phi);
+ fPythia->OmegaDirect();
+ fPythia->OmegaDalitz();
+ fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
+ }
+ else if(idpart == 331){
+ fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 0);
+ Lu1Ent(0, idpart, energy, theta, phi);
+ fPythia->EtaprimeDalitz();
+ fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 1);
+ }
+ else if(idpart == 333){
+ fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 0);
+ Lu1Ent(0, idpart, energy, theta, phi);
+ fPythia->PhiDirect();
+ fPythia->PhiDalitz();
+ fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 1);
+ }
+ else if(idpart == 443){
+ Lu1Ent(0, idpart, energy, theta, phi);
+ fPythia->JPsiDirect();
+ }
+ }
+
+ fPythia->GetPrimaries();
}
Int_t AliDecayerPythia::ImportParticles(TClonesArray *particles)
{
// Import the decay products
//
- return fPythia->ImportParticles(particles, "All");
+ 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 = fPythia->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 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;
//
// select mode
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+/-
+ 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
+ break;
case kSemiMuonic:
ForceParticleDecay( 411,13,1); // D+/-
ForceParticleDecay( 421,13,1); // D0
ForceParticleDecay( 5122,13,1); // Lambda_b
ForceParticleDecay( 5132,13,1); // Xsi_b
ForceParticleDecay( 5232,13,1); // Sigma_b
- ForceParticleDecay( 5332,13,1); // Omega_b
- ForceParticleDecay( 24,13,1); // W+/- -> mu+/-
+ 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(100553,13,2);// Upsilon'
ForceParticleDecay(200553,13,2);// Upsilon''
break;
+ case kJpsiDiMuon:
+ ForceParticleDecay( 443,13,2); // J/Psi
+ 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
+ break;
case kSemiElectronic:
ForceParticleDecay( 411,11,1); // D+/-
ForceParticleDecay( 421,11,1); // D0
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(100553,11,2);// Upsilon'
ForceParticleDecay(200553,11,2);// Upsilon''
break;
+ case kPsiPrimeJpsiDiElectron:
+ products[0] = 443;
+ products[1] = 211;
+ mult[0] = 1;
+ mult[1] = 2;
+ ForceParticleDecay( 100443, products, mult, 2, 1);
+ ForceParticleDecay( 443,11,2);
+ break;
case kBJpsiDiMuon:
products[0] = 443;
ForceParticleDecay( 531,443,1); // B_s
ForceParticleDecay( 5122,443,1); // Lambda_b
break;
+ case kBJpsiUndecayed:
+ ForceParticleDecay( 511,443,1); // B0
+ ForceParticleDecay( 521,443,1); // B+/-
+ ForceParticleDecay( 531,443,1); // B_s
+ ForceParticleDecay( 5122,443,1); // Lambda_b
+ fPythia->SetMDCY(fPythia->Pycomp(443),1,0); // switch-off J/psi
+ break;
case kBPsiPrimeDiElectron:
ForceParticleDecay( 511,100443,1); // B0
ForceParticleDecay( 521,100443,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 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 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-
break;
case kOmega:
ForceOmega();
+ break;
+ case kLambda:
+ ForceLambda();
+ break;
case kAll:
break;
case kNoDecay:
fPythia->SetMSTJ(21,0);
break;
case kNoDecayHeavy:
+ case kNeutralPion:
+ break;
+ case kNoDecayBeauty:
+ SwitchOffBDecay();
+ break;
+ case kElectronEM:
+ ForceParticleDecay( 111,11,1); // pi^0
+ ForceParticleDecay( 221,11,1); // eta
+ ForceParticleDecay( 113,11,1); // rho
+ ForceParticleDecay( 223,11,1); // omega
+ ForceParticleDecay( 331,11,1); // etaprime
+ ForceParticleDecay( 333,11,1); // phi
+ ForceParticleDecay( 443,11,1); // jpsi
+ break;
+ case kDiElectronEM:
+ ForceParticleDecay( 111,11,2); // pi^0
+ ForceParticleDecay( 221,11,2); // eta
+ ForceParticleDecay( 113,11,2); // rho
+ ForceParticleDecay( 223,11,2); // omega
+ ForceParticleDecay( 331,11,2); // etaprime
+ ForceParticleDecay( 333,11,2); // phi
+ ForceParticleDecay( 443,11,2); // jpsi
break;
+ case kGammaEM:
+ ForceParticleDecay( 111,22,1); // pi^0
+ ForceParticleDecay( 221,22,1); // eta
+ ForceParticleDecay( 113,22,1); // rho
+ ForceParticleDecay( 223,22,1); // omega
+ ForceParticleDecay( 331,22,1); // etaprime
+ ForceParticleDecay( 333,22,1); // phi
+ ForceParticleDecay( 443,22,1); // jpsi
+ break;
+ case kBeautyUpgrade:
+ ForceBeautyUpgrade();
+ 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::ForceBeautyUpgrade()
+{
+ //
+ // Force dedicated decay channels of signals ineresting
+ // for the ITS upgrade (Lb, Lc, Xi_c, B)
+ //
+
+ // Lb: 50% of them in Lc pi+ and the rest with a Lc in final state
+ if(gRandom->Rndm()<0.50) {
+ const Int_t prod3[3]={4122,211,0};
+ Int_t mult3[3]={1,1,1};
+ ForceParticleDecay(5122,prod3,mult3,3,1);
+ } else {
+ ForceParticleDecay( 5122, 4122, 1);
+ }
+ // Xi_c
+ ForceParticleDecay( 4232, 3312, 1);
+ // B+ -> D0pi+
+ const Int_t prod[2]={421,211};
+ Int_t mult[2]={1,1};
+ ForceParticleDecay(521,prod,mult,2,1);
+ // B0 -> D*-pi+
+ const Int_t prod2[2]={413,211};
+ ForceParticleDecay(511,prod2,mult,2,1);
+ // force charm hadronic decays (D mesons and Lc)
+ ForceHadronicD(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};
- Int_t decayP[kNHadrons][3] =
- {
- {kKMinus, kPiPlus, kPiPlus},
- {kKMinus, kPiPlus, 0 },
- {-1 , -1 , -1 },
- {-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);
+ // for Lambda_c -> Lambda(1520) pi+ -> p K- pi+
+ Int_t iLambda1520 = 3124;
+ Int_t productsL[2] = {kProton, kKMinus}, multL[2] = {1, 1};
+ ForceParticleDecay(iLambda1520, productsL, multL, 2);
+ // for Lambda_c -> Lambda pi+
+ Int_t iLambda=3122;
+ //for Lambda_c -> antiK0 p
+ Int_t iK0bar=-311;
+
+
+ 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] =
+ {
+ {kPiPlus , iPhi , 0 , 0},
+ {kKMinus , kPiPlus, iRho0 , 0 },
+ {-1 , -1 , -1 , -1},
+ {-1 , -1 , -1 , -1},
+ {kProton , kKMinus, kPiPlus , 0}
+ };
+ // for Lambda_c -> Lambda_1520 pi+ -> p K- pi+, D0-> K*0 pi+ pi- -> K3pi
+ Int_t decayP4[kNHadrons][4] =
{
- 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
+ {iKstarbar0 , kKPlus , 0 , 0},
+ {iKstarbar0 , kPiPlus , kPiMinus, 0},
+ {-1 , -1 , -1 , -1},
+ {-1 , -1 , -1 , -1},
+ {iLambda1520 , kPiPlus , 0 , 0}
+ };
+ // for Lambda_c -> Lambda pi+
+ Int_t decayP5[kNHadrons][4] =
+ {
+ {-1 , -1 , -1 , -1},
+ {-1 , -1 , -1 , -1},
+ {-1 , -1 , -1 , -1},
+ {-1 , -1 , -1 , -1},
+ {iLambda , kPiPlus, 0 , 0}
+ };
+
+ // for Lambda_c -> K0bar p
+ Int_t decayP6[kNHadrons][4] =
+ {
+ {-1 , -1 , -1 , -1},
+ {-1 , -1 , -1 , -1},
+ {-1 , -1 , -1 , -1},
+ {-1 , -1 , -1 , -1},
+ {kProton , iK0bar, 0 , 0}
+ };
+
+ if(optUse4Bodies==0){
+ for(Int_t iDau=0;iDau<4;iDau++){
+ decayP2[1][iDau]=-1;
+ decayP3[1][iDau]=-1;
+ decayP4[1][iDau]=-1;
+ }
+ }
+
+ for (Int_t ihadron = 0; ihadron < kNHadrons; ihadron++)
+ {
+ Int_t kc = fPythia->Pycomp(hadron[ihadron]);
+ Int_t ifirst = fPythia->GetMDCY(kc,2);
+ Int_t ilast = ifirst + fPythia->GetMDCY(kc,3)-1;
+ Double_t norm = 0.;
+ for (channel=ifirst; channel<=ilast;channel++) norm+=fPythia->GetBRAT(channel);
+ if (norm < 1.-1.e-12 || norm > 1.+1.e-12) {
+ char pName[16];
+ fPythia->Pyname(hadron[ihadron],pName);
+ AliWarning(Form("Total branching ratio of %s (PDG code = %d) not equal to 1 (= %f)",pName,hadron[ihadron],norm));
+ }
+ fBraPart[kc] = norm;
+ fPythia->SetMDCY(kc,1,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->GetKFDP(channel,1) == decayP4[ihadron][0] &&
+ fPythia->GetKFDP(channel,2) == decayP4[ihadron][1] &&
+ fPythia->GetKFDP(channel,3) == decayP4[ihadron][2] &&
+ fPythia->GetKFDP(channel,4) == decayP4[ihadron][3] &&
+ fPythia->GetKFDP(channel,5) == 0
+ ) || (
+ fPythia->GetKFDP(channel,1) == decayP5[ihadron][0] &&
+ fPythia->GetKFDP(channel,2) == decayP5[ihadron][1] &&
+ fPythia->GetKFDP(channel,3) == decayP5[ihadron][2] &&
+ fPythia->GetKFDP(channel,4) == decayP5[ihadron][3] &&
+ fPythia->GetKFDP(channel,5) == 0
+ ) || (
+ fPythia->GetKFDP(channel,1) == decayP6[ihadron][0] &&
+ fPythia->GetKFDP(channel,2) == decayP6[ihadron][1] &&
+ fPythia->GetKFDP(channel,3) == decayP6[ihadron][2] &&
+ fPythia->GetKFDP(channel,4) == decayP6[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
+ if (norm > 0) fBraPart[kc] /= norm;
} // hadrons
}
+
+
void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t mult)
{
//
fPythia->SetMDCY(kc,1,1);
Int_t ifirst=fPythia->GetMDCY(kc,2);
Int_t ilast=ifirst+fPythia->GetMDCY(kc,3)-1;
- fBraPart[kc] = 1;
+ Double_t norm = 0.;
+ for (Int_t channel=ifirst; channel<=ilast;channel++){
+ norm+=fPythia->GetBRAT(channel);
+ }
+ if (norm < 1.-1.e-12 || norm > 1.+1.e-12) {
+ char pName[16];
+ fPythia->Pyname(particle,pName);
+ AliWarning(Form("Total branching ratio of %s (PDG code = %d) not equal to 1 (= %f)",pName,particle,norm));
+ }
+ fBraPart[kc] = norm;
//
// Loop over decay channels
for (Int_t channel=ifirst; channel<=ilast;channel++) {
fPythia->SetMDME(channel,1,1);
} else {
fPythia->SetMDME(channel,1,0);
- fBraPart[kc]-=fPythia->GetBRAT(channel);
+ fBraPart[kc] -= fPythia->GetBRAT(channel);
}
}
+ if (norm > 0.) fBraPart[kc] /= norm;
}
-void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t* products, Int_t* mult, Int_t npart)
+void AliDecayerPythia::ForceParticleDecay(Int_t particle, const Int_t* products, Int_t* mult, Int_t npart, Bool_t flag)
{
//
// Force decay of particle into products with multiplicity mult
fPythia->SetMDCY(kc,1,1);
Int_t ifirst=fPythia->GetMDCY(kc,2);
Int_t ilast=ifirst+fPythia->GetMDCY(kc,3)-1;
- fBraPart[kc] = 1;
+ Double_t norm = 0.;
+ for (Int_t channel=ifirst; channel<=ilast;channel++) norm+=fPythia->GetBRAT(channel);
+ if (norm < 1.-1.e-12 || norm > 1.+1.e-12) {
+ char pName[16];
+ fPythia->Pyname(particle,pName);
+ AliWarning(Form("Total branching ratio of %s (PDG code = %d) not equal to 1 (= %f)",pName,particle,norm));
+ }
+ fBraPart[kc] = norm;
//
// Loop over decay channels
for (Int_t channel = ifirst; channel <= ilast; channel++) {
for (Int_t i = 0; i < npart; i++) {
nprod += (CountProducts(channel, products[i]) >= mult[i]);
}
- if (nprod) {
+ if ((nprod && !flag) || ((nprod == npart) && flag)) {
fPythia->SetMDME(channel,1,1);
- } else {
+ } else { //
fPythia->SetMDME(channel,1,0);
- fBraPart[kc]-=fPythia->GetBRAT(channel);
+ fBraPart[kc] -= fPythia->GetBRAT(channel);
}
}
+ if (norm > 0.) fBraPart[kc] /= norm;
+
+
+
}
void AliDecayerPythia::DefineParticles()
} // decay channels
}
+void AliDecayerPythia::ForceLambda()
+{
+ // Force Lambda -> p pi-
+ Int_t kc=fPythia->Pycomp(3122);
+ 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) == kProton &&
+ fPythia->GetKFDP(channel,2) == kPiMinus &&
+ fPythia->GetKFDP(channel,3) == 0
+ )
+ {
+ fPythia->SetMDME(channel,1,1);
+ } else {
+ fPythia->SetMDME(channel,1,0);
+ } // selected channel ?
+ } // decay channels
+}
+
+void AliDecayerPythia::SwitchOffBDecay()
+{
+// Switch off B-decays
+ Int_t heavyB[]={511,521,531,5122,5132,5232,5332};
+ for(int i=0;i<4;i++)
+ {
+ fPythia->SetMDCY(fPythia->Pycomp(heavyB[i]),1,0);
+ }
+}
Float_t AliDecayerPythia::GetPartialBranchingRatio(Int_t kf)
{
}
-#ifdef never
-void AliDecayerPythia::Streamer(TBuffer &R__b)
-{
- // Stream an object of class AliDecayerPythia.
-
- if (R__b.IsReading()) {
- Version_t R__v = R__b.ReadVersion(); if (R__v) { }
- AliDecayer::Streamer(R__b);
- (AliPythia::Instance())->Streamer(R__b);
- R__b >> (Int_t&)fDecay;
- R__b.ReadStaticArray(fBraPart);
- } else {
- R__b.WriteVersion(AliDecayerPythia::IsA());
- AliDecayer::Streamer(R__b);
- R__b << fPythia;
- R__b << (Int_t)fDecay;
- R__b.WriteArray(fBraPart, 501);
- }
-}
-#endif
-
void AliDecayerPythia::Copy(TObject &) const
{
//
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff