X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PYTHIA6%2FAliDecayerPythia.cxx;h=7579cfecfecc7da0816197d96084c79398cf3e55;hb=fb44eb8e13f9490abf4660b881b154cb919ac30a;hp=36a24b1b08f4c5a1b718982197a7d5247d31756d;hpb=2f29952613c5764368dda5a624df48a3171aa797;p=u%2Fmrichter%2FAliRoot.git diff --git a/PYTHIA6/AliDecayerPythia.cxx b/PYTHIA6/AliDecayerPythia.cxx index 36a24b1b08f..7579cfecfec 100644 --- a/PYTHIA6/AliDecayerPythia.cxx +++ b/PYTHIA6/AliDecayerPythia.cxx @@ -54,14 +54,26 @@ extern "C" void type_of_call Bool_t AliDecayerPythia::fgInit = kFALSE; -AliDecayerPythia::AliDecayerPythia() +AliDecayerPythia::AliDecayerPythia(): + fPythia(AliPythia::Instance()), + fDecay(kAll), + fHeavyFlavour(kTRUE) { // 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) +{ + // Copy Constructor + decayer.Copy(*this); +} + void AliDecayerPythia::Init() { // Initialisation @@ -89,7 +101,25 @@ void AliDecayerPythia::Init() } } } +//...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); + } + + + 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(); } @@ -103,8 +133,6 @@ void AliDecayerPythia::Decay(Int_t idpart, TLorentzVector* p) Lu1Ent(0, idpart, energy, theta, phi); fPythia->GetPrimaries(); -// fPythia->Pylist(1); - } Int_t AliDecayerPythia::ImportParticles(TClonesArray *particles) @@ -118,6 +146,9 @@ Int_t AliDecayerPythia::ImportParticles(TClonesArray *particles) 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; @@ -146,8 +177,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 @@ -156,6 +186,34 @@ 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+/- + 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 @@ -183,6 +241,15 @@ void AliDecayerPythia::ForceDecay() ForceParticleDecay(100553,13,2);// Upsilon' ForceParticleDecay(200553,13,2);// Upsilon'' 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 @@ -257,6 +324,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; @@ -276,8 +347,14 @@ 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- @@ -290,10 +367,23 @@ void AliDecayerPythia::ForceDecay() fPythia->SetMSTJ(21,0); break; case kNoDecayHeavy: + case kNeutralPion: 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::Lu1Ent(Int_t flag, Int_t idpart, Double_t mom, Double_t theta, Double_t phi) { @@ -316,47 +406,103 @@ 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}; - 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); + + + 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] = + { + {-1 , -1 , -1 , -1}, + {kKMinus , kPiPlus, iRho0 , 0 }, + {-1 , -1 , -1 , -1}, + {-1 , -1 , -1 , -1}, + {kProton , kKMinus, kPiPlus , 0} + }; + if(optUse4Bodies==0){ + for(Int_t iDau=0;iDau<4;iDau++){ + decayP2[1][iDau]=-1; + decayP3[1][iDau]=-1; + } + } + + 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 + 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) == 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->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) { //