#include <G4HEKaonZeroShortInelastic.hh>
Inheritance diagram for G4HEKaonZeroShortInelastic:
Public Member Functions | |
G4HEKaonZeroShortInelastic () | |
~G4HEKaonZeroShortInelastic () | |
virtual void | ModelDescription (std::ostream &) const |
G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) |
G4int | GetNumberOfSecondaries () |
void | FirstIntInCasKaonZero (G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, const G4double atomicWeight) |
void | FirstIntInCasAntiKaonZero (G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle) |
Data Fields | |
G4int | vecLength |
Definition at line 53 of file G4HEKaonZeroShortInelastic.hh.
G4HEKaonZeroShortInelastic::G4HEKaonZeroShortInelastic | ( | ) | [inline] |
Definition at line 56 of file G4HEKaonZeroShortInelastic.hh.
References G4cout, G4endl, G4HEInelastic::MAXPART, G4HadronicInteraction::theMaxEnergy, G4HadronicInteraction::theMinEnergy, and G4HEInelastic::verboseLevel.
00056 : G4HEInelastic("G4HEKaonZeroShortInelastic") 00057 { 00058 theMinEnergy = 20*CLHEP::GeV; 00059 theMaxEnergy = 10*CLHEP::TeV; 00060 MAXPART = 2048; 00061 verboseLevel = 0; 00062 G4cout << "WARNING: model G4HEKaonZeroShortInelastic is being deprecated and will\n" 00063 << "disappear in Geant4 version 10.0" << G4endl; 00064 }
G4HEKaonZeroShortInelastic::~G4HEKaonZeroShortInelastic | ( | ) | [inline] |
G4HadFinalState * G4HEKaonZeroShortInelastic::ApplyYourself | ( | const G4HadProjectile & | aTrack, | |
G4Nucleus & | targetNucleus | |||
) | [virtual] |
Implements G4HadronicInteraction.
Definition at line 58 of file G4HEKaonZeroShortInelastic.cc.
References G4HEInelastic::AntiKaonZero, G4HEInelastic::ElasticScattering(), G4HEInelastic::FillParticleChange(), FirstIntInCasAntiKaonZero(), FirstIntInCasKaonZero(), G4cout, G4endl, G4UniformRand, G4Nucleus::GetA_asInt(), G4HEVector::getCode(), G4HEVector::getEnergy(), G4HEVector::getMass(), G4HEVector::getName(), G4Nucleus::GetZ_asInt(), G4HEInelastic::HighEnergyCascading(), G4HEInelastic::HighEnergyClusterProduction(), G4HEInelastic::KaonZero, G4HEInelastic::KaonZeroLong, G4HEInelastic::KaonZeroShort, G4HEInelastic::MAXPART, G4HEInelastic::MediumEnergyCascading(), G4HEInelastic::MediumEnergyClusterProduction(), G4HEInelastic::NuclearExcitation(), G4HEInelastic::NuclearInelasticity(), G4HEInelastic::QuasiElasticScattering(), G4HEVector::setDefinition(), G4HadFinalState::SetStatusChange(), stopAndKill, G4HEInelastic::StrangeParticlePairProduction(), G4HadronicInteraction::theParticleChange, vecLength, and G4HEInelastic::verboseLevel.
00060 { 00061 G4HEVector* pv = new G4HEVector[MAXPART]; 00062 const G4HadProjectile* aParticle = &aTrack; 00063 const G4double atomicWeight = targetNucleus.GetA_asInt(); 00064 const G4double atomicNumber = targetNucleus.GetZ_asInt(); 00065 G4HEVector incidentParticle(aParticle); 00066 00067 G4int incidentCode = incidentParticle.getCode(); 00068 G4double incidentMass = incidentParticle.getMass(); 00069 G4double incidentTotalEnergy = incidentParticle.getEnergy(); 00070 00071 // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 00072 // DHW 19 May 2011: variable set but not used 00073 00074 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass; 00075 00076 if(incidentKineticEnergy < 1) 00077 G4cout << "GHEKaonZeroShortInelastic: incident energy < 1 GeV" << G4endl; 00078 00079 if(verboseLevel > 1) { 00080 G4cout << "G4HEKaonZeroShortInelastic::ApplyYourself" << G4endl; 00081 G4cout << "incident particle " << incidentParticle.getName() 00082 << "mass " << incidentMass 00083 << "kinetic energy " << incidentKineticEnergy 00084 << G4endl; 00085 G4cout << "target material with (A,Z) = (" 00086 << atomicWeight << "," << atomicNumber << ")" << G4endl; 00087 } 00088 00089 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 00090 atomicWeight, atomicNumber); 00091 if(verboseLevel > 1) 00092 G4cout << "nuclear inelasticity = " << inelasticity << G4endl; 00093 00094 incidentKineticEnergy -= inelasticity; 00095 00096 G4double excitationEnergyGNP = 0.; 00097 G4double excitationEnergyDTA = 0.; 00098 00099 G4double excitation = NuclearExcitation(incidentKineticEnergy, 00100 atomicWeight, atomicNumber, 00101 excitationEnergyGNP, 00102 excitationEnergyDTA); 00103 if(verboseLevel > 1) 00104 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 00105 << excitationEnergyDTA << G4endl; 00106 00107 incidentKineticEnergy -= excitation; 00108 incidentTotalEnergy = incidentKineticEnergy + incidentMass; 00109 // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass) 00110 // *(incidentTotalEnergy+incidentMass)); 00111 // DHW 19 May 2011: variable set but not used 00112 00113 G4HEVector targetParticle; 00114 if(G4UniformRand() < atomicNumber/atomicWeight) { 00115 targetParticle.setDefinition("Proton"); 00116 } else { 00117 targetParticle.setDefinition("Neutron"); 00118 } 00119 00120 G4double targetMass = targetParticle.getMass(); 00121 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass 00122 + targetMass*targetMass 00123 + 2.0*targetMass*incidentTotalEnergy); 00124 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass; 00125 00126 G4bool inElastic = true; 00127 vecLength = 0; 00128 00129 if(verboseLevel > 1) 00130 G4cout << "ApplyYourself: CallFirstIntInCascade for particle " 00131 << incidentCode << G4endl; 00132 00133 G4bool successful = false; 00134 00135 // Split K0L into K0 and K0bar 00136 if (G4UniformRand() < 0.5) 00137 FirstIntInCasAntiKaonZero(inElastic, availableEnergy, pv, vecLength, 00138 incidentParticle, targetParticle ); 00139 else 00140 FirstIntInCasKaonZero(inElastic, availableEnergy, pv, vecLength, 00141 incidentParticle, targetParticle, atomicWeight ); 00142 00143 // Do nuclear interaction with either K0 or K0bar 00144 if ((vecLength > 0) && (availableEnergy > 1.)) 00145 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy, 00146 pv, vecLength, 00147 incidentParticle, targetParticle); 00148 00149 HighEnergyCascading(successful, pv, vecLength, 00150 excitationEnergyGNP, excitationEnergyDTA, 00151 incidentParticle, targetParticle, 00152 atomicWeight, atomicNumber); 00153 if (!successful) 00154 HighEnergyClusterProduction(successful, pv, vecLength, 00155 excitationEnergyGNP, excitationEnergyDTA, 00156 incidentParticle, targetParticle, 00157 atomicWeight, atomicNumber); 00158 if (!successful) 00159 MediumEnergyCascading(successful, pv, vecLength, 00160 excitationEnergyGNP, excitationEnergyDTA, 00161 incidentParticle, targetParticle, 00162 atomicWeight, atomicNumber); 00163 00164 if (!successful) 00165 MediumEnergyClusterProduction(successful, pv, vecLength, 00166 excitationEnergyGNP, excitationEnergyDTA, 00167 incidentParticle, targetParticle, 00168 atomicWeight, atomicNumber); 00169 if (!successful) 00170 QuasiElasticScattering(successful, pv, vecLength, 00171 excitationEnergyGNP, excitationEnergyDTA, 00172 incidentParticle, targetParticle, 00173 atomicWeight, atomicNumber); 00174 00175 if (!successful) 00176 ElasticScattering(successful, pv, vecLength, 00177 incidentParticle, 00178 atomicWeight, atomicNumber); 00179 00180 if (!successful) 00181 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 00182 << G4endl; 00183 00184 // Check for K0, K0bar and change particle types to K0L, K0S if necessary 00185 G4int kcode; 00186 for (G4int i = 0; i < vecLength; i++) { 00187 kcode = pv[i].getCode(); 00188 if (kcode == KaonZero.getCode() || kcode == AntiKaonZero.getCode()) { 00189 if (G4UniformRand() < 0.5) 00190 pv[i] = KaonZeroShort; 00191 else 00192 pv[i] = KaonZeroLong; 00193 } 00194 } 00195 00196 // ................ 00197 00198 FillParticleChange(pv, vecLength); 00199 delete [] pv; 00200 theParticleChange.SetStatusChange(stopAndKill); 00201 return &theParticleChange; 00202 }
void G4HEKaonZeroShortInelastic::FirstIntInCasAntiKaonZero | ( | G4bool & | inElastic, | |
const G4double | availableEnergy, | |||
G4HEVector | pv[], | |||
G4int & | vecLen, | |||
const G4HEVector & | incidentParticle, | |||
const G4HEVector & | targetParticle | |||
) |
Definition at line 528 of file G4HEKaonZeroShortInelastic.cc.
References G4HEInelastic::AntiKaonZero, G4cout, G4endl, G4UniformRand, G4HEVector::getCode(), G4HEVector::getMass(), G4HEVector::getName(), G4HEVector::getTotalMomentum(), G4HEInelastic::KaonMinus, G4HEInelastic::KaonZero, G4HEInelastic::Lambda, CLHEP::detail::n, G4HEInelastic::Neutron, G4INCL::Math::pi, G4HEInelastic::PionMinus, G4HEInelastic::PionPlus, G4HEInelastic::PionZero, G4HEInelastic::pmltpc(), G4HEInelastic::Proton, G4HEInelastic::SigmaMinus, G4HEInelastic::SigmaPlus, G4HEInelastic::SigmaZero, and G4HEInelastic::verboseLevel.
Referenced by ApplyYourself().
00542 { 00543 static const G4double expxu = 82.; // upper bound for arg. of exp 00544 static const G4double expxl = -expxu; // lower bound for arg. of exp 00545 00546 static const G4double protb = 0.7; 00547 static const G4double neutb = 0.7; 00548 static const G4double c = 1.25; 00549 00550 static const G4int numMul = 1200; 00551 static const G4int numSec = 60; 00552 00553 G4int neutronCode = Neutron.getCode(); 00554 G4int protonCode = Proton.getCode(); 00555 G4int kaonMinusCode = KaonMinus.getCode(); 00556 G4int kaonZeroCode = KaonZero.getCode(); 00557 G4int antiKaonZeroCode = AntiKaonZero.getCode(); 00558 00559 G4int targetCode = targetParticle.getCode(); 00560 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 00561 00562 static G4bool first = true; 00563 static G4double protmul[numMul], protnorm[numSec]; // proton constants 00564 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 00565 00566 // misc. local variables 00567 // npos = number of pi+, nneg = number of pi-, nzero = number of pi0 00568 00569 G4int i, counter, nt, npos, nneg, nzero; 00570 00571 if (first) { 00572 // compute normalization constants, this will only be done once 00573 first = false; 00574 for( i=0; i<numMul; i++ )protmul[i] = 0.0; 00575 for( i=0; i<numSec; i++ )protnorm[i] = 0.0; 00576 counter = -1; 00577 for(npos=0; npos<(numSec/3); npos++) { 00578 for(nneg=std::max(0,npos-2); nneg<=npos; nneg++) { 00579 for(nzero=0; nzero<numSec/3; nzero++) { 00580 if(++counter < numMul) { 00581 nt = npos+nneg+nzero; 00582 if( (nt>0) && (nt<=numSec) ) { 00583 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c) ; 00584 protnorm[nt-1] += protmul[counter]; 00585 } 00586 } 00587 } 00588 } 00589 } 00590 00591 for( i=0; i<numMul; i++ )neutmul[i] = 0.0; 00592 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0; 00593 counter = -1; 00594 for(npos=0; npos<numSec/3; npos++) { 00595 for(nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) { 00596 for(nzero=0; nzero<numSec/3; nzero++) { 00597 if(++counter < numMul) { 00598 nt = npos+nneg+nzero; 00599 if( (nt>0) && (nt<=numSec) ) { 00600 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c); 00601 neutnorm[nt-1] += neutmul[counter]; 00602 } 00603 } 00604 } 00605 } 00606 } 00607 00608 for(i=0; i<numSec; i++) { 00609 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i]; 00610 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i]; 00611 } 00612 } // end of initialization 00613 00614 // initialize the first two particles 00615 // the same as beam and target 00616 pv[0] = incidentParticle; 00617 pv[1] = targetParticle; 00618 vecLen = 2; 00619 00620 if (!inElastic || (availableEnergy <= PionPlus.getMass())) 00621 return; 00622 00623 // Inelastic scattering 00624 00625 npos = 0, nneg = 0, nzero = 0; 00626 G4double cech[] = { 1., 1., 1., 0.70, 0.60, 0.55, 0.35, 0.25, 0.18, 0.15}; 00627 G4int iplab = G4int( incidentTotalMomentum*5.); 00628 if( (iplab < 10) && (G4UniformRand() < cech[iplab]) ) { 00629 G4int ipl = std::min(19, G4int( incidentTotalMomentum*5.)); 00630 G4double cnk0[] = {0.17, 0.18, 0.17, 0.24, 0.26, 0.20, 0.22, 0.21, 0.34, 0.45, 00631 0.58, 0.55, 0.36, 0.29, 0.29, 0.32, 0.32, 0.33, 0.33, 0.33}; 00632 if(G4UniformRand() < cnk0[ipl]) { 00633 if(targetCode == protonCode) { 00634 return; 00635 } else { 00636 pv[0] = KaonMinus; 00637 pv[1] = Proton; 00638 return; 00639 } 00640 } 00641 00642 G4double ran = G4UniformRand(); 00643 if(targetCode == protonCode) { 00644 00645 // target is a proton 00646 if( ran < 0.25 ) { 00647 ; 00648 } else if (ran < 0.50) { 00649 pv[0] = PionPlus; 00650 pv[1] = SigmaZero; 00651 } else if (ran < 0.75) { 00652 ; 00653 } else { 00654 pv[0] = PionPlus; 00655 pv[1] = Lambda; 00656 } 00657 } else { 00658 00659 // target is a neutron 00660 if( ran < 0.25 ) { 00661 pv[0] = PionMinus; 00662 pv[1] = SigmaPlus; 00663 } else if (ran < 0.50) { 00664 pv[0] = PionZero; 00665 pv[1] = SigmaZero; 00666 } else if (ran < 0.75) { 00667 pv[0] = PionPlus; 00668 pv[1] = SigmaMinus; 00669 } else { 00670 pv[0] = PionZero; 00671 pv[1] = Lambda; 00672 } 00673 } 00674 return; 00675 00676 } else { 00677 // number of total particles vs. centre of mass Energy - 2*proton mass 00678 00679 G4double aleab = std::log(availableEnergy); 00680 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514 00681 + aleab*(0.117712+0.0136912*aleab))) - 2.0; 00682 00683 // Normalization constant for kno-distribution. 00684 // Calculate first the sum of all constants, check for numerical problems. 00685 G4double test, dum, anpn = 0.0; 00686 00687 for (nt=1; nt<=numSec; nt++) { 00688 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00689 dum = pi*nt/(2.0*n*n); 00690 if (std::fabs(dum) < 1.0) { 00691 if( test >= 1.0e-10 )anpn += dum*test; 00692 } else { 00693 anpn += dum*test; 00694 } 00695 } 00696 00697 G4double ran = G4UniformRand(); 00698 G4double excs = 0.0; 00699 if (targetCode == protonCode) { 00700 counter = -1; 00701 for (npos=0; npos<numSec/3; npos++) { 00702 for (nneg=std::max(0,npos-2); nneg<=npos; nneg++) { 00703 for (nzero=0; nzero<numSec/3; nzero++) { 00704 if (++counter < numMul) { 00705 nt = npos+nneg+nzero; 00706 if( (nt>0) && (nt<=numSec) ) { 00707 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00708 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n); 00709 00710 if (std::fabs(dum) < 1.0) { 00711 if( test >= 1.0e-10 )excs += dum*test; 00712 } else { 00713 excs += dum*test; 00714 } 00715 00716 if (ran < excs) goto outOfLoop; //-----------------------> 00717 } 00718 } 00719 } 00720 } 00721 } 00722 // 3 previous loops continued to the end 00723 inElastic = false; // quasi-elastic scattering 00724 return; 00725 00726 } else { // target must be a neutron 00727 counter = -1; 00728 for (npos=0; npos<numSec/3; npos++) { 00729 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) { 00730 for (nzero=0; nzero<numSec/3; nzero++) { 00731 if (++counter < numMul) { 00732 nt = npos+nneg+nzero; 00733 if( (nt>=1) && (nt<=numSec) ) { 00734 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00735 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n); 00736 00737 if (std::fabs(dum) < 1.0) { 00738 if( test >= 1.0e-10 )excs += dum*test; 00739 } else { 00740 excs += dum*test; 00741 } 00742 00743 if (ran < excs) goto outOfLoop; // --------------------------> 00744 } 00745 } 00746 } 00747 } 00748 } 00749 // 3 previous loops continued to the end 00750 inElastic = false; // quasi-elastic scattering. 00751 return; 00752 } 00753 } 00754 outOfLoop: // <------------------------------------------------------------------------ 00755 00756 if( targetCode == protonCode) 00757 { 00758 if( npos == nneg) 00759 { 00760 } 00761 else if (npos == (1+nneg)) 00762 { 00763 if( G4UniformRand() < 0.5) 00764 { 00765 pv[0] = KaonMinus; 00766 } 00767 else 00768 { 00769 pv[1] = Neutron; 00770 } 00771 } 00772 else 00773 { 00774 pv[0] = KaonMinus; 00775 pv[1] = Neutron; 00776 } 00777 } 00778 else 00779 { 00780 if( npos == nneg) 00781 { 00782 if( G4UniformRand() < 0.75) 00783 { 00784 } 00785 else 00786 { 00787 pv[0] = KaonMinus; 00788 pv[1] = Proton; 00789 } 00790 } 00791 else if ( npos == (1+nneg)) 00792 { 00793 pv[0] = KaonMinus; 00794 } 00795 else 00796 { 00797 pv[1] = Proton; 00798 } 00799 } 00800 00801 00802 if( G4UniformRand() < 0.5 ) 00803 { 00804 if( ( (pv[0].getCode() == kaonMinusCode) 00805 && (pv[1].getCode() == neutronCode) ) 00806 || ( (pv[0].getCode() == kaonZeroCode) 00807 && (pv[1].getCode() == protonCode) ) 00808 || ( (pv[0].getCode() == antiKaonZeroCode) 00809 && (pv[1].getCode() == protonCode) ) ) 00810 { 00811 G4double ran = G4UniformRand(); 00812 if( pv[1].getCode() == protonCode) 00813 { 00814 if(ran < 0.68) 00815 { 00816 pv[0] = PionPlus; 00817 pv[1] = Lambda; 00818 } 00819 else if (ran < 0.84) 00820 { 00821 pv[0] = PionZero; 00822 pv[1] = SigmaPlus; 00823 } 00824 else 00825 { 00826 pv[0] = PionPlus; 00827 pv[1] = SigmaZero; 00828 } 00829 } 00830 else 00831 { 00832 if(ran < 0.68) 00833 { 00834 pv[0] = PionMinus; 00835 pv[1] = Lambda; 00836 } 00837 else if (ran < 0.84) 00838 { 00839 pv[0] = PionMinus; 00840 pv[1] = SigmaZero; 00841 } 00842 else 00843 { 00844 pv[0] = PionZero; 00845 pv[1] = SigmaMinus; 00846 } 00847 } 00848 } 00849 else 00850 { 00851 G4double ran = G4UniformRand(); 00852 if (ran < 0.67) 00853 { 00854 pv[0] = PionZero; 00855 pv[1] = Lambda; 00856 } 00857 else if (ran < 0.78) 00858 { 00859 pv[0] = PionMinus; 00860 pv[1] = SigmaPlus; 00861 } 00862 else if (ran < 0.89) 00863 { 00864 pv[0] = PionZero; 00865 pv[1] = SigmaZero; 00866 } 00867 else 00868 { 00869 pv[0] = PionPlus; 00870 pv[1] = SigmaMinus; 00871 } 00872 } 00873 } 00874 00875 nt = npos + nneg + nzero; 00876 while ( nt > 0) { 00877 G4double ran = G4UniformRand(); 00878 if ( ran < (G4double)npos/nt) { 00879 if( npos > 0 ) { 00880 pv[vecLen++] = PionPlus; 00881 npos--; 00882 } 00883 } else if (ran < (G4double)(npos+nneg)/nt) { 00884 if( nneg > 0 ) { 00885 pv[vecLen++] = PionMinus; 00886 nneg--; 00887 } 00888 } else { 00889 if( nzero > 0 ) { 00890 pv[vecLen++] = PionZero; 00891 nzero--; 00892 } 00893 } 00894 nt = npos + nneg + nzero; 00895 } 00896 00897 if (verboseLevel > 1) { 00898 G4cout << "Particles produced: " ; 00899 G4cout << pv[0].getName() << " " ; 00900 G4cout << pv[1].getName() << " " ; 00901 for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ; 00902 G4cout << G4endl; 00903 } 00904 00905 return; 00906 }
void G4HEKaonZeroShortInelastic::FirstIntInCasKaonZero | ( | G4bool & | inElastic, | |
const G4double | availableEnergy, | |||
G4HEVector | pv[], | |||
G4int & | vecLen, | |||
const G4HEVector & | incidentParticle, | |||
const G4HEVector & | targetParticle, | |||
const G4double | atomicWeight | |||
) |
Definition at line 206 of file G4HEKaonZeroShortInelastic.cc.
References G4cout, G4endl, G4UniformRand, G4HEVector::getCode(), G4HEVector::getMass(), G4HEVector::getName(), G4HEVector::getTotalMomentum(), G4HEInelastic::KaonPlus, CLHEP::detail::n, G4HEInelastic::Neutron, neutronCode, G4INCL::Math::pi, G4HEInelastic::PionMinus, G4HEInelastic::PionPlus, G4HEInelastic::PionZero, G4HEInelastic::pmltpc(), G4HEInelastic::Proton, sqr(), and G4HEInelastic::verboseLevel.
Referenced by ApplyYourself().
00221 { 00222 static const G4double expxu = 82.; // upper bound for arg. of exp 00223 static const G4double expxl = -expxu; // lower bound for arg. of exp 00224 00225 static const G4double protb = 0.7; 00226 static const G4double neutb = 0.7; 00227 static const G4double c = 1.25; 00228 00229 static const G4int numMul = 1200; 00230 static const G4int numSec = 60; 00231 00232 G4int neutronCode = Neutron.getCode(); 00233 G4int protonCode = Proton.getCode(); 00234 00235 G4int targetCode = targetParticle.getCode(); 00236 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 00237 00238 static G4bool first = true; 00239 static G4double protmul[numMul], protnorm[numSec]; // proton constants 00240 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 00241 00242 // misc. local variables 00243 // npos = number of pi+, nneg = number of pi-, nzero = number of pi0 00244 00245 G4int i, counter, nt, npos, nneg, nzero; 00246 00247 if (first) { 00248 // compute normalization constants, this will only be done once 00249 first = false; 00250 for( i=0; i<numMul; i++ )protmul[i] = 0.0; 00251 for( i=0; i<numSec; i++ )protnorm[i] = 0.0; 00252 counter = -1; 00253 for (npos=0; npos<(numSec/3); npos++) { 00254 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) { 00255 for (nzero=0; nzero<numSec/3; nzero++) { 00256 if (++counter < numMul) { 00257 nt = npos+nneg+nzero; 00258 if( (nt>0) && (nt<=numSec) ) { 00259 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c) ; 00260 protnorm[nt-1] += protmul[counter]; 00261 } 00262 } 00263 } 00264 } 00265 } 00266 00267 for( i=0; i<numMul; i++ )neutmul[i] = 0.0; 00268 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0; 00269 counter = -1; 00270 for (npos=0; npos<numSec/3; npos++) { 00271 for (nneg=npos; nneg<=(npos+2); nneg++) { 00272 for (nzero=0; nzero<numSec/3; nzero++) { 00273 if (++counter < numMul) { 00274 nt = npos+nneg+nzero; 00275 if( (nt>0) && (nt<=numSec) ) { 00276 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c); 00277 neutnorm[nt-1] += neutmul[counter]; 00278 } 00279 } 00280 } 00281 } 00282 } 00283 00284 for (i=0; i<numSec; i++) { 00285 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i]; 00286 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i]; 00287 } 00288 } // end of initialization 00289 00290 00291 // Initialize the first two particles 00292 // the same as beam and target 00293 pv[0] = incidentParticle; 00294 pv[1] = targetParticle; 00295 vecLen = 2; 00296 00297 if( !inElastic ) { 00298 // quasi-elastic scattering, no pions produced 00299 if( targetCode == protonCode) { 00300 G4double cech[] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.10,0.09,0.07}; 00301 G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*5. ) ); 00302 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42)) { 00303 // charge exchange K+ n -> K0 p 00304 pv[0] = KaonPlus; 00305 pv[1] = Neutron; 00306 } 00307 } 00308 return; 00309 } else if (availableEnergy <= PionPlus.getMass()) { 00310 return; 00311 } 00312 00313 // Inelastic scattering 00314 00315 npos = 0, nneg = 0, nzero = 0; 00316 G4double eab = availableEnergy; 00317 G4int ieab = G4int( eab*5.0 ); 00318 00319 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98}; 00320 if( (ieab <= 9) && (G4UniformRand() >= supp[ieab])) { 00321 // Suppress high multiplicity events at low momentum 00322 // only one additional pion will be produced 00323 G4double w0, wp, wm, wt, ran; 00324 if (targetCode == neutronCode) { 00325 // target is a neutron 00326 w0 = - sqr(1.+protb)/(2.*c*c); 00327 w0 = std::exp(w0); 00328 wm = - sqr(-1.+protb)/(2.*c*c); 00329 wm = std::exp(wm); 00330 w0 = w0/2.; 00331 wm = wm*1.5; 00332 if (G4UniformRand() < w0/(w0+wm) ) { 00333 npos = 0; 00334 nneg = 0; 00335 nzero = 1; 00336 } else { 00337 npos = 0; 00338 nneg = 1; 00339 nzero = 0; 00340 } 00341 00342 } else { 00343 // target is a proton 00344 w0 = -sqr(1.+neutb)/(2.*c*c); 00345 wp = w0 = std::exp(w0); 00346 wm = -sqr(-1.+neutb)/(2.*c*c); 00347 wm = std::exp(wm); 00348 wt = w0+wp+wm; 00349 wp = w0+wp; 00350 ran = G4UniformRand(); 00351 if ( ran < w0/wt) { 00352 npos = 0; 00353 nneg = 0; 00354 nzero = 1; 00355 } else if (ran < wp/wt) { 00356 npos = 1; 00357 nneg = 0; 00358 nzero = 0; 00359 } else { 00360 npos = 0; 00361 nneg = 1; 00362 nzero = 0; 00363 } 00364 } 00365 } else { 00366 // number of total particles vs. centre of mass Energy - 2*proton mass 00367 00368 G4double aleab = std::log(availableEnergy); 00369 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514 00370 + aleab*(0.117712+0.0136912*aleab))) - 2.0; 00371 00372 // Normalization constant for kno-distribution. 00373 // Calculate first the sum of all constants, check for numerical problems. 00374 G4double test, dum, anpn = 0.0; 00375 00376 for (nt=1; nt<=numSec; nt++) { 00377 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00378 dum = pi*nt/(2.0*n*n); 00379 if (std::fabs(dum) < 1.0) { 00380 if( test >= 1.0e-10 )anpn += dum*test; 00381 } else { 00382 anpn += dum*test; 00383 } 00384 } 00385 00386 G4double ran = G4UniformRand(); 00387 G4double excs = 0.0; 00388 if( targetCode == protonCode ) 00389 { 00390 counter = -1; 00391 for( npos=0; npos<numSec/3; npos++ ) 00392 { 00393 for( nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++ ) 00394 { 00395 for (nzero=0; nzero<numSec/3; nzero++) { 00396 if (++counter < numMul) { 00397 nt = npos+nneg+nzero; 00398 if ( (nt>0) && (nt<=numSec) ) { 00399 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00400 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n); 00401 if (std::fabs(dum) < 1.0) { 00402 if( test >= 1.0e-10 )excs += dum*test; 00403 } else { 00404 excs += dum*test; 00405 } 00406 if (ran < excs) goto outOfLoop; //-----------------------> 00407 } 00408 } 00409 } 00410 } 00411 } 00412 00413 // 3 previous loops continued to the end 00414 inElastic = false; // quasi-elastic scattering 00415 return; 00416 } 00417 else 00418 { // target must be a neutron 00419 counter = -1; 00420 for( npos=0; npos<numSec/3; npos++ ) 00421 { 00422 for( nneg=npos; nneg<=(npos+2); nneg++ ) 00423 { 00424 for (nzero=0; nzero<numSec/3; nzero++) { 00425 if (++counter < numMul) { 00426 nt = npos+nneg+nzero; 00427 if ( (nt>=1) && (nt<=numSec) ) { 00428 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) ); 00429 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n); 00430 if (std::fabs(dum) < 1.0) { 00431 if( test >= 1.0e-10 )excs += dum*test; 00432 } else { 00433 excs += dum*test; 00434 } 00435 if (ran < excs) goto outOfLoop; // --------------------------> 00436 } 00437 } 00438 } 00439 } 00440 } 00441 // 3 previous loops continued to the end 00442 inElastic = false; // quasi-elastic scattering. 00443 return; 00444 } 00445 } 00446 outOfLoop: // <----------------------------------------------- 00447 00448 if( targetCode == neutronCode) 00449 { 00450 if( npos == nneg) 00451 { 00452 } 00453 else if (npos == (nneg-1)) 00454 { 00455 if( G4UniformRand() < 0.5) 00456 { 00457 pv[0] = KaonPlus; 00458 } 00459 else 00460 { 00461 pv[1] = Proton; 00462 } 00463 } 00464 else 00465 { 00466 pv[0] = KaonPlus; 00467 pv[1] = Proton; 00468 } 00469 } 00470 else 00471 { 00472 if( npos == nneg ) 00473 { 00474 if( G4UniformRand() < 0.25) 00475 { 00476 pv[0] = KaonPlus; 00477 pv[1] = Neutron; 00478 } 00479 else 00480 { 00481 } 00482 } 00483 else if ( npos == (nneg+1)) 00484 { 00485 pv[1] = Neutron; 00486 } 00487 else 00488 { 00489 pv[0] = KaonPlus; 00490 } 00491 } 00492 00493 nt = npos + nneg + nzero; 00494 while (nt > 0) { 00495 G4double ran = G4UniformRand(); 00496 if (ran < (G4double)npos/nt) { 00497 if (npos > 0) { 00498 pv[vecLen++] = PionPlus; 00499 npos--; 00500 } 00501 } else if ( ran < (G4double)(npos+nneg)/nt) { 00502 if (nneg > 0) { 00503 pv[vecLen++] = PionMinus; 00504 nneg--; 00505 } 00506 } else { 00507 if (nzero > 0) { 00508 pv[vecLen++] = PionZero; 00509 nzero--; 00510 } 00511 } 00512 nt = npos + nneg + nzero; 00513 } 00514 00515 if (verboseLevel > 1) { 00516 G4cout << "Particles produced: " ; 00517 G4cout << pv[0].getName() << " " ; 00518 G4cout << pv[1].getName() << " " ; 00519 for (i=2; i < vecLen; i++) G4cout << pv[i].getName() << " " ; 00520 G4cout << G4endl; 00521 } 00522 00523 return; 00524 }
G4int G4HEKaonZeroShortInelastic::GetNumberOfSecondaries | ( | ) | [inline] |
Definition at line 75 of file G4HEKaonZeroShortInelastic.hh.
References vecLength.
00076 { return vecLength;};
void G4HEKaonZeroShortInelastic::ModelDescription | ( | std::ostream & | ) | const [virtual] |
Reimplemented from G4HadronicInteraction.
Definition at line 43 of file G4HEKaonZeroShortInelastic.cc.
00044 { 00045 outFile << "G4HEKaonZeroShortInelastic is one of the High Energy\n" 00046 << "Parameterized (HEP) models used to implement inelastic\n" 00047 << "K0S scattering from nuclei. It is a re-engineered\n" 00048 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n" 00049 << "initial collision products into backward- and forward-going\n" 00050 << "clusters which are then decayed into final state hadrons.\n" 00051 << "The model does not conserve energy on an event-by-event\n" 00052 << "basis. It may be applied to K0S with initial energies\n" 00053 << "above 20 GeV.\n"; 00054 }
Definition at line 70 of file G4HEKaonZeroShortInelastic.hh.
Referenced by ApplyYourself(), and GetNumberOfSecondaries().