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