Changeset 1347 for trunk/source/processes/hadronic/models/high_energy/src/G4HEAntiXiMinusInelastic.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/G4HEAntiXiMinusInelastic.cc
r1340 r1347 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4HEAntiXiMinusInelastic.cc,v 1.15 2008/03/17 20:49:17 dennis Exp $ 28 // GEANT4 tag $Name: geant4-09-03-ref-09 $ 29 // 26 // $Id: G4HEAntiXiMinusInelastic.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 "G4HEAntiXiMinusInelastic.hh" 46 43 47 G4HadFinalState * G4HEAntiXiMinusInelastic::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 54 55 44 G4HadFinalState* 45 G4HEAntiXiMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack, 46 G4Nucleus& targetNucleus) 47 { 48 G4HEVector* pv = new G4HEVector[MAXPART]; 49 const G4HadProjectile* aParticle = &aTrack; 50 const G4double A = targetNucleus.GetN(); 51 const G4double Z = targetNucleus.GetZ(); 52 G4HEVector incidentParticle(aParticle); 56 53 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 { 68 G4cout << "GHEAntiXiMinusInelastic: incident energy < 1 GeV" << G4endl; 69 } 70 if(verboseLevel > 1) 71 { 72 G4cout << "G4HEAntiXiMinusInelastic::ApplyYourself" << G4endl; 73 G4cout << "incident particle " << incidentParticle.getName() 74 << "mass " << incidentMass 75 << "kinetic energy " << incidentKineticEnergy 76 << G4endl; 77 G4cout << "target material with (A,Z) = (" 78 << atomicWeight << "," << atomicNumber << ")" << G4endl; 79 } 80 81 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 82 atomicWeight, atomicNumber); 83 if(verboseLevel > 1) 84 G4cout << "nuclear inelasticity = " << inelasticity << G4endl; 54 G4double atomicNumber = Z; 55 G4double atomicWeight = A; 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 G4cout << "GHEAntiXiMinusInelastic: incident energy < 1 GeV" << G4endl; 65 66 if (verboseLevel > 1) { 67 G4cout << "G4HEAntiXiMinusInelastic::ApplyYourself" << G4endl; 68 G4cout << "incident particle " << incidentParticle.getName() 69 << "mass " << incidentMass 70 << "kinetic energy " << incidentKineticEnergy 71 << G4endl; 72 G4cout << "target material with (A,Z) = (" 73 << atomicWeight << "," << atomicNumber << ")" << G4endl; 74 } 75 76 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 77 atomicWeight, atomicNumber); 78 if (verboseLevel > 1) 79 G4cout << "nuclear inelasticity = " << inelasticity << G4endl; 85 80 86 81 incidentKineticEnergy -= inelasticity; 87 82 88 G4double excitationEnergyGNP = 0.; 89 G4double excitationEnergyDTA = 0.; 90 91 G4double excitation = NuclearExcitation(incidentKineticEnergy, 92 atomicWeight, atomicNumber, 93 excitationEnergyGNP, 94 excitationEnergyDTA); 95 if(verboseLevel > 1) 96 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 97 << excitationEnergyDTA << G4endl; 98 99 100 incidentKineticEnergy -= excitation; 101 incidentTotalEnergy = incidentKineticEnergy + incidentMass; 102 incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass) 103 *(incidentTotalEnergy+incidentMass)); 104 105 106 G4HEVector targetParticle; 107 if(G4UniformRand() < atomicNumber/atomicWeight) 108 { 109 targetParticle.setDefinition("Proton"); 110 } 111 else 112 { 113 targetParticle.setDefinition("Neutron"); 114 } 115 116 G4double targetMass = targetParticle.getMass(); 117 G4double centerOfMassEnergy = std::sqrt( incidentMass*incidentMass + targetMass*targetMass 118 + 2.0*targetMass*incidentTotalEnergy); 119 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass; 120 121 // this was the meaning of inElastic in the 122 // original Gheisha stand-alone version. 123 // G4bool inElastic = InElasticCrossSectionInFirstInt 124 // (availableEnergy, incidentCode, incidentTotalMomentum); 125 // by unknown reasons, it has been replaced 126 // to the following code in Geant??? 127 G4bool inElastic = true; 128 // if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false; 129 130 vecLength = 0; 83 G4double excitationEnergyGNP = 0.; 84 G4double excitationEnergyDTA = 0.; 85 86 G4double excitation = NuclearExcitation(incidentKineticEnergy, 87 atomicWeight, atomicNumber, 88 excitationEnergyGNP, 89 excitationEnergyDTA); 90 if (verboseLevel > 1) 91 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 92 << excitationEnergyDTA << G4endl; 93 94 incidentKineticEnergy -= excitation; 95 incidentTotalEnergy = incidentKineticEnergy + incidentMass; 96 incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass) 97 *(incidentTotalEnergy+incidentMass)); 98 99 G4HEVector targetParticle; 100 if (G4UniformRand() < atomicNumber/atomicWeight) { 101 targetParticle.setDefinition("Proton"); 102 } else { 103 targetParticle.setDefinition("Neutron"); 104 } 105 106 G4double targetMass = targetParticle.getMass(); 107 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass 108 + targetMass*targetMass 109 + 2.0*targetMass*incidentTotalEnergy); 110 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass; 111 112 G4bool inElastic = true; 113 vecLength = 0; 131 114 132 if(verboseLevel > 1)133 115 if (verboseLevel > 1) 116 G4cout << "ApplyYourself: CallFirstIntInCascade for particle " 134 117 << incidentCode << G4endl; 135 118 136 119 G4bool successful = false; 137 120 138 if(inElastic || (!inElastic && atomicWeight < 1.5)) 139 { 140 FirstIntInCasAntiXiMinus(inElastic, availableEnergy, pv, vecLength, 141 incidentParticle, targetParticle, atomicWeight); 142 143 if(verboseLevel > 1) 144 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl; 145 146 147 if ((vecLength > 0) && (availableEnergy > 1.)) 148 StrangeParticlePairProduction( availableEnergy, centerOfMassEnergy, 149 pv, vecLength, 150 incidentParticle, targetParticle); 151 HighEnergyCascading( successful, pv, vecLength, 152 excitationEnergyGNP, excitationEnergyDTA, 153 incidentParticle, targetParticle, 154 atomicWeight, atomicNumber); 155 if (!successful) 156 HighEnergyClusterProduction( successful, pv, vecLength, 157 excitationEnergyGNP, excitationEnergyDTA, 158 incidentParticle, targetParticle, 159 atomicWeight, atomicNumber); 160 if (!successful) 161 MediumEnergyCascading( successful, pv, vecLength, 162 excitationEnergyGNP, excitationEnergyDTA, 163 incidentParticle, targetParticle, 164 atomicWeight, atomicNumber); 165 166 if (!successful) 167 MediumEnergyClusterProduction( successful, pv, vecLength, 168 excitationEnergyGNP, excitationEnergyDTA, 169 incidentParticle, targetParticle, 170 atomicWeight, atomicNumber); 171 if (!successful) 172 QuasiElasticScattering( successful, pv, vecLength, 173 excitationEnergyGNP, excitationEnergyDTA, 174 incidentParticle, targetParticle, 175 atomicWeight, atomicNumber); 176 } 177 if (!successful) 178 { 179 ElasticScattering( successful, pv, vecLength, 180 incidentParticle, 181 atomicWeight, atomicNumber); 182 } 183 184 if (!successful) 185 { 186 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"; 187 } 188 FillParticleChange(pv, vecLength); 189 delete [] pv; 190 theParticleChange.SetStatusChange(stopAndKill); 191 return & theParticleChange; 192 } 121 FirstIntInCasAntiXiMinus(inElastic, availableEnergy, pv, vecLength, 122 incidentParticle, targetParticle, atomicWeight); 123 124 if (verboseLevel > 1) 125 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl; 126 127 if ((vecLength > 0) && (availableEnergy > 1.)) 128 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy, 129 pv, vecLength, 130 incidentParticle, targetParticle); 131 HighEnergyCascading(successful, pv, vecLength, 132 excitationEnergyGNP, excitationEnergyDTA, 133 incidentParticle, targetParticle, 134 atomicWeight, atomicNumber); 135 if (!successful) 136 HighEnergyClusterProduction(successful, pv, vecLength, 137 excitationEnergyGNP, excitationEnergyDTA, 138 incidentParticle, targetParticle, 139 atomicWeight, atomicNumber); 140 if (!successful) 141 MediumEnergyCascading(successful, pv, vecLength, 142 excitationEnergyGNP, excitationEnergyDTA, 143 incidentParticle, targetParticle, 144 atomicWeight, atomicNumber); 145 146 if (!successful) 147 MediumEnergyClusterProduction(successful, pv, vecLength, 148 excitationEnergyGNP, excitationEnergyDTA, 149 incidentParticle, targetParticle, 150 atomicWeight, atomicNumber); 151 if (!successful) 152 QuasiElasticScattering(successful, pv, vecLength, 153 excitationEnergyGNP, excitationEnergyDTA, 154 incidentParticle, targetParticle, 155 atomicWeight, atomicNumber); 156 if (!successful) 157 ElasticScattering(successful, pv, vecLength, 158 incidentParticle, 159 atomicWeight, atomicNumber); 160 161 if (!successful) 162 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 163 << G4endl; 164 165 FillParticleChange(pv, vecLength); 166 delete [] pv; 167 theParticleChange.SetStatusChange(stopAndKill); 168 return &theParticleChange; 169 } 170 193 171 194 172 void 195 G4HEAntiXiMinusInelastic::FirstIntInCasAntiXiMinus( G4bool &inElastic,196 197 198 G4int &vecLen,199 G4HEVectorincidentParticle,200 G4HEVectortargetParticle,201 173 G4HEAntiXiMinusInelastic::FirstIntInCasAntiXiMinus(G4bool& inElastic, 174 const G4double availableEnergy, 175 G4HEVector pv[], 176 G4int& vecLen, 177 const G4HEVector& incidentParticle, 178 const G4HEVector& targetParticle, 179 const G4double atomicWeight) 202 180 203 181 // AntiXi- undergoes interaction with nucleon within a nucleus. 204 182 // As in Geant3, we think that this routine has absolutely no influence 205 183 // on the whole performance of the program. Take AntiLambda instaed. 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 protonCode = Proton.getCode(); 197 198 G4int targetCode = targetParticle.getCode(); 199 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 200 201 static G4bool first = true; 202 static G4double protmul[numMul], protnorm[numSec]; // proton constants 203 static G4double protmulAn[numMulAn],protnormAn[numSec]; 204 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 205 static G4double neutmulAn[numMulAn],neutnormAn[numSec]; 206 207 // misc. local variables 208 // np = number of pi+, nm = number of pi-, nz = number of pi0 209 210 G4int i, counter, nt, np, nm, nz; 237 211 238 212 if( first ) 239 { 213 { // compute normalization constants, this will only be done once 240 214 first = false; 241 215 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
Note: See TracChangeset
for help on using the changeset viewer.
