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