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 // 26 // 27 // 27 // 28 // Physics model class G4NuclNuclDiffuseElasti 28 // Physics model class G4NuclNuclDiffuseElastic 29 // 29 // 30 // 30 // 31 // G4 Model: optical diffuse elastic scatterin 31 // G4 Model: optical diffuse elastic scattering with 4-momentum balance 32 // 32 // 33 // 24-May-07 V. Grichine 33 // 24-May-07 V. Grichine 34 // 34 // 35 35 36 #include "G4NuclNuclDiffuseElastic.hh" 36 #include "G4NuclNuclDiffuseElastic.hh" 37 #include "G4ParticleTable.hh" 37 #include "G4ParticleTable.hh" 38 #include "G4ParticleDefinition.hh" 38 #include "G4ParticleDefinition.hh" 39 #include "G4IonTable.hh" 39 #include "G4IonTable.hh" 40 #include "G4NucleiProperties.hh" 40 #include "G4NucleiProperties.hh" 41 41 42 #include "Randomize.hh" 42 #include "Randomize.hh" 43 #include "G4Integrator.hh" 43 #include "G4Integrator.hh" 44 #include "globals.hh" 44 #include "globals.hh" 45 #include "G4PhysicalConstants.hh" 45 #include "G4PhysicalConstants.hh" 46 #include "G4SystemOfUnits.hh" 46 #include "G4SystemOfUnits.hh" 47 47 48 #include "G4Proton.hh" 48 #include "G4Proton.hh" 49 #include "G4Neutron.hh" 49 #include "G4Neutron.hh" 50 #include "G4Deuteron.hh" 50 #include "G4Deuteron.hh" 51 #include "G4Alpha.hh" 51 #include "G4Alpha.hh" 52 #include "G4PionPlus.hh" 52 #include "G4PionPlus.hh" 53 #include "G4PionMinus.hh" 53 #include "G4PionMinus.hh" 54 54 55 #include "G4Element.hh" 55 #include "G4Element.hh" 56 #include "G4ElementTable.hh" 56 #include "G4ElementTable.hh" 57 #include "G4NistManager.hh" 57 #include "G4NistManager.hh" 58 #include "G4PhysicsTable.hh" 58 #include "G4PhysicsTable.hh" 59 #include "G4PhysicsLogVector.hh" 59 #include "G4PhysicsLogVector.hh" 60 #include "G4PhysicsFreeVector.hh" 60 #include "G4PhysicsFreeVector.hh" 61 #include "G4HadronicParameters.hh" 61 #include "G4HadronicParameters.hh" 62 62 63 ////////////////////////////////////////////// 63 ///////////////////////////////////////////////////////////////////////// 64 // 64 // 65 // Test Constructor. Just to check xsc 65 // Test Constructor. Just to check xsc 66 66 67 67 68 G4NuclNuclDiffuseElastic::G4NuclNuclDiffuseEla 68 G4NuclNuclDiffuseElastic::G4NuclNuclDiffuseElastic() 69 : G4HadronElastic("NNDiffuseElastic"), fPart 69 : G4HadronElastic("NNDiffuseElastic"), fParticle(0) 70 { 70 { 71 SetMinEnergy( 50*MeV ); 71 SetMinEnergy( 50*MeV ); 72 SetMaxEnergy( G4HadronicParameters::Instance 72 SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() ); 73 verboseLevel = 0; 73 verboseLevel = 0; 74 lowEnergyRecoilLimit = 100.*keV; 74 lowEnergyRecoilLimit = 100.*keV; 75 lowEnergyLimitQ = 0.0*GeV; 75 lowEnergyLimitQ = 0.0*GeV; 76 lowEnergyLimitHE = 0.0*GeV; 76 lowEnergyLimitHE = 0.0*GeV; 77 lowestEnergyLimit= 0.0*keV; 77 lowestEnergyLimit= 0.0*keV; 78 plabLowLimit = 20.0*MeV; 78 plabLowLimit = 20.0*MeV; 79 79 80 theProton = G4Proton::Proton(); 80 theProton = G4Proton::Proton(); 81 theNeutron = G4Neutron::Neutron(); 81 theNeutron = G4Neutron::Neutron(); 82 theDeuteron = G4Deuteron::Deuteron(); 82 theDeuteron = G4Deuteron::Deuteron(); 83 theAlpha = G4Alpha::Alpha(); 83 theAlpha = G4Alpha::Alpha(); 84 thePionPlus = G4PionPlus::PionPlus(); 84 thePionPlus = G4PionPlus::PionPlus(); 85 thePionMinus= G4PionMinus::PionMinus(); 85 thePionMinus= G4PionMinus::PionMinus(); 86 86 87 fEnergyBin = 300; // Increased from the ori 87 fEnergyBin = 300; // Increased from the original 200 to have no wider log-energy-bins up to 10 PeV 88 fAngleBin = 200; 88 fAngleBin = 200; 89 89 90 fEnergyVector = new G4PhysicsLogVector( the 90 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin ); 91 fAngleTable = 0; 91 fAngleTable = 0; 92 92 93 fParticle = 0; 93 fParticle = 0; 94 fWaveVector = 0.; 94 fWaveVector = 0.; 95 fAtomicWeight = 0.; 95 fAtomicWeight = 0.; 96 fAtomicNumber = 0.; 96 fAtomicNumber = 0.; 97 fNuclearRadius = 0.; 97 fNuclearRadius = 0.; 98 fBeta = 0.; 98 fBeta = 0.; 99 fZommerfeld = 0.; 99 fZommerfeld = 0.; 100 fAm = 0.; 100 fAm = 0.; 101 fAddCoulomb = false; 101 fAddCoulomb = false; 102 // Ranges of angle table relative to current 102 // Ranges of angle table relative to current Rutherford (Coulomb grazing) angle 103 103 104 // Empirical parameters 104 // Empirical parameters 105 105 106 fCofAlphaMax = 1.5; 106 fCofAlphaMax = 1.5; 107 fCofAlphaCoulomb = 0.5; 107 fCofAlphaCoulomb = 0.5; 108 108 109 fProfileDelta = 1.; 109 fProfileDelta = 1.; 110 fProfileAlpha = 0.5; 110 fProfileAlpha = 0.5; 111 111 112 fCofLambda = 1.0; 112 fCofLambda = 1.0; 113 fCofDelta = 0.04; 113 fCofDelta = 0.04; 114 fCofAlpha = 0.095; 114 fCofAlpha = 0.095; 115 115 116 fNuclearRadius1 = fNuclearRadius2 = fNuclear 116 fNuclearRadius1 = fNuclearRadius2 = fNuclearRadiusSquare 117 = fRutherfordRatio = fCoulombPhase0 = fHal 117 = fRutherfordRatio = fCoulombPhase0 = fHalfRutThetaTg = fHalfRutThetaTg2 118 = fRutherfordTheta = fProfileLambda = fCof 118 = fRutherfordTheta = fProfileLambda = fCofPhase = fCofFar 119 = fSumSigma = fEtaRatio = fReZ = 0.0; 119 = fSumSigma = fEtaRatio = fReZ = 0.0; 120 fMaxL = 0; 120 fMaxL = 0; 121 121 122 fNuclearRadiusCof = 1.0; 122 fNuclearRadiusCof = 1.0; 123 fCoulombMuC = 0.0; 123 fCoulombMuC = 0.0; 124 } 124 } 125 125 126 ////////////////////////////////////////////// 126 ////////////////////////////////////////////////////////////////////////////// 127 // 127 // 128 // Destructor 128 // Destructor 129 129 130 G4NuclNuclDiffuseElastic::~G4NuclNuclDiffuseEl 130 G4NuclNuclDiffuseElastic::~G4NuclNuclDiffuseElastic() 131 { 131 { 132 if ( fEnergyVector ) { 132 if ( fEnergyVector ) { 133 delete fEnergyVector; 133 delete fEnergyVector; 134 fEnergyVector = 0; 134 fEnergyVector = 0; 135 } 135 } 136 136 137 for ( std::vector<G4PhysicsTable*>::iterator 137 for ( std::vector<G4PhysicsTable*>::iterator it = fAngleBank.begin(); 138 it != fAngleBank.end(); ++it ) { 138 it != fAngleBank.end(); ++it ) { 139 if ( (*it) ) (*it)->clearAndDestroy(); 139 if ( (*it) ) (*it)->clearAndDestroy(); 140 delete *it; 140 delete *it; 141 *it = 0; 141 *it = 0; 142 } 142 } 143 fAngleTable = 0; 143 fAngleTable = 0; 144 } 144 } 145 145 146 ////////////////////////////////////////////// 146 ////////////////////////////////////////////////////////////////////////////// 147 // 147 // 148 // Initialisation for given particle using ele 148 // Initialisation for given particle using element table of application 149 149 150 void G4NuclNuclDiffuseElastic::Initialise() 150 void G4NuclNuclDiffuseElastic::Initialise() 151 { 151 { 152 152 153 // fEnergyVector = new G4PhysicsLogVector( t 153 // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin ); 154 154 155 const G4ElementTable* theElementTable = G4El 155 const G4ElementTable* theElementTable = G4Element::GetElementTable(); 156 std::size_t jEl, numOfEl = G4Element::GetNum << 156 size_t jEl, numOfEl = G4Element::GetNumberOfElements(); 157 157 158 // projectile radius 158 // projectile radius 159 159 160 G4double A1 = G4double( fParticle->GetBaryon 160 G4double A1 = G4double( fParticle->GetBaryonNumber() ); 161 G4double R1 = CalculateNuclearRad(A1); 161 G4double R1 = CalculateNuclearRad(A1); 162 162 163 for(jEl = 0 ; jEl < numOfEl; ++jEl) // appli 163 for(jEl = 0 ; jEl < numOfEl; ++jEl) // application element loop 164 { 164 { 165 fAtomicNumber = (*theElementTable)[jEl]->G 165 fAtomicNumber = (*theElementTable)[jEl]->GetZ(); // atomic number 166 fAtomicWeight = G4NistManager::Instance()- 166 fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( fAtomicNumber ) ); 167 167 168 fNuclearRadius = CalculateNuclearRad(fAtom 168 fNuclearRadius = CalculateNuclearRad(fAtomicWeight); 169 fNuclearRadius += R1; 169 fNuclearRadius += R1; 170 170 171 if(verboseLevel > 0) 171 if(verboseLevel > 0) 172 { 172 { 173 G4cout<<"G4NuclNuclDiffuseElastic::Initi 173 G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element: " 174 <<(*theElementTable)[jEl]->GetName()<<G4 174 <<(*theElementTable)[jEl]->GetName()<<G4endl; 175 } 175 } 176 fElementNumberVector.push_back(fAtomicNumb 176 fElementNumberVector.push_back(fAtomicNumber); 177 fElementNameVector.push_back((*theElementT 177 fElementNameVector.push_back((*theElementTable)[jEl]->GetName()); 178 178 179 BuildAngleTable(); 179 BuildAngleTable(); 180 fAngleBank.push_back(fAngleTable); 180 fAngleBank.push_back(fAngleTable); 181 } 181 } 182 } 182 } 183 183 184 184 185 ////////////////////////////////////////////// 185 //////////////////////////////////////////////////////////////////////////// 186 // 186 // 187 // return differential elastic cross section d 187 // return differential elastic cross section d(sigma)/d(omega) 188 188 189 G4double 189 G4double 190 G4NuclNuclDiffuseElastic::GetDiffuseElasticXsc 190 G4NuclNuclDiffuseElastic::GetDiffuseElasticXsc( const G4ParticleDefinition* particle, 191 G4doub 191 G4double theta, 192 G4double momentum, 192 G4double momentum, 193 G4doub 193 G4double A ) 194 { 194 { 195 fParticle = particle; 195 fParticle = particle; 196 fWaveVector = momentum/hbarc; 196 fWaveVector = momentum/hbarc; 197 fAtomicWeight = A; 197 fAtomicWeight = A; 198 fAddCoulomb = false; 198 fAddCoulomb = false; 199 fNuclearRadius = CalculateNuclearRad(A); 199 fNuclearRadius = CalculateNuclearRad(A); 200 200 201 G4double sigma = fNuclearRadius*fNuclearRadi 201 G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta); 202 202 203 return sigma; 203 return sigma; 204 } 204 } 205 205 206 ////////////////////////////////////////////// 206 //////////////////////////////////////////////////////////////////////////// 207 // 207 // 208 // return invariant differential elastic cross 208 // return invariant differential elastic cross section d(sigma)/d(tMand) 209 209 210 G4double 210 G4double 211 G4NuclNuclDiffuseElastic::GetInvElasticXsc( co 211 G4NuclNuclDiffuseElastic::GetInvElasticXsc( const G4ParticleDefinition* particle, 212 G4doub 212 G4double tMand, 213 G4double plab, 213 G4double plab, 214 G4doub 214 G4double A, G4double Z ) 215 { 215 { 216 G4double m1 = particle->GetPDGMass(); 216 G4double m1 = particle->GetPDGMass(); 217 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla 217 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1)); 218 218 219 G4int iZ = static_cast<G4int>(Z+0.5); 219 G4int iZ = static_cast<G4int>(Z+0.5); 220 G4int iA = static_cast<G4int>(A+0.5); 220 G4int iA = static_cast<G4int>(A+0.5); 221 G4ParticleDefinition * theDef = 0; 221 G4ParticleDefinition * theDef = 0; 222 222 223 if (iZ == 1 && iA == 1) theDef = thePro 223 if (iZ == 1 && iA == 1) theDef = theProton; 224 else if (iZ == 1 && iA == 2) theDef = theDeu 224 else if (iZ == 1 && iA == 2) theDef = theDeuteron; 225 else if (iZ == 1 && iA == 3) theDef = G4Trit 225 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton(); 226 else if (iZ == 2 && iA == 3) theDef = G4He3: 226 else if (iZ == 2 && iA == 3) theDef = G4He3::He3(); 227 else if (iZ == 2 && iA == 4) theDef = theAlp 227 else if (iZ == 2 && iA == 4) theDef = theAlpha; 228 else theDef = G4ParticleTable::GetParticleTa 228 else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0); 229 229 230 G4double tmass = theDef->GetPDGMass(); 230 G4double tmass = theDef->GetPDGMass(); 231 231 232 G4LorentzVector lv(0.0,0.0,0.0,tmass); 232 G4LorentzVector lv(0.0,0.0,0.0,tmass); 233 lv += lv1; 233 lv += lv1; 234 234 235 G4ThreeVector bst = lv.boostVector(); 235 G4ThreeVector bst = lv.boostVector(); 236 lv1.boost(-bst); 236 lv1.boost(-bst); 237 237 238 G4ThreeVector p1 = lv1.vect(); 238 G4ThreeVector p1 = lv1.vect(); 239 G4double ptot = p1.mag(); 239 G4double ptot = p1.mag(); 240 G4double ptot2 = ptot*ptot; 240 G4double ptot2 = ptot*ptot; 241 G4double cost = 1 - 0.5*std::fabs(tMand)/pto 241 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2; 242 242 243 if( cost >= 1.0 ) cost = 1.0; 243 if( cost >= 1.0 ) cost = 1.0; 244 else if( cost <= -1.0) cost = -1.0; 244 else if( cost <= -1.0) cost = -1.0; 245 245 246 G4double thetaCMS = std::acos(cost); 246 G4double thetaCMS = std::acos(cost); 247 247 248 G4double sigma = GetDiffuseElasticXsc( parti 248 G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A); 249 249 250 sigma *= pi/ptot2; 250 sigma *= pi/ptot2; 251 251 252 return sigma; 252 return sigma; 253 } 253 } 254 254 255 ////////////////////////////////////////////// 255 //////////////////////////////////////////////////////////////////////////// 256 // 256 // 257 // return differential elastic cross section d 257 // return differential elastic cross section d(sigma)/d(omega) with Coulomb 258 // correction 258 // correction 259 259 260 G4double 260 G4double 261 G4NuclNuclDiffuseElastic::GetDiffuseElasticSum 261 G4NuclNuclDiffuseElastic::GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle, 262 G4doub 262 G4double theta, 263 G4double momentum, 263 G4double momentum, 264 G4doub 264 G4double A, G4double Z ) 265 { 265 { 266 fParticle = particle; 266 fParticle = particle; 267 fWaveVector = momentum/hbarc; 267 fWaveVector = momentum/hbarc; 268 fAtomicWeight = A; 268 fAtomicWeight = A; 269 fAtomicNumber = Z; 269 fAtomicNumber = Z; 270 fNuclearRadius = CalculateNuclearRad(A); 270 fNuclearRadius = CalculateNuclearRad(A); 271 fAddCoulomb = false; 271 fAddCoulomb = false; 272 272 273 G4double z = particle->GetPDGCharge(); 273 G4double z = particle->GetPDGCharge(); 274 274 275 G4double kRt = fWaveVector*fNuclearRadius* 275 G4double kRt = fWaveVector*fNuclearRadius*theta; 276 G4double kRtC = 1.9; 276 G4double kRtC = 1.9; 277 277 278 if( z && (kRt > kRtC) ) 278 if( z && (kRt > kRtC) ) 279 { 279 { 280 fAddCoulomb = true; 280 fAddCoulomb = true; 281 fBeta = CalculateParticleBeta( parti 281 fBeta = CalculateParticleBeta( particle, momentum); 282 fZommerfeld = CalculateZommerfeld( fBeta, 282 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber); 283 fAm = CalculateAm( momentum, fZomm 283 fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber); 284 } 284 } 285 G4double sigma = fNuclearRadius*fNuclearRadi 285 G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta); 286 286 287 return sigma; 287 return sigma; 288 } 288 } 289 289 290 ////////////////////////////////////////////// 290 //////////////////////////////////////////////////////////////////////////// 291 // 291 // 292 // return invariant differential elastic cross 292 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb 293 // correction 293 // correction 294 294 295 G4double 295 G4double 296 G4NuclNuclDiffuseElastic::GetInvElasticSumXsc( 296 G4NuclNuclDiffuseElastic::GetInvElasticSumXsc( const G4ParticleDefinition* particle, 297 G4doub 297 G4double tMand, 298 G4double plab, 298 G4double plab, 299 G4doub 299 G4double A, G4double Z ) 300 { 300 { 301 G4double m1 = particle->GetPDGMass(); 301 G4double m1 = particle->GetPDGMass(); 302 302 303 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla 303 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1)); 304 304 305 G4int iZ = static_cast<G4int>(Z+0.5); 305 G4int iZ = static_cast<G4int>(Z+0.5); 306 G4int iA = static_cast<G4int>(A+0.5); 306 G4int iA = static_cast<G4int>(A+0.5); 307 307 308 G4ParticleDefinition* theDef = 0; 308 G4ParticleDefinition* theDef = 0; 309 309 310 if (iZ == 1 && iA == 1) theDef = thePro 310 if (iZ == 1 && iA == 1) theDef = theProton; 311 else if (iZ == 1 && iA == 2) theDef = theDeu 311 else if (iZ == 1 && iA == 2) theDef = theDeuteron; 312 else if (iZ == 1 && iA == 3) theDef = G4Trit 312 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton(); 313 else if (iZ == 2 && iA == 3) theDef = G4He3: 313 else if (iZ == 2 && iA == 3) theDef = G4He3::He3(); 314 else if (iZ == 2 && iA == 4) theDef = theAlp 314 else if (iZ == 2 && iA == 4) theDef = theAlpha; 315 else theDef = G4ParticleTable::GetParticleTa 315 else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0); 316 316 317 G4double tmass = theDef->GetPDGMass(); 317 G4double tmass = theDef->GetPDGMass(); 318 318 319 G4LorentzVector lv(0.0,0.0,0.0,tmass); 319 G4LorentzVector lv(0.0,0.0,0.0,tmass); 320 lv += lv1; 320 lv += lv1; 321 321 322 G4ThreeVector bst = lv.boostVector(); 322 G4ThreeVector bst = lv.boostVector(); 323 lv1.boost(-bst); 323 lv1.boost(-bst); 324 324 325 G4ThreeVector p1 = lv1.vect(); 325 G4ThreeVector p1 = lv1.vect(); 326 G4double ptot = p1.mag(); 326 G4double ptot = p1.mag(); 327 G4double ptot2 = ptot*ptot; 327 G4double ptot2 = ptot*ptot; 328 G4double cost = 1 - 0.5*std::fabs(tMand)/ 328 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2; 329 329 330 if( cost >= 1.0 ) cost = 1.0; 330 if( cost >= 1.0 ) cost = 1.0; 331 else if( cost <= -1.0) cost = -1.0; 331 else if( cost <= -1.0) cost = -1.0; 332 332 333 G4double thetaCMS = std::acos(cost); 333 G4double thetaCMS = std::acos(cost); 334 334 335 G4double sigma = GetDiffuseElasticSumXsc( pa 335 G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z ); 336 336 337 sigma *= pi/ptot2; 337 sigma *= pi/ptot2; 338 338 339 return sigma; 339 return sigma; 340 } 340 } 341 341 342 ////////////////////////////////////////////// 342 //////////////////////////////////////////////////////////////////////////// 343 // 343 // 344 // return invariant differential elastic cross 344 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb 345 // correction 345 // correction 346 346 347 G4double 347 G4double 348 G4NuclNuclDiffuseElastic::GetInvCoulombElastic 348 G4NuclNuclDiffuseElastic::GetInvCoulombElasticXsc( const G4ParticleDefinition* particle, 349 G4doub 349 G4double tMand, 350 G4double plab, 350 G4double plab, 351 G4doub 351 G4double A, G4double Z ) 352 { 352 { 353 G4double m1 = particle->GetPDGMass(); 353 G4double m1 = particle->GetPDGMass(); 354 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla 354 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1)); 355 355 356 G4int iZ = static_cast<G4int>(Z+0.5); 356 G4int iZ = static_cast<G4int>(Z+0.5); 357 G4int iA = static_cast<G4int>(A+0.5); 357 G4int iA = static_cast<G4int>(A+0.5); 358 G4ParticleDefinition * theDef = 0; 358 G4ParticleDefinition * theDef = 0; 359 359 360 if (iZ == 1 && iA == 1) theDef = thePro 360 if (iZ == 1 && iA == 1) theDef = theProton; 361 else if (iZ == 1 && iA == 2) theDef = theDeu 361 else if (iZ == 1 && iA == 2) theDef = theDeuteron; 362 else if (iZ == 1 && iA == 3) theDef = G4Trit 362 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton(); 363 else if (iZ == 2 && iA == 3) theDef = G4He3: 363 else if (iZ == 2 && iA == 3) theDef = G4He3::He3(); 364 else if (iZ == 2 && iA == 4) theDef = theAlp 364 else if (iZ == 2 && iA == 4) theDef = theAlpha; 365 else theDef = G4ParticleTable::GetParticleTa 365 else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0); 366 366 367 G4double tmass = theDef->GetPDGMass(); 367 G4double tmass = theDef->GetPDGMass(); 368 368 369 G4LorentzVector lv(0.0,0.0,0.0,tmass); 369 G4LorentzVector lv(0.0,0.0,0.0,tmass); 370 lv += lv1; 370 lv += lv1; 371 371 372 G4ThreeVector bst = lv.boostVector(); 372 G4ThreeVector bst = lv.boostVector(); 373 lv1.boost(-bst); 373 lv1.boost(-bst); 374 374 375 G4ThreeVector p1 = lv1.vect(); 375 G4ThreeVector p1 = lv1.vect(); 376 G4double ptot = p1.mag(); 376 G4double ptot = p1.mag(); 377 G4double ptot2 = ptot*ptot; 377 G4double ptot2 = ptot*ptot; 378 G4double cost = 1 - 0.5*std::fabs(tMand)/pto 378 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2; 379 379 380 if( cost >= 1.0 ) cost = 1.0; 380 if( cost >= 1.0 ) cost = 1.0; 381 else if( cost <= -1.0) cost = -1.0; 381 else if( cost <= -1.0) cost = -1.0; 382 382 383 G4double thetaCMS = std::acos(cost); 383 G4double thetaCMS = std::acos(cost); 384 384 385 G4double sigma = GetCoulombElasticXsc( parti 385 G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z ); 386 386 387 sigma *= pi/ptot2; 387 sigma *= pi/ptot2; 388 388 389 return sigma; 389 return sigma; 390 } 390 } 391 391 392 ////////////////////////////////////////////// 392 //////////////////////////////////////////////////////////////////////////// 393 // 393 // 394 // return differential elastic probability d(p 394 // return differential elastic probability d(probability)/d(omega) 395 395 396 G4double 396 G4double 397 G4NuclNuclDiffuseElastic::GetDiffElasticProb( 397 G4NuclNuclDiffuseElastic::GetDiffElasticProb( // G4ParticleDefinition* particle, 398 G4doub 398 G4double theta 399 // G4double momentum, 399 // G4double momentum, 400 // G4double A 400 // G4double A 401 ) 401 ) 402 { 402 { 403 G4double sigma, bzero, bzero2, bonebyarg, bo 403 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2; 404 G4double delta, diffuse, gamma; 404 G4double delta, diffuse, gamma; 405 G4double e1, e2, bone, bone2; 405 G4double e1, e2, bone, bone2; 406 406 407 // G4double wavek = momentum/hbarc; // wave 407 // G4double wavek = momentum/hbarc; // wave vector 408 // G4double r0 = 1.08*fermi; 408 // G4double r0 = 1.08*fermi; 409 // G4double rad = r0*G4Pow::GetInstance()- 409 // G4double rad = r0*G4Pow::GetInstance()->A13(A); 410 410 411 G4double kr = fWaveVector*fNuclearRadius; 411 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad; 412 G4double kr2 = kr*kr; 412 G4double kr2 = kr*kr; 413 G4double krt = kr*theta; 413 G4double krt = kr*theta; 414 414 415 bzero = BesselJzero(krt); 415 bzero = BesselJzero(krt); 416 bzero2 = bzero*bzero; 416 bzero2 = bzero*bzero; 417 bone = BesselJone(krt); 417 bone = BesselJone(krt); 418 bone2 = bone*bone; 418 bone2 = bone*bone; 419 bonebyarg = BesselOneByArg(krt); 419 bonebyarg = BesselOneByArg(krt); 420 bonebyarg2 = bonebyarg*bonebyarg; 420 bonebyarg2 = bonebyarg*bonebyarg; 421 421 422 // VI - Coverity complains 422 // VI - Coverity complains 423 /* 423 /* 424 if (fParticle == theProton) 424 if (fParticle == theProton) 425 { 425 { 426 diffuse = 0.63*fermi; 426 diffuse = 0.63*fermi; 427 gamma = 0.3*fermi; 427 gamma = 0.3*fermi; 428 delta = 0.1*fermi*fermi; 428 delta = 0.1*fermi*fermi; 429 e1 = 0.3*fermi; 429 e1 = 0.3*fermi; 430 e2 = 0.35*fermi; 430 e2 = 0.35*fermi; 431 } 431 } 432 else // as proton, if were not defined 432 else // as proton, if were not defined 433 { 433 { 434 */ 434 */ 435 diffuse = 0.63*fermi; 435 diffuse = 0.63*fermi; 436 gamma = 0.3*fermi; 436 gamma = 0.3*fermi; 437 delta = 0.1*fermi*fermi; 437 delta = 0.1*fermi*fermi; 438 e1 = 0.3*fermi; 438 e1 = 0.3*fermi; 439 e2 = 0.35*fermi; 439 e2 = 0.35*fermi; 440 //} 440 //} 441 G4double lambda = 15.; // 15 ok 441 G4double lambda = 15.; // 15 ok 442 442 443 // G4double kgamma = fWaveVector*gamma; 443 // G4double kgamma = fWaveVector*gamma; // wavek*delta; 444 444 445 G4double kgamma = lambda*(1.-G4Exp(-fWave 445 G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta; 446 G4double kgamma2 = kgamma*kgamma; 446 G4double kgamma2 = kgamma*kgamma; 447 447 448 // G4double dk2t = delta*fWaveVector*fWaveV 448 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta; 449 // G4double dk2t2 = dk2t*dk2t; 449 // G4double dk2t2 = dk2t*dk2t; 450 // G4double pikdt = pi*fWaveVector*diffuse*t 450 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta; 451 451 452 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa 452 G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta; 453 453 454 damp = DampFactor(pikdt); 454 damp = DampFactor(pikdt); 455 damp2 = damp*damp; 455 damp2 = damp*damp; 456 456 457 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector 457 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector; 458 G4double e2dk3t = -2.*e2*delta*fWaveVector* 458 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta; 459 459 460 460 461 sigma = kgamma2; 461 sigma = kgamma2; 462 // sigma += dk2t2; 462 // sigma += dk2t2; 463 sigma *= bzero2; 463 sigma *= bzero2; 464 sigma += mode2k2*bone2 + e2dk3t*bzero*bone; 464 sigma += mode2k2*bone2 + e2dk3t*bzero*bone; 465 sigma += kr2*bonebyarg2; 465 sigma += kr2*bonebyarg2; 466 sigma *= damp2; // *rad*rad; 466 sigma *= damp2; // *rad*rad; 467 467 468 return sigma; 468 return sigma; 469 } 469 } 470 470 471 ////////////////////////////////////////////// 471 //////////////////////////////////////////////////////////////////////////// 472 // 472 // 473 // return differential elastic probability d(p 473 // return differential elastic probability d(probability)/d(omega) with 474 // Coulomb correction 474 // Coulomb correction 475 475 476 G4double 476 G4double 477 G4NuclNuclDiffuseElastic::GetDiffElasticSumPro 477 G4NuclNuclDiffuseElastic::GetDiffElasticSumProb( // G4ParticleDefinition* particle, 478 G4doub 478 G4double theta 479 // G4double momentum, 479 // G4double momentum, 480 // G4double A 480 // G4double A 481 ) 481 ) 482 { 482 { 483 G4double sigma, bzero, bzero2, bonebyarg, bo 483 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2; 484 G4double delta, diffuse, gamma; 484 G4double delta, diffuse, gamma; 485 G4double e1, e2, bone, bone2; 485 G4double e1, e2, bone, bone2; 486 486 487 // G4double wavek = momentum/hbarc; // wave 487 // G4double wavek = momentum/hbarc; // wave vector 488 // G4double r0 = 1.08*fermi; 488 // G4double r0 = 1.08*fermi; 489 // G4double rad = r0*G4Pow::GetInstance()- 489 // G4double rad = r0*G4Pow::GetInstance()->A13(A); 490 490 491 G4double kr = fWaveVector*fNuclearRadius; 491 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad; 492 G4double kr2 = kr*kr; 492 G4double kr2 = kr*kr; 493 G4double krt = kr*theta; 493 G4double krt = kr*theta; 494 494 495 bzero = BesselJzero(krt); 495 bzero = BesselJzero(krt); 496 bzero2 = bzero*bzero; 496 bzero2 = bzero*bzero; 497 bone = BesselJone(krt); 497 bone = BesselJone(krt); 498 bone2 = bone*bone; 498 bone2 = bone*bone; 499 bonebyarg = BesselOneByArg(krt); 499 bonebyarg = BesselOneByArg(krt); 500 bonebyarg2 = bonebyarg*bonebyarg; 500 bonebyarg2 = bonebyarg*bonebyarg; 501 501 502 if (fParticle == theProton) 502 if (fParticle == theProton) 503 { 503 { 504 diffuse = 0.63*fermi; 504 diffuse = 0.63*fermi; 505 // diffuse = 0.6*fermi; 505 // diffuse = 0.6*fermi; 506 gamma = 0.3*fermi; 506 gamma = 0.3*fermi; 507 delta = 0.1*fermi*fermi; 507 delta = 0.1*fermi*fermi; 508 e1 = 0.3*fermi; 508 e1 = 0.3*fermi; 509 e2 = 0.35*fermi; 509 e2 = 0.35*fermi; 510 } 510 } 511 else // as proton, if were not defined 511 else // as proton, if were not defined 512 { 512 { 513 diffuse = 0.63*fermi; 513 diffuse = 0.63*fermi; 514 gamma = 0.3*fermi; 514 gamma = 0.3*fermi; 515 delta = 0.1*fermi*fermi; 515 delta = 0.1*fermi*fermi; 516 e1 = 0.3*fermi; 516 e1 = 0.3*fermi; 517 e2 = 0.35*fermi; 517 e2 = 0.35*fermi; 518 } 518 } 519 G4double lambda = 15.; // 15 ok 519 G4double lambda = 15.; // 15 ok 520 // G4double kgamma = fWaveVector*gamma; 520 // G4double kgamma = fWaveVector*gamma; // wavek*delta; 521 G4double kgamma = lambda*(1.-G4Exp(-fWave 521 G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta; 522 522 523 // G4cout<<"kgamma = "<<kgamma<<G4endl; 523 // G4cout<<"kgamma = "<<kgamma<<G4endl; 524 524 525 if(fAddCoulomb) // add Coulomb correction 525 if(fAddCoulomb) // add Coulomb correction 526 { 526 { 527 G4double sinHalfTheta = std::sin(0.5*thet 527 G4double sinHalfTheta = std::sin(0.5*theta); 528 G4double sinHalfTheta2 = sinHalfTheta*sinH 528 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta; 529 529 530 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta 530 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0() 531 // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe 531 // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0() 532 } 532 } 533 533 534 G4double kgamma2 = kgamma*kgamma; 534 G4double kgamma2 = kgamma*kgamma; 535 535 536 // G4double dk2t = delta*fWaveVector*fWaveV 536 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta; 537 // G4cout<<"dk2t = "<<dk2t<<G4endl; 537 // G4cout<<"dk2t = "<<dk2t<<G4endl; 538 // G4double dk2t2 = dk2t*dk2t; 538 // G4double dk2t2 = dk2t*dk2t; 539 // G4double pikdt = pi*fWaveVector*diffuse*t 539 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta; 540 540 541 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa 541 G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta; 542 542 543 // G4cout<<"pikdt = "<<pikdt<<G4endl; 543 // G4cout<<"pikdt = "<<pikdt<<G4endl; 544 544 545 damp = DampFactor(pikdt); 545 damp = DampFactor(pikdt); 546 damp2 = damp*damp; 546 damp2 = damp*damp; 547 547 548 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector 548 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector; 549 G4double e2dk3t = -2.*e2*delta*fWaveVector* 549 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta; 550 550 551 sigma = kgamma2; 551 sigma = kgamma2; 552 // sigma += dk2t2; 552 // sigma += dk2t2; 553 sigma *= bzero2; 553 sigma *= bzero2; 554 sigma += mode2k2*bone2; 554 sigma += mode2k2*bone2; 555 sigma += e2dk3t*bzero*bone; 555 sigma += e2dk3t*bzero*bone; 556 556 557 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerf 557 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/() 558 sigma += kr2*bonebyarg2; // correction at J 558 sigma += kr2*bonebyarg2; // correction at J1()/() 559 559 560 sigma *= damp2; // *rad*rad; 560 sigma *= damp2; // *rad*rad; 561 561 562 return sigma; 562 return sigma; 563 } 563 } 564 564 565 565 566 ////////////////////////////////////////////// 566 //////////////////////////////////////////////////////////////////////////// 567 // 567 // 568 // return differential elastic probability d(p 568 // return differential elastic probability d(probability)/d(t) with 569 // Coulomb correction 569 // Coulomb correction 570 570 571 G4double 571 G4double 572 G4NuclNuclDiffuseElastic::GetDiffElasticSumPro 572 G4NuclNuclDiffuseElastic::GetDiffElasticSumProbA( G4double alpha ) 573 { 573 { 574 G4double theta; 574 G4double theta; 575 575 576 theta = std::sqrt(alpha); 576 theta = std::sqrt(alpha); 577 577 578 // theta = std::acos( 1 - alpha/2. ); 578 // theta = std::acos( 1 - alpha/2. ); 579 579 580 G4double sigma, bzero, bzero2, bonebyarg, bo 580 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2; 581 G4double delta, diffuse, gamma; 581 G4double delta, diffuse, gamma; 582 G4double e1, e2, bone, bone2; 582 G4double e1, e2, bone, bone2; 583 583 584 // G4double wavek = momentum/hbarc; // wave 584 // G4double wavek = momentum/hbarc; // wave vector 585 // G4double r0 = 1.08*fermi; 585 // G4double r0 = 1.08*fermi; 586 // G4double rad = r0*G4Pow::GetInstance()- 586 // G4double rad = r0*G4Pow::GetInstance()->A13(A); 587 587 588 G4double kr = fWaveVector*fNuclearRadius; 588 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad; 589 G4double kr2 = kr*kr; 589 G4double kr2 = kr*kr; 590 G4double krt = kr*theta; 590 G4double krt = kr*theta; 591 591 592 bzero = BesselJzero(krt); 592 bzero = BesselJzero(krt); 593 bzero2 = bzero*bzero; 593 bzero2 = bzero*bzero; 594 bone = BesselJone(krt); 594 bone = BesselJone(krt); 595 bone2 = bone*bone; 595 bone2 = bone*bone; 596 bonebyarg = BesselOneByArg(krt); 596 bonebyarg = BesselOneByArg(krt); 597 bonebyarg2 = bonebyarg*bonebyarg; 597 bonebyarg2 = bonebyarg*bonebyarg; 598 598 599 if (fParticle == theProton) 599 if (fParticle == theProton) 600 { 600 { 601 diffuse = 0.63*fermi; 601 diffuse = 0.63*fermi; 602 // diffuse = 0.6*fermi; 602 // diffuse = 0.6*fermi; 603 gamma = 0.3*fermi; 603 gamma = 0.3*fermi; 604 delta = 0.1*fermi*fermi; 604 delta = 0.1*fermi*fermi; 605 e1 = 0.3*fermi; 605 e1 = 0.3*fermi; 606 e2 = 0.35*fermi; 606 e2 = 0.35*fermi; 607 } 607 } 608 else // as proton, if were not defined 608 else // as proton, if were not defined 609 { 609 { 610 diffuse = 0.63*fermi; 610 diffuse = 0.63*fermi; 611 gamma = 0.3*fermi; 611 gamma = 0.3*fermi; 612 delta = 0.1*fermi*fermi; 612 delta = 0.1*fermi*fermi; 613 e1 = 0.3*fermi; 613 e1 = 0.3*fermi; 614 e2 = 0.35*fermi; 614 e2 = 0.35*fermi; 615 } 615 } 616 G4double lambda = 15.; // 15 ok 616 G4double lambda = 15.; // 15 ok 617 // G4double kgamma = fWaveVector*gamma; 617 // G4double kgamma = fWaveVector*gamma; // wavek*delta; 618 G4double kgamma = lambda*(1.-G4Exp(-fWave 618 G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta; 619 619 620 // G4cout<<"kgamma = "<<kgamma<<G4endl; 620 // G4cout<<"kgamma = "<<kgamma<<G4endl; 621 621 622 if(fAddCoulomb) // add Coulomb correction 622 if(fAddCoulomb) // add Coulomb correction 623 { 623 { 624 G4double sinHalfTheta = theta*0.5; // std 624 G4double sinHalfTheta = theta*0.5; // std::sin(0.5*theta); 625 G4double sinHalfTheta2 = sinHalfTheta*sinH 625 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta; 626 626 627 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta 627 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0() 628 // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe 628 // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0() 629 } 629 } 630 630 631 G4double kgamma2 = kgamma*kgamma; 631 G4double kgamma2 = kgamma*kgamma; 632 632 633 // G4double dk2t = delta*fWaveVector*fWaveV 633 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta; 634 // G4cout<<"dk2t = "<<dk2t<<G4endl; 634 // G4cout<<"dk2t = "<<dk2t<<G4endl; 635 // G4double dk2t2 = dk2t*dk2t; 635 // G4double dk2t2 = dk2t*dk2t; 636 // G4double pikdt = pi*fWaveVector*diffuse*t 636 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta; 637 637 638 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa 638 G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta; 639 639 640 // G4cout<<"pikdt = "<<pikdt<<G4endl; 640 // G4cout<<"pikdt = "<<pikdt<<G4endl; 641 641 642 damp = DampFactor(pikdt); 642 damp = DampFactor(pikdt); 643 damp2 = damp*damp; 643 damp2 = damp*damp; 644 644 645 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector 645 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector; 646 G4double e2dk3t = -2.*e2*delta*fWaveVector* 646 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta; 647 647 648 sigma = kgamma2; 648 sigma = kgamma2; 649 // sigma += dk2t2; 649 // sigma += dk2t2; 650 sigma *= bzero2; 650 sigma *= bzero2; 651 sigma += mode2k2*bone2; 651 sigma += mode2k2*bone2; 652 sigma += e2dk3t*bzero*bone; 652 sigma += e2dk3t*bzero*bone; 653 653 654 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerf 654 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/() 655 sigma += kr2*bonebyarg2; // correction at J 655 sigma += kr2*bonebyarg2; // correction at J1()/() 656 656 657 sigma *= damp2; // *rad*rad; 657 sigma *= damp2; // *rad*rad; 658 658 659 return sigma; 659 return sigma; 660 } 660 } 661 661 662 662 663 ////////////////////////////////////////////// 663 //////////////////////////////////////////////////////////////////////////// 664 // 664 // 665 // return differential elastic probability 2*p 665 // return differential elastic probability 2*pi*sin(theta)*d(probability)/d(omega) 666 666 667 G4double 667 G4double 668 G4NuclNuclDiffuseElastic::GetIntegrandFunction 668 G4NuclNuclDiffuseElastic::GetIntegrandFunction( G4double alpha ) 669 { 669 { 670 G4double result; 670 G4double result; 671 671 672 result = GetDiffElasticSumProbA(alpha); 672 result = GetDiffElasticSumProbA(alpha); 673 673 674 // result *= 2*pi*std::sin(theta); 674 // result *= 2*pi*std::sin(theta); 675 675 676 return result; 676 return result; 677 } 677 } 678 678 679 ////////////////////////////////////////////// 679 //////////////////////////////////////////////////////////////////////////// 680 // 680 // 681 // return integral elastic cross section d(sig 681 // return integral elastic cross section d(sigma)/d(omega) integrated 0 - theta 682 682 683 G4double 683 G4double 684 G4NuclNuclDiffuseElastic::IntegralElasticProb( 684 G4NuclNuclDiffuseElastic::IntegralElasticProb( const G4ParticleDefinition* particle, 685 G4doub 685 G4double theta, 686 G4double momentum, 686 G4double momentum, 687 G4doub 687 G4double A ) 688 { 688 { 689 G4double result; 689 G4double result; 690 fParticle = particle; 690 fParticle = particle; 691 fWaveVector = momentum/hbarc; 691 fWaveVector = momentum/hbarc; 692 fAtomicWeight = A; 692 fAtomicWeight = A; 693 693 694 fNuclearRadius = CalculateNuclearRad(A); 694 fNuclearRadius = CalculateNuclearRad(A); 695 695 696 696 697 G4Integrator<G4NuclNuclDiffuseElastic,G4doub 697 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral; 698 698 699 // result = integral.Legendre10(this,&G4Nucl 699 // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta ); 700 result = integral.Legendre96(this,&G4NuclNuc 700 result = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta ); 701 701 702 return result; 702 return result; 703 } 703 } 704 704 705 ////////////////////////////////////////////// 705 //////////////////////////////////////////////////////////////////////////// 706 // 706 // 707 // Return inv momentum transfer -t > 0 707 // Return inv momentum transfer -t > 0 708 708 709 G4double G4NuclNuclDiffuseElastic::SampleT( co 709 G4double G4NuclNuclDiffuseElastic::SampleT( const G4ParticleDefinition* aParticle, 710 G4double p, G4double A) 710 G4double p, G4double A) 711 { 711 { 712 G4double theta = SampleThetaCMS( aParticle, 712 G4double theta = SampleThetaCMS( aParticle, p, A); // sample theta in cms 713 G4double t = 2*p*p*( 1 - std::cos(theta) 713 G4double t = 2*p*p*( 1 - std::cos(theta) ); // -t !!! 714 return t; 714 return t; 715 } 715 } 716 716 717 ////////////////////////////////////////////// 717 //////////////////////////////////////////////////////////////////////////// 718 // 718 // 719 // Return scattering angle sampled in cms 719 // Return scattering angle sampled in cms 720 720 721 721 722 G4double 722 G4double 723 G4NuclNuclDiffuseElastic::SampleThetaCMS(const 723 G4NuclNuclDiffuseElastic::SampleThetaCMS(const G4ParticleDefinition* particle, 724 G4doubl 724 G4double momentum, G4double A) 725 { 725 { 726 G4int i, iMax = 100; 726 G4int i, iMax = 100; 727 G4double norm, theta1, theta2, thetaMax; << 727 G4double norm, result, theta1, theta2, thetaMax, sum = 0.; 728 G4double result = 0., sum = 0.; << 729 728 730 fParticle = particle; 729 fParticle = particle; 731 fWaveVector = momentum/hbarc; 730 fWaveVector = momentum/hbarc; 732 fAtomicWeight = A; 731 fAtomicWeight = A; 733 732 734 fNuclearRadius = CalculateNuclearRad(A); 733 fNuclearRadius = CalculateNuclearRad(A); 735 734 736 thetaMax = 10.174/fWaveVector/fNuclearRadius 735 thetaMax = 10.174/fWaveVector/fNuclearRadius; 737 736 738 if (thetaMax > pi) thetaMax = pi; 737 if (thetaMax > pi) thetaMax = pi; 739 738 740 G4Integrator<G4NuclNuclDiffuseElastic,G4doub 739 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral; 741 740 742 // result = integral.Legendre10(this,&G4Nucl 741 // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta ); 743 norm = integral.Legendre96(this,&G4NuclNuclD 742 norm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., thetaMax ); 744 743 745 norm *= G4UniformRand(); 744 norm *= G4UniformRand(); 746 745 747 for(i = 1; i <= iMax; i++) 746 for(i = 1; i <= iMax; i++) 748 { 747 { 749 theta1 = (i-1)*thetaMax/iMax; 748 theta1 = (i-1)*thetaMax/iMax; 750 theta2 = i*thetaMax/iMax; 749 theta2 = i*thetaMax/iMax; 751 sum += integral.Legendre10(this,&G4NuclN 750 sum += integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, theta1, theta2); 752 751 753 if ( sum >= norm ) 752 if ( sum >= norm ) 754 { 753 { 755 result = 0.5*(theta1 + theta2); 754 result = 0.5*(theta1 + theta2); 756 break; 755 break; 757 } 756 } 758 } 757 } 759 if (i > iMax ) result = 0.5*(theta1 + theta2 758 if (i > iMax ) result = 0.5*(theta1 + theta2); 760 759 761 G4double sigma = pi*thetaMax/iMax; 760 G4double sigma = pi*thetaMax/iMax; 762 761 763 result += G4RandGauss::shoot(0.,sigma); 762 result += G4RandGauss::shoot(0.,sigma); 764 763 765 if(result < 0.) result = 0.; 764 if(result < 0.) result = 0.; 766 if(result > thetaMax) result = thetaMax; 765 if(result > thetaMax) result = thetaMax; 767 766 768 return result; 767 return result; 769 } 768 } 770 769 771 ////////////////////////////////////////////// 770 ///////////////////////////////////////////////////////////////////////////// 772 ///////////////////// Table preparation and r 771 ///////////////////// Table preparation and reading //////////////////////// 773 772 774 ////////////////////////////////////////////// 773 //////////////////////////////////////////////////////////////////////////// 775 // 774 // 776 // Interface function. Return inv momentum tra 775 // Interface function. Return inv momentum transfer -t > 0 from initialisation table 777 776 778 G4double G4NuclNuclDiffuseElastic::SampleInvar 777 G4double G4NuclNuclDiffuseElastic::SampleInvariantT( const G4ParticleDefinition* aParticle, G4double p, 779 778 G4int Z, G4int A) 780 { 779 { 781 fParticle = aParticle; 780 fParticle = aParticle; 782 fAtomicNumber = G4double(Z); 781 fAtomicNumber = G4double(Z); 783 fAtomicWeight = G4double(A); 782 fAtomicWeight = G4double(A); 784 783 785 G4double t(0.), m1 = fParticle->GetPDGMass() 784 G4double t(0.), m1 = fParticle->GetPDGMass(); 786 G4double totElab = std::sqrt(m1*m1+p*p); 785 G4double totElab = std::sqrt(m1*m1+p*p); 787 G4double mass2 = G4NucleiProperties::GetNucl 786 G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z); 788 G4LorentzVector lv1(p,0.0,0.0,totElab); 787 G4LorentzVector lv1(p,0.0,0.0,totElab); 789 G4LorentzVector lv(0.0,0.0,0.0,mass2); 788 G4LorentzVector lv(0.0,0.0,0.0,mass2); 790 lv += lv1; 789 lv += lv1; 791 790 792 G4ThreeVector bst = lv.boostVector(); 791 G4ThreeVector bst = lv.boostVector(); 793 lv1.boost(-bst); 792 lv1.boost(-bst); 794 793 795 G4ThreeVector p1 = lv1.vect(); 794 G4ThreeVector p1 = lv1.vect(); 796 G4double momentumCMS = p1.mag(); 795 G4double momentumCMS = p1.mag(); 797 796 798 // t = SampleTableT( aParticle, momentumCMS 797 // t = SampleTableT( aParticle, momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms 799 798 800 t = SampleCoulombMuCMS( aParticle, momentum 799 t = SampleCoulombMuCMS( aParticle, momentumCMS); // sample theta2 in cms 801 800 802 801 803 return t; 802 return t; 804 } 803 } 805 804 806 ////////////////////////////////////////////// 805 //////////////////////////////////////////////////////////////////////////// 807 // 806 // 808 // Return inv momentum transfer -t > 0 as Coul 807 // Return inv momentum transfer -t > 0 as Coulomb scattering <= thetaC 809 808 810 G4double G4NuclNuclDiffuseElastic::SampleCoulo 809 G4double G4NuclNuclDiffuseElastic::SampleCoulombMuCMS( const G4ParticleDefinition* aParticle, G4double p) 811 { 810 { 812 G4double t(0.), rand(0.), mu(0.); 811 G4double t(0.), rand(0.), mu(0.); 813 812 814 G4double A1 = G4double( aParticle->GetBaryon 813 G4double A1 = G4double( aParticle->GetBaryonNumber() ); 815 G4double R1 = CalculateNuclearRad(A1); 814 G4double R1 = CalculateNuclearRad(A1); 816 815 817 fNuclearRadius = CalculateNuclearRad(fAtomi 816 fNuclearRadius = CalculateNuclearRad(fAtomicWeight); 818 fNuclearRadius += R1; 817 fNuclearRadius += R1; 819 818 820 InitDynParameters(fParticle, p); 819 InitDynParameters(fParticle, p); 821 820 822 fCoulombMuC = fHalfRutThetaTg2/(1.+fHalfRutT 821 fCoulombMuC = fHalfRutThetaTg2/(1.+fHalfRutThetaTg2); 823 822 824 rand = G4UniformRand(); 823 rand = G4UniformRand(); 825 824 826 // sample (1-cosTheta) in 0 < Theta < ThetaC 825 // sample (1-cosTheta) in 0 < Theta < ThetaC as Coulomb scattering 827 826 828 mu = fCoulombMuC*rand*fAm; 827 mu = fCoulombMuC*rand*fAm; 829 mu /= fAm + fCoulombMuC*( 1. - rand ); 828 mu /= fAm + fCoulombMuC*( 1. - rand ); 830 829 831 t = 4.*p*p*mu; 830 t = 4.*p*p*mu; 832 831 833 return t; 832 return t; 834 } 833 } 835 834 836 835 837 ////////////////////////////////////////////// 836 //////////////////////////////////////////////////////////////////////////// 838 // 837 // 839 // Return inv momentum transfer -t > 0 from in 838 // Return inv momentum transfer -t > 0 from initialisation table 840 839 841 G4double G4NuclNuclDiffuseElastic::SampleTable 840 G4double G4NuclNuclDiffuseElastic::SampleTableT( const G4ParticleDefinition* aParticle, G4double p, 842 841 G4double Z, G4double A) 843 { 842 { 844 G4double alpha = SampleTableThetaCMS( aParti 843 G4double alpha = SampleTableThetaCMS( aParticle, p, Z, A); // sample theta2 in cms 845 // G4double t = 2*p*p*( 1 - std::cos(std 844 // G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) ); // -t !!! 846 G4double t = p*p*alpha; // - 845 G4double t = p*p*alpha; // -t !!! 847 return t; 846 return t; 848 } 847 } 849 848 850 ////////////////////////////////////////////// 849 //////////////////////////////////////////////////////////////////////////// 851 // 850 // 852 // Return scattering angle2 sampled in cms acc 851 // Return scattering angle2 sampled in cms according to precalculated table. 853 852 854 853 855 G4double 854 G4double 856 G4NuclNuclDiffuseElastic::SampleTableThetaCMS( 855 G4NuclNuclDiffuseElastic::SampleTableThetaCMS(const G4ParticleDefinition* particle, 857 G4doubl 856 G4double momentum, G4double Z, G4double A) 858 { 857 { 859 std::size_t iElement; << 858 size_t iElement; 860 G4int iMomentum, iAngle; 859 G4int iMomentum, iAngle; 861 G4double randAngle, position, theta1, theta2 860 G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W; 862 G4double m1 = particle->GetPDGMass(); 861 G4double m1 = particle->GetPDGMass(); 863 862 864 for(iElement = 0; iElement < fElementNumberV 863 for(iElement = 0; iElement < fElementNumberVector.size(); iElement++) 865 { 864 { 866 if( std::fabs(Z - fElementNumberVector[iEl 865 if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break; 867 } 866 } 868 if ( iElement == fElementNumberVector.size() 867 if ( iElement == fElementNumberVector.size() ) 869 { 868 { 870 InitialiseOnFly(Z,A); // table preparation 869 InitialiseOnFly(Z,A); // table preparation, if needed 871 870 872 // iElement--; 871 // iElement--; 873 872 874 // G4cout << "G4NuclNuclDiffuseElastic: El 873 // G4cout << "G4NuclNuclDiffuseElastic: Element with atomic number " << Z 875 // << " is not found, return zero angle" < 874 // << " is not found, return zero angle" << G4endl; 876 // return 0.; // no table for this element 875 // return 0.; // no table for this element 877 } 876 } 878 // G4cout<<"iElement = "<<iElement<<G4endl; 877 // G4cout<<"iElement = "<<iElement<<G4endl; 879 878 880 fAngleTable = fAngleBank[iElement]; 879 fAngleTable = fAngleBank[iElement]; 881 880 882 G4double kinE = std::sqrt(momentum*momentum 881 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1; 883 882 884 for( iMomentum = 0; iMomentum < fEnergyBin; 883 for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++) 885 { 884 { 886 // G4cout<<iMomentum<<", kinE = "<<kinE/Me 885 // G4cout<<iMomentum<<", kinE = "<<kinE/MeV<<", vectorE = "<<fEnergyVector->GetLowEdgeEnergy(iMomentum)/MeV<<G4endl; 887 886 888 if( kinE < fEnergyVector->GetLowEdgeEnergy 887 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break; 889 } 888 } 890 // G4cout<<"iMomentum = "<<iMomentum<<", fEn 889 // G4cout<<"iMomentum = "<<iMomentum<<", fEnergyBin -1 = "<<fEnergyBin -1<<G4endl; 891 890 892 891 893 if ( iMomentum >= fEnergyBin ) iMomentum = f 892 if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1; // kinE is more then theMaxEnergy 894 if ( iMomentum < 0 ) iMomentum = 0 893 if ( iMomentum < 0 ) iMomentum = 0; // against negative index, kinE < theMinEnergy 895 894 896 895 897 if (iMomentum == fEnergyBin -1 || iMomentum 896 if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges 898 { 897 { 899 position = (*(*fAngleTable)(iMomentum))(fA 898 position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand(); 900 899 901 // G4cout<<"position = "<<position<<G4endl 900 // G4cout<<"position = "<<position<<G4endl; 902 901 903 for(iAngle = 0; iAngle < fAngleBin-1; iAng 902 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++) 904 { 903 { 905 if( position < (*(*fAngleTable)(iMomentu 904 if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break; 906 } 905 } 907 if (iAngle >= fAngleBin-1) iAngle = fAngle 906 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2; 908 907 909 // G4cout<<"iAngle = "<<iAngle<<G4endl; 908 // G4cout<<"iAngle = "<<iAngle<<G4endl; 910 909 911 randAngle = GetScatteringAngle(iMomentum, 910 randAngle = GetScatteringAngle(iMomentum, iAngle, position); 912 911 913 // G4cout<<"randAngle = "<<randAngle<<G4en 912 // G4cout<<"randAngle = "<<randAngle<<G4endl; 914 } 913 } 915 else // kinE inside between energy table ed 914 else // kinE inside between energy table edges 916 { 915 { 917 // position = (*(*fAngleTable)(iMomentum)) 916 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand(); 918 position = (*(*fAngleTable)(iMomentum))(0) 917 position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand(); 919 918 920 // G4cout<<"position = "<<position<<G4endl 919 // G4cout<<"position = "<<position<<G4endl; 921 920 922 for(iAngle = 0; iAngle < fAngleBin-1; iAng 921 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++) 923 { 922 { 924 // if( position < (*(*fAngleTable)(iMome 923 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break; 925 if( position > (*(*fAngleTable)(iMomentu 924 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break; 926 } 925 } 927 if (iAngle >= fAngleBin-1) iAngle = fAngle 926 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2; 928 927 929 // G4cout<<"iAngle = "<<iAngle<<G4endl; 928 // G4cout<<"iAngle = "<<iAngle<<G4endl; 930 929 931 theta2 = GetScatteringAngle(iMomentum, iA 930 theta2 = GetScatteringAngle(iMomentum, iAngle, position); 932 931 933 // G4cout<<"theta2 = "<<theta2<<G4endl; 932 // G4cout<<"theta2 = "<<theta2<<G4endl; 934 933 935 E2 = fEnergyVector->GetLowEdgeEnergy(iMome 934 E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum); 936 935 937 // G4cout<<"E2 = "<<E2<<G4endl; 936 // G4cout<<"E2 = "<<E2<<G4endl; 938 937 939 iMomentum--; 938 iMomentum--; 940 939 941 // position = (*(*fAngleTable)(iMomentum)) 940 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand(); 942 941 943 // G4cout<<"position = "<<position<<G4endl 942 // G4cout<<"position = "<<position<<G4endl; 944 943 945 for(iAngle = 0; iAngle < fAngleBin-1; iAng 944 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++) 946 { 945 { 947 // if( position < (*(*fAngleTable)(iMome 946 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break; 948 if( position > (*(*fAngleTable)(iMomentu 947 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break; 949 } 948 } 950 if (iAngle >= fAngleBin-1) iAngle = fAngle 949 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2; 951 950 952 theta1 = GetScatteringAngle(iMomentum, iA 951 theta1 = GetScatteringAngle(iMomentum, iAngle, position); 953 952 954 // G4cout<<"theta1 = "<<theta1<<G4endl; 953 // G4cout<<"theta1 = "<<theta1<<G4endl; 955 954 956 E1 = fEnergyVector->GetLowEdgeEnergy(iMome 955 E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum); 957 956 958 // G4cout<<"E1 = "<<E1<<G4endl; 957 // G4cout<<"E1 = "<<E1<<G4endl; 959 958 960 W = 1.0/(E2 - E1); 959 W = 1.0/(E2 - E1); 961 W1 = (E2 - kinE)*W; 960 W1 = (E2 - kinE)*W; 962 W2 = (kinE - E1)*W; 961 W2 = (kinE - E1)*W; 963 962 964 randAngle = W1*theta1 + W2*theta2; 963 randAngle = W1*theta1 + W2*theta2; 965 964 966 // randAngle = theta2; 965 // randAngle = theta2; 967 } 966 } 968 // G4double angle = randAngle; 967 // G4double angle = randAngle; 969 // if (randAngle > 0.) randAngle /= 2*pi*std 968 // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle); 970 // G4cout<<"randAngle = "<<randAngle/degree< 969 // G4cout<<"randAngle = "<<randAngle/degree<<G4endl; 971 970 972 return randAngle; 971 return randAngle; 973 } 972 } 974 973 975 ////////////////////////////////////////////// 974 ////////////////////////////////////////////////////////////////////////////// 976 // 975 // 977 // Initialisation for given particle on fly us 976 // Initialisation for given particle on fly using new element number 978 977 979 void G4NuclNuclDiffuseElastic::InitialiseOnFly 978 void G4NuclNuclDiffuseElastic::InitialiseOnFly(G4double Z, G4double A) 980 { 979 { 981 fAtomicNumber = Z; // atomic number 980 fAtomicNumber = Z; // atomic number 982 fAtomicWeight = G4NistManager::Instance()-> 981 fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( Z ) ); 983 982 984 G4double A1 = G4double( fParticle->GetBaryon 983 G4double A1 = G4double( fParticle->GetBaryonNumber() ); 985 G4double R1 = CalculateNuclearRad(A1); 984 G4double R1 = CalculateNuclearRad(A1); 986 985 987 fNuclearRadius = CalculateNuclearRad(fAtomic 986 fNuclearRadius = CalculateNuclearRad(fAtomicWeight) + R1; 988 987 989 if( verboseLevel > 0 ) 988 if( verboseLevel > 0 ) 990 { 989 { 991 G4cout<<"G4NuclNuclDiffuseElastic::Initial 990 G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element with Z = " 992 <<Z<<"; and A = "<<A<<G4endl; 991 <<Z<<"; and A = "<<A<<G4endl; 993 } 992 } 994 fElementNumberVector.push_back(fAtomicNumber 993 fElementNumberVector.push_back(fAtomicNumber); 995 994 996 BuildAngleTable(); 995 BuildAngleTable(); 997 996 998 fAngleBank.push_back(fAngleTable); 997 fAngleBank.push_back(fAngleTable); 999 998 1000 return; 999 return; 1001 } 1000 } 1002 1001 1003 ///////////////////////////////////////////// 1002 /////////////////////////////////////////////////////////////////////////////// 1004 // 1003 // 1005 // Build for given particle and element table 1004 // Build for given particle and element table of momentum, angle probability. 1006 // For the moment in lab system. 1005 // For the moment in lab system. 1007 1006 1008 void G4NuclNuclDiffuseElastic::BuildAngleTabl 1007 void G4NuclNuclDiffuseElastic::BuildAngleTable() 1009 { 1008 { 1010 G4int i, j; 1009 G4int i, j; 1011 G4double partMom, kinE, m1 = fParticle->Get 1010 G4double partMom, kinE, m1 = fParticle->GetPDGMass(); 1012 G4double alpha1, alpha2, alphaMax, alphaCou 1011 G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.; 1013 1012 1014 // G4cout<<"particle z = "<<z<<"; particle 1013 // G4cout<<"particle z = "<<z<<"; particle m1 = "<<m1/GeV<<" GeV"<<G4endl; 1015 1014 1016 G4Integrator<G4NuclNuclDiffuseElastic,G4dou 1015 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral; 1017 1016 1018 fAngleTable = new G4PhysicsTable(fEnergyBin 1017 fAngleTable = new G4PhysicsTable(fEnergyBin); 1019 1018 1020 for( i = 0; i < fEnergyBin; i++) 1019 for( i = 0; i < fEnergyBin; i++) 1021 { 1020 { 1022 kinE = fEnergyVector->GetLowEdgeEn 1021 kinE = fEnergyVector->GetLowEdgeEnergy(i); 1023 1022 1024 // G4cout<<G4endl; 1023 // G4cout<<G4endl; 1025 // G4cout<<"kinE = "<<kinE/MeV<<" MeV"<<G 1024 // G4cout<<"kinE = "<<kinE/MeV<<" MeV"<<G4endl; 1026 1025 1027 partMom = std::sqrt( kinE*(kinE + 2*m 1026 partMom = std::sqrt( kinE*(kinE + 2*m1) ); 1028 1027 1029 InitDynParameters(fParticle, partMom); 1028 InitDynParameters(fParticle, partMom); 1030 1029 1031 alphaMax = fRutherfordTheta*fCofAlphaMax; 1030 alphaMax = fRutherfordTheta*fCofAlphaMax; 1032 1031 1033 if(alphaMax > pi) alphaMax = pi; 1032 if(alphaMax > pi) alphaMax = pi; 1034 1033 1035 // VI: Coverity complain 1034 // VI: Coverity complain 1036 //alphaMax = pi2; 1035 //alphaMax = pi2; 1037 1036 1038 alphaCoulomb = fRutherfordTheta*fCofAlpha 1037 alphaCoulomb = fRutherfordTheta*fCofAlphaCoulomb; 1039 1038 1040 // G4cout<<"alphaCoulomb = "<<alphaCoulom 1039 // G4cout<<"alphaCoulomb = "<<alphaCoulomb/degree<<"; alphaMax = "<<alphaMax/degree<<G4endl; 1041 1040 1042 G4PhysicsFreeVector* angleVector = new G4 1041 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1); 1043 1042 1044 // G4PhysicsLogVector* angleBins = new G 1043 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin ); 1045 1044 1046 G4double delth = (alphaMax-alphaCoulomb)/ 1045 G4double delth = (alphaMax-alphaCoulomb)/fAngleBin; 1047 1046 1048 sum = 0.; 1047 sum = 0.; 1049 1048 1050 // fAddCoulomb = false; 1049 // fAddCoulomb = false; 1051 fAddCoulomb = true; 1050 fAddCoulomb = true; 1052 1051 1053 // for(j = 1; j < fAngleBin; j++) 1052 // for(j = 1; j < fAngleBin; j++) 1054 for(j = fAngleBin-1; j >= 1; j--) 1053 for(j = fAngleBin-1; j >= 1; j--) 1055 { 1054 { 1056 // alpha1 = angleBins->GetLowEdgeEnergy 1055 // alpha1 = angleBins->GetLowEdgeEnergy(j-1); 1057 // alpha2 = angleBins->GetLowEdgeEnergy 1056 // alpha2 = angleBins->GetLowEdgeEnergy(j); 1058 1057 1059 // alpha1 = alphaMax*(j-1)/fAngleBin; 1058 // alpha1 = alphaMax*(j-1)/fAngleBin; 1060 // alpha2 = alphaMax*( j )/fAngleBin; 1059 // alpha2 = alphaMax*( j )/fAngleBin; 1061 1060 1062 alpha1 = alphaCoulomb + delth*(j-1); 1061 alpha1 = alphaCoulomb + delth*(j-1); 1063 // if(alpha1 < kRlim2) alpha1 = kRlim2; 1062 // if(alpha1 < kRlim2) alpha1 = kRlim2; 1064 alpha2 = alpha1 + delth; 1063 alpha2 = alpha1 + delth; 1065 1064 1066 delta = integral.Legendre10(this, &G4Nu 1065 delta = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc, alpha1, alpha2); 1067 // delta = integral.Legendre96(this, &G 1066 // delta = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2); 1068 1067 1069 sum += delta; 1068 sum += delta; 1070 1069 1071 angleVector->PutValue( j-1 , alpha1, su 1070 angleVector->PutValue( j-1 , alpha1, sum ); // alpha2 1072 // G4cout<<"j-1 = "<<j-1<<"; alpha2 = " 1071 // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2/degree<<"; sum = "<<sum<<G4endl; 1073 } 1072 } 1074 fAngleTable->insertAt(i,angleVector); 1073 fAngleTable->insertAt(i,angleVector); 1075 1074 1076 // delete[] angleVector; 1075 // delete[] angleVector; 1077 // delete[] angleBins; 1076 // delete[] angleBins; 1078 } 1077 } 1079 return; 1078 return; 1080 } 1079 } 1081 1080 1082 ///////////////////////////////////////////// 1081 ///////////////////////////////////////////////////////////////////////////////// 1083 // 1082 // 1084 // 1083 // 1085 1084 1086 G4double 1085 G4double 1087 G4NuclNuclDiffuseElastic:: GetScatteringAngle 1086 G4NuclNuclDiffuseElastic:: GetScatteringAngle( G4int iMomentum, G4int iAngle, G4double position ) 1088 { 1087 { 1089 G4double x1, x2, y1, y2, randAngle; 1088 G4double x1, x2, y1, y2, randAngle; 1090 1089 1091 if( iAngle == 0 ) 1090 if( iAngle == 0 ) 1092 { 1091 { 1093 randAngle = (*fAngleTable)(iMomentum)->Ge 1092 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle); 1094 // iAngle++; 1093 // iAngle++; 1095 } 1094 } 1096 else 1095 else 1097 { 1096 { 1098 if ( iAngle >= G4int((*fAngleTable)(iMome 1097 if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) ) 1099 { 1098 { 1100 iAngle = G4int((*fAngleTable)(iMomentum << 1099 iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1; 1101 } 1100 } 1102 y1 = (*(*fAngleTable)(iMomentum))(iAngle- 1101 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1); 1103 y2 = (*(*fAngleTable)(iMomentum))(iAngle) 1102 y2 = (*(*fAngleTable)(iMomentum))(iAngle); 1104 1103 1105 x1 = (*fAngleTable)(iMomentum)->GetLowEdg 1104 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1); 1106 x2 = (*fAngleTable)(iMomentum)->GetLowEdg 1105 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle); 1107 1106 1108 if ( x1 == x2 ) randAngle = x2; 1107 if ( x1 == x2 ) randAngle = x2; 1109 else 1108 else 1110 { 1109 { 1111 if ( y1 == y2 ) randAngle = x1 + ( x2 - 1110 if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand(); 1112 else 1111 else 1113 { 1112 { 1114 randAngle = x1 + ( position - y1 )*( 1113 randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 ); 1115 } 1114 } 1116 } 1115 } 1117 } 1116 } 1118 return randAngle; 1117 return randAngle; 1119 } 1118 } 1120 1119 1121 1120 1122 1121 1123 ///////////////////////////////////////////// 1122 //////////////////////////////////////////////////////////////////////////// 1124 // 1123 // 1125 // Return scattering angle sampled in lab sys 1124 // Return scattering angle sampled in lab system (target at rest) 1126 1125 1127 1126 1128 1127 1129 G4double 1128 G4double 1130 G4NuclNuclDiffuseElastic::SampleThetaLab( con 1129 G4NuclNuclDiffuseElastic::SampleThetaLab( const G4HadProjectile* aParticle, 1131 G4dou 1130 G4double tmass, G4double A) 1132 { 1131 { 1133 const G4ParticleDefinition* theParticle = a 1132 const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); 1134 G4double m1 = theParticle->GetPDGMass(); 1133 G4double m1 = theParticle->GetPDGMass(); 1135 G4double plab = aParticle->GetTotalMomentum 1134 G4double plab = aParticle->GetTotalMomentum(); 1136 G4LorentzVector lv1 = aParticle->Get4Moment 1135 G4LorentzVector lv1 = aParticle->Get4Momentum(); 1137 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1136 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1138 lv += lv1; 1137 lv += lv1; 1139 1138 1140 G4ThreeVector bst = lv.boostVector(); 1139 G4ThreeVector bst = lv.boostVector(); 1141 lv1.boost(-bst); 1140 lv1.boost(-bst); 1142 1141 1143 G4ThreeVector p1 = lv1.vect(); 1142 G4ThreeVector p1 = lv1.vect(); 1144 G4double ptot = p1.mag(); 1143 G4double ptot = p1.mag(); 1145 G4double tmax = 4.0*ptot*ptot; 1144 G4double tmax = 4.0*ptot*ptot; 1146 G4double t = 0.0; 1145 G4double t = 0.0; 1147 1146 1148 1147 1149 // 1148 // 1150 // Sample t 1149 // Sample t 1151 // 1150 // 1152 1151 1153 t = SampleT( theParticle, ptot, A); 1152 t = SampleT( theParticle, ptot, A); 1154 1153 1155 // NaN finder 1154 // NaN finder 1156 if(!(t < 0.0 || t >= 0.0)) 1155 if(!(t < 0.0 || t >= 0.0)) 1157 { 1156 { 1158 if (verboseLevel > 0) 1157 if (verboseLevel > 0) 1159 { 1158 { 1160 G4cout << "G4NuclNuclDiffuseElastic:WAR 1159 G4cout << "G4NuclNuclDiffuseElastic:WARNING: A = " << A 1161 << " mom(GeV)= " << plab/GeV 1160 << " mom(GeV)= " << plab/GeV 1162 << " S-wave will be sampled" 1161 << " S-wave will be sampled" 1163 << G4endl; 1162 << G4endl; 1164 } 1163 } 1165 t = G4UniformRand()*tmax; 1164 t = G4UniformRand()*tmax; 1166 } 1165 } 1167 if(verboseLevel>1) 1166 if(verboseLevel>1) 1168 { 1167 { 1169 G4cout <<" t= " << t << " tmax= " << tmax 1168 G4cout <<" t= " << t << " tmax= " << tmax 1170 << " ptot= " << ptot << G4endl; 1169 << " ptot= " << ptot << G4endl; 1171 } 1170 } 1172 // Sampling of angles in CM system 1171 // Sampling of angles in CM system 1173 1172 1174 G4double phi = G4UniformRand()*twopi; 1173 G4double phi = G4UniformRand()*twopi; 1175 G4double cost = 1. - 2.0*t/tmax; 1174 G4double cost = 1. - 2.0*t/tmax; 1176 G4double sint; 1175 G4double sint; 1177 1176 1178 if( cost >= 1.0 ) 1177 if( cost >= 1.0 ) 1179 { 1178 { 1180 cost = 1.0; 1179 cost = 1.0; 1181 sint = 0.0; 1180 sint = 0.0; 1182 } 1181 } 1183 else if( cost <= -1.0) 1182 else if( cost <= -1.0) 1184 { 1183 { 1185 cost = -1.0; 1184 cost = -1.0; 1186 sint = 0.0; 1185 sint = 0.0; 1187 } 1186 } 1188 else 1187 else 1189 { 1188 { 1190 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1189 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1191 } 1190 } 1192 if (verboseLevel>1) 1191 if (verboseLevel>1) 1193 { 1192 { 1194 G4cout << "cos(t)=" << cost << " std::sin 1193 G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl; 1195 } 1194 } 1196 G4ThreeVector v1(sint*std::cos(phi),sint*st 1195 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost); 1197 v1 *= ptot; 1196 v1 *= ptot; 1198 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s 1197 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1)); 1199 1198 1200 nlv1.boost(bst); 1199 nlv1.boost(bst); 1201 1200 1202 G4ThreeVector np1 = nlv1.vect(); 1201 G4ThreeVector np1 = nlv1.vect(); 1203 1202 1204 // G4double theta = std::acos( np1.z()/np 1203 // G4double theta = std::acos( np1.z()/np1.mag() ); // degree; 1205 1204 1206 G4double theta = np1.theta(); 1205 G4double theta = np1.theta(); 1207 1206 1208 return theta; 1207 return theta; 1209 } 1208 } 1210 1209 1211 ///////////////////////////////////////////// 1210 //////////////////////////////////////////////////////////////////////////// 1212 // 1211 // 1213 // Return scattering angle in lab system (tar 1212 // Return scattering angle in lab system (target at rest) knowing theta in CMS 1214 1213 1215 1214 1216 1215 1217 G4double 1216 G4double 1218 G4NuclNuclDiffuseElastic::ThetaCMStoThetaLab( 1217 G4NuclNuclDiffuseElastic::ThetaCMStoThetaLab( const G4DynamicParticle* aParticle, 1219 G4dou 1218 G4double tmass, G4double thetaCMS) 1220 { 1219 { 1221 const G4ParticleDefinition* theParticle = a 1220 const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); 1222 G4double m1 = theParticle->GetPDGMass(); 1221 G4double m1 = theParticle->GetPDGMass(); 1223 // G4double plab = aParticle->GetTotalMomen 1222 // G4double plab = aParticle->GetTotalMomentum(); 1224 G4LorentzVector lv1 = aParticle->Get4Moment 1223 G4LorentzVector lv1 = aParticle->Get4Momentum(); 1225 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1224 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1226 1225 1227 lv += lv1; 1226 lv += lv1; 1228 1227 1229 G4ThreeVector bst = lv.boostVector(); 1228 G4ThreeVector bst = lv.boostVector(); 1230 1229 1231 lv1.boost(-bst); 1230 lv1.boost(-bst); 1232 1231 1233 G4ThreeVector p1 = lv1.vect(); 1232 G4ThreeVector p1 = lv1.vect(); 1234 G4double ptot = p1.mag(); 1233 G4double ptot = p1.mag(); 1235 1234 1236 G4double phi = G4UniformRand()*twopi; 1235 G4double phi = G4UniformRand()*twopi; 1237 G4double cost = std::cos(thetaCMS); 1236 G4double cost = std::cos(thetaCMS); 1238 G4double sint; 1237 G4double sint; 1239 1238 1240 if( cost >= 1.0 ) 1239 if( cost >= 1.0 ) 1241 { 1240 { 1242 cost = 1.0; 1241 cost = 1.0; 1243 sint = 0.0; 1242 sint = 0.0; 1244 } 1243 } 1245 else if( cost <= -1.0) 1244 else if( cost <= -1.0) 1246 { 1245 { 1247 cost = -1.0; 1246 cost = -1.0; 1248 sint = 0.0; 1247 sint = 0.0; 1249 } 1248 } 1250 else 1249 else 1251 { 1250 { 1252 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1251 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1253 } 1252 } 1254 if (verboseLevel>1) 1253 if (verboseLevel>1) 1255 { 1254 { 1256 G4cout << "cos(tcms)=" << cost << " std:: 1255 G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl; 1257 } 1256 } 1258 G4ThreeVector v1(sint*std::cos(phi),sint*st 1257 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost); 1259 v1 *= ptot; 1258 v1 *= ptot; 1260 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s 1259 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1)); 1261 1260 1262 nlv1.boost(bst); 1261 nlv1.boost(bst); 1263 1262 1264 G4ThreeVector np1 = nlv1.vect(); 1263 G4ThreeVector np1 = nlv1.vect(); 1265 1264 1266 1265 1267 G4double thetaLab = np1.theta(); 1266 G4double thetaLab = np1.theta(); 1268 1267 1269 return thetaLab; 1268 return thetaLab; 1270 } 1269 } 1271 1270 1272 ///////////////////////////////////////////// 1271 //////////////////////////////////////////////////////////////////////////// 1273 // 1272 // 1274 // Return scattering angle in CMS system (tar 1273 // Return scattering angle in CMS system (target at rest) knowing theta in Lab 1275 1274 1276 1275 1277 1276 1278 G4double 1277 G4double 1279 G4NuclNuclDiffuseElastic::ThetaLabToThetaCMS( 1278 G4NuclNuclDiffuseElastic::ThetaLabToThetaCMS( const G4DynamicParticle* aParticle, 1280 G4dou 1279 G4double tmass, G4double thetaLab) 1281 { 1280 { 1282 const G4ParticleDefinition* theParticle = a 1281 const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); 1283 G4double m1 = theParticle->GetPDGMass(); 1282 G4double m1 = theParticle->GetPDGMass(); 1284 G4double plab = aParticle->GetTotalMomentum 1283 G4double plab = aParticle->GetTotalMomentum(); 1285 G4LorentzVector lv1 = aParticle->Get4Moment 1284 G4LorentzVector lv1 = aParticle->Get4Momentum(); 1286 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1285 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1287 1286 1288 lv += lv1; 1287 lv += lv1; 1289 1288 1290 G4ThreeVector bst = lv.boostVector(); 1289 G4ThreeVector bst = lv.boostVector(); 1291 1290 1292 // lv1.boost(-bst); 1291 // lv1.boost(-bst); 1293 1292 1294 // G4ThreeVector p1 = lv1.vect(); 1293 // G4ThreeVector p1 = lv1.vect(); 1295 // G4double ptot = p1.mag(); 1294 // G4double ptot = p1.mag(); 1296 1295 1297 G4double phi = G4UniformRand()*twopi; 1296 G4double phi = G4UniformRand()*twopi; 1298 G4double cost = std::cos(thetaLab); 1297 G4double cost = std::cos(thetaLab); 1299 G4double sint; 1298 G4double sint; 1300 1299 1301 if( cost >= 1.0 ) 1300 if( cost >= 1.0 ) 1302 { 1301 { 1303 cost = 1.0; 1302 cost = 1.0; 1304 sint = 0.0; 1303 sint = 0.0; 1305 } 1304 } 1306 else if( cost <= -1.0) 1305 else if( cost <= -1.0) 1307 { 1306 { 1308 cost = -1.0; 1307 cost = -1.0; 1309 sint = 0.0; 1308 sint = 0.0; 1310 } 1309 } 1311 else 1310 else 1312 { 1311 { 1313 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1312 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1314 } 1313 } 1315 if (verboseLevel>1) 1314 if (verboseLevel>1) 1316 { 1315 { 1317 G4cout << "cos(tlab)=" << cost << " std:: 1316 G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl; 1318 } 1317 } 1319 G4ThreeVector v1(sint*std::cos(phi),sint*st 1318 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost); 1320 v1 *= plab; 1319 v1 *= plab; 1321 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s 1320 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1)); 1322 1321 1323 nlv1.boost(-bst); 1322 nlv1.boost(-bst); 1324 1323 1325 G4ThreeVector np1 = nlv1.vect(); 1324 G4ThreeVector np1 = nlv1.vect(); 1326 1325 1327 1326 1328 G4double thetaCMS = np1.theta(); 1327 G4double thetaCMS = np1.theta(); 1329 1328 1330 return thetaCMS; 1329 return thetaCMS; 1331 } 1330 } 1332 1331 1333 ///////////////////////////////////////////// 1332 /////////////////////////////////////////////////////////////////////////////// 1334 // 1333 // 1335 // Test for given particle and element table 1334 // Test for given particle and element table of momentum, angle probability. 1336 // For the moment in lab system. 1335 // For the moment in lab system. 1337 1336 1338 void G4NuclNuclDiffuseElastic::TestAngleTable 1337 void G4NuclNuclDiffuseElastic::TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom, 1339 G 1338 G4double Z, G4double A) 1340 { 1339 { 1341 fAtomicNumber = Z; // atomic number 1340 fAtomicNumber = Z; // atomic number 1342 fAtomicWeight = A; // number of nucleo 1341 fAtomicWeight = A; // number of nucleons 1343 fNuclearRadius = CalculateNuclearRad(fAtomi 1342 fNuclearRadius = CalculateNuclearRad(fAtomicWeight); 1344 1343 1345 1344 1346 1345 1347 G4cout<<"G4NuclNuclDiffuseElastic::TestAngl 1346 G4cout<<"G4NuclNuclDiffuseElastic::TestAngleTable() init the element with Z = " 1348 <<Z<<"; and A = "<<A<<G4endl; 1347 <<Z<<"; and A = "<<A<<G4endl; 1349 1348 1350 fElementNumberVector.push_back(fAtomicNumbe 1349 fElementNumberVector.push_back(fAtomicNumber); 1351 1350 1352 1351 1353 1352 1354 1353 1355 G4int i=0, j; 1354 G4int i=0, j; 1356 G4double a = 0., z = theParticle->GetPDGCha 1355 G4double a = 0., z = theParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass(); 1357 G4double alpha1=0., alpha2=0., alphaMax=0., 1356 G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.; 1358 G4double deltaL10 = 0., deltaL96 = 0., delt 1357 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.; 1359 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0. 1358 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.; 1360 G4double epsilon = 0.001; 1359 G4double epsilon = 0.001; 1361 1360 1362 G4Integrator<G4NuclNuclDiffuseElastic,G4dou 1361 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral; 1363 1362 1364 fAngleTable = new G4PhysicsTable(fEnergyBin 1363 fAngleTable = new G4PhysicsTable(fEnergyBin); 1365 1364 1366 fWaveVector = partMom/hbarc; 1365 fWaveVector = partMom/hbarc; 1367 1366 1368 G4double kR = fWaveVector*fNuclearRadiu 1367 G4double kR = fWaveVector*fNuclearRadius; 1369 G4double kR2 = kR*kR; 1368 G4double kR2 = kR*kR; 1370 G4double kRmax = 10.6; // 10.6, 18, 10.174 1369 G4double kRmax = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25. 1371 G4double kRcoul = 1.2; // 1.4, 2.5; // on t 1370 G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1 1372 1371 1373 alphaMax = kRmax*kRmax/kR2; 1372 alphaMax = kRmax*kRmax/kR2; 1374 1373 1375 if (alphaMax > 4.) alphaMax = 4.; // vmg05 1374 if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2 1376 1375 1377 alphaCoulomb = kRcoul*kRcoul/kR2; 1376 alphaCoulomb = kRcoul*kRcoul/kR2; 1378 1377 1379 if( z ) 1378 if( z ) 1380 { 1379 { 1381 a = partMom/m1; // beta*gamma 1380 a = partMom/m1; // beta*gamma for m1 1382 fBeta = a/std::sqrt(1+a*a); 1381 fBeta = a/std::sqrt(1+a*a); 1383 fZommerfeld = CalculateZommerfeld( fBet 1382 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber); 1384 fAm = CalculateAm( partMom, fZo 1383 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber); 1385 } 1384 } 1386 G4PhysicsFreeVector* angleVector = new G4Ph 1385 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1); 1387 1386 1388 // G4PhysicsLogVector* angleBins = new G4P 1387 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin ); 1389 1388 1390 1389 1391 fAddCoulomb = false; 1390 fAddCoulomb = false; 1392 1391 1393 for(j = 1; j < fAngleBin; j++) 1392 for(j = 1; j < fAngleBin; j++) 1394 { 1393 { 1395 // alpha1 = angleBins->GetLowEdgeEnergy 1394 // alpha1 = angleBins->GetLowEdgeEnergy(j-1); 1396 // alpha2 = angleBins->GetLowEdgeEnergy 1395 // alpha2 = angleBins->GetLowEdgeEnergy(j); 1397 1396 1398 alpha1 = alphaMax*(j-1)/fAngleBin; 1397 alpha1 = alphaMax*(j-1)/fAngleBin; 1399 alpha2 = alphaMax*( j )/fAngleBin; 1398 alpha2 = alphaMax*( j )/fAngleBin; 1400 1399 1401 if( ( alpha2 > alphaCoulomb ) && z ) fAdd 1400 if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true; 1402 1401 1403 deltaL10 = integral.Legendre10(this, &G4N 1402 deltaL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2); 1404 deltaL96 = integral.Legendre96(this, &G4N 1403 deltaL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2); 1405 deltaAG = integral.AdaptiveGauss(this, & 1404 deltaAG = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 1406 alpha1 1405 alpha1, alpha2,epsilon); 1407 1406 1408 // G4cout<<alpha1<<"\t"<<std::sqrt(alph 1407 // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t" 1409 // <<deltaL10<<"\t"<<deltaL96<<"\t" 1408 // <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl; 1410 1409 1411 sumL10 += deltaL10; 1410 sumL10 += deltaL10; 1412 sumL96 += deltaL96; 1411 sumL96 += deltaL96; 1413 sumAG += deltaAG; 1412 sumAG += deltaAG; 1414 1413 1415 G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/d 1414 G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t" 1416 <<sumL10<<"\t"<<sumL96<<"\t"<<sum 1415 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl; 1417 1416 1418 angleVector->PutValue( j-1 , alpha1, sumL 1417 angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2 1419 } 1418 } 1420 fAngleTable->insertAt(i,angleVector); 1419 fAngleTable->insertAt(i,angleVector); 1421 fAngleBank.push_back(fAngleTable); 1420 fAngleBank.push_back(fAngleTable); 1422 1421 1423 /* 1422 /* 1424 // Integral over all angle range - Bad accu 1423 // Integral over all angle range - Bad accuracy !!! 1425 1424 1426 sumL10 = integral.Legendre10(this, &G4NuclN 1425 sumL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2); 1427 sumL96 = integral.Legendre96(this, &G4NuclN 1426 sumL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2); 1428 sumAG = integral.AdaptiveGauss(this, &G4Nu 1427 sumAG = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 1429 0., al 1428 0., alpha2,epsilon); 1430 G4cout<<G4endl; 1429 G4cout<<G4endl; 1431 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/deg 1430 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t" 1432 <<sumL10<<"\t"<<sumL96<<"\t"<<sum 1431 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl; 1433 */ 1432 */ 1434 return; 1433 return; 1435 } 1434 } 1436 1435 1437 ///////////////////////////////////////////// 1436 ///////////////////////////////////////////////////////////////// 1438 // 1437 // 1439 // 1438 // 1440 1439 1441 G4double G4NuclNuclDiffuseElastic::GetLegend 1440 G4double G4NuclNuclDiffuseElastic::GetLegendrePol(G4int n, G4double theta) 1442 { 1441 { 1443 G4double legPol, epsilon = 1.e-6; 1442 G4double legPol, epsilon = 1.e-6; 1444 G4double x = std::cos(theta); 1443 G4double x = std::cos(theta); 1445 1444 1446 if ( n < 0 ) legPol = 0.; 1445 if ( n < 0 ) legPol = 0.; 1447 else if( n == 0 ) legPol = 1.; 1446 else if( n == 0 ) legPol = 1.; 1448 else if( n == 1 ) legPol = x; 1447 else if( n == 1 ) legPol = x; 1449 else if( n == 2 ) legPol = (3.*x*x-1.)/2.; 1448 else if( n == 2 ) legPol = (3.*x*x-1.)/2.; 1450 else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/ 1449 else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/2.; 1451 else if( n == 4 ) legPol = (35.*x*x*x*x-30. 1450 else if( n == 4 ) legPol = (35.*x*x*x*x-30.*x*x+3.)/8.; 1452 else if( n == 5 ) legPol = (63.*x*x*x*x*x-7 1451 else if( n == 5 ) legPol = (63.*x*x*x*x*x-70.*x*x*x+15.*x)/8.; 1453 else if( n == 6 ) legPol = (231.*x*x*x*x*x* 1452 else if( n == 6 ) legPol = (231.*x*x*x*x*x*x-315.*x*x*x*x+105.*x*x-5.)/16.; 1454 else 1453 else 1455 { 1454 { 1456 // legPol = ( (2*n-1)*x*GetLegendrePol(n- 1455 // legPol = ( (2*n-1)*x*GetLegendrePol(n-1,x) - (n-1)*GetLegendrePol(n-2,x) )/n; 1457 1456 1458 legPol = std::sqrt( 2./(n*CLHEP::pi*std:: 1457 legPol = std::sqrt( 2./(n*CLHEP::pi*std::sin(theta+epsilon)) )*std::sin( (n+0.5)*theta+0.25*CLHEP::pi ); 1459 } 1458 } 1460 return legPol; 1459 return legPol; 1461 } 1460 } 1462 1461 1463 ///////////////////////////////////////////// 1462 ///////////////////////////////////////////////////////////////// 1464 // 1463 // 1465 // 1464 // 1466 1465 1467 G4complex G4NuclNuclDiffuseElastic::GetErfCom 1466 G4complex G4NuclNuclDiffuseElastic::GetErfComp(G4complex z, G4int nMax) 1468 { 1467 { 1469 G4int n; 1468 G4int n; 1470 G4double n2, cofn, shny, chny, fn, gn; 1469 G4double n2, cofn, shny, chny, fn, gn; 1471 1470 1472 G4double x = z.real(); 1471 G4double x = z.real(); 1473 G4double y = z.imag(); 1472 G4double y = z.imag(); 1474 1473 1475 G4double outRe = 0., outIm = 0.; 1474 G4double outRe = 0., outIm = 0.; 1476 1475 1477 G4double twox = 2.*x; 1476 G4double twox = 2.*x; 1478 G4double twoxy = twox*y; 1477 G4double twoxy = twox*y; 1479 G4double twox2 = twox*twox; 1478 G4double twox2 = twox*twox; 1480 1479 1481 G4double cof1 = G4Exp(-x*x)/CLHEP::pi; 1480 G4double cof1 = G4Exp(-x*x)/CLHEP::pi; 1482 1481 1483 G4double cos2xy = std::cos(twoxy); 1482 G4double cos2xy = std::cos(twoxy); 1484 G4double sin2xy = std::sin(twoxy); 1483 G4double sin2xy = std::sin(twoxy); 1485 1484 1486 G4double twoxcos2xy = twox*cos2xy; 1485 G4double twoxcos2xy = twox*cos2xy; 1487 G4double twoxsin2xy = twox*sin2xy; 1486 G4double twoxsin2xy = twox*sin2xy; 1488 1487 1489 for( n = 1; n <= nMax; n++) 1488 for( n = 1; n <= nMax; n++) 1490 { 1489 { 1491 n2 = n*n; 1490 n2 = n*n; 1492 1491 1493 cofn = G4Exp(-0.5*n2)/(n2+twox2); // /(n 1492 cofn = G4Exp(-0.5*n2)/(n2+twox2); // /(n2+0.5*twox2); 1494 1493 1495 chny = std::cosh(n*y); 1494 chny = std::cosh(n*y); 1496 shny = std::sinh(n*y); 1495 shny = std::sinh(n*y); 1497 1496 1498 fn = twox - twoxcos2xy*chny + n*sin2xy*s 1497 fn = twox - twoxcos2xy*chny + n*sin2xy*shny; 1499 gn = twoxsin2xy*chny + n*cos2xy*s 1498 gn = twoxsin2xy*chny + n*cos2xy*shny; 1500 1499 1501 fn *= cofn; 1500 fn *= cofn; 1502 gn *= cofn; 1501 gn *= cofn; 1503 1502 1504 outRe += fn; 1503 outRe += fn; 1505 outIm += gn; 1504 outIm += gn; 1506 } 1505 } 1507 outRe *= 2*cof1; 1506 outRe *= 2*cof1; 1508 outIm *= 2*cof1; 1507 outIm *= 2*cof1; 1509 1508 1510 if(std::abs(x) < 0.0001) 1509 if(std::abs(x) < 0.0001) 1511 { 1510 { 1512 outRe += GetErf(x); 1511 outRe += GetErf(x); 1513 outIm += cof1*y; 1512 outIm += cof1*y; 1514 } 1513 } 1515 else 1514 else 1516 { 1515 { 1517 outRe += GetErf(x) + cof1*(1-cos2xy)/twox 1516 outRe += GetErf(x) + cof1*(1-cos2xy)/twox; 1518 outIm += cof1*sin2xy/twox; 1517 outIm += cof1*sin2xy/twox; 1519 } 1518 } 1520 return G4complex(outRe, outIm); 1519 return G4complex(outRe, outIm); 1521 } 1520 } 1522 1521 1523 1522 1524 ///////////////////////////////////////////// 1523 ///////////////////////////////////////////////////////////////// 1525 // 1524 // 1526 // 1525 // 1527 1526 1528 G4complex G4NuclNuclDiffuseElastic::GetErfInt 1527 G4complex G4NuclNuclDiffuseElastic::GetErfInt(G4complex z) // , G4int nMax) 1529 { 1528 { 1530 G4double outRe, outIm; 1529 G4double outRe, outIm; 1531 1530 1532 G4double x = z.real(); 1531 G4double x = z.real(); 1533 G4double y = z.imag(); 1532 G4double y = z.imag(); 1534 fReZ = x; 1533 fReZ = x; 1535 1534 1536 G4Integrator<G4NuclNuclDiffuseElastic,G4dou 1535 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral; 1537 1536 1538 outRe = integral.Legendre96(this,&G4NuclNuc 1537 outRe = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpSin, 0., y ); 1539 outIm = integral.Legendre96(this,&G4NuclNuc 1538 outIm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpCos, 0., y ); 1540 1539 1541 outRe *= 2./std::sqrt(CLHEP::pi); 1540 outRe *= 2./std::sqrt(CLHEP::pi); 1542 outIm *= 2./std::sqrt(CLHEP::pi); 1541 outIm *= 2./std::sqrt(CLHEP::pi); 1543 1542 1544 outRe += GetErf(x); 1543 outRe += GetErf(x); 1545 1544 1546 return G4complex(outRe, outIm); 1545 return G4complex(outRe, outIm); 1547 } 1546 } 1548 1547 1549 ///////////////////////////////////////////// 1548 ///////////////////////////////////////////////////////////////// 1550 // 1549 // 1551 // 1550 // 1552 1551 1553 1552 1554 G4complex G4NuclNuclDiffuseElastic::GammaLess 1553 G4complex G4NuclNuclDiffuseElastic::GammaLess(G4double theta) 1555 { 1554 { 1556 G4double sinThetaR = 2.*fHalfRutThetaT 1555 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2); 1557 G4double cosHalfThetaR2 = 1./(1. + fHalfRut 1556 G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2); 1558 1557 1559 G4double u = std::sqrt(0.5*fPr 1558 G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR); 1560 G4double kappa = u/std::sqrt(CLHEP 1559 G4double kappa = u/std::sqrt(CLHEP::pi); 1561 G4double dTheta = theta - fRutherfo 1560 G4double dTheta = theta - fRutherfordTheta; 1562 u *= dTheta; 1561 u *= dTheta; 1563 G4double u2 = u*u; 1562 G4double u2 = u*u; 1564 G4double u2m2p3 = u2*2./3.; 1563 G4double u2m2p3 = u2*2./3.; 1565 1564 1566 G4complex im = G4complex(0.,1.); 1565 G4complex im = G4complex(0.,1.); 1567 G4complex order = G4complex(u,u); 1566 G4complex order = G4complex(u,u); 1568 order /= std::sqrt(2.); 1567 order /= std::sqrt(2.); 1569 1568 1570 G4complex gamma = CLHEP::pi*kappa*G 1569 G4complex gamma = CLHEP::pi*kappa*GetErfcInt(-order)*std::exp(im*(u*u+0.25*CLHEP::pi)); 1571 G4complex a0 = 0.5*(1. + 4.*(1.+ 1570 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR; 1572 G4complex a1 = 0.5*(1. + 2.*(1.+ 1571 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR; 1573 G4complex out = gamma*(1. - a1*dT 1572 G4complex out = gamma*(1. - a1*dTheta) - a0; 1574 1573 1575 return out; 1574 return out; 1576 } 1575 } 1577 1576 1578 ///////////////////////////////////////////// 1577 ///////////////////////////////////////////////////////////////// 1579 // 1578 // 1580 // 1579 // 1581 1580 1582 G4complex G4NuclNuclDiffuseElastic::GammaMore 1581 G4complex G4NuclNuclDiffuseElastic::GammaMore(G4double theta) 1583 { 1582 { 1584 G4double sinThetaR = 2.*fHalfRutThetaT 1583 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2); 1585 G4double cosHalfThetaR2 = 1./(1. + fHalfRut 1584 G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2); 1586 1585 1587 G4double u = std::sqrt(0.5*fPr 1586 G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR); 1588 G4double kappa = u/std::sqrt(CLHEP 1587 G4double kappa = u/std::sqrt(CLHEP::pi); 1589 G4double dTheta = theta - fRutherfo 1588 G4double dTheta = theta - fRutherfordTheta; 1590 u *= dTheta; 1589 u *= dTheta; 1591 G4double u2 = u*u; 1590 G4double u2 = u*u; 1592 G4double u2m2p3 = u2*2./3.; 1591 G4double u2m2p3 = u2*2./3.; 1593 1592 1594 G4complex im = G4complex(0.,1.); 1593 G4complex im = G4complex(0.,1.); 1595 G4complex order = G4complex(u,u); 1594 G4complex order = G4complex(u,u); 1596 order /= std::sqrt(2.); 1595 order /= std::sqrt(2.); 1597 G4complex gamma = CLHEP::pi*kappa*G 1596 G4complex gamma = CLHEP::pi*kappa*GetErfcInt(order)*std::exp(im*(u*u+0.25*CLHEP::pi)); 1598 G4complex a0 = 0.5*(1. + 4.*(1.+ 1597 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR; 1599 G4complex a1 = 0.5*(1. + 2.*(1.+ 1598 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR; 1600 G4complex out = -gamma*(1. - a1*d 1599 G4complex out = -gamma*(1. - a1*dTheta) - a0; 1601 1600 1602 return out; 1601 return out; 1603 } 1602 } 1604 1603 1605 ///////////////////////////////////////////// 1604 ///////////////////////////////////////////////////////////////// 1606 // 1605 // 1607 // 1606 // 1608 1607 1609 G4complex G4NuclNuclDiffuseElastic::Amplitud 1608 G4complex G4NuclNuclDiffuseElastic::AmplitudeNear(G4double theta) 1610 { 1609 { 1611 G4double kappa = std::sqrt(0.5*fProfileLamb 1610 G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi); 1612 G4complex out = G4complex(kappa/fWaveVector 1611 G4complex out = G4complex(kappa/fWaveVector,0.); 1613 1612 1614 out *= PhaseNear(theta); 1613 out *= PhaseNear(theta); 1615 1614 1616 if( theta <= fRutherfordTheta ) 1615 if( theta <= fRutherfordTheta ) 1617 { 1616 { 1618 out *= GammaLess(theta) + ProfileNear(the 1617 out *= GammaLess(theta) + ProfileNear(theta); 1619 // out *= GammaMore(theta) + ProfileNear( 1618 // out *= GammaMore(theta) + ProfileNear(theta); 1620 out += CoulombAmplitude(theta); 1619 out += CoulombAmplitude(theta); 1621 } 1620 } 1622 else 1621 else 1623 { 1622 { 1624 out *= GammaMore(theta) + ProfileNear(the 1623 out *= GammaMore(theta) + ProfileNear(theta); 1625 // out *= GammaLess(theta) + ProfileNear( 1624 // out *= GammaLess(theta) + ProfileNear(theta); 1626 } 1625 } 1627 return out; 1626 return out; 1628 } 1627 } 1629 1628 1630 ///////////////////////////////////////////// 1629 ///////////////////////////////////////////////////////////////// 1631 // 1630 // 1632 // 1631 // 1633 1632 1634 G4complex G4NuclNuclDiffuseElastic::Amplitude 1633 G4complex G4NuclNuclDiffuseElastic::AmplitudeSim(G4double theta) 1635 { 1634 { 1636 G4double sinThetaR = 2.*fHalfRutThetaTg/(1 1635 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2); 1637 G4double dTheta = 0.5*(theta - fRutherf 1636 G4double dTheta = 0.5*(theta - fRutherfordTheta); 1638 G4double sindTheta = std::sin(dTheta); 1637 G4double sindTheta = std::sin(dTheta); 1639 G4double persqrt2 = std::sqrt(0.5); 1638 G4double persqrt2 = std::sqrt(0.5); 1640 1639 1641 G4complex order = G4complex(persqrt2,pe 1640 G4complex order = G4complex(persqrt2,persqrt2); 1642 order *= std::sqrt(0.5*fProfil 1641 order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*sindTheta; 1643 // order *= std::sqrt(0.5*fPro 1642 // order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*dTheta; 1644 1643 1645 G4complex out; 1644 G4complex out; 1646 1645 1647 if ( theta <= fRutherfordTheta ) 1646 if ( theta <= fRutherfordTheta ) 1648 { 1647 { 1649 out = 1. - 0.5*GetErfcInt(-order)*Profile 1648 out = 1. - 0.5*GetErfcInt(-order)*ProfileNear(theta); 1650 } 1649 } 1651 else 1650 else 1652 { 1651 { 1653 out = 0.5*GetErfcInt(order)*ProfileNear(t 1652 out = 0.5*GetErfcInt(order)*ProfileNear(theta); 1654 } 1653 } 1655 1654 1656 out *= CoulombAmplitude(theta); 1655 out *= CoulombAmplitude(theta); 1657 1656 1658 return out; 1657 return out; 1659 } 1658 } 1660 1659 1661 ///////////////////////////////////////////// 1660 ///////////////////////////////////////////////////////////////// 1662 // 1661 // 1663 // 1662 // 1664 1663 1665 G4complex G4NuclNuclDiffuseElastic::Amplitud 1664 G4complex G4NuclNuclDiffuseElastic::AmplitudeGla(G4double theta) 1666 { 1665 { 1667 G4int n; 1666 G4int n; 1668 G4double T12b, b, b2; // cosTheta = std::co 1667 G4double T12b, b, b2; // cosTheta = std::cos(theta); 1669 G4complex out = G4complex(0.,0.), shiftC, s 1668 G4complex out = G4complex(0.,0.), shiftC, shiftN; 1670 G4complex im = G4complex(0.,1.); 1669 G4complex im = G4complex(0.,1.); 1671 1670 1672 for( n = 0; n < fMaxL; n++) 1671 for( n = 0; n < fMaxL; n++) 1673 { 1672 { 1674 shiftC = std::exp( im*2.*CalculateCoulomb 1673 shiftC = std::exp( im*2.*CalculateCoulombPhase(n) ); 1675 // b = ( fZommerfeld + std::sqrt( fZommer 1674 // b = ( fZommerfeld + std::sqrt( fZommerfeld*fZommerfeld + n*(n+1) ) )/fWaveVector; 1676 b = ( std::sqrt( G4double(n*(n+1)) ) )/fW 1675 b = ( std::sqrt( G4double(n*(n+1)) ) )/fWaveVector; 1677 b2 = b*b; 1676 b2 = b*b; 1678 T12b = fSumSigma*G4Exp(-b2/fNuclearRadius 1677 T12b = fSumSigma*G4Exp(-b2/fNuclearRadiusSquare)/CLHEP::pi/fNuclearRadiusSquare; 1679 shiftN = std::exp( -0.5*(1.-im*fEtaRatio) 1678 shiftN = std::exp( -0.5*(1.-im*fEtaRatio)*T12b ) - 1.; 1680 out += (2.*n+1.)*shiftC*shiftN*GetLegend 1679 out += (2.*n+1.)*shiftC*shiftN*GetLegendrePol(n, theta); 1681 } 1680 } 1682 out /= 2.*im*fWaveVector; 1681 out /= 2.*im*fWaveVector; 1683 out += CoulombAmplitude(theta); 1682 out += CoulombAmplitude(theta); 1684 return out; 1683 return out; 1685 } 1684 } 1686 1685 1687 1686 1688 ///////////////////////////////////////////// 1687 ///////////////////////////////////////////////////////////////// 1689 // 1688 // 1690 // 1689 // 1691 1690 1692 G4complex G4NuclNuclDiffuseElastic::Amplitud 1691 G4complex G4NuclNuclDiffuseElastic::AmplitudeGG(G4double theta) 1693 { 1692 { 1694 G4int n; 1693 G4int n; 1695 G4double T12b, a, aTemp, b2, sinThetaH = st 1694 G4double T12b, a, aTemp, b2, sinThetaH = std::sin(0.5*theta); 1696 G4double sinThetaH2 = sinThetaH*sinThetaH; 1695 G4double sinThetaH2 = sinThetaH*sinThetaH; 1697 G4complex out = G4complex(0.,0.); 1696 G4complex out = G4complex(0.,0.); 1698 G4complex im = G4complex(0.,1.); 1697 G4complex im = G4complex(0.,1.); 1699 1698 1700 a = -fSumSigma/CLHEP::twopi/fNuclearRadius 1699 a = -fSumSigma/CLHEP::twopi/fNuclearRadiusSquare; 1701 b2 = fWaveVector*fWaveVector*fNuclearRadius 1700 b2 = fWaveVector*fWaveVector*fNuclearRadiusSquare*sinThetaH2; 1702 1701 1703 aTemp = a; 1702 aTemp = a; 1704 1703 1705 for( n = 1; n < fMaxL; n++) 1704 for( n = 1; n < fMaxL; n++) 1706 { 1705 { 1707 T12b = aTemp*G4Exp(-b2/n)/n; 1706 T12b = aTemp*G4Exp(-b2/n)/n; 1708 aTemp *= a; 1707 aTemp *= a; 1709 out += T12b; 1708 out += T12b; 1710 G4cout<<"out = "<<out<<G4endl; 1709 G4cout<<"out = "<<out<<G4endl; 1711 } 1710 } 1712 out *= -4.*im*fWaveVector/CLHEP::pi; 1711 out *= -4.*im*fWaveVector/CLHEP::pi; 1713 out += CoulombAmplitude(theta); 1712 out += CoulombAmplitude(theta); 1714 return out; 1713 return out; 1715 } 1714 } 1716 1715 1717 1716 1718 ///////////////////////////////////////////// 1717 /////////////////////////////////////////////////////////////////////////////// 1719 // 1718 // 1720 // Test for given particle and element table 1719 // Test for given particle and element table of momentum, angle probability. 1721 // For the partMom in CMS. 1720 // For the partMom in CMS. 1722 1721 1723 void G4NuclNuclDiffuseElastic::InitParameters 1722 void G4NuclNuclDiffuseElastic::InitParameters(const G4ParticleDefinition* theParticle, 1724 G4d 1723 G4double partMom, G4double Z, G4double A) 1725 { 1724 { 1726 fAtomicNumber = Z; // atomic number 1725 fAtomicNumber = Z; // atomic number 1727 fAtomicWeight = A; // number of nucleo 1726 fAtomicWeight = A; // number of nucleons 1728 1727 1729 fNuclearRadius2 = CalculateNuclearRad(fAtom 1728 fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight); 1730 G4double A1 = G4double( theParticle->Ge 1729 G4double A1 = G4double( theParticle->GetBaryonNumber() ); 1731 fNuclearRadius1 = CalculateNuclearRad(A1); 1730 fNuclearRadius1 = CalculateNuclearRad(A1); 1732 // fNuclearRadius = std::sqrt(fNuclearRadiu 1731 // fNuclearRadius = std::sqrt(fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2); 1733 fNuclearRadius = fNuclearRadius1 + fNuclear 1732 fNuclearRadius = fNuclearRadius1 + fNuclearRadius2; 1734 1733 1735 G4double a = 0.; 1734 G4double a = 0.; 1736 G4double z = theParticle->GetPDGCharge(); 1735 G4double z = theParticle->GetPDGCharge(); 1737 G4double m1 = theParticle->GetPDGMass(); 1736 G4double m1 = theParticle->GetPDGMass(); 1738 1737 1739 fWaveVector = partMom/CLHEP::hbarc; 1738 fWaveVector = partMom/CLHEP::hbarc; 1740 1739 1741 G4double lambda = fCofLambda*fWaveVector*fN 1740 G4double lambda = fCofLambda*fWaveVector*fNuclearRadius; 1742 G4cout<<"kR = "<<lambda<<G4endl; 1741 G4cout<<"kR = "<<lambda<<G4endl; 1743 1742 1744 if( z ) 1743 if( z ) 1745 { 1744 { 1746 a = partMom/m1; // beta*gamma f 1745 a = partMom/m1; // beta*gamma for m1 1747 fBeta = a/std::sqrt(1+a*a); 1746 fBeta = a/std::sqrt(1+a*a); 1748 fZommerfeld = CalculateZommerfeld( fBeta, 1747 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber); 1749 fRutherfordRatio = fZommerfeld/fWaveVecto 1748 fRutherfordRatio = fZommerfeld/fWaveVector; 1750 fAm = CalculateAm( partMom, fZomm 1749 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber); 1751 } 1750 } 1752 G4cout<<"fZommerfeld = "<<fZommerfeld<<G4en 1751 G4cout<<"fZommerfeld = "<<fZommerfeld<<G4endl; 1753 fProfileLambda = lambda; // *std::sqrt(1.-2 1752 fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda); 1754 G4cout<<"fProfileLambda = "<<fProfileLambda 1753 G4cout<<"fProfileLambda = "<<fProfileLambda<<G4endl; 1755 fProfileDelta = fCofDelta*fProfileLambda; 1754 fProfileDelta = fCofDelta*fProfileLambda; 1756 fProfileAlpha = fCofAlpha*fProfileLambda; 1755 fProfileAlpha = fCofAlpha*fProfileLambda; 1757 1756 1758 CalculateCoulombPhaseZero(); 1757 CalculateCoulombPhaseZero(); 1759 CalculateRutherfordAnglePar(); 1758 CalculateRutherfordAnglePar(); 1760 1759 1761 return; 1760 return; 1762 } 1761 } 1763 ///////////////////////////////////////////// 1762 /////////////////////////////////////////////////////////////////////////////// 1764 // 1763 // 1765 // Test for given particle and element table 1764 // Test for given particle and element table of momentum, angle probability. 1766 // For the partMom in CMS. 1765 // For the partMom in CMS. 1767 1766 1768 void G4NuclNuclDiffuseElastic::InitDynParamet 1767 void G4NuclNuclDiffuseElastic::InitDynParameters(const G4ParticleDefinition* theParticle, 1769 G4d 1768 G4double partMom) 1770 { 1769 { 1771 G4double a = 0.; 1770 G4double a = 0.; 1772 G4double z = theParticle->GetPDGCharge(); 1771 G4double z = theParticle->GetPDGCharge(); 1773 G4double m1 = theParticle->GetPDGMass(); 1772 G4double m1 = theParticle->GetPDGMass(); 1774 1773 1775 fWaveVector = partMom/CLHEP::hbarc; 1774 fWaveVector = partMom/CLHEP::hbarc; 1776 1775 1777 G4double lambda = fCofLambda*fWaveVector*fN 1776 G4double lambda = fCofLambda*fWaveVector*fNuclearRadius; 1778 1777 1779 if( z ) 1778 if( z ) 1780 { 1779 { 1781 a = partMom/m1; // beta*gamma f 1780 a = partMom/m1; // beta*gamma for m1 1782 fBeta = a/std::sqrt(1+a*a); 1781 fBeta = a/std::sqrt(1+a*a); 1783 fZommerfeld = CalculateZommerfeld( fBeta, 1782 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber); 1784 fRutherfordRatio = fZommerfeld/fWaveVecto 1783 fRutherfordRatio = fZommerfeld/fWaveVector; 1785 fAm = CalculateAm( partMom, fZomm 1784 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber); 1786 } 1785 } 1787 fProfileLambda = lambda; // *std::sqrt(1.-2 1786 fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda); 1788 fProfileDelta = fCofDelta*fProfileLambda; 1787 fProfileDelta = fCofDelta*fProfileLambda; 1789 fProfileAlpha = fCofAlpha*fProfileLambda; 1788 fProfileAlpha = fCofAlpha*fProfileLambda; 1790 1789 1791 CalculateCoulombPhaseZero(); 1790 CalculateCoulombPhaseZero(); 1792 CalculateRutherfordAnglePar(); 1791 CalculateRutherfordAnglePar(); 1793 1792 1794 return; 1793 return; 1795 } 1794 } 1796 1795 1797 1796 1798 ///////////////////////////////////////////// 1797 /////////////////////////////////////////////////////////////////////////////// 1799 // 1798 // 1800 // Test for given particle and element table 1799 // Test for given particle and element table of momentum, angle probability. 1801 // For the partMom in CMS. 1800 // For the partMom in CMS. 1802 1801 1803 void G4NuclNuclDiffuseElastic::InitParameters 1802 void G4NuclNuclDiffuseElastic::InitParametersGla(const G4DynamicParticle* aParticle, 1804 G4d 1803 G4double partMom, G4double Z, G4double A) 1805 { 1804 { 1806 fAtomicNumber = Z; // target atomic nu 1805 fAtomicNumber = Z; // target atomic number 1807 fAtomicWeight = A; // target number of 1806 fAtomicWeight = A; // target number of nucleons 1808 1807 1809 fNuclearRadius2 = CalculateNuclearRad(fAtom 1808 fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight); // target nucleus radius 1810 G4double A1 = G4double( aParticle->GetD 1809 G4double A1 = G4double( aParticle->GetDefinition()->GetBaryonNumber() ); 1811 fNuclearRadius1 = CalculateNuclearRad(A1); 1810 fNuclearRadius1 = CalculateNuclearRad(A1); // projectile nucleus radius 1812 fNuclearRadiusSquare = fNuclearRadius1*fNuc 1811 fNuclearRadiusSquare = fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2; 1813 1812 1814 1813 1815 G4double a = 0., kR12; 1814 G4double a = 0., kR12; 1816 G4double z = aParticle->GetDefinition()->G 1815 G4double z = aParticle->GetDefinition()->GetPDGCharge(); 1817 G4double m1 = aParticle->GetDefinition()->G 1816 G4double m1 = aParticle->GetDefinition()->GetPDGMass(); 1818 1817 1819 fWaveVector = partMom/CLHEP::hbarc; 1818 fWaveVector = partMom/CLHEP::hbarc; 1820 1819 1821 G4double pN = A1 - z; 1820 G4double pN = A1 - z; 1822 if( pN < 0. ) pN = 0.; 1821 if( pN < 0. ) pN = 0.; 1823 1822 1824 G4double tN = A - Z; 1823 G4double tN = A - Z; 1825 if( tN < 0. ) tN = 0.; 1824 if( tN < 0. ) tN = 0.; 1826 1825 1827 G4double pTkin = aParticle->GetKineticEnerg 1826 G4double pTkin = aParticle->GetKineticEnergy(); 1828 pTkin /= A1; 1827 pTkin /= A1; 1829 1828 1830 1829 1831 fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXsc 1830 fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) + 1832 (z*tN+pN*Z)*GetHadronNucleonXsc 1831 (z*tN+pN*Z)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron); 1833 1832 1834 G4cout<<"fSumSigma = "<<fSumSigma/CLHEP::mi 1833 G4cout<<"fSumSigma = "<<fSumSigma/CLHEP::millibarn<<" mb"<<G4endl; 1835 G4cout<<"pi*R2 = "<<CLHEP::pi*fNuclearRadiu 1834 G4cout<<"pi*R2 = "<<CLHEP::pi*fNuclearRadiusSquare/CLHEP::millibarn<<" mb"<<G4endl; 1836 kR12 = fWaveVector*std::sqrt(fNuclearRadius 1835 kR12 = fWaveVector*std::sqrt(fNuclearRadiusSquare); 1837 G4cout<<"k*sqrt(R2) = "<<kR12<<" "<<G4endl; 1836 G4cout<<"k*sqrt(R2) = "<<kR12<<" "<<G4endl; 1838 fMaxL = (G4int(kR12)+1)*4; 1837 fMaxL = (G4int(kR12)+1)*4; 1839 G4cout<<"fMaxL = "<<fMaxL<<" "<<G4endl; 1838 G4cout<<"fMaxL = "<<fMaxL<<" "<<G4endl; 1840 1839 1841 if( z ) 1840 if( z ) 1842 { 1841 { 1843 a = partMom/m1; // beta*gamma 1842 a = partMom/m1; // beta*gamma for m1 1844 fBeta = a/std::sqrt(1+a*a); 1843 fBeta = a/std::sqrt(1+a*a); 1845 fZommerfeld = CalculateZommerfeld( fBet 1844 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber); 1846 fAm = CalculateAm( partMom, fZo 1845 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber); 1847 } 1846 } 1848 1847 1849 CalculateCoulombPhaseZero(); 1848 CalculateCoulombPhaseZero(); 1850 1849 1851 1850 1852 return; 1851 return; 1853 } 1852 } 1854 1853 1855 1854 1856 ///////////////////////////////////////////// 1855 ///////////////////////////////////////////////////////////////////////////////////// 1857 // 1856 // 1858 // Returns nucleon-nucleon cross-section base 1857 // Returns nucleon-nucleon cross-section based on N. Starkov parametrisation of 1859 // data from mainly http://wwwppds.ihep.su:80 1858 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database 1860 // projectile nucleon is pParticle with pTkin 1859 // projectile nucleon is pParticle with pTkin shooting target nucleon tParticle 1861 1860 1862 G4double 1861 G4double 1863 G4NuclNuclDiffuseElastic::GetHadronNucleonXsc 1862 G4NuclNuclDiffuseElastic::GetHadronNucleonXscNS( G4ParticleDefinition* pParticle, 1864 1863 G4double pTkin, 1865 1864 G4ParticleDefinition* tParticle) 1866 { 1865 { 1867 G4double xsection(0), /*Delta,*/ A0, B0; 1866 G4double xsection(0), /*Delta,*/ A0, B0; 1868 G4double hpXsc(0); 1867 G4double hpXsc(0); 1869 G4double hnXsc(0); 1868 G4double hnXsc(0); 1870 1869 1871 1870 1872 G4double targ_mass = tParticle->GetPDGM 1871 G4double targ_mass = tParticle->GetPDGMass(); 1873 G4double proj_mass = pParticle->GetPDGM 1872 G4double proj_mass = pParticle->GetPDGMass(); 1874 1873 1875 G4double proj_energy = proj_mass + pTkin; 1874 G4double proj_energy = proj_mass + pTkin; 1876 G4double proj_momentum = std::sqrt(pTkin*(p 1875 G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass)); 1877 1876 1878 G4double sMand = CalcMandelstamS ( proj_mas 1877 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum ); 1879 1878 1880 sMand /= CLHEP::GeV*CLHEP::GeV; // 1879 sMand /= CLHEP::GeV*CLHEP::GeV; // in GeV for parametrisation 1881 proj_momentum /= CLHEP::GeV; 1880 proj_momentum /= CLHEP::GeV; 1882 proj_energy /= CLHEP::GeV; 1881 proj_energy /= CLHEP::GeV; 1883 proj_mass /= CLHEP::GeV; 1882 proj_mass /= CLHEP::GeV; 1884 G4double logS = G4Log(sMand); 1883 G4double logS = G4Log(sMand); 1885 1884 1886 // General PDG fit constants 1885 // General PDG fit constants 1887 1886 1888 1887 1889 // fEtaRatio=Re[f(0)]/Im[f(0)] 1888 // fEtaRatio=Re[f(0)]/Im[f(0)] 1890 1889 1891 if( proj_momentum >= 1.2 ) 1890 if( proj_momentum >= 1.2 ) 1892 { 1891 { 1893 fEtaRatio = 0.13*(logS - 5.8579332)*G4Po 1892 fEtaRatio = 0.13*(logS - 5.8579332)*G4Pow::GetInstance()->powA(sMand,-0.18); 1894 } 1893 } 1895 else if( proj_momentum >= 0.6 ) 1894 else if( proj_momentum >= 0.6 ) 1896 { 1895 { 1897 fEtaRatio = -75.5*(G4Pow::GetInstance()-> 1896 fEtaRatio = -75.5*(G4Pow::GetInstance()->powA(proj_momentum,0.25)-0.95)/ 1898 (G4Pow::GetInstance()->powA(3*proj_moment 1897 (G4Pow::GetInstance()->powA(3*proj_momentum,2.2)+1); 1899 } 1898 } 1900 else 1899 else 1901 { 1900 { 1902 fEtaRatio = 15.5*proj_momentum/(27*proj_m 1901 fEtaRatio = 15.5*proj_momentum/(27*proj_momentum*proj_momentum*proj_momentum+2); 1903 } 1902 } 1904 G4cout<<"fEtaRatio = "<<fEtaRatio<<G4endl; 1903 G4cout<<"fEtaRatio = "<<fEtaRatio<<G4endl; 1905 1904 1906 // xsc 1905 // xsc 1907 1906 1908 if( proj_momentum >= 10. ) // high energy: 1907 if( proj_momentum >= 10. ) // high energy: pp = nn = np 1909 // if( proj_momentum >= 2.) 1908 // if( proj_momentum >= 2.) 1910 { 1909 { 1911 //Delta = 1.; 1910 //Delta = 1.; 1912 1911 1913 //if( proj_energy < 40. ) Delta = 0.916+0 1912 //if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy; 1914 1913 1915 //AR-12Aug2016 if( proj_momentum >= 10.) 1914 //AR-12Aug2016 if( proj_momentum >= 10.) 1916 { 1915 { 1917 B0 = 7.5; 1916 B0 = 7.5; 1918 A0 = 100. - B0*G4Log(3.0e7); 1917 A0 = 100. - B0*G4Log(3.0e7); 1919 1918 1920 xsection = A0 + B0*G4Log(proj_energy) 1919 xsection = A0 + B0*G4Log(proj_energy) - 11 1921 + 103*G4Pow::GetInstance()- 1920 + 103*G4Pow::GetInstance()->powA(2*0.93827*proj_energy + proj_mass*proj_mass+ 1922 0.93827*0.93827,-0.165); 1921 0.93827*0.93827,-0.165); // mb 1923 } 1922 } 1924 } 1923 } 1925 else // low energy pp = nn != np 1924 else // low energy pp = nn != np 1926 { 1925 { 1927 if(pParticle == tParticle) // pp or nn 1926 if(pParticle == tParticle) // pp or nn // nn to be pp 1928 { 1927 { 1929 if( proj_momentum < 0.73 ) 1928 if( proj_momentum < 0.73 ) 1930 { 1929 { 1931 hnXsc = 23 + 50*( G4Pow::GetInstanc 1930 hnXsc = 23 + 50*( G4Pow::GetInstance()->powA( G4Log(0.73/proj_momentum), 3.5 ) ); 1932 } 1931 } 1933 else if( proj_momentum < 1.05 ) 1932 else if( proj_momentum < 1.05 ) 1934 { 1933 { 1935 hnXsc = 23 + 40*(G4Log(proj_momentu 1934 hnXsc = 23 + 40*(G4Log(proj_momentum/0.73))* 1936 (G4Log(proj_momentum 1935 (G4Log(proj_momentum/0.73)); 1937 } 1936 } 1938 else // if( proj_momentum < 10. ) 1937 else // if( proj_momentum < 10. ) 1939 { 1938 { 1940 hnXsc = 39.0 + 1939 hnXsc = 39.0 + 1941 75*(proj_momentum - 1.2)/(G4Pow 1940 75*(proj_momentum - 1.2)/(G4Pow::GetInstance()->powA(proj_momentum,3.0) + 0.15); 1942 } 1941 } 1943 xsection = hnXsc; 1942 xsection = hnXsc; 1944 } 1943 } 1945 else // pn to be np 1944 else // pn to be np 1946 { 1945 { 1947 if( proj_momentum < 0.8 ) 1946 if( proj_momentum < 0.8 ) 1948 { 1947 { 1949 hpXsc = 33+30*G4Pow::GetInstance()- 1948 hpXsc = 33+30*G4Pow::GetInstance()->powA(G4Log(proj_momentum/1.3),4.0); 1950 } 1949 } 1951 else if( proj_momentum < 1.4 ) 1950 else if( proj_momentum < 1.4 ) 1952 { 1951 { 1953 hpXsc = 33+30*G4Pow::GetInstance()- 1952 hpXsc = 33+30*G4Pow::GetInstance()->powA(G4Log(proj_momentum/0.95),2.0); 1954 } 1953 } 1955 else // if( proj_momentum < 10. ) 1954 else // if( proj_momentum < 10. ) 1956 { 1955 { 1957 hpXsc = 33.3+ 1956 hpXsc = 33.3+ 1958 20.8*(G4Pow::GetInstance()->pow 1957 20.8*(G4Pow::GetInstance()->powA(proj_momentum,2.0)-1.35)/ 1959 (G4Pow::GetInstance()->powA( 1958 (G4Pow::GetInstance()->powA(proj_momentum,2.50)+0.95); 1960 } 1959 } 1961 xsection = hpXsc; 1960 xsection = hpXsc; 1962 } 1961 } 1963 } 1962 } 1964 xsection *= CLHEP::millibarn; // parametris 1963 xsection *= CLHEP::millibarn; // parametrised in mb 1965 G4cout<<"xsection = "<<xsection/CLHEP::mill 1964 G4cout<<"xsection = "<<xsection/CLHEP::millibarn<<" mb"<<G4endl; 1966 return xsection; 1965 return xsection; 1967 } 1966 } 1968 1967 1969 ///////////////////////////////////////////// 1968 ///////////////////////////////////////////////////////////////// 1970 // 1969 // 1971 // The ratio el/ruth for Fresnel smooth nucle 1970 // The ratio el/ruth for Fresnel smooth nucleus profile 1972 1971 1973 G4double G4NuclNuclDiffuseElastic::GetRatioGe 1972 G4double G4NuclNuclDiffuseElastic::GetRatioGen(G4double theta) 1974 { 1973 { 1975 G4double sinThetaR = 2.*fHalfRutThetaTg/(1 1974 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2); 1976 G4double dTheta = 0.5*(theta - fRutherf 1975 G4double dTheta = 0.5*(theta - fRutherfordTheta); 1977 G4double sindTheta = std::sin(dTheta); 1976 G4double sindTheta = std::sin(dTheta); 1978 1977 1979 G4double prof = Profile(theta); 1978 G4double prof = Profile(theta); 1980 G4double prof2 = prof*prof; 1979 G4double prof2 = prof*prof; 1981 // G4double profmod = std::abs(prof); 1980 // G4double profmod = std::abs(prof); 1982 G4double order = std::sqrt(fProfileLam 1981 G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta; 1983 1982 1984 order = std::abs(order); // since sin chang 1983 order = std::abs(order); // since sin changes sign! 1985 // G4cout<<"order = "<<order<<G4endl; 1984 // G4cout<<"order = "<<order<<G4endl; 1986 1985 1987 G4double cosFresnel = GetCint(order); 1986 G4double cosFresnel = GetCint(order); 1988 G4double sinFresnel = GetSint(order); 1987 G4double sinFresnel = GetSint(order); 1989 1988 1990 G4double out; 1989 G4double out; 1991 1990 1992 if ( theta <= fRutherfordTheta ) 1991 if ( theta <= fRutherfordTheta ) 1993 { 1992 { 1994 out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-c 1993 out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2; 1995 out += ( cosFresnel + sinFresnel - 1. )*p 1994 out += ( cosFresnel + sinFresnel - 1. )*prof; 1996 } 1995 } 1997 else 1996 else 1998 { 1997 { 1999 out = 0.5*( (0.5-cosFresnel)*(0.5-cosFres 1998 out = 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2; 2000 } 1999 } 2001 2000 2002 return out; 2001 return out; 2003 } 2002 } 2004 2003 2005 ///////////////////////////////////////////// 2004 /////////////////////////////////////////////////////////////////// 2006 // 2005 // 2007 // For the calculation of arg Gamma(z) one ne 2006 // For the calculation of arg Gamma(z) one needs complex extension 2008 // of ln(Gamma(z)) 2007 // of ln(Gamma(z)) 2009 2008 2010 G4complex G4NuclNuclDiffuseElastic::GammaLoga 2009 G4complex G4NuclNuclDiffuseElastic::GammaLogarithm(G4complex zz) 2011 { 2010 { 2012 const G4double cof[6] = { 76.18009172947146 2011 const G4double cof[6] = { 76.18009172947146, -86.50532032941677, 2013 24.0140982408309 2012 24.01409824083091, -1.231739572450155, 2014 0.1208650973866 2013 0.1208650973866179e-2, -0.5395239384953e-5 } ; 2015 G4int j; 2014 G4int j; 2016 G4complex z = zz - 1.0; 2015 G4complex z = zz - 1.0; 2017 G4complex tmp = z + 5.5; 2016 G4complex tmp = z + 5.5; 2018 tmp -= (z + 0.5) * std::log(tmp); 2017 tmp -= (z + 0.5) * std::log(tmp); 2019 G4complex ser = G4complex(1.000000000190015 2018 G4complex ser = G4complex(1.000000000190015,0.); 2020 2019 2021 for ( j = 0; j <= 5; j++ ) 2020 for ( j = 0; j <= 5; j++ ) 2022 { 2021 { 2023 z += 1.0; 2022 z += 1.0; 2024 ser += cof[j]/z; 2023 ser += cof[j]/z; 2025 } 2024 } 2026 return -tmp + std::log(2.5066282746310005*s 2025 return -tmp + std::log(2.5066282746310005*ser); 2027 } 2026 } 2028 2027 2029 ///////////////////////////////////////////// 2028 ///////////////////////////////////////////////////////////// 2030 // 2029 // 2031 // Bessel J0 function based on rational appro 2030 // Bessel J0 function based on rational approximation from 2032 // J.F. Hart, Computer Approximations, New Yo 2031 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141 2033 2032 2034 G4double G4NuclNuclDiffuseElastic::BesselJzer 2033 G4double G4NuclNuclDiffuseElastic::BesselJzero(G4double value) 2035 { 2034 { 2036 G4double modvalue, value2, fact1, fact2, ar 2035 G4double modvalue, value2, fact1, fact2, arg, shift, bessel; 2037 2036 2038 modvalue = std::fabs(value); 2037 modvalue = std::fabs(value); 2039 2038 2040 if ( value < 8.0 && value > -8.0 ) 2039 if ( value < 8.0 && value > -8.0 ) 2041 { 2040 { 2042 value2 = value*value; 2041 value2 = value*value; 2043 2042 2044 fact1 = 57568490574.0 + value2*(-1336259 2043 fact1 = 57568490574.0 + value2*(-13362590354.0 2045 + value2*( 6516196 2044 + value2*( 651619640.7 2046 + value2*(-1121442 2045 + value2*(-11214424.18 2047 + value2*( 77392.3 2046 + value2*( 77392.33017 2048 + value2*(-184.905 2047 + value2*(-184.9052456 ) ) ) ) ); 2049 2048 2050 fact2 = 57568490411.0 + value2*( 1029532 2049 fact2 = 57568490411.0 + value2*( 1029532985.0 2051 + value2*( 9494680 2050 + value2*( 9494680.718 2052 + value2*(59272.64 2051 + value2*(59272.64853 2053 + value2*(267.8532 2052 + value2*(267.8532712 2054 + value2*1.0 2053 + value2*1.0 ) ) ) ); 2055 2054 2056 bessel = fact1/fact2; 2055 bessel = fact1/fact2; 2057 } 2056 } 2058 else 2057 else 2059 { 2058 { 2060 arg = 8.0/modvalue; 2059 arg = 8.0/modvalue; 2061 2060 2062 value2 = arg*arg; 2061 value2 = arg*arg; 2063 2062 2064 shift = modvalue-0.785398164; 2063 shift = modvalue-0.785398164; 2065 2064 2066 fact1 = 1.0 + value2*(-0.1098628627e-2 2065 fact1 = 1.0 + value2*(-0.1098628627e-2 2067 + value2*(0.2734510407e-4 2066 + value2*(0.2734510407e-4 2068 + value2*(-0.2073370639e-5 2067 + value2*(-0.2073370639e-5 2069 + value2*0.2093887211e-6 2068 + value2*0.2093887211e-6 ) ) ); 2070 2069 2071 fact2 = -0.1562499995e-1 + value2*(0.143 2070 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3 2072 + value2*(-0.69 2071 + value2*(-0.6911147651e-5 2073 + value2*(0.762 2072 + value2*(0.7621095161e-6 2074 - value2*0.9349 2073 - value2*0.934945152e-7 ) ) ); 2075 2074 2076 bessel = std::sqrt(0.636619772/modvalue)* 2075 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 ); 2077 } 2076 } 2078 return bessel; 2077 return bessel; 2079 } 2078 } 2080 2079 2081 ///////////////////////////////////////////// 2080 ///////////////////////////////////////////////////////////// 2082 // 2081 // 2083 // Bessel J1 function based on rational appro 2082 // Bessel J1 function based on rational approximation from 2084 // J.F. Hart, Computer Approximations, New Yo 2083 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141 2085 2084 2086 G4double G4NuclNuclDiffuseElastic::BesselJone 2085 G4double G4NuclNuclDiffuseElastic::BesselJone(G4double value) 2087 { 2086 { 2088 G4double modvalue, value2, fact1, fact2, ar 2087 G4double modvalue, value2, fact1, fact2, arg, shift, bessel; 2089 2088 2090 modvalue = std::fabs(value); 2089 modvalue = std::fabs(value); 2091 2090 2092 if ( modvalue < 8.0 ) 2091 if ( modvalue < 8.0 ) 2093 { 2092 { 2094 value2 = value*value; 2093 value2 = value*value; 2095 2094 2096 fact1 = value*(72362614232.0 + value2*(- 2095 fact1 = value*(72362614232.0 + value2*(-7895059235.0 2097 + value2*( 2096 + value2*( 242396853.1 2098 + value2*(- 2097 + value2*(-2972611.439 2099 + value2*( 2098 + value2*( 15704.48260 2100 + value2*(- 2099 + value2*(-30.16036606 ) ) ) ) ) ); 2101 2100 2102 fact2 = 144725228442.0 + value2*(2300535 2101 fact2 = 144725228442.0 + value2*(2300535178.0 2103 + value2*(1858330 2102 + value2*(18583304.74 2104 + value2*(99447.4 2103 + value2*(99447.43394 2105 + value2*(376.999 2104 + value2*(376.9991397 2106 + value2*1.0 2105 + value2*1.0 ) ) ) ); 2107 bessel = fact1/fact2; 2106 bessel = fact1/fact2; 2108 } 2107 } 2109 else 2108 else 2110 { 2109 { 2111 arg = 8.0/modvalue; 2110 arg = 8.0/modvalue; 2112 2111 2113 value2 = arg*arg; 2112 value2 = arg*arg; 2114 2113 2115 shift = modvalue - 2.356194491; 2114 shift = modvalue - 2.356194491; 2116 2115 2117 fact1 = 1.0 + value2*( 0.183105e-2 2116 fact1 = 1.0 + value2*( 0.183105e-2 2118 + value2*(-0.3516396496e-4 2117 + value2*(-0.3516396496e-4 2119 + value2*(0.2457520174e-5 2118 + value2*(0.2457520174e-5 2120 + value2*(-0.240337019e-6 2119 + value2*(-0.240337019e-6 ) ) ) ); 2121 2120 2122 fact2 = 0.04687499995 + value2*(-0.200269 2121 fact2 = 0.04687499995 + value2*(-0.2002690873e-3 2123 + value2*( 0.844919 2122 + value2*( 0.8449199096e-5 2124 + value2*(-0.882289 2123 + value2*(-0.88228987e-6 2125 + value2*0.10578741 2124 + value2*0.105787412e-6 ) ) ); 2126 2125 2127 bessel = std::sqrt( 0.636619772/modvalue) 2126 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2); 2128 2127 2129 if (value < 0.0) bessel = -bessel; 2128 if (value < 0.0) bessel = -bessel; 2130 } 2129 } 2131 return bessel; 2130 return bessel; 2132 } 2131 } 2133 2132 2134 // 2133 // 2135 // 2134 // 2136 ///////////////////////////////////////////// 2135 ///////////////////////////////////////////////////////////////////////////////// 2137 2136