Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 24.11.08 V. Grichine - first implementation 27 // 28 // 04.09.18 V. Ivantchenko Major revision of i 29 // 27.05.19 V. Ivantchenko Removed obsolete me 30 31 #include "G4ComponentGGNuclNuclXsc.hh" 32 33 #include "G4PhysicalConstants.hh" 34 #include "G4SystemOfUnits.hh" 35 #include "G4NucleiProperties.hh" 36 #include "G4ParticleDefinition.hh" 37 #include "G4HadronNucleonXsc.hh" 38 #include "G4ComponentGGHadronNucleusXsc.hh" 39 #include "G4NuclearRadii.hh" 40 #include "G4Pow.hh" 41 42 namespace 43 { 44 const G4double inve = 1./CLHEP::eplus; 45 } 46 47 G4ComponentGGNuclNuclXsc::G4ComponentGGNuclNuc 48 : G4VComponentCrossSection("Glauber-Gribov Nu 49 { 50 theProton = G4Proton::Proton(); 51 theNeutron = G4Neutron::Neutron(); 52 theLambda = G4Lambda::Lambda(); 53 fHNXsc = new G4HadronNucleonXsc(); 54 fHadrNucl = new G4ComponentGGHadronNucleusXs 55 } 56 57 G4ComponentGGNuclNuclXsc::~G4ComponentGGNuclNu 58 { 59 delete fHNXsc; 60 } 61 62 ////////////////////////////////////////////// 63 64 G4double G4ComponentGGNuclNuclXsc::GetTotalEle 65 const G4ParticleDefinition* aParticle 66 G4int Z, G4double A) 67 { 68 ComputeCrossSections(aParticle, kinEnergy, Z 69 return fTotalXsc; 70 } 71 72 ////////////////////////////////////////////// 73 74 G4double G4ComponentGGNuclNuclXsc::GetTotalIso 75 const G4ParticleDefinition* aParticle 76 G4int Z, G4int A) 77 { 78 ComputeCrossSections(aParticle, kinEnergy, Z 79 return fTotalXsc; 80 } 81 82 ////////////////////////////////////////////// 83 84 G4double G4ComponentGGNuclNuclXsc::GetInelasti 85 const G4ParticleDefinition* aParticle 86 G4int Z, G4double A) 87 { 88 ComputeCrossSections(aParticle, kinEnergy, Z 89 return fInelasticXsc; 90 } 91 92 ////////////////////////////////////////////// 93 94 G4double G4ComponentGGNuclNuclXsc::GetInelasti 95 const G4ParticleDefinition* aParticle 96 G4int Z, G4int A) 97 { 98 ComputeCrossSections(aParticle, kinEnergy, Z 99 return fInelasticXsc; 100 } 101 102 ////////////////////////////////////////////// 103 104 G4double G4ComponentGGNuclNuclXsc::GetElasticE 105 const G4ParticleDefinition* aParticle 106 G4int Z, G4double A) 107 { 108 ComputeCrossSections(aParticle, kinEnergy, Z 109 return fElasticXsc; 110 } 111 112 ////////////////////////////////////////////// 113 114 G4double G4ComponentGGNuclNuclXsc::GetElasticI 115 const G4ParticleDefinition* aParticle 116 G4int Z, G4int A) 117 { 118 ComputeCrossSections(aParticle, kinEnergy, Z 119 return fElasticXsc; 120 } 121 122 ////////////////////////////////////////////// 123 124 G4double G4ComponentGGNuclNuclXsc::ComputeQuas 125 const G4ParticleDefinition* aParticle 126 G4int Z, G4int A) 127 { 128 ComputeCrossSections(aParticle, kinEnergy, Z 129 return (fInelasticXsc > fProductionXsc) 130 ? (fInelasticXsc - fProductionXsc)/fInelas 131 } 132 133 ////////////////////////////////////////////// 134 135 void G4ComponentGGNuclNuclXsc::BuildPhysicsTab 136 {} 137 138 ////////////////////////////////////////////// 139 140 void G4ComponentGGNuclNuclXsc::DumpPhysicsTabl 141 { 142 G4cout << "G4ComponentGGNuclNuclXsc: uses Gl 143 } 144 145 ////////////////////////////////////////////// 146 147 void G4ComponentGGNuclNuclXsc::Description(std 148 { 149 outFile << "G4ComponentGGNuclNuclXsc calcula 150 << "elastic cross sections for nucle 151 << "the Glauber model with Gribov co 152 << "all incident energies above 100 153 << "For the hydrogen target G4HadronNucleo 154 } 155 156 ////////////////////////////////////////////// 157 // 158 // Calculates total and inelastic Xsc, derives 159 // accordong to Glauber model with Gribov corr 160 // approximation on light cone. Gaussian densi 161 // to calculate rest integrals of the model. [ 162 // nucl-th/0306044 + simplification above 163 164 void G4ComponentGGNuclNuclXsc::ComputeCrossSec 165 const G4ParticleDefinition* aParticle, G4 166 G4int Z, G4int A) 167 { 168 // check cache 169 if(aParticle == fParticle && fZ == Z && fA = 170 { return; } 171 fParticle = aParticle; 172 fZ = Z; 173 fA = A; 174 fEnergy = kinEnergy; 175 G4Pow* pG4Pow = G4Pow::GetInstance(); 176 177 G4int pZ = G4lrint(aParticle->GetPDGCharge() 178 G4int pA = aParticle->GetBaryonNumber(); 179 G4int pL = aParticle->GetNumberOfLambdasInHy 180 G4bool pHN = aParticle->IsHypernucleus(); 181 G4double cHN(0.88); 182 183 // hydrogen 184 if(1 == Z && 1 == A) { 185 G4double e = kinEnergy*CLHEP::proton_mass_ 186 fHadrNucl->ComputeCrossSections( theProton 187 fTotalXsc = fHadrNucl->GetTotalGlauberGrib 188 fElasticXsc = fHadrNucl->GetElasticGlauber 189 fInelasticXsc = fHadrNucl->GetInelasticGla 190 fProductionXsc = fHadrNucl->GetProductionG 191 fDiffractionXsc = fHadrNucl->GetDiffractio 192 return; 193 } 194 static const G4double cofInelastic = 2.4; 195 static const G4double cofTotal = 2.0; 196 197 G4double pTkin = kinEnergy/(G4double)pA; 198 199 G4int pN = pA - pZ; 200 G4int tN = A - Z; 201 202 G4double tR = G4NuclearRadii::Radius(Z, A); 203 G4double pR = G4NuclearRadii::Radius(pZ, pA) 204 205 if(pHN) 206 pR *= std::sqrt( pG4Pow->Z23( pA - pL ) + 207 208 G4double cB = ComputeCoulombBarier(aParticle 209 210 if ( cB > 0. ) 211 { 212 G4double sigma = (pZ*Z+pN*tN)*fHNXsc->Hadr 213 if(pHN) sigma += pL*A*fHNXsc->HadronNucleo 214 G4double ppInXsc = fHNXsc->GetInelasticHad 215 216 sigma += (pZ*tN+pN*Z)*fHNXsc->HadronNucleo 217 G4double npInXsc = fHNXsc->GetInelasticHad 218 219 G4double nucleusSquare = cofTotal*CLHEP::p 220 221 G4double ratio= sigma/nucleusSquare; 222 fTotalXsc = nucleusSquare*G4Log( 1. + 223 fInelasticXsc = nucleusSquare*G4Log( 1. + 224 fElasticXsc = std::max(fTotalXsc - fInel 225 226 G4double difratio = ratio/(1.+ratio); 227 fDiffractionXsc = 0.5*nucleusSquare*( difr 228 229 G4double xratio= ((pZ*Z+pN*tN)*ppInXsc + ( 230 fProductionXsc = nucleusSquare*G4Log( 1. + 231 fProductionXsc = std::min(fProductionXsc, 232 } 233 else 234 { 235 fInelasticXsc = fTotalXsc = fElasticXsc = 236 } 237 } 238 239 ////////////////////////////////////////////// 240 241 G4double G4ComponentGGNuclNuclXsc::ComputeCoul 242 const G4ParticleDefinitio 243 G4double pTkin, G4int Z, G4int A, 244 G4double pR, G4double tR) 245 { 246 G4int pZ = aParticle->GetPDGCharge()*inve; 247 G4double pM = aParticle->GetPDGMass(); 248 G4double tM = G4NucleiProperties::GetNuclear 249 G4double pElab = pTkin + pM; 250 G4double totEcm = std::sqrt(pM*pM + tM*tM + 251 G4double totTcm = totEcm - pM -tM; 252 253 // 0.5 defines shape of Cross section correc 254 // at cB = totTcm it become zero 255 static const G4double qfact = 0.5*CLHEP::elm 256 G4double bC = qfact*pZ*Z/(pR + tR); 257 258 G4double ratio = (totTcm <= bC) ? 0. : 1. - 259 260 #ifdef G4VERBOSE 261 if (GetVerboseLevel() > 1) { 262 G4cout << "G4ComponentGGNuclNuclXsc::Compu 263 << "; pTkin(GeV)=" << pTkin/CLHEP::MeV 264 << " totTcm= " << totTcm/CLHEP::MeV<< "; 265 << G4endl; 266 } 267 #endif 268 return ratio; 269 } 270 271 ////////////////////////////////////////////// 272 // 273 // Return single-diffraction/inelastic cross-s 274 275 G4double G4ComponentGGNuclNuclXsc::GetRatioSD( 276 const G4DynamicParticle* aParticle, G 277 { 278 ComputeCrossSections(aParticle->GetDefinitio 279 aParticle->GetKineticEn 280 G4lrint(tZ), G4lrint(tA)); 281 282 return (fInelasticXsc > 0.0) ? fDiffractionX 283 } 284 285 ////////////////////////////////////////////// 286 // 287 // Return quasi-elastic/inelastic cross-sectio 288 289 G4double G4ComponentGGNuclNuclXsc::GetRatioQE( 290 const G4DynamicParticle* aParticle, G 291 { 292 ComputeCrossSections(aParticle->GetDefinitio 293 aParticle->GetKineticEn 294 G4lrint(tZ), G4lrint(tA)); 295 296 return (fInelasticXsc > 0.0) ? 1.0 - fProduc 297 } 298 299 ////////////////////////////////////////////// 300