Changeset 1347 for trunk/source/processes/hadronic/models/high_energy/src/G4HEAntiOmegaMinusInelastic.cc
- 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/G4HEAntiOmegaMinusInelastic.cc
r1340 r1347 24 24 // ******************************************************************** 25 25 // 26 // 27 // $Id: G4HEAntiOmegaMinusInelastic.cc,v 1.15 2008/03/17 20:49:17 dennis Exp $ 28 // GEANT4 tag $Name: geant4-09-03-ref-09 $ 29 // 26 // $Id: G4HEAntiOmegaMinusInelastic.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 is 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 "G4HEAntiOmegaMinusInelastic.hh" 46 43 47 G4HadFinalState * G4HEAntiOmegaMinusInelastic:: 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 << "GHEAntiOmegaMinusInelastic: incident energy < 1 GeV" << G4endl; 66 } 67 if(verboseLevel > 1) 68 { 44 G4HadFinalState* 45 G4HEAntiOmegaMinusInelastic::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 << "GHEAntiOmegaMinusInelastic: incident energy < 1 GeV" << G4endl; 62 63 if (verboseLevel > 1) { 69 64 G4cout << "G4HEAntiOmegaMinusInelastic::ApplyYourself" << G4endl; 70 65 G4cout << "incident particle " << incidentParticle.getName() … … 74 69 G4cout << "target material with (A,Z) = (" 75 70 << atomicWeight << "," << atomicNumber << ")" << G4endl; 76 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 86 87 88 G4double excitation= NuclearExcitation(incidentKineticEnergy,89 90 91 92 if(verboseLevel > 1)93 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 94 89 << excitationEnergyDTA << G4endl; 95 90 96 91 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; 92 incidentKineticEnergy -= excitation; 93 incidentTotalEnergy = incidentKineticEnergy + incidentMass; 94 incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass) 95 *(incidentTotalEnergy+incidentMass)); 96 97 G4HEVector targetParticle; 98 if (G4UniformRand() < atomicNumber/atomicWeight) { 99 targetParticle.setDefinition("Proton"); 100 } else { 101 targetParticle.setDefinition("Neutron"); 102 } 103 104 G4double targetMass = targetParticle.getMass(); 105 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass 106 + targetMass*targetMass 107 + 2.0*targetMass*incidentTotalEnergy); 108 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass; 109 110 G4bool inElastic = true; 111 vecLength = 0; 128 112 129 if(verboseLevel > 1)130 113 if (verboseLevel > 1) 114 G4cout << "ApplyYourself: CallFirstIntInCascade for particle " 131 115 << incidentCode << G4endl; 132 116 133 117 G4bool successful = false; 134 118 135 if(inElastic || (!inElastic && atomicWeight < 1.5)) 136 { 137 FirstIntInCasAntiOmegaMinus(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 } 119 FirstIntInCasAntiOmegaMinus(inElastic, availableEnergy, pv, vecLength, 120 incidentParticle, targetParticle, atomicWeight); 121 122 if (verboseLevel > 1) 123 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl; 124 125 if ((vecLength > 0) && (availableEnergy > 1.)) 126 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy, 127 pv, vecLength, 128 incidentParticle, targetParticle); 129 HighEnergyCascading(successful, pv, vecLength, 130 excitationEnergyGNP, excitationEnergyDTA, 131 incidentParticle, targetParticle, 132 atomicWeight, atomicNumber); 133 if (!successful) 134 HighEnergyClusterProduction(successful, pv, vecLength, 135 excitationEnergyGNP, excitationEnergyDTA, 136 incidentParticle, targetParticle, 137 atomicWeight, atomicNumber); 138 if (!successful) 139 MediumEnergyCascading(successful, pv, vecLength, 140 excitationEnergyGNP, excitationEnergyDTA, 141 incidentParticle, targetParticle, 142 atomicWeight, atomicNumber); 143 144 if (!successful) 145 MediumEnergyClusterProduction(successful, pv, vecLength, 146 excitationEnergyGNP, excitationEnergyDTA, 147 incidentParticle, targetParticle, 148 atomicWeight, atomicNumber); 149 if (!successful) 150 QuasiElasticScattering(successful, pv, vecLength, 151 excitationEnergyGNP, excitationEnergyDTA, 152 incidentParticle, targetParticle, 153 atomicWeight, atomicNumber); 154 if (!successful) 155 ElasticScattering(successful, pv, vecLength, 156 incidentParticle, 157 atomicWeight, atomicNumber); 158 159 if (!successful) 160 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 161 << G4endl; 162 163 FillParticleChange(pv, vecLength); 164 delete [] pv; 165 theParticleChange.SetStatusChange(stopAndKill); 166 return &theParticleChange; 167 } 168 190 169 191 170 void 192 G4HEAntiOmegaMinusInelastic::FirstIntInCasAntiOmegaMinus( G4bool &inElastic,193 194 195 G4int &vecLen,196 G4HEVectorincidentParticle,197 G4HEVectortargetParticle,198 const G4double atomicWeight)171 G4HEAntiOmegaMinusInelastic::FirstIntInCasAntiOmegaMinus(G4bool& inElastic, 172 const G4double availableEnergy, 173 G4HEVector pv[], 174 G4int& vecLen, 175 const G4HEVector& incidentParticle, 176 const G4HEVector& targetParticle, 177 const G4double atomicWeight) 199 178 200 179 // AntiOmega undergoes interaction with nucleon within a nucleus. 201 180 // As in Geant3, we think that this routine has absolutely no influence 202 181 // on the whole performance of the program. Take AntiLambda instaed. 203 204 { 205 static const G4double expxu = std::log(MAXFLOAT); // upper bound for arg. of exp 206 static const G4double expxl = -expxu; // lower bound for arg. of exp 207 208 static const G4double protb = 0.7; 209 static const G4double neutb = 0.7; 210 static const G4double c = 1.25; 211 212 static const G4int numMul = 1200; 213 static const G4int numMulAn = 400; 214 static const G4int numSec = 60; 215 216 // G4int neutronCode = Neutron.getCode(); 217 G4int protonCode = Proton.getCode(); 218 219 G4int targetCode = targetParticle.getCode(); 220 // G4double incidentMass = incidentParticle.getMass(); 221 // G4double incidentEnergy = incidentParticle.getEnergy(); 222 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 223 224 static G4bool first = true; 225 static G4double protmul[numMul], protnorm[numSec]; // proton constants 226 static G4double protmulAn[numMulAn],protnormAn[numSec]; 227 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 228 static G4double neutmulAn[numMulAn],neutnormAn[numSec]; 229 230 // misc. local variables 231 // np = number of pi+, nm = number of pi-, nz = number of pi0 182 { 183 static const G4double expxu = std::log(MAXFLOAT); // upper bound for arg. of exp 184 static const G4double expxl = -expxu; // lower bound for arg. of exp 185 186 static const G4double protb = 0.7; 187 static const G4double neutb = 0.7; 188 static const G4double c = 1.25; 189 190 static const G4int numMul = 1200; 191 static const G4int numMulAn = 400; 192 static const G4int numSec = 60; 193 194 G4int protonCode = Proton.getCode(); 195 196 G4int targetCode = targetParticle.getCode(); 197 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum(); 198 199 static G4bool first = true; 200 static G4double protmul[numMul], protnorm[numSec]; // proton constants 201 static G4double protmulAn[numMulAn],protnormAn[numSec]; 202 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants 203 static G4double neutmulAn[numMulAn],neutnormAn[numSec]; 204 205 // misc. local variables 206 // np = number of pi+, nm = number of pi-, nz = number of pi0 232 207 233 208 G4int i, counter, nt, np, nm, nz; 234 209 235 210 if( first ) 236 { 211 { // compute normalization constants, this will only be done once 237 212 first = false; 238 213 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
Note: See TracChangeset
for help on using the changeset viewer.