X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PYTHIA6%2FAliDecayerPythia.cxx;h=1465fdf9865ed16cd5b5f016fe828eec06a0a412;hb=04f85d205e80e82527ddf666908b3e47678d39c4;hp=75bde8aaa3b551350d38ae900f96083e70ac242d;hpb=e345c61934f23672489b32609c87d970e8a5a67a;p=u%2Fmrichter%2FAliRoot.git diff --git a/PYTHIA6/AliDecayerPythia.cxx b/PYTHIA6/AliDecayerPythia.cxx index 75bde8aaa3b..1465fdf9865 100644 --- a/PYTHIA6/AliDecayerPythia.cxx +++ b/PYTHIA6/AliDecayerPythia.cxx @@ -24,9 +24,11 @@ #include "AliDecayerPythia.h" #include "AliPythia.h" +#include "AliLog.h" #include #include #include +#include ClassImp(AliDecayerPythia) @@ -56,7 +58,10 @@ Bool_t AliDecayerPythia::fgInit = kFALSE; AliDecayerPythia::AliDecayerPythia(): fPythia(AliPythia::Instance()), - fDecay(kAll) + fDecay(kAll), + fHeavyFlavour(kTRUE), + fLongLived(kFALSE), + fPatchOmegaDalitz(0) { // Default Constructor for (Int_t i=0; i< 501; i++) fBraPart[i]= 1.; @@ -66,14 +71,19 @@ AliDecayerPythia::AliDecayerPythia(): AliDecayerPythia::AliDecayerPythia(const AliDecayerPythia &decayer): AliDecayer(decayer), fPythia(0), - fDecay(kAll) + fDecay(kAll), + fHeavyFlavour(kTRUE), + fLongLived(kFALSE), + fPatchOmegaDalitz(0) { // Copy Constructor decayer.Copy(*this); + for (Int_t i = 0; i < 501; i++) fBraPart[i] = 0.; } void AliDecayerPythia::Init() { + // Initialisation // if (!fgInit) { @@ -100,15 +110,28 @@ void AliDecayerPythia::Init() } } //...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); + +/* + if (fDecay != kNeutralPion) { + fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0); + } else { + fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1); + } +*/ + + 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(); } @@ -120,27 +143,43 @@ void AliDecayerPythia::Decay(Int_t idpart, TLorentzVector* p) Float_t energy = p->Energy(); Float_t theta = p->Theta(); Float_t phi = p->Phi(); - - Lu1Ent(0, idpart, energy, theta, phi); + if (!fPatchOmegaDalitz) { + Lu1Ent(0, idpart, energy, theta, phi); + } else { + fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0); + Lu1Ent(0, idpart, energy, theta, phi); + fPythia->DalitzDecays(); + fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1); + fPythia->Pyexec(); + } fPythia->GetPrimaries(); -// fPythia->Pylist(1); - } 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;ipAt(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 @@ -166,8 +205,7 @@ void AliDecayerPythia::ForceDecay() 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 @@ -176,6 +214,25 @@ void AliDecayerPythia::ForceDecay() 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+/- @@ -212,6 +269,9 @@ void AliDecayerPythia::ForceDecay() 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+/- @@ -248,6 +308,14 @@ void AliDecayerPythia::ForceDecay() 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; @@ -282,6 +350,13 @@ void AliDecayerPythia::ForceDecay() 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+/- @@ -295,6 +370,10 @@ void AliDecayerPythia::ForceDecay() 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; @@ -314,24 +393,92 @@ void AliDecayerPythia::ForceDecay() 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 + 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 + 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 + 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) + // + + ForceParticleDecay( 5122, 4122, 1); + ForceParticleDecay( 4232, 3312, 1); + const Int_t prod[2]={421,211}; + Int_t mult[2]={1,1}; + ForceParticleDecay(521,prod,mult,2,1); + ForceHadronicD(1); +} + void AliDecayerPythia::Lu1Ent(Int_t flag, Int_t idpart, Double_t mom, Double_t theta, Double_t phi) { @@ -354,70 +501,169 @@ Int_t AliDecayerPythia::CountProducts(Int_t channel, Int_t particle) } -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); + // 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] = + { + {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] = { - 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 + {-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) { // @@ -427,7 +673,14 @@ void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t m 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++) { @@ -435,12 +688,13 @@ void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t product, Int_t m 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 @@ -449,7 +703,14 @@ void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t* products, Int_t 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++) { @@ -457,13 +718,14 @@ void AliDecayerPythia::ForceParticleDecay(Int_t particle, Int_t* products, Int_t 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() @@ -547,7 +809,37 @@ void AliDecayerPythia::ForceOmega() } // 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) { @@ -581,27 +873,6 @@ void AliDecayerPythia::ReadDecayTable() } -#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 { //