- Timestamp:
- Dec 22, 2010, 3:52:27 PM (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/high_energy/src/G4HEPionMinusInelastic.cc
r1340 r1347 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4HEPionMinusInelastic.cc,v 1.17 2010/11/29 05:44:44 dennis Exp $ 27 // GEANT4 tag $Name: geant4-09-04-ref-00 $ 26 28 // 27 // $Id: G4HEPionMinusInelastic.cc,v 1.15 2008/03/17 20:49:17 dennis Exp $28 // GEANT4 tag $Name: geant4-09-03-ref-09 $29 30 29 // 11-OCT-2007 F.W. Jones: fixed incorrect Imax (should be Imin) in 31 30 // sampling of charge exchange. … … 35 34 #include "G4ios.hh" 36 35 37 //38 36 // G4 Process: Gheisha High Energy Collision model. 39 37 // This includes the high energy cascading model, the two-body-resonance model 40 // and the low energy two-body model. Not included are the low energy stuff like41 // nuclear reactions, nuclear fission without any cascading and all processes for42 // p articles at rest.38 // and the low energy two-body model. Not included are the low energy stuff 39 // like nuclear reactions, nuclear fission without any cascading and all 40 // processes for particles at rest. 43 41 // First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96. 44 42 // H. Fesefeldt, RWTH-Aachen, 23-October-1996 … … 47 45 #include "G4HEPionMinusInelastic.hh" 48 46 49 G4HadFinalState * G4HEPionMinusInelastic::50 ApplyYourself( const G4HadProjectile &aTrack, G4Nucleus &targetNucleus ) 51 {52 G4HEVector * pv = new G4HEVector[MAXPART]; 53 const G4HadProjectile *aParticle = &aTrack;54 // G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();55 56 57 47 G4HadFinalState* 48 G4HEPionMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack, 49 G4Nucleus& targetNucleus) 50 { 51 G4HEVector* pv = new G4HEVector[MAXPART]; 52 const G4HadProjectile* aParticle = &aTrack; 53 const G4double A = targetNucleus.GetN(); 54 const G4double Z = targetNucleus.GetZ(); 55 G4HEVector incidentParticle(aParticle); 58 56 59 G4double atomicNumber = Z; 60 G4double atomicWeight = A; 61 62 G4int incidentCode = incidentParticle.getCode(); 63 G4double incidentMass = incidentParticle.getMass(); 64 G4double incidentTotalEnergy = incidentParticle.getEnergy(); 65 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 66 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass; 67 68 if(incidentKineticEnergy < 1.) 69 { 70 G4cout << "GHEPionMinusInelastic: incident energy < 1 GeV" << G4endl ; 71 } 72 if(verboseLevel > 1) 73 { 74 G4cout << "G4HEPionMinusInelastic::ApplyYourself" << G4endl; 75 G4cout << "incident particle " << incidentParticle.getName() 76 << "mass " << incidentMass 77 << "kinetic energy " << incidentKineticEnergy 78 << G4endl; 79 G4cout << "target material with (A,Z) = (" 80 << atomicWeight << "," << atomicNumber << ")" << G4endl; 81 } 82 83 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 84 atomicWeight, atomicNumber); 85 if(verboseLevel > 1) 86 G4cout << "nuclear inelasticity = " << inelasticity << G4endl; 87 88 incidentKineticEnergy -= inelasticity; 57 G4double atomicNumber = Z; 58 G4double atomicWeight = A; 59 60 G4int incidentCode = incidentParticle.getCode(); 61 G4double incidentMass = incidentParticle.getMass(); 62 G4double incidentTotalEnergy = incidentParticle.getEnergy(); 63 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 64 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass; 65 66 if (incidentKineticEnergy < 1.) 67 G4cout << "GHEPionMinusInelastic: incident energy < 1 GeV" << G4endl; 68 69 if (verboseLevel > 1) { 70 G4cout << "G4HEPionMinusInelastic::ApplyYourself" << G4endl; 71 G4cout << "incident particle " << incidentParticle.getName() 72 << "mass " << incidentMass 73 << "kinetic energy " << incidentKineticEnergy 74 << G4endl; 75 G4cout << "target material with (A,Z) = (" 76 << atomicWeight << "," << atomicNumber << ")" << G4endl; 77 } 78 79 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 80 atomicWeight, atomicNumber); 81 if (verboseLevel > 1) 82 G4cout << "nuclear inelasticity = " << inelasticity << G4endl; 83 84 incidentKineticEnergy -= inelasticity; 89 85 90 G4double excitationEnergyGNP = 0.; 91 G4double excitationEnergyDTA = 0.; 92 93 G4double excitation = NuclearExcitation(incidentKineticEnergy, 94 atomicWeight, atomicNumber, 95 excitationEnergyGNP, 96 excitationEnergyDTA); 97 if(verboseLevel > 1) 98 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 99 << excitationEnergyDTA << G4endl; 100 101 102 incidentKineticEnergy -= excitation; 103 incidentTotalEnergy = incidentKineticEnergy + incidentMass; 104 incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass) 105 *(incidentTotalEnergy+incidentMass)); 106 107 G4HEVector targetParticle; 86 G4double excitationEnergyGNP = 0.; 87 G4double excitationEnergyDTA = 0.; 88 89 G4double excitation = NuclearExcitation(incidentKineticEnergy, 90 atomicWeight, atomicNumber, 91 excitationEnergyGNP, 92 excitationEnergyDTA); 93 if (verboseLevel > 1) 94 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 95 << excitationEnergyDTA << G4endl; 96 97 incidentKineticEnergy -= excitation; 98 incidentTotalEnergy = incidentKineticEnergy + incidentMass; 99 incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass) 100 *(incidentTotalEnergy+incidentMass)); 101 102 G4HEVector targetParticle; 108 103 109 if(G4UniformRand() < atomicNumber/atomicWeight) 110 { 111 targetParticle.setDefinition("Proton"); 112 } 113 else 114 { 115 targetParticle.setDefinition("Neutron"); 116 } 117 118 G4double targetMass = targetParticle.getMass(); 119 G4double centerOfMassEnergy = std::sqrt( incidentMass*incidentMass + targetMass*targetMass 120 + 2.0*targetMass*incidentTotalEnergy); 121 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass; 122 123 // this was the meaning of inElastic in the 124 // original Gheisha stand-alone version. 125 // G4bool inElastic = InElasticCrossSectionInFirstInt 126 // (availableEnergy, incidentCode, incidentTotalMomentum); 127 // by unknown reasons, it has been replaced 128 // to the following code in Geant??? 129 G4bool inElastic = true; 130 // if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false; 131 132 vecLength = 0; 104 if (G4UniformRand() < atomicNumber/atomicWeight) { 105 targetParticle.setDefinition("Proton"); 106 } else { 107 targetParticle.setDefinition("Neutron"); 108 } 109 110 G4double targetMass = targetParticle.getMass(); 111 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass 112 + targetMass*targetMass 113 + 2.0*targetMass*incidentTotalEnergy); 114 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass; 115 116 // The value of the inElastic flag was originally defined in the Gheisha 117 // stand-alone code as follows: 118 // G4bool inElastic = InElasticCrossSectionInFirstInt 119 // (availableEnergy, incidentCode, incidentTotalMomentum); 120 // For unknown reasons, it was replaced by the following code in Geant 121 122 G4bool inElastic = true; 123 // if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false; 124 125 vecLength = 0; 133 126 134 if(verboseLevel > 1)135 127 if (verboseLevel > 1) 128 G4cout << "ApplyYourself: CallFirstIntInCascade for particle " 136 129 << incidentCode << G4endl; 137 130 138 131 G4bool successful = false; 139 132 140 if(inElastic || (!inElastic && atomicWeight < 1.5)) 141 { 142 FirstIntInCasPionMinus(inElastic, availableEnergy, pv, vecLength, 143 incidentParticle, targetParticle); 144 145 if(verboseLevel > 1) 146 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl; 147 148 149 if ((vecLength > 0) && (availableEnergy > 1.)) 150 StrangeParticlePairProduction( availableEnergy, centerOfMassEnergy, 151 pv, vecLength, 152 incidentParticle, targetParticle); 153 HighEnergyCascading( successful, pv, vecLength, 154 excitationEnergyGNP, excitationEnergyDTA, 155 incidentParticle, targetParticle, 156 atomicWeight, atomicNumber); 157 if (!successful) 158 HighEnergyClusterProduction( successful, pv, vecLength, 159 excitationEnergyGNP, excitationEnergyDTA, 160 incidentParticle, targetParticle, 161 atomicWeight, atomicNumber); 162 if (!successful) 163 MediumEnergyCascading( successful, pv, vecLength, 164 excitationEnergyGNP, excitationEnergyDTA, 165 incidentParticle, targetParticle, 166 atomicWeight, atomicNumber); 167 168 if (!successful) 169 MediumEnergyClusterProduction( successful, pv, vecLength, 170 excitationEnergyGNP, excitationEnergyDTA, 171 incidentParticle, targetParticle, 172 atomicWeight, atomicNumber); 173 if (!successful) 174 QuasiElasticScattering( successful, pv, vecLength, 175 excitationEnergyGNP, excitationEnergyDTA, 176 incidentParticle, targetParticle, 177 atomicWeight, atomicNumber); 178 } 179 if (!successful) 180 { 181 ElasticScattering( successful, pv, vecLength, 182 incidentParticle, 183 atomicWeight, atomicNumber); 184 } 185 186 if (!successful) 187 { 188 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" << G4endl; 189 } 190 FillParticleChange(pv, vecLength); 191 delete [] pv; 192 theParticleChange.SetStatusChange(stopAndKill); 193 return & theParticleChange; 194 } 133 FirstIntInCasPionMinus(inElastic, availableEnergy, pv, vecLength, 134 incidentParticle, targetParticle); 135 136 if (verboseLevel > 1) 137 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl; 138 139 if ((vecLength > 0) && (availableEnergy > 1.)) 140 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy, 141 pv, vecLength, 142 incidentParticle, targetParticle); 143 144 HighEnergyCascading(successful, pv, vecLength, 145 excitationEnergyGNP, excitationEnergyDTA, 146 incidentParticle, targetParticle, 147 atomicWeight, atomicNumber); 148 if (!successful) 149 HighEnergyClusterProduction(successful, pv, vecLength, 150 excitationEnergyGNP, excitationEnergyDTA, 151 incidentParticle, targetParticle, 152 atomicWeight, atomicNumber); 153 if (!successful) 154 MediumEnergyCascading(successful, pv, vecLength, 155 excitationEnergyGNP, excitationEnergyDTA, 156 incidentParticle, targetParticle, 157 atomicWeight, atomicNumber); 158 159 if (!successful) 160 MediumEnergyClusterProduction(successful, pv, vecLength, 161 excitationEnergyGNP, excitationEnergyDTA, 162 incidentParticle, targetParticle, 163 atomicWeight, atomicNumber); 164 if (!successful) 165 QuasiElasticScattering(successful, pv, vecLength, 166 excitationEnergyGNP, excitationEnergyDTA, 167 incidentParticle, targetParticle, 168 atomicWeight, atomicNumber); 169 if (!successful) 170 ElasticScattering(successful, pv, vecLength, 171 incidentParticle, 172 atomicWeight, atomicNumber); 173 174 if (!successful) 175 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 176 << G4endl; 177 178 FillParticleChange(pv, vecLength); 179 delete [] pv; 180 theParticleChange.SetStatusChange(stopAndKill); 181 return &theParticleChange; 182 } 183 195 184 196 185 void 197 G4HEPionMinusInelastic::FirstIntInCasPionMinus( G4bool &inElastic,198 199 200 G4int &vecLen,201 G4HEVectorincidentParticle,202 G4HEVectortargetParticle)186 G4HEPionMinusInelastic::FirstIntInCasPionMinus(G4bool& inElastic, 187 const G4double availableEnergy, 188 G4HEVector pv[], 189 G4int& vecLen, 190 const G4HEVector& incidentParticle, 191 const G4HEVector& targetParticle) 203 192 204 193 // Pion- undergoes interaction with nucleon within a nucleus. Check if it is … … 209 198 // protons/neutrons by kaons or strange baryons according to the average 210 199 // multiplicity per inelastic reaction. 211 212 { 213 static const G4double expxu = std::log(MAXFLOAT); // upper bound for arg. of exp 214 static const G4double expxl = -expxu; // lower bound for arg. of exp 215 216 static const G4double protb = 0.7; 217 static const G4double neutb = 0.7; 218 static const G4double c = 1.25; 219 220 static const G4int numMul = 1200; 221 static const G4int numSec = 60; 222 223 // G4int neutronCode = Neutron.getCode(); 224 G4int protonCode = Proton.getCode(); 225 226 G4int targetCode = targetParticle.getCode(); 227 // G4double incidentMass = incidentParticle.getMass(); 228 // G4double incidentEnergy = incidentParticle.getEnergy(); 229 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 230 231 static G4bool first = true; 232 static G4double protmul[numMul], protnorm[numSec]; // proton constants 233 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 234 235 // misc. local variables 236 // np = number of pi+, nm = number of pi-, nz = number of pi0 237 238 G4int i, counter, nt, np, nm, nz; 200 { 201 static const G4double expxu = std::log(MAXFLOAT); // upper bound for arg. of exp 202 static const G4double expxl = -expxu; // lower bound for arg. of exp 203 204 static const G4double protb = 0.7; 205 static const G4double neutb = 0.7; 206 static const G4double c = 1.25; 207 208 static const G4int numMul = 1200; 209 static const G4int numSec = 60; 210 211 G4int protonCode = Proton.getCode(); 212 213 G4int targetCode = targetParticle.getCode(); 214 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 215 216 static G4bool first = true; 217 static G4double protmul[numMul], protnorm[numSec]; // proton constants 218 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 219 220 // misc. local variables 221 // np = number of pi+, nm = number of pi-, nz = number of pi0 222 223 G4int i, counter, nt, np, nm, nz; 239 224 240 225 if( first ) 241 { 226 { // compute normalization constants, this will only be done once 242 227 first = false; 243 228 for( i=0; i<numMul; i++ )protmul[i] = 0.0;
Note: See TracChangeset
for help on using the changeset viewer.
