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