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 // $Id: G4DiffuseElastic.cc 93440 2015-10-22 14:11:41Z gcosmo $ 26 // 27 // 27 // 28 // 28 // Physics model class G4DiffuseElastic 29 // Physics model class G4DiffuseElastic 29 // 30 // 30 // 31 // 31 // G4 Model: optical diffuse elastic scatterin 32 // G4 Model: optical diffuse elastic scattering with 4-momentum balance 32 // 33 // 33 // 24-May-07 V. Grichine 34 // 24-May-07 V. Grichine 34 // 35 // 35 // 21.10.15 V. Grichine 36 // 21.10.15 V. Grichine 36 // Bug fixed in BuildAngleTable, i 37 // Bug fixed in BuildAngleTable, improving accuracy for 37 // angle bins at high energies > 5 38 // angle bins at high energies > 50 GeV for pions. 38 // 39 // 39 40 40 #include "G4DiffuseElastic.hh" 41 #include "G4DiffuseElastic.hh" 41 #include "G4ParticleTable.hh" 42 #include "G4ParticleTable.hh" 42 #include "G4ParticleDefinition.hh" 43 #include "G4ParticleDefinition.hh" 43 #include "G4IonTable.hh" 44 #include "G4IonTable.hh" 44 #include "G4NucleiProperties.hh" 45 #include "G4NucleiProperties.hh" 45 46 46 #include "Randomize.hh" 47 #include "Randomize.hh" 47 #include "G4Integrator.hh" 48 #include "G4Integrator.hh" 48 #include "globals.hh" 49 #include "globals.hh" 49 #include "G4PhysicalConstants.hh" 50 #include "G4PhysicalConstants.hh" 50 #include "G4SystemOfUnits.hh" 51 #include "G4SystemOfUnits.hh" 51 52 52 #include "G4Proton.hh" 53 #include "G4Proton.hh" 53 #include "G4Neutron.hh" 54 #include "G4Neutron.hh" 54 #include "G4Deuteron.hh" 55 #include "G4Deuteron.hh" 55 #include "G4Alpha.hh" 56 #include "G4Alpha.hh" 56 #include "G4PionPlus.hh" 57 #include "G4PionPlus.hh" 57 #include "G4PionMinus.hh" 58 #include "G4PionMinus.hh" 58 59 59 #include "G4Element.hh" 60 #include "G4Element.hh" 60 #include "G4ElementTable.hh" 61 #include "G4ElementTable.hh" 61 #include "G4NistManager.hh" 62 #include "G4NistManager.hh" 62 #include "G4PhysicsTable.hh" 63 #include "G4PhysicsTable.hh" 63 #include "G4PhysicsLogVector.hh" 64 #include "G4PhysicsLogVector.hh" 64 #include "G4PhysicsFreeVector.hh" 65 #include "G4PhysicsFreeVector.hh" 65 66 66 #include "G4Exp.hh" 67 #include "G4Exp.hh" 67 68 68 #include "G4HadronicParameters.hh" << 69 << 70 ////////////////////////////////////////////// 69 ///////////////////////////////////////////////////////////////////////// 71 // 70 // 72 // Test Constructor. Just to check xsc 71 // Test Constructor. Just to check xsc 73 72 74 73 75 G4DiffuseElastic::G4DiffuseElastic() 74 G4DiffuseElastic::G4DiffuseElastic() 76 : G4HadronElastic("DiffuseElastic"), fPartic 75 : G4HadronElastic("DiffuseElastic"), fParticle(0) 77 { 76 { 78 SetMinEnergy( 0.01*MeV ); << 77 SetMinEnergy( 0.01*MeV ); // 0.01*GeV ); 79 SetMaxEnergy( G4HadronicParameters::Instance << 78 SetMaxEnergy( 1.*TeV ); 80 79 81 verboseLevel = 0; 80 verboseLevel = 0; 82 lowEnergyRecoilLimit = 100.*keV; 81 lowEnergyRecoilLimit = 100.*keV; 83 lowEnergyLimitQ = 0.0*GeV; 82 lowEnergyLimitQ = 0.0*GeV; 84 lowEnergyLimitHE = 0.0*GeV; 83 lowEnergyLimitHE = 0.0*GeV; 85 lowestEnergyLimit = 0.0*keV; 84 lowestEnergyLimit = 0.0*keV; 86 plabLowLimit = 20.0*MeV; 85 plabLowLimit = 20.0*MeV; 87 86 88 theProton = G4Proton::Proton(); 87 theProton = G4Proton::Proton(); 89 theNeutron = G4Neutron::Neutron(); 88 theNeutron = G4Neutron::Neutron(); 90 theDeuteron = G4Deuteron::Deuteron(); 89 theDeuteron = G4Deuteron::Deuteron(); 91 theAlpha = G4Alpha::Alpha(); 90 theAlpha = G4Alpha::Alpha(); 92 thePionPlus = G4PionPlus::PionPlus(); 91 thePionPlus = G4PionPlus::PionPlus(); 93 thePionMinus = G4PionMinus::PionMinus(); 92 thePionMinus = G4PionMinus::PionMinus(); 94 93 95 fEnergyBin = 300; // Increased from the ori << 94 fEnergyBin = 200; 96 fAngleBin = 200; 95 fAngleBin = 200; 97 96 98 fEnergyVector = new G4PhysicsLogVector( the 97 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin ); 99 98 100 fAngleTable = 0; 99 fAngleTable = 0; 101 100 102 fParticle = 0; 101 fParticle = 0; 103 fWaveVector = 0.; 102 fWaveVector = 0.; 104 fAtomicWeight = 0.; 103 fAtomicWeight = 0.; 105 fAtomicNumber = 0.; 104 fAtomicNumber = 0.; 106 fNuclearRadius = 0.; 105 fNuclearRadius = 0.; 107 fBeta = 0.; 106 fBeta = 0.; 108 fZommerfeld = 0.; 107 fZommerfeld = 0.; 109 fAm = 0.; 108 fAm = 0.; 110 fAddCoulomb = false; 109 fAddCoulomb = false; 111 } 110 } 112 111 113 ////////////////////////////////////////////// 112 ////////////////////////////////////////////////////////////////////////////// 114 // 113 // 115 // Destructor 114 // Destructor 116 115 117 G4DiffuseElastic::~G4DiffuseElastic() 116 G4DiffuseElastic::~G4DiffuseElastic() 118 { 117 { 119 if ( fEnergyVector ) 118 if ( fEnergyVector ) 120 { 119 { 121 delete fEnergyVector; 120 delete fEnergyVector; 122 fEnergyVector = 0; 121 fEnergyVector = 0; 123 } 122 } 124 for ( std::vector<G4PhysicsTable*>::iterator 123 for ( std::vector<G4PhysicsTable*>::iterator it = fAngleBank.begin(); 125 it != fAngleBank.end(); ++it ) 124 it != fAngleBank.end(); ++it ) 126 { 125 { 127 if ( (*it) ) (*it)->clearAndDestroy(); 126 if ( (*it) ) (*it)->clearAndDestroy(); 128 127 129 delete *it; 128 delete *it; 130 *it = 0; 129 *it = 0; 131 } 130 } 132 fAngleTable = 0; 131 fAngleTable = 0; 133 } 132 } 134 133 135 ////////////////////////////////////////////// 134 ////////////////////////////////////////////////////////////////////////////// 136 // 135 // 137 // Initialisation for given particle using ele 136 // Initialisation for given particle using element table of application 138 137 139 void G4DiffuseElastic::Initialise() 138 void G4DiffuseElastic::Initialise() 140 { 139 { 141 140 142 // fEnergyVector = new G4PhysicsLogVector( t 141 // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin ); 143 142 144 const G4ElementTable* theElementTable = G4El 143 const G4ElementTable* theElementTable = G4Element::GetElementTable(); 145 144 146 std::size_t jEl, numOfEl = G4Element::GetNum << 145 size_t jEl, numOfEl = G4Element::GetNumberOfElements(); 147 146 148 for( jEl = 0; jEl < numOfEl; ++jEl) // appli 147 for( jEl = 0; jEl < numOfEl; ++jEl) // application element loop 149 { 148 { 150 fAtomicNumber = (*theElementTable)[jEl]->G 149 fAtomicNumber = (*theElementTable)[jEl]->GetZ(); // atomic number 151 fAtomicWeight = G4NistManager::Instance()- 150 fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( fAtomicNumber ) ); 152 fNuclearRadius = CalculateNuclearRad(fAtom 151 fNuclearRadius = CalculateNuclearRad(fAtomicWeight); 153 152 154 if( verboseLevel > 0 ) 153 if( verboseLevel > 0 ) 155 { 154 { 156 G4cout<<"G4DiffuseElastic::Initialise() 155 G4cout<<"G4DiffuseElastic::Initialise() the element: " 157 <<(*theElementTable)[jEl]->GetName()<<G4 156 <<(*theElementTable)[jEl]->GetName()<<G4endl; 158 } 157 } 159 fElementNumberVector.push_back(fAtomicNumb 158 fElementNumberVector.push_back(fAtomicNumber); 160 fElementNameVector.push_back((*theElementT 159 fElementNameVector.push_back((*theElementTable)[jEl]->GetName()); 161 160 162 BuildAngleTable(); 161 BuildAngleTable(); 163 fAngleBank.push_back(fAngleTable); 162 fAngleBank.push_back(fAngleTable); 164 } 163 } 165 return; 164 return; 166 } 165 } 167 166 168 ////////////////////////////////////////////// 167 //////////////////////////////////////////////////////////////////////////// 169 // 168 // 170 // return differential elastic cross section d 169 // return differential elastic cross section d(sigma)/d(omega) 171 170 172 G4double 171 G4double 173 G4DiffuseElastic::GetDiffuseElasticXsc( const 172 G4DiffuseElastic::GetDiffuseElasticXsc( const G4ParticleDefinition* particle, 174 G4doub 173 G4double theta, 175 G4double momentum, 174 G4double momentum, 176 G4doub 175 G4double A ) 177 { 176 { 178 fParticle = particle; 177 fParticle = particle; 179 fWaveVector = momentum/hbarc; 178 fWaveVector = momentum/hbarc; 180 fAtomicWeight = A; 179 fAtomicWeight = A; 181 fAddCoulomb = false; 180 fAddCoulomb = false; 182 fNuclearRadius = CalculateNuclearRad(A); 181 fNuclearRadius = CalculateNuclearRad(A); 183 182 184 G4double sigma = fNuclearRadius*fNuclearRadi 183 G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta); 185 184 186 return sigma; 185 return sigma; 187 } 186 } 188 187 189 ////////////////////////////////////////////// 188 //////////////////////////////////////////////////////////////////////////// 190 // 189 // 191 // return invariant differential elastic cross 190 // return invariant differential elastic cross section d(sigma)/d(tMand) 192 191 193 G4double 192 G4double 194 G4DiffuseElastic::GetInvElasticXsc( const G4Pa 193 G4DiffuseElastic::GetInvElasticXsc( const G4ParticleDefinition* particle, 195 G4doub 194 G4double tMand, 196 G4double plab, 195 G4double plab, 197 G4doub 196 G4double A, G4double Z ) 198 { 197 { 199 G4double m1 = particle->GetPDGMass(); 198 G4double m1 = particle->GetPDGMass(); 200 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla 199 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1)); 201 200 202 G4int iZ = static_cast<G4int>(Z+0.5); 201 G4int iZ = static_cast<G4int>(Z+0.5); 203 G4int iA = static_cast<G4int>(A+0.5); 202 G4int iA = static_cast<G4int>(A+0.5); 204 G4ParticleDefinition * theDef = 0; 203 G4ParticleDefinition * theDef = 0; 205 204 206 if (iZ == 1 && iA == 1) theDef = thePro 205 if (iZ == 1 && iA == 1) theDef = theProton; 207 else if (iZ == 1 && iA == 2) theDef = theDeu 206 else if (iZ == 1 && iA == 2) theDef = theDeuteron; 208 else if (iZ == 1 && iA == 3) theDef = G4Trit 207 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton(); 209 else if (iZ == 2 && iA == 3) theDef = G4He3: 208 else if (iZ == 2 && iA == 3) theDef = G4He3::He3(); 210 else if (iZ == 2 && iA == 4) theDef = theAlp 209 else if (iZ == 2 && iA == 4) theDef = theAlpha; 211 else theDef = G4ParticleTable::GetParticleTa 210 else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0); 212 211 213 G4double tmass = theDef->GetPDGMass(); 212 G4double tmass = theDef->GetPDGMass(); 214 213 215 G4LorentzVector lv(0.0,0.0,0.0,tmass); 214 G4LorentzVector lv(0.0,0.0,0.0,tmass); 216 lv += lv1; 215 lv += lv1; 217 216 218 G4ThreeVector bst = lv.boostVector(); 217 G4ThreeVector bst = lv.boostVector(); 219 lv1.boost(-bst); 218 lv1.boost(-bst); 220 219 221 G4ThreeVector p1 = lv1.vect(); 220 G4ThreeVector p1 = lv1.vect(); 222 G4double ptot = p1.mag(); 221 G4double ptot = p1.mag(); 223 G4double ptot2 = ptot*ptot; 222 G4double ptot2 = ptot*ptot; 224 G4double cost = 1 - 0.5*std::fabs(tMand)/pto 223 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2; 225 224 226 if( cost >= 1.0 ) cost = 1.0; 225 if( cost >= 1.0 ) cost = 1.0; 227 else if( cost <= -1.0) cost = -1.0; 226 else if( cost <= -1.0) cost = -1.0; 228 227 229 G4double thetaCMS = std::acos(cost); 228 G4double thetaCMS = std::acos(cost); 230 229 231 G4double sigma = GetDiffuseElasticXsc( parti 230 G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A); 232 231 233 sigma *= pi/ptot2; 232 sigma *= pi/ptot2; 234 233 235 return sigma; 234 return sigma; 236 } 235 } 237 236 238 ////////////////////////////////////////////// 237 //////////////////////////////////////////////////////////////////////////// 239 // 238 // 240 // return differential elastic cross section d 239 // return differential elastic cross section d(sigma)/d(omega) with Coulomb 241 // correction 240 // correction 242 241 243 G4double 242 G4double 244 G4DiffuseElastic::GetDiffuseElasticSumXsc( con 243 G4DiffuseElastic::GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle, 245 G4doub 244 G4double theta, 246 G4double momentum, 245 G4double momentum, 247 G4doub 246 G4double A, G4double Z ) 248 { 247 { 249 fParticle = particle; 248 fParticle = particle; 250 fWaveVector = momentum/hbarc; 249 fWaveVector = momentum/hbarc; 251 fAtomicWeight = A; 250 fAtomicWeight = A; 252 fAtomicNumber = Z; 251 fAtomicNumber = Z; 253 fNuclearRadius = CalculateNuclearRad(A); 252 fNuclearRadius = CalculateNuclearRad(A); 254 fAddCoulomb = false; 253 fAddCoulomb = false; 255 254 256 G4double z = particle->GetPDGCharge(); 255 G4double z = particle->GetPDGCharge(); 257 256 258 G4double kRt = fWaveVector*fNuclearRadius* 257 G4double kRt = fWaveVector*fNuclearRadius*theta; 259 G4double kRtC = 1.9; 258 G4double kRtC = 1.9; 260 259 261 if( z && (kRt > kRtC) ) 260 if( z && (kRt > kRtC) ) 262 { 261 { 263 fAddCoulomb = true; 262 fAddCoulomb = true; 264 fBeta = CalculateParticleBeta( parti 263 fBeta = CalculateParticleBeta( particle, momentum); 265 fZommerfeld = CalculateZommerfeld( fBeta, 264 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber); 266 fAm = CalculateAm( momentum, fZomm 265 fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber); 267 } 266 } 268 G4double sigma = fNuclearRadius*fNuclearRadi 267 G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta); 269 268 270 return sigma; 269 return sigma; 271 } 270 } 272 271 273 ////////////////////////////////////////////// 272 //////////////////////////////////////////////////////////////////////////// 274 // 273 // 275 // return invariant differential elastic cross 274 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb 276 // correction 275 // correction 277 276 278 G4double 277 G4double 279 G4DiffuseElastic::GetInvElasticSumXsc( const G 278 G4DiffuseElastic::GetInvElasticSumXsc( const G4ParticleDefinition* particle, 280 G4doub 279 G4double tMand, 281 G4double plab, 280 G4double plab, 282 G4doub 281 G4double A, G4double Z ) 283 { 282 { 284 G4double m1 = particle->GetPDGMass(); 283 G4double m1 = particle->GetPDGMass(); 285 284 286 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla 285 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1)); 287 286 288 G4int iZ = static_cast<G4int>(Z+0.5); 287 G4int iZ = static_cast<G4int>(Z+0.5); 289 G4int iA = static_cast<G4int>(A+0.5); 288 G4int iA = static_cast<G4int>(A+0.5); 290 289 291 G4ParticleDefinition* theDef = 0; 290 G4ParticleDefinition* theDef = 0; 292 291 293 if (iZ == 1 && iA == 1) theDef = thePro 292 if (iZ == 1 && iA == 1) theDef = theProton; 294 else if (iZ == 1 && iA == 2) theDef = theDeu 293 else if (iZ == 1 && iA == 2) theDef = theDeuteron; 295 else if (iZ == 1 && iA == 3) theDef = G4Trit 294 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton(); 296 else if (iZ == 2 && iA == 3) theDef = G4He3: 295 else if (iZ == 2 && iA == 3) theDef = G4He3::He3(); 297 else if (iZ == 2 && iA == 4) theDef = theAlp 296 else if (iZ == 2 && iA == 4) theDef = theAlpha; 298 else theDef = G4ParticleTable::GetParticleTa 297 else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0); 299 298 300 G4double tmass = theDef->GetPDGMass(); 299 G4double tmass = theDef->GetPDGMass(); 301 300 302 G4LorentzVector lv(0.0,0.0,0.0,tmass); 301 G4LorentzVector lv(0.0,0.0,0.0,tmass); 303 lv += lv1; 302 lv += lv1; 304 303 305 G4ThreeVector bst = lv.boostVector(); 304 G4ThreeVector bst = lv.boostVector(); 306 lv1.boost(-bst); 305 lv1.boost(-bst); 307 306 308 G4ThreeVector p1 = lv1.vect(); 307 G4ThreeVector p1 = lv1.vect(); 309 G4double ptot = p1.mag(); 308 G4double ptot = p1.mag(); 310 G4double ptot2 = ptot*ptot; 309 G4double ptot2 = ptot*ptot; 311 G4double cost = 1 - 0.5*std::fabs(tMand)/ 310 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2; 312 311 313 if( cost >= 1.0 ) cost = 1.0; 312 if( cost >= 1.0 ) cost = 1.0; 314 else if( cost <= -1.0) cost = -1.0; 313 else if( cost <= -1.0) cost = -1.0; 315 314 316 G4double thetaCMS = std::acos(cost); 315 G4double thetaCMS = std::acos(cost); 317 316 318 G4double sigma = GetDiffuseElasticSumXsc( pa 317 G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z ); 319 318 320 sigma *= pi/ptot2; 319 sigma *= pi/ptot2; 321 320 322 return sigma; 321 return sigma; 323 } 322 } 324 323 325 ////////////////////////////////////////////// 324 //////////////////////////////////////////////////////////////////////////// 326 // 325 // 327 // return invariant differential elastic cross 326 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb 328 // correction 327 // correction 329 328 330 G4double 329 G4double 331 G4DiffuseElastic::GetInvCoulombElasticXsc( con 330 G4DiffuseElastic::GetInvCoulombElasticXsc( const G4ParticleDefinition* particle, 332 G4doub 331 G4double tMand, 333 G4double plab, 332 G4double plab, 334 G4doub 333 G4double A, G4double Z ) 335 { 334 { 336 G4double m1 = particle->GetPDGMass(); 335 G4double m1 = particle->GetPDGMass(); 337 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla 336 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1)); 338 337 339 G4int iZ = static_cast<G4int>(Z+0.5); 338 G4int iZ = static_cast<G4int>(Z+0.5); 340 G4int iA = static_cast<G4int>(A+0.5); 339 G4int iA = static_cast<G4int>(A+0.5); 341 G4ParticleDefinition * theDef = 0; 340 G4ParticleDefinition * theDef = 0; 342 341 343 if (iZ == 1 && iA == 1) theDef = thePro 342 if (iZ == 1 && iA == 1) theDef = theProton; 344 else if (iZ == 1 && iA == 2) theDef = theDeu 343 else if (iZ == 1 && iA == 2) theDef = theDeuteron; 345 else if (iZ == 1 && iA == 3) theDef = G4Trit 344 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton(); 346 else if (iZ == 2 && iA == 3) theDef = G4He3: 345 else if (iZ == 2 && iA == 3) theDef = G4He3::He3(); 347 else if (iZ == 2 && iA == 4) theDef = theAlp 346 else if (iZ == 2 && iA == 4) theDef = theAlpha; 348 else theDef = G4ParticleTable::GetParticleTa 347 else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0); 349 348 350 G4double tmass = theDef->GetPDGMass(); 349 G4double tmass = theDef->GetPDGMass(); 351 350 352 G4LorentzVector lv(0.0,0.0,0.0,tmass); 351 G4LorentzVector lv(0.0,0.0,0.0,tmass); 353 lv += lv1; 352 lv += lv1; 354 353 355 G4ThreeVector bst = lv.boostVector(); 354 G4ThreeVector bst = lv.boostVector(); 356 lv1.boost(-bst); 355 lv1.boost(-bst); 357 356 358 G4ThreeVector p1 = lv1.vect(); 357 G4ThreeVector p1 = lv1.vect(); 359 G4double ptot = p1.mag(); 358 G4double ptot = p1.mag(); 360 G4double ptot2 = ptot*ptot; 359 G4double ptot2 = ptot*ptot; 361 G4double cost = 1 - 0.5*std::fabs(tMand)/pto 360 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2; 362 361 363 if( cost >= 1.0 ) cost = 1.0; 362 if( cost >= 1.0 ) cost = 1.0; 364 else if( cost <= -1.0) cost = -1.0; 363 else if( cost <= -1.0) cost = -1.0; 365 364 366 G4double thetaCMS = std::acos(cost); 365 G4double thetaCMS = std::acos(cost); 367 366 368 G4double sigma = GetCoulombElasticXsc( parti 367 G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z ); 369 368 370 sigma *= pi/ptot2; 369 sigma *= pi/ptot2; 371 370 372 return sigma; 371 return sigma; 373 } 372 } 374 373 375 ////////////////////////////////////////////// 374 //////////////////////////////////////////////////////////////////////////// 376 // 375 // 377 // return differential elastic probability d(p 376 // return differential elastic probability d(probability)/d(omega) 378 377 379 G4double 378 G4double 380 G4DiffuseElastic::GetDiffElasticProb( // G4Par 379 G4DiffuseElastic::GetDiffElasticProb( // G4ParticleDefinition* particle, 381 G4doub 380 G4double theta 382 // G4double momentum, 381 // G4double momentum, 383 // G4double A 382 // G4double A 384 ) 383 ) 385 { 384 { 386 G4double sigma, bzero, bzero2, bonebyarg, bo 385 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2; 387 G4double delta, diffuse, gamma; 386 G4double delta, diffuse, gamma; 388 G4double e1, e2, bone, bone2; 387 G4double e1, e2, bone, bone2; 389 388 390 // G4double wavek = momentum/hbarc; // wave 389 // G4double wavek = momentum/hbarc; // wave vector 391 // G4double r0 = 1.08*fermi; 390 // G4double r0 = 1.08*fermi; 392 // G4double rad = r0*G4Pow::GetInstance()- 391 // G4double rad = r0*G4Pow::GetInstance()->A13(A); 393 392 394 if (fParticle == theProton) 393 if (fParticle == theProton) 395 { 394 { 396 diffuse = 0.63*fermi; 395 diffuse = 0.63*fermi; 397 gamma = 0.3*fermi; 396 gamma = 0.3*fermi; 398 delta = 0.1*fermi*fermi; 397 delta = 0.1*fermi*fermi; 399 e1 = 0.3*fermi; 398 e1 = 0.3*fermi; 400 e2 = 0.35*fermi; 399 e2 = 0.35*fermi; 401 } 400 } 402 else if (fParticle == theNeutron) 401 else if (fParticle == theNeutron) 403 { 402 { 404 diffuse = 0.63*fermi; // 1.63*fermi; // 403 diffuse = 0.63*fermi; // 1.63*fermi; // 405 G4double k0 = 1*GeV/hbarc; 404 G4double k0 = 1*GeV/hbarc; 406 diffuse *= k0/fWaveVector; 405 diffuse *= k0/fWaveVector; 407 406 408 gamma = 0.3*fermi; 407 gamma = 0.3*fermi; 409 delta = 0.1*fermi*fermi; 408 delta = 0.1*fermi*fermi; 410 e1 = 0.3*fermi; 409 e1 = 0.3*fermi; 411 e2 = 0.35*fermi; 410 e2 = 0.35*fermi; 412 } 411 } 413 else // as proton, if were not defined 412 else // as proton, if were not defined 414 { 413 { 415 diffuse = 0.63*fermi; 414 diffuse = 0.63*fermi; 416 gamma = 0.3*fermi; 415 gamma = 0.3*fermi; 417 delta = 0.1*fermi*fermi; 416 delta = 0.1*fermi*fermi; 418 e1 = 0.3*fermi; 417 e1 = 0.3*fermi; 419 e2 = 0.35*fermi; 418 e2 = 0.35*fermi; 420 } 419 } 421 G4double kr = fWaveVector*fNuclearRadius; 420 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad; 422 G4double kr2 = kr*kr; 421 G4double kr2 = kr*kr; 423 G4double krt = kr*theta; 422 G4double krt = kr*theta; 424 423 425 bzero = BesselJzero(krt); 424 bzero = BesselJzero(krt); 426 bzero2 = bzero*bzero; 425 bzero2 = bzero*bzero; 427 bone = BesselJone(krt); 426 bone = BesselJone(krt); 428 bone2 = bone*bone; 427 bone2 = bone*bone; 429 bonebyarg = BesselOneByArg(krt); 428 bonebyarg = BesselOneByArg(krt); 430 bonebyarg2 = bonebyarg*bonebyarg; 429 bonebyarg2 = bonebyarg*bonebyarg; 431 430 432 G4double lambda = 15.; // 15 ok 431 G4double lambda = 15.; // 15 ok 433 432 434 // G4double kgamma = fWaveVector*gamma; 433 // G4double kgamma = fWaveVector*gamma; // wavek*delta; 435 434 436 G4double kgamma = lambda*(1.-G4Exp(-fWave 435 G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta; 437 G4double kgamma2 = kgamma*kgamma; 436 G4double kgamma2 = kgamma*kgamma; 438 437 439 // G4double dk2t = delta*fWaveVector*fWaveV 438 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta; 440 // G4double dk2t2 = dk2t*dk2t; 439 // G4double dk2t2 = dk2t*dk2t; 441 // G4double pikdt = pi*fWaveVector*diffuse*t 440 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta; 442 441 443 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa 442 G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta; 444 443 445 damp = DampFactor(pikdt); 444 damp = DampFactor(pikdt); 446 damp2 = damp*damp; 445 damp2 = damp*damp; 447 446 448 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector 447 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector; 449 G4double e2dk3t = -2.*e2*delta*fWaveVector* 448 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta; 450 449 451 450 452 sigma = kgamma2; 451 sigma = kgamma2; 453 // sigma += dk2t2; 452 // sigma += dk2t2; 454 sigma *= bzero2; 453 sigma *= bzero2; 455 sigma += mode2k2*bone2 + e2dk3t*bzero*bone; 454 sigma += mode2k2*bone2 + e2dk3t*bzero*bone; 456 sigma += kr2*bonebyarg2; 455 sigma += kr2*bonebyarg2; 457 sigma *= damp2; // *rad*rad; 456 sigma *= damp2; // *rad*rad; 458 457 459 return sigma; 458 return sigma; 460 } 459 } 461 460 462 ////////////////////////////////////////////// 461 //////////////////////////////////////////////////////////////////////////// 463 // 462 // 464 // return differential elastic probability d(p 463 // return differential elastic probability d(probability)/d(omega) with 465 // Coulomb correction 464 // Coulomb correction 466 465 467 G4double 466 G4double 468 G4DiffuseElastic::GetDiffElasticSumProb( // G4 467 G4DiffuseElastic::GetDiffElasticSumProb( // G4ParticleDefinition* particle, 469 G4doub 468 G4double theta 470 // G4double momentum, 469 // G4double momentum, 471 // G4double A 470 // G4double A 472 ) 471 ) 473 { 472 { 474 G4double sigma, bzero, bzero2, bonebyarg, bo 473 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2; 475 G4double delta, diffuse, gamma; 474 G4double delta, diffuse, gamma; 476 G4double e1, e2, bone, bone2; 475 G4double e1, e2, bone, bone2; 477 476 478 // G4double wavek = momentum/hbarc; // wave 477 // G4double wavek = momentum/hbarc; // wave vector 479 // G4double r0 = 1.08*fermi; 478 // G4double r0 = 1.08*fermi; 480 // G4double rad = r0*G4Pow::GetInstance()- 479 // G4double rad = r0*G4Pow::GetInstance()->A13(A); 481 480 482 G4double kr = fWaveVector*fNuclearRadius; 481 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad; 483 G4double kr2 = kr*kr; 482 G4double kr2 = kr*kr; 484 G4double krt = kr*theta; 483 G4double krt = kr*theta; 485 484 486 bzero = BesselJzero(krt); 485 bzero = BesselJzero(krt); 487 bzero2 = bzero*bzero; 486 bzero2 = bzero*bzero; 488 bone = BesselJone(krt); 487 bone = BesselJone(krt); 489 bone2 = bone*bone; 488 bone2 = bone*bone; 490 bonebyarg = BesselOneByArg(krt); 489 bonebyarg = BesselOneByArg(krt); 491 bonebyarg2 = bonebyarg*bonebyarg; 490 bonebyarg2 = bonebyarg*bonebyarg; 492 491 493 if (fParticle == theProton) 492 if (fParticle == theProton) 494 { 493 { 495 diffuse = 0.63*fermi; 494 diffuse = 0.63*fermi; 496 // diffuse = 0.6*fermi; 495 // diffuse = 0.6*fermi; 497 gamma = 0.3*fermi; 496 gamma = 0.3*fermi; 498 delta = 0.1*fermi*fermi; 497 delta = 0.1*fermi*fermi; 499 e1 = 0.3*fermi; 498 e1 = 0.3*fermi; 500 e2 = 0.35*fermi; 499 e2 = 0.35*fermi; 501 } 500 } 502 else if (fParticle == theNeutron) 501 else if (fParticle == theNeutron) 503 { 502 { 504 diffuse = 0.63*fermi; 503 diffuse = 0.63*fermi; 505 // diffuse = 0.6*fermi; 504 // diffuse = 0.6*fermi; 506 G4double k0 = 1*GeV/hbarc; 505 G4double k0 = 1*GeV/hbarc; 507 diffuse *= k0/fWaveVector; 506 diffuse *= k0/fWaveVector; 508 gamma = 0.3*fermi; 507 gamma = 0.3*fermi; 509 delta = 0.1*fermi*fermi; 508 delta = 0.1*fermi*fermi; 510 e1 = 0.3*fermi; 509 e1 = 0.3*fermi; 511 e2 = 0.35*fermi; 510 e2 = 0.35*fermi; 512 } 511 } 513 else // as proton, if were not defined 512 else // as proton, if were not defined 514 { 513 { 515 diffuse = 0.63*fermi; 514 diffuse = 0.63*fermi; 516 gamma = 0.3*fermi; 515 gamma = 0.3*fermi; 517 delta = 0.1*fermi*fermi; 516 delta = 0.1*fermi*fermi; 518 e1 = 0.3*fermi; 517 e1 = 0.3*fermi; 519 e2 = 0.35*fermi; 518 e2 = 0.35*fermi; 520 } 519 } 521 G4double lambda = 15.; // 15 ok 520 G4double lambda = 15.; // 15 ok 522 // G4double kgamma = fWaveVector*gamma; 521 // G4double kgamma = fWaveVector*gamma; // wavek*delta; 523 G4double kgamma = lambda*(1.-G4Exp(-fWave 522 G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta; 524 523 525 // G4cout<<"kgamma = "<<kgamma<<G4endl; 524 // G4cout<<"kgamma = "<<kgamma<<G4endl; 526 525 527 if(fAddCoulomb) // add Coulomb correction 526 if(fAddCoulomb) // add Coulomb correction 528 { 527 { 529 G4double sinHalfTheta = std::sin(0.5*thet 528 G4double sinHalfTheta = std::sin(0.5*theta); 530 G4double sinHalfTheta2 = sinHalfTheta*sinH 529 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta; 531 530 532 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta 531 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0() 533 // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe 532 // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0() 534 } 533 } 535 534 536 G4double kgamma2 = kgamma*kgamma; 535 G4double kgamma2 = kgamma*kgamma; 537 536 538 // G4double dk2t = delta*fWaveVector*fWaveV 537 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta; 539 // G4cout<<"dk2t = "<<dk2t<<G4endl; 538 // G4cout<<"dk2t = "<<dk2t<<G4endl; 540 // G4double dk2t2 = dk2t*dk2t; 539 // G4double dk2t2 = dk2t*dk2t; 541 // G4double pikdt = pi*fWaveVector*diffuse*t 540 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta; 542 541 543 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa 542 G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta; 544 543 545 // G4cout<<"pikdt = "<<pikdt<<G4endl; 544 // G4cout<<"pikdt = "<<pikdt<<G4endl; 546 545 547 damp = DampFactor(pikdt); 546 damp = DampFactor(pikdt); 548 damp2 = damp*damp; 547 damp2 = damp*damp; 549 548 550 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector 549 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector; 551 G4double e2dk3t = -2.*e2*delta*fWaveVector* 550 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta; 552 551 553 sigma = kgamma2; 552 sigma = kgamma2; 554 // sigma += dk2t2; 553 // sigma += dk2t2; 555 sigma *= bzero2; 554 sigma *= bzero2; 556 sigma += mode2k2*bone2; 555 sigma += mode2k2*bone2; 557 sigma += e2dk3t*bzero*bone; 556 sigma += e2dk3t*bzero*bone; 558 557 559 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerf 558 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/() 560 sigma += kr2*bonebyarg2; // correction at J 559 sigma += kr2*bonebyarg2; // correction at J1()/() 561 560 562 sigma *= damp2; // *rad*rad; 561 sigma *= damp2; // *rad*rad; 563 562 564 return sigma; 563 return sigma; 565 } 564 } 566 565 567 566 568 ////////////////////////////////////////////// 567 //////////////////////////////////////////////////////////////////////////// 569 // 568 // 570 // return differential elastic probability d(p 569 // return differential elastic probability d(probability)/d(t) with 571 // Coulomb correction. It is called from Build 570 // Coulomb correction. It is called from BuildAngleTable() 572 571 573 G4double 572 G4double 574 G4DiffuseElastic::GetDiffElasticSumProbA( G4do 573 G4DiffuseElastic::GetDiffElasticSumProbA( G4double alpha ) 575 { 574 { 576 G4double theta; 575 G4double theta; 577 576 578 theta = std::sqrt(alpha); 577 theta = std::sqrt(alpha); 579 578 580 // theta = std::acos( 1 - alpha/2. ); 579 // theta = std::acos( 1 - alpha/2. ); 581 580 582 G4double sigma, bzero, bzero2, bonebyarg, bo 581 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2; 583 G4double delta, diffuse, gamma; 582 G4double delta, diffuse, gamma; 584 G4double e1, e2, bone, bone2; 583 G4double e1, e2, bone, bone2; 585 584 586 // G4double wavek = momentum/hbarc; // wave 585 // G4double wavek = momentum/hbarc; // wave vector 587 // G4double r0 = 1.08*fermi; 586 // G4double r0 = 1.08*fermi; 588 // G4double rad = r0*G4Pow::GetInstance()- 587 // G4double rad = r0*G4Pow::GetInstance()->A13(A); 589 588 590 G4double kr = fWaveVector*fNuclearRadius; 589 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad; 591 G4double kr2 = kr*kr; 590 G4double kr2 = kr*kr; 592 G4double krt = kr*theta; 591 G4double krt = kr*theta; 593 592 594 bzero = BesselJzero(krt); 593 bzero = BesselJzero(krt); 595 bzero2 = bzero*bzero; 594 bzero2 = bzero*bzero; 596 bone = BesselJone(krt); 595 bone = BesselJone(krt); 597 bone2 = bone*bone; 596 bone2 = bone*bone; 598 bonebyarg = BesselOneByArg(krt); 597 bonebyarg = BesselOneByArg(krt); 599 bonebyarg2 = bonebyarg*bonebyarg; 598 bonebyarg2 = bonebyarg*bonebyarg; 600 599 601 if ( fParticle == theProton ) 600 if ( fParticle == theProton ) 602 { 601 { 603 diffuse = 0.63*fermi; 602 diffuse = 0.63*fermi; 604 // diffuse = 0.6*fermi; 603 // diffuse = 0.6*fermi; 605 gamma = 0.3*fermi; 604 gamma = 0.3*fermi; 606 delta = 0.1*fermi*fermi; 605 delta = 0.1*fermi*fermi; 607 e1 = 0.3*fermi; 606 e1 = 0.3*fermi; 608 e2 = 0.35*fermi; 607 e2 = 0.35*fermi; 609 } 608 } 610 else if ( fParticle == theNeutron ) 609 else if ( fParticle == theNeutron ) 611 { 610 { 612 diffuse = 0.63*fermi; 611 diffuse = 0.63*fermi; 613 // diffuse = 0.6*fermi; 612 // diffuse = 0.6*fermi; 614 // G4double k0 = 0.8*GeV/hbarc; 613 // G4double k0 = 0.8*GeV/hbarc; 615 // diffuse *= k0/fWaveVector; 614 // diffuse *= k0/fWaveVector; 616 gamma = 0.3*fermi; 615 gamma = 0.3*fermi; 617 delta = 0.1*fermi*fermi; 616 delta = 0.1*fermi*fermi; 618 e1 = 0.3*fermi; 617 e1 = 0.3*fermi; 619 e2 = 0.35*fermi; 618 e2 = 0.35*fermi; 620 } 619 } 621 else // as proton, if were not defined 620 else // as proton, if were not defined 622 { 621 { 623 diffuse = 0.63*fermi; 622 diffuse = 0.63*fermi; 624 gamma = 0.3*fermi; 623 gamma = 0.3*fermi; 625 delta = 0.1*fermi*fermi; 624 delta = 0.1*fermi*fermi; 626 e1 = 0.3*fermi; 625 e1 = 0.3*fermi; 627 e2 = 0.35*fermi; 626 e2 = 0.35*fermi; 628 } 627 } 629 G4double lambda = 15.; // 15 ok 628 G4double lambda = 15.; // 15 ok 630 // G4double kgamma = fWaveVector*gamma; 629 // G4double kgamma = fWaveVector*gamma; // wavek*delta; 631 G4double kgamma = lambda*(1.-G4Exp(-fWave 630 G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta; 632 631 633 // G4cout<<"kgamma = "<<kgamma<<G4endl; 632 // G4cout<<"kgamma = "<<kgamma<<G4endl; 634 633 635 if( fAddCoulomb ) // add Coulomb correction 634 if( fAddCoulomb ) // add Coulomb correction 636 { 635 { 637 G4double sinHalfTheta = theta*0.5; // std 636 G4double sinHalfTheta = theta*0.5; // std::sin(0.5*theta); 638 G4double sinHalfTheta2 = sinHalfTheta*sinH 637 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta; 639 638 640 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta 639 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0() 641 // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe 640 // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0() 642 } 641 } 643 G4double kgamma2 = kgamma*kgamma; 642 G4double kgamma2 = kgamma*kgamma; 644 643 645 // G4double dk2t = delta*fWaveVector*fWaveV 644 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta; 646 // G4cout<<"dk2t = "<<dk2t<<G4endl; 645 // G4cout<<"dk2t = "<<dk2t<<G4endl; 647 // G4double dk2t2 = dk2t*dk2t; 646 // G4double dk2t2 = dk2t*dk2t; 648 // G4double pikdt = pi*fWaveVector*diffuse*t 647 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta; 649 648 650 G4double pikdt = lambda*(1. - G4Exp( -pi* 649 G4double pikdt = lambda*(1. - G4Exp( -pi*fWaveVector*diffuse*theta/lambda ) ); // wavek*delta; 651 650 652 // G4cout<<"pikdt = "<<pikdt<<G4endl; 651 // G4cout<<"pikdt = "<<pikdt<<G4endl; 653 652 654 damp = DampFactor( pikdt ); 653 damp = DampFactor( pikdt ); 655 damp2 = damp*damp; 654 damp2 = damp*damp; 656 655 657 G4double mode2k2 = ( e1*e1 + e2*e2 )*fWaveVe 656 G4double mode2k2 = ( e1*e1 + e2*e2 )*fWaveVector*fWaveVector; 658 G4double e2dk3t = -2.*e2*delta*fWaveVector* 657 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta; 659 658 660 sigma = kgamma2; 659 sigma = kgamma2; 661 // sigma += dk2t2; 660 // sigma += dk2t2; 662 sigma *= bzero2; 661 sigma *= bzero2; 663 sigma += mode2k2*bone2; 662 sigma += mode2k2*bone2; 664 sigma += e2dk3t*bzero*bone; 663 sigma += e2dk3t*bzero*bone; 665 664 666 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerf 665 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/() 667 sigma += kr2*bonebyarg2; // correction at J 666 sigma += kr2*bonebyarg2; // correction at J1()/() 668 667 669 sigma *= damp2; // *rad*rad; 668 sigma *= damp2; // *rad*rad; 670 669 671 return sigma; 670 return sigma; 672 } 671 } 673 672 674 673 675 ////////////////////////////////////////////// 674 //////////////////////////////////////////////////////////////////////////// 676 // 675 // 677 // return differential elastic probability 2*p 676 // return differential elastic probability 2*pi*sin(theta)*d(probability)/d(omega) 678 677 679 G4double 678 G4double 680 G4DiffuseElastic::GetIntegrandFunction( G4doub 679 G4DiffuseElastic::GetIntegrandFunction( G4double alpha ) 681 { 680 { 682 G4double result; 681 G4double result; 683 682 684 result = GetDiffElasticSumProbA(alpha); 683 result = GetDiffElasticSumProbA(alpha); 685 684 686 // result *= 2*pi*std::sin(theta); 685 // result *= 2*pi*std::sin(theta); 687 686 688 return result; 687 return result; 689 } 688 } 690 689 691 ////////////////////////////////////////////// 690 //////////////////////////////////////////////////////////////////////////// 692 // 691 // 693 // return integral elastic cross section d(sig 692 // return integral elastic cross section d(sigma)/d(omega) integrated 0 - theta 694 693 695 G4double 694 G4double 696 G4DiffuseElastic::IntegralElasticProb( const 695 G4DiffuseElastic::IntegralElasticProb( const G4ParticleDefinition* particle, 697 G4doub 696 G4double theta, 698 G4double momentum, 697 G4double momentum, 699 G4doub 698 G4double A ) 700 { 699 { 701 G4double result; 700 G4double result; 702 fParticle = particle; 701 fParticle = particle; 703 fWaveVector = momentum/hbarc; 702 fWaveVector = momentum/hbarc; 704 fAtomicWeight = A; 703 fAtomicWeight = A; 705 704 706 fNuclearRadius = CalculateNuclearRad(A); 705 fNuclearRadius = CalculateNuclearRad(A); 707 706 708 707 709 G4Integrator<G4DiffuseElastic,G4double(G4Dif 708 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral; 710 709 711 // result = integral.Legendre10(this,&G4Diff 710 // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta ); 712 result = integral.Legendre96(this,&G4Diffuse 711 result = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta ); 713 712 714 return result; 713 return result; 715 } 714 } 716 715 717 ////////////////////////////////////////////// 716 //////////////////////////////////////////////////////////////////////////// 718 // 717 // 719 // Return inv momentum transfer -t > 0 718 // Return inv momentum transfer -t > 0 720 719 721 G4double G4DiffuseElastic::SampleT( const G4Pa 720 G4double G4DiffuseElastic::SampleT( const G4ParticleDefinition* aParticle, G4double p, G4double A) 722 { 721 { 723 G4double theta = SampleThetaCMS( aParticle, 722 G4double theta = SampleThetaCMS( aParticle, p, A); // sample theta in cms 724 G4double t = 2*p*p*( 1 - std::cos(theta) 723 G4double t = 2*p*p*( 1 - std::cos(theta) ); // -t !!! 725 return t; 724 return t; 726 } 725 } 727 726 728 ////////////////////////////////////////////// 727 //////////////////////////////////////////////////////////////////////////// 729 // 728 // 730 // Return scattering angle sampled in cms 729 // Return scattering angle sampled in cms 731 730 732 731 733 G4double 732 G4double 734 G4DiffuseElastic::SampleThetaCMS(const G4Parti 733 G4DiffuseElastic::SampleThetaCMS(const G4ParticleDefinition* particle, 735 G4doubl 734 G4double momentum, G4double A) 736 { 735 { 737 G4int i, iMax = 100; 736 G4int i, iMax = 100; 738 G4double norm, theta1, theta2, thetaMax; << 737 G4double norm, result, theta1, theta2, thetaMax, sum = 0.; 739 G4double result = 0., sum = 0.; << 740 738 741 fParticle = particle; 739 fParticle = particle; 742 fWaveVector = momentum/hbarc; 740 fWaveVector = momentum/hbarc; 743 fAtomicWeight = A; 741 fAtomicWeight = A; 744 742 745 fNuclearRadius = CalculateNuclearRad(A); 743 fNuclearRadius = CalculateNuclearRad(A); 746 744 747 thetaMax = 10.174/fWaveVector/fNuclearRadius 745 thetaMax = 10.174/fWaveVector/fNuclearRadius; 748 746 749 if (thetaMax > pi) thetaMax = pi; 747 if (thetaMax > pi) thetaMax = pi; 750 748 751 G4Integrator<G4DiffuseElastic,G4double(G4Dif 749 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral; 752 750 753 // result = integral.Legendre10(this,&G4Diff 751 // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta ); 754 norm = integral.Legendre96(this,&G4DiffuseEl 752 norm = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., thetaMax ); 755 753 756 norm *= G4UniformRand(); 754 norm *= G4UniformRand(); 757 755 758 for(i = 1; i <= iMax; i++) 756 for(i = 1; i <= iMax; i++) 759 { 757 { 760 theta1 = (i-1)*thetaMax/iMax; 758 theta1 = (i-1)*thetaMax/iMax; 761 theta2 = i*thetaMax/iMax; 759 theta2 = i*thetaMax/iMax; 762 sum += integral.Legendre10(this,&G4Diffu 760 sum += integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, theta1, theta2); 763 761 764 if ( sum >= norm ) 762 if ( sum >= norm ) 765 { 763 { 766 result = 0.5*(theta1 + theta2); 764 result = 0.5*(theta1 + theta2); 767 break; 765 break; 768 } 766 } 769 } 767 } 770 if (i > iMax ) result = 0.5*(theta1 + theta2 768 if (i > iMax ) result = 0.5*(theta1 + theta2); 771 769 772 G4double sigma = pi*thetaMax/iMax; 770 G4double sigma = pi*thetaMax/iMax; 773 771 774 result += G4RandGauss::shoot(0.,sigma); 772 result += G4RandGauss::shoot(0.,sigma); 775 773 776 if(result < 0.) result = 0.; 774 if(result < 0.) result = 0.; 777 if(result > thetaMax) result = thetaMax; 775 if(result > thetaMax) result = thetaMax; 778 776 779 return result; 777 return result; 780 } 778 } 781 779 782 ////////////////////////////////////////////// 780 ///////////////////////////////////////////////////////////////////////////// 783 ///////////////////// Table preparation and r 781 ///////////////////// Table preparation and reading //////////////////////// 784 ////////////////////////////////////////////// 782 //////////////////////////////////////////////////////////////////////////// 785 // 783 // 786 // Return inv momentum transfer -t > 0 from in 784 // Return inv momentum transfer -t > 0 from initialisation table 787 785 788 G4double G4DiffuseElastic::SampleInvariantT( c 786 G4double G4DiffuseElastic::SampleInvariantT( const G4ParticleDefinition* aParticle, G4double p, 789 787 G4int Z, G4int A) 790 { 788 { 791 fParticle = aParticle; 789 fParticle = aParticle; 792 G4double m1 = fParticle->GetPDGMass(), t; 790 G4double m1 = fParticle->GetPDGMass(), t; 793 G4double totElab = std::sqrt(m1*m1+p*p); 791 G4double totElab = std::sqrt(m1*m1+p*p); 794 G4double mass2 = G4NucleiProperties::GetNucl 792 G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z); 795 G4LorentzVector lv1(p,0.0,0.0,totElab); 793 G4LorentzVector lv1(p,0.0,0.0,totElab); 796 G4LorentzVector lv(0.0,0.0,0.0,mass2); 794 G4LorentzVector lv(0.0,0.0,0.0,mass2); 797 lv += lv1; 795 lv += lv1; 798 796 799 G4ThreeVector bst = lv.boostVector(); 797 G4ThreeVector bst = lv.boostVector(); 800 lv1.boost(-bst); 798 lv1.boost(-bst); 801 799 802 G4ThreeVector p1 = lv1.vect(); 800 G4ThreeVector p1 = lv1.vect(); 803 G4double momentumCMS = p1.mag(); 801 G4double momentumCMS = p1.mag(); 804 802 805 if( aParticle == theNeutron) 803 if( aParticle == theNeutron) 806 { 804 { 807 G4double Tmax = NeutronTuniform( Z ); 805 G4double Tmax = NeutronTuniform( Z ); 808 G4double pCMS2 = momentumCMS*momentumCMS; 806 G4double pCMS2 = momentumCMS*momentumCMS; 809 G4double Tkin = std::sqrt(pCMS2+m1*m1)-m1; 807 G4double Tkin = std::sqrt(pCMS2+m1*m1)-m1; 810 808 811 if( Tkin <= Tmax ) 809 if( Tkin <= Tmax ) 812 { 810 { 813 t = 4.*pCMS2*G4UniformRand(); 811 t = 4.*pCMS2*G4UniformRand(); 814 // G4cout<<Tkin<<", "<<Tmax<<", "<<std:: 812 // G4cout<<Tkin<<", "<<Tmax<<", "<<std::sqrt(t)<<"; "; 815 return t; 813 return t; 816 } 814 } 817 } 815 } 818 816 819 t = SampleTableT( aParticle, momentumCMS, G 817 t = SampleTableT( aParticle, momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms 820 818 821 return t; 819 return t; 822 } 820 } 823 821 824 ////////////////////////////////////////////// 822 /////////////////////////////////////////////////////// 825 823 826 G4double G4DiffuseElastic::NeutronTuniform(G4i 824 G4double G4DiffuseElastic::NeutronTuniform(G4int Z) 827 { 825 { 828 G4double elZ = G4double(Z); 826 G4double elZ = G4double(Z); 829 elZ -= 1.; 827 elZ -= 1.; 830 // G4double Tkin = 20.*G4Exp(-elZ/10.) + 1.; 828 // G4double Tkin = 20.*G4Exp(-elZ/10.) + 1.; 831 G4double Tkin = 12.*G4Exp(-elZ/10.) + 1.; 829 G4double Tkin = 12.*G4Exp(-elZ/10.) + 1.; 832 return Tkin; 830 return Tkin; 833 } 831 } 834 832 835 833 836 ////////////////////////////////////////////// 834 //////////////////////////////////////////////////////////////////////////// 837 // 835 // 838 // Return inv momentum transfer -t > 0 from in 836 // Return inv momentum transfer -t > 0 from initialisation table 839 837 840 G4double G4DiffuseElastic::SampleTableT( const 838 G4double G4DiffuseElastic::SampleTableT( const G4ParticleDefinition* aParticle, G4double p, 841 839 G4double Z, G4double A) 842 { 840 { 843 G4double alpha = SampleTableThetaCMS( aParti 841 G4double alpha = SampleTableThetaCMS( aParticle, p, Z, A); // sample theta2 in cms 844 G4double t = 2*p*p*( 1 - std::cos(std::s 842 G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) ); // -t !!! 845 // G4double t = p*p*alpha; / 843 // G4double t = p*p*alpha; // -t !!! 846 return t; 844 return t; 847 } 845 } 848 846 849 ////////////////////////////////////////////// 847 //////////////////////////////////////////////////////////////////////////// 850 // 848 // 851 // Return scattering angle2 sampled in cms acc 849 // Return scattering angle2 sampled in cms according to precalculated table. 852 850 853 851 854 G4double 852 G4double 855 G4DiffuseElastic::SampleTableThetaCMS(const G4 853 G4DiffuseElastic::SampleTableThetaCMS(const G4ParticleDefinition* particle, 856 G4doubl 854 G4double momentum, G4double Z, G4double A) 857 { 855 { 858 std::size_t iElement; << 856 size_t iElement; 859 G4int iMomentum, iAngle; 857 G4int iMomentum, iAngle; 860 G4double randAngle, position, theta1, theta2 858 G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W; 861 G4double m1 = particle->GetPDGMass(); 859 G4double m1 = particle->GetPDGMass(); 862 860 863 for(iElement = 0; iElement < fElementNumberV 861 for(iElement = 0; iElement < fElementNumberVector.size(); iElement++) 864 { 862 { 865 if( std::fabs(Z - fElementNumberVector[iEl 863 if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break; 866 } 864 } 867 if ( iElement == fElementNumberVector.size() 865 if ( iElement == fElementNumberVector.size() ) 868 { 866 { 869 InitialiseOnFly(Z,A); // table preparation 867 InitialiseOnFly(Z,A); // table preparation, if needed 870 868 871 // iElement--; 869 // iElement--; 872 870 873 // G4cout << "G4DiffuseElastic: Element wi 871 // G4cout << "G4DiffuseElastic: Element with atomic number " << Z 874 // << " is not found, return zero angle" < 872 // << " is not found, return zero angle" << G4endl; 875 // return 0.; // no table for this element 873 // return 0.; // no table for this element 876 } 874 } 877 // G4cout<<"iElement = "<<iElement<<G4endl; 875 // G4cout<<"iElement = "<<iElement<<G4endl; 878 876 879 fAngleTable = fAngleBank[iElement]; 877 fAngleTable = fAngleBank[iElement]; 880 878 881 G4double kinE = std::sqrt(momentum*momentum 879 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1; 882 880 883 for( iMomentum = 0; iMomentum < fEnergyBin; 881 for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++) 884 { 882 { 885 if( kinE < fEnergyVector->GetLowEdgeEnergy 883 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break; 886 } 884 } 887 if ( iMomentum >= fEnergyBin ) iMomentum = f 885 if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1; // kinE is more then theMaxEnergy 888 if ( iMomentum < 0 ) iMomentum = 0 886 if ( iMomentum < 0 ) iMomentum = 0; // against negative index, kinE < theMinEnergy 889 887 890 // G4cout<<"iMomentum = "<<iMomentum<<G4endl 888 // G4cout<<"iMomentum = "<<iMomentum<<G4endl; 891 889 892 if (iMomentum == fEnergyBin -1 || iMomentum 890 if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges 893 { 891 { 894 position = (*(*fAngleTable)(iMomentum))(fA 892 position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand(); 895 893 896 // G4cout<<"position = "<<position<<G4endl 894 // G4cout<<"position = "<<position<<G4endl; 897 895 898 for(iAngle = 0; iAngle < fAngleBin-1; iAng 896 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++) 899 { 897 { 900 if( position > (*(*fAngleTable)(iMomentu 898 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break; 901 } 899 } 902 if (iAngle >= fAngleBin-1) iAngle = fAngle 900 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2; 903 901 904 // G4cout<<"iAngle = "<<iAngle<<G4endl; 902 // G4cout<<"iAngle = "<<iAngle<<G4endl; 905 903 906 randAngle = GetScatteringAngle(iMomentum, 904 randAngle = GetScatteringAngle(iMomentum, iAngle, position); 907 905 908 // G4cout<<"randAngle = "<<randAngle<<G4en 906 // G4cout<<"randAngle = "<<randAngle<<G4endl; 909 } 907 } 910 else // kinE inside between energy table ed 908 else // kinE inside between energy table edges 911 { 909 { 912 // position = (*(*fAngleTable)(iMomentum)) 910 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand(); 913 position = (*(*fAngleTable)(iMomentum))(0) 911 position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand(); 914 912 915 // G4cout<<"position = "<<position<<G4endl 913 // G4cout<<"position = "<<position<<G4endl; 916 914 917 for(iAngle = 0; iAngle < fAngleBin-1; iAng 915 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++) 918 { 916 { 919 // if( position < (*(*fAngleTable)(iMome 917 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break; 920 if( position > (*(*fAngleTable)(iMomentu 918 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break; 921 } 919 } 922 if (iAngle >= fAngleBin-1) iAngle = fAngle 920 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2; 923 921 924 // G4cout<<"iAngle = "<<iAngle<<G4endl; 922 // G4cout<<"iAngle = "<<iAngle<<G4endl; 925 923 926 theta2 = GetScatteringAngle(iMomentum, iA 924 theta2 = GetScatteringAngle(iMomentum, iAngle, position); 927 925 928 // G4cout<<"theta2 = "<<theta2<<G4endl; 926 // G4cout<<"theta2 = "<<theta2<<G4endl; 929 E2 = fEnergyVector->GetLowEdgeEnergy(iMome 927 E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum); 930 928 931 // G4cout<<"E2 = "<<E2<<G4endl; 929 // G4cout<<"E2 = "<<E2<<G4endl; 932 930 933 iMomentum--; 931 iMomentum--; 934 932 935 // position = (*(*fAngleTable)(iMomentum)) 933 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand(); 936 934 937 // G4cout<<"position = "<<position<<G4endl 935 // G4cout<<"position = "<<position<<G4endl; 938 936 939 for(iAngle = 0; iAngle < fAngleBin-1; iAng 937 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++) 940 { 938 { 941 // if( position < (*(*fAngleTable)(iMome 939 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break; 942 if( position > (*(*fAngleTable)(iMomentu 940 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break; 943 } 941 } 944 if (iAngle >= fAngleBin-1) iAngle = fAngle 942 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2; 945 943 946 theta1 = GetScatteringAngle(iMomentum, iA 944 theta1 = GetScatteringAngle(iMomentum, iAngle, position); 947 945 948 // G4cout<<"theta1 = "<<theta1<<G4endl; 946 // G4cout<<"theta1 = "<<theta1<<G4endl; 949 947 950 E1 = fEnergyVector->GetLowEdgeEnergy(iMome 948 E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum); 951 949 952 // G4cout<<"E1 = "<<E1<<G4endl; 950 // G4cout<<"E1 = "<<E1<<G4endl; 953 951 954 W = 1.0/(E2 - E1); 952 W = 1.0/(E2 - E1); 955 W1 = (E2 - kinE)*W; 953 W1 = (E2 - kinE)*W; 956 W2 = (kinE - E1)*W; 954 W2 = (kinE - E1)*W; 957 955 958 randAngle = W1*theta1 + W2*theta2; 956 randAngle = W1*theta1 + W2*theta2; 959 957 960 // randAngle = theta2; 958 // randAngle = theta2; 961 // G4cout<<"randAngle = "<<randAngle<<G4en 959 // G4cout<<"randAngle = "<<randAngle<<G4endl; 962 } 960 } 963 // G4double angle = randAngle; 961 // G4double angle = randAngle; 964 // if (randAngle > 0.) randAngle /= 2*pi*std 962 // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle); 965 963 966 if(randAngle < 0.) randAngle = 0.; 964 if(randAngle < 0.) randAngle = 0.; 967 965 968 return randAngle; 966 return randAngle; 969 } 967 } 970 968 971 ////////////////////////////////////////////// 969 ////////////////////////////////////////////////////////////////////////////// 972 // 970 // 973 // Initialisation for given particle on fly us 971 // Initialisation for given particle on fly using new element number 974 972 975 void G4DiffuseElastic::InitialiseOnFly(G4doubl 973 void G4DiffuseElastic::InitialiseOnFly(G4double Z, G4double A) 976 { 974 { 977 fAtomicNumber = Z; // atomic number 975 fAtomicNumber = Z; // atomic number 978 fAtomicWeight = G4NistManager::Instance()-> 976 fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( Z ) ); 979 977 980 fNuclearRadius = CalculateNuclearRad(fAtomic 978 fNuclearRadius = CalculateNuclearRad(fAtomicWeight); 981 979 982 if( verboseLevel > 0 ) 980 if( verboseLevel > 0 ) 983 { 981 { 984 G4cout<<"G4DiffuseElastic::InitialiseOnFly 982 G4cout<<"G4DiffuseElastic::InitialiseOnFly() the element with Z = " 985 <<Z<<"; and A = "<<A<<G4endl; 983 <<Z<<"; and A = "<<A<<G4endl; 986 } 984 } 987 fElementNumberVector.push_back(fAtomicNumber 985 fElementNumberVector.push_back(fAtomicNumber); 988 986 989 BuildAngleTable(); 987 BuildAngleTable(); 990 988 991 fAngleBank.push_back(fAngleTable); 989 fAngleBank.push_back(fAngleTable); 992 990 993 return; 991 return; 994 } 992 } 995 993 996 ////////////////////////////////////////////// 994 /////////////////////////////////////////////////////////////////////////////// 997 // 995 // 998 // Build for given particle and element table 996 // Build for given particle and element table of momentum, angle probability. 999 // For the moment in lab system. 997 // For the moment in lab system. 1000 998 1001 void G4DiffuseElastic::BuildAngleTable() 999 void G4DiffuseElastic::BuildAngleTable() 1002 { 1000 { 1003 G4int i, j; 1001 G4int i, j; 1004 G4double partMom, kinE, a = 0., z = fPartic 1002 G4double partMom, kinE, a = 0., z = fParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass(); 1005 G4double alpha1, alpha2, alphaMax, alphaCou 1003 G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.; 1006 1004 1007 G4Integrator<G4DiffuseElastic,G4double(G4Di 1005 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral; 1008 1006 1009 fAngleTable = new G4PhysicsTable( fEnergyBi 1007 fAngleTable = new G4PhysicsTable( fEnergyBin ); 1010 1008 1011 for( i = 0; i < fEnergyBin; i++) 1009 for( i = 0; i < fEnergyBin; i++) 1012 { 1010 { 1013 kinE = fEnergyVector->GetLowEdgeEn 1011 kinE = fEnergyVector->GetLowEdgeEnergy(i); 1014 partMom = std::sqrt( kinE*(kinE + 2*m 1012 partMom = std::sqrt( kinE*(kinE + 2*m1) ); 1015 1013 1016 fWaveVector = partMom/hbarc; 1014 fWaveVector = partMom/hbarc; 1017 1015 1018 G4double kR = fWaveVector*fNuclearRad 1016 G4double kR = fWaveVector*fNuclearRadius; 1019 G4double kR2 = kR*kR; 1017 G4double kR2 = kR*kR; 1020 G4double kRmax = 18.6; // 10.6; 10.6, 18 1018 G4double kRmax = 18.6; // 10.6; 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25. 1021 G4double kRcoul = 1.9; // 1.2; 1.4, 2.5; 1019 G4double kRcoul = 1.9; // 1.2; 1.4, 2.5; // on the first slope of J1 1022 // G4double kRlim = 1.2; 1020 // G4double kRlim = 1.2; 1023 // G4double kRlim2 = kRlim*kRlim/kR2; 1021 // G4double kRlim2 = kRlim*kRlim/kR2; 1024 1022 1025 alphaMax = kRmax*kRmax/kR2; 1023 alphaMax = kRmax*kRmax/kR2; 1026 1024 1027 1025 1028 // if (alphaMax > 4.) alphaMax = 4.; // 1026 // if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2 1029 // if ( alphaMax > 4. || alphaMax < 1. ) 1027 // if ( alphaMax > 4. || alphaMax < 1. ) alphaMax = 15.; // vmg27.11.14 1030 1028 1031 // if ( alphaMax > 4. || alphaMax < 1. ) 1029 // if ( alphaMax > 4. || alphaMax < 1. ) alphaMax = CLHEP::pi*CLHEP::pi; // vmg06.01.15 1032 1030 1033 // G4cout<<"alphaMax = "<<alphaMax<<", "; 1031 // G4cout<<"alphaMax = "<<alphaMax<<", "; 1034 1032 1035 if ( alphaMax >= CLHEP::pi*CLHEP::pi ) al 1033 if ( alphaMax >= CLHEP::pi*CLHEP::pi ) alphaMax = CLHEP::pi*CLHEP::pi; // vmg21.10.15 1036 1034 1037 alphaCoulomb = kRcoul*kRcoul/kR2; 1035 alphaCoulomb = kRcoul*kRcoul/kR2; 1038 1036 1039 if( z ) 1037 if( z ) 1040 { 1038 { 1041 a = partMom/m1; // be 1039 a = partMom/m1; // beta*gamma for m1 1042 fBeta = a/std::sqrt(1+a*a); 1040 fBeta = a/std::sqrt(1+a*a); 1043 fZommerfeld = CalculateZommerfeld( fBet 1041 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber); 1044 fAm = CalculateAm( partMom, fZo 1042 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber); 1045 } 1043 } 1046 G4PhysicsFreeVector* angleVector = new G4 1044 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1); 1047 1045 1048 // G4PhysicsLogVector* angleBins = new G 1046 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin ); 1049 1047 1050 G4double delth = alphaMax/fAngleBin; 1048 G4double delth = alphaMax/fAngleBin; 1051 1049 1052 sum = 0.; 1050 sum = 0.; 1053 1051 1054 // fAddCoulomb = false; 1052 // fAddCoulomb = false; 1055 fAddCoulomb = true; 1053 fAddCoulomb = true; 1056 1054 1057 // for(j = 1; j < fAngleBin; j++) 1055 // for(j = 1; j < fAngleBin; j++) 1058 for(j = fAngleBin-1; j >= 1; j--) 1056 for(j = fAngleBin-1; j >= 1; j--) 1059 { 1057 { 1060 // alpha1 = angleBins->GetLowEdgeEnergy 1058 // alpha1 = angleBins->GetLowEdgeEnergy(j-1); 1061 // alpha2 = angleBins->GetLowEdgeEnergy 1059 // alpha2 = angleBins->GetLowEdgeEnergy(j); 1062 1060 1063 // alpha1 = alphaMax*(j-1)/fAngleBin; 1061 // alpha1 = alphaMax*(j-1)/fAngleBin; 1064 // alpha2 = alphaMax*( j )/fAngleBin; 1062 // alpha2 = alphaMax*( j )/fAngleBin; 1065 1063 1066 alpha1 = delth*(j-1); 1064 alpha1 = delth*(j-1); 1067 // if(alpha1 < kRlim2) alpha1 = kRlim2; 1065 // if(alpha1 < kRlim2) alpha1 = kRlim2; 1068 alpha2 = alpha1 + delth; 1066 alpha2 = alpha1 + delth; 1069 1067 1070 // if( ( alpha2 > alphaCoulomb ) && z ) 1068 // if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true; 1071 if( ( alpha1 < alphaCoulomb ) && z ) fA 1069 if( ( alpha1 < alphaCoulomb ) && z ) fAddCoulomb = false; 1072 1070 1073 delta = integral.Legendre10(this, &G4Di 1071 delta = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2); 1074 // delta = integral.Legendre96(this, &G 1072 // delta = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2); 1075 1073 1076 sum += delta; 1074 sum += delta; 1077 1075 1078 angleVector->PutValue( j-1 , alpha1, su 1076 angleVector->PutValue( j-1 , alpha1, sum ); // alpha2 1079 // G4cout<<"j-1 = "<<j-1<<"; alpha << 1077 // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2<<"; sum = "<<sum<<G4endl; 1080 } 1078 } 1081 fAngleTable->insertAt(i, angleVector); 1079 fAngleTable->insertAt(i, angleVector); 1082 1080 1083 // delete[] angleVector; 1081 // delete[] angleVector; 1084 // delete[] angleBins; 1082 // delete[] angleBins; 1085 } 1083 } 1086 return; 1084 return; 1087 } 1085 } 1088 1086 1089 ///////////////////////////////////////////// 1087 ///////////////////////////////////////////////////////////////////////////////// 1090 // 1088 // 1091 // 1089 // 1092 1090 1093 G4double 1091 G4double 1094 G4DiffuseElastic:: GetScatteringAngle( G4int 1092 G4DiffuseElastic:: GetScatteringAngle( G4int iMomentum, G4int iAngle, G4double position ) 1095 { 1093 { 1096 G4double x1, x2, y1, y2, randAngle; 1094 G4double x1, x2, y1, y2, randAngle; 1097 1095 1098 if( iAngle == 0 ) 1096 if( iAngle == 0 ) 1099 { 1097 { 1100 randAngle = (*fAngleTable)(iMomentum)->Ge 1098 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle); 1101 // iAngle++; 1099 // iAngle++; 1102 } 1100 } 1103 else 1101 else 1104 { 1102 { 1105 if ( iAngle >= G4int((*fAngleTable)(iMome 1103 if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) ) 1106 { 1104 { 1107 iAngle = G4int((*fAngleTable)(iMomentum << 1105 iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1; 1108 } 1106 } 1109 y1 = (*(*fAngleTable)(iMomentum))(iAngle- 1107 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1); 1110 y2 = (*(*fAngleTable)(iMomentum))(iAngle) 1108 y2 = (*(*fAngleTable)(iMomentum))(iAngle); 1111 1109 1112 x1 = (*fAngleTable)(iMomentum)->GetLowEdg 1110 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1); 1113 x2 = (*fAngleTable)(iMomentum)->GetLowEdg 1111 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle); 1114 1112 1115 if ( x1 == x2 ) randAngle = x2; 1113 if ( x1 == x2 ) randAngle = x2; 1116 else 1114 else 1117 { 1115 { 1118 if ( y1 == y2 ) randAngle = x1 + ( x2 - 1116 if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand(); 1119 else 1117 else 1120 { 1118 { 1121 randAngle = x1 + ( position - y1 )*( 1119 randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 ); 1122 } 1120 } 1123 } 1121 } 1124 } 1122 } 1125 return randAngle; 1123 return randAngle; 1126 } 1124 } 1127 1125 1128 1126 1129 1127 1130 ///////////////////////////////////////////// 1128 //////////////////////////////////////////////////////////////////////////// 1131 // 1129 // 1132 // Return scattering angle sampled in lab sys 1130 // Return scattering angle sampled in lab system (target at rest) 1133 1131 1134 1132 1135 1133 1136 G4double 1134 G4double 1137 G4DiffuseElastic::SampleThetaLab( const G4Had 1135 G4DiffuseElastic::SampleThetaLab( const G4HadProjectile* aParticle, 1138 G4dou 1136 G4double tmass, G4double A) 1139 { 1137 { 1140 const G4ParticleDefinition* theParticle = a 1138 const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); 1141 G4double m1 = theParticle->GetPDGMass(); 1139 G4double m1 = theParticle->GetPDGMass(); 1142 G4double plab = aParticle->GetTotalMomentum 1140 G4double plab = aParticle->GetTotalMomentum(); 1143 G4LorentzVector lv1 = aParticle->Get4Moment 1141 G4LorentzVector lv1 = aParticle->Get4Momentum(); 1144 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1142 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1145 lv += lv1; 1143 lv += lv1; 1146 1144 1147 G4ThreeVector bst = lv.boostVector(); 1145 G4ThreeVector bst = lv.boostVector(); 1148 lv1.boost(-bst); 1146 lv1.boost(-bst); 1149 1147 1150 G4ThreeVector p1 = lv1.vect(); 1148 G4ThreeVector p1 = lv1.vect(); 1151 G4double ptot = p1.mag(); 1149 G4double ptot = p1.mag(); 1152 G4double tmax = 4.0*ptot*ptot; 1150 G4double tmax = 4.0*ptot*ptot; 1153 G4double t = 0.0; 1151 G4double t = 0.0; 1154 1152 1155 1153 1156 // 1154 // 1157 // Sample t 1155 // Sample t 1158 // 1156 // 1159 1157 1160 t = SampleT( theParticle, ptot, A); 1158 t = SampleT( theParticle, ptot, A); 1161 1159 1162 // NaN finder 1160 // NaN finder 1163 if(!(t < 0.0 || t >= 0.0)) 1161 if(!(t < 0.0 || t >= 0.0)) 1164 { 1162 { 1165 if (verboseLevel > 0) 1163 if (verboseLevel > 0) 1166 { 1164 { 1167 G4cout << "G4DiffuseElastic:WARNING: A 1165 G4cout << "G4DiffuseElastic:WARNING: A = " << A 1168 << " mom(GeV)= " << plab/GeV 1166 << " mom(GeV)= " << plab/GeV 1169 << " S-wave will be sampled" 1167 << " S-wave will be sampled" 1170 << G4endl; 1168 << G4endl; 1171 } 1169 } 1172 t = G4UniformRand()*tmax; 1170 t = G4UniformRand()*tmax; 1173 } 1171 } 1174 if(verboseLevel>1) 1172 if(verboseLevel>1) 1175 { 1173 { 1176 G4cout <<" t= " << t << " tmax= " << tmax 1174 G4cout <<" t= " << t << " tmax= " << tmax 1177 << " ptot= " << ptot << G4endl; 1175 << " ptot= " << ptot << G4endl; 1178 } 1176 } 1179 // Sampling of angles in CM system 1177 // Sampling of angles in CM system 1180 1178 1181 G4double phi = G4UniformRand()*twopi; 1179 G4double phi = G4UniformRand()*twopi; 1182 G4double cost = 1. - 2.0*t/tmax; 1180 G4double cost = 1. - 2.0*t/tmax; 1183 G4double sint; 1181 G4double sint; 1184 1182 1185 if( cost >= 1.0 ) 1183 if( cost >= 1.0 ) 1186 { 1184 { 1187 cost = 1.0; 1185 cost = 1.0; 1188 sint = 0.0; 1186 sint = 0.0; 1189 } 1187 } 1190 else if( cost <= -1.0) 1188 else if( cost <= -1.0) 1191 { 1189 { 1192 cost = -1.0; 1190 cost = -1.0; 1193 sint = 0.0; 1191 sint = 0.0; 1194 } 1192 } 1195 else 1193 else 1196 { 1194 { 1197 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1195 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1198 } 1196 } 1199 if (verboseLevel>1) 1197 if (verboseLevel>1) 1200 { 1198 { 1201 G4cout << "cos(t)=" << cost << " std::sin 1199 G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl; 1202 } 1200 } 1203 G4ThreeVector v1(sint*std::cos(phi),sint*st 1201 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost); 1204 v1 *= ptot; 1202 v1 *= ptot; 1205 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s 1203 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1)); 1206 1204 1207 nlv1.boost(bst); 1205 nlv1.boost(bst); 1208 1206 1209 G4ThreeVector np1 = nlv1.vect(); 1207 G4ThreeVector np1 = nlv1.vect(); 1210 1208 1211 // G4double theta = std::acos( np1.z()/np 1209 // G4double theta = std::acos( np1.z()/np1.mag() ); // degree; 1212 1210 1213 G4double theta = np1.theta(); 1211 G4double theta = np1.theta(); 1214 1212 1215 return theta; 1213 return theta; 1216 } 1214 } 1217 1215 1218 ///////////////////////////////////////////// 1216 //////////////////////////////////////////////////////////////////////////// 1219 // 1217 // 1220 // Return scattering angle in lab system (tar 1218 // Return scattering angle in lab system (target at rest) knowing theta in CMS 1221 1219 1222 1220 1223 1221 1224 G4double 1222 G4double 1225 G4DiffuseElastic::ThetaCMStoThetaLab( const G 1223 G4DiffuseElastic::ThetaCMStoThetaLab( const G4DynamicParticle* aParticle, 1226 G4dou 1224 G4double tmass, G4double thetaCMS) 1227 { 1225 { 1228 const G4ParticleDefinition* theParticle = a 1226 const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); 1229 G4double m1 = theParticle->GetPDGMass(); 1227 G4double m1 = theParticle->GetPDGMass(); 1230 // G4double plab = aParticle->GetTotalMomen 1228 // G4double plab = aParticle->GetTotalMomentum(); 1231 G4LorentzVector lv1 = aParticle->Get4Moment 1229 G4LorentzVector lv1 = aParticle->Get4Momentum(); 1232 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1230 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1233 1231 1234 lv += lv1; 1232 lv += lv1; 1235 1233 1236 G4ThreeVector bst = lv.boostVector(); 1234 G4ThreeVector bst = lv.boostVector(); 1237 1235 1238 lv1.boost(-bst); 1236 lv1.boost(-bst); 1239 1237 1240 G4ThreeVector p1 = lv1.vect(); 1238 G4ThreeVector p1 = lv1.vect(); 1241 G4double ptot = p1.mag(); 1239 G4double ptot = p1.mag(); 1242 1240 1243 G4double phi = G4UniformRand()*twopi; 1241 G4double phi = G4UniformRand()*twopi; 1244 G4double cost = std::cos(thetaCMS); 1242 G4double cost = std::cos(thetaCMS); 1245 G4double sint; 1243 G4double sint; 1246 1244 1247 if( cost >= 1.0 ) 1245 if( cost >= 1.0 ) 1248 { 1246 { 1249 cost = 1.0; 1247 cost = 1.0; 1250 sint = 0.0; 1248 sint = 0.0; 1251 } 1249 } 1252 else if( cost <= -1.0) 1250 else if( cost <= -1.0) 1253 { 1251 { 1254 cost = -1.0; 1252 cost = -1.0; 1255 sint = 0.0; 1253 sint = 0.0; 1256 } 1254 } 1257 else 1255 else 1258 { 1256 { 1259 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1257 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1260 } 1258 } 1261 if (verboseLevel>1) 1259 if (verboseLevel>1) 1262 { 1260 { 1263 G4cout << "cos(tcms)=" << cost << " std:: 1261 G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl; 1264 } 1262 } 1265 G4ThreeVector v1(sint*std::cos(phi),sint*st 1263 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost); 1266 v1 *= ptot; 1264 v1 *= ptot; 1267 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s 1265 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1)); 1268 1266 1269 nlv1.boost(bst); 1267 nlv1.boost(bst); 1270 1268 1271 G4ThreeVector np1 = nlv1.vect(); 1269 G4ThreeVector np1 = nlv1.vect(); 1272 1270 1273 1271 1274 G4double thetaLab = np1.theta(); 1272 G4double thetaLab = np1.theta(); 1275 1273 1276 return thetaLab; 1274 return thetaLab; 1277 } 1275 } 1278 ///////////////////////////////////////////// 1276 //////////////////////////////////////////////////////////////////////////// 1279 // 1277 // 1280 // Return scattering angle in CMS system (tar 1278 // Return scattering angle in CMS system (target at rest) knowing theta in Lab 1281 1279 1282 1280 1283 1281 1284 G4double 1282 G4double 1285 G4DiffuseElastic::ThetaLabToThetaCMS( const G 1283 G4DiffuseElastic::ThetaLabToThetaCMS( const G4DynamicParticle* aParticle, 1286 G4dou 1284 G4double tmass, G4double thetaLab) 1287 { 1285 { 1288 const G4ParticleDefinition* theParticle = a 1286 const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); 1289 G4double m1 = theParticle->GetPDGMass(); 1287 G4double m1 = theParticle->GetPDGMass(); 1290 G4double plab = aParticle->GetTotalMomentum 1288 G4double plab = aParticle->GetTotalMomentum(); 1291 G4LorentzVector lv1 = aParticle->Get4Moment 1289 G4LorentzVector lv1 = aParticle->Get4Momentum(); 1292 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1290 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1293 1291 1294 lv += lv1; 1292 lv += lv1; 1295 1293 1296 G4ThreeVector bst = lv.boostVector(); 1294 G4ThreeVector bst = lv.boostVector(); 1297 1295 1298 // lv1.boost(-bst); 1296 // lv1.boost(-bst); 1299 1297 1300 // G4ThreeVector p1 = lv1.vect(); 1298 // G4ThreeVector p1 = lv1.vect(); 1301 // G4double ptot = p1.mag(); 1299 // G4double ptot = p1.mag(); 1302 1300 1303 G4double phi = G4UniformRand()*twopi; 1301 G4double phi = G4UniformRand()*twopi; 1304 G4double cost = std::cos(thetaLab); 1302 G4double cost = std::cos(thetaLab); 1305 G4double sint; 1303 G4double sint; 1306 1304 1307 if( cost >= 1.0 ) 1305 if( cost >= 1.0 ) 1308 { 1306 { 1309 cost = 1.0; 1307 cost = 1.0; 1310 sint = 0.0; 1308 sint = 0.0; 1311 } 1309 } 1312 else if( cost <= -1.0) 1310 else if( cost <= -1.0) 1313 { 1311 { 1314 cost = -1.0; 1312 cost = -1.0; 1315 sint = 0.0; 1313 sint = 0.0; 1316 } 1314 } 1317 else 1315 else 1318 { 1316 { 1319 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1317 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1320 } 1318 } 1321 if (verboseLevel>1) 1319 if (verboseLevel>1) 1322 { 1320 { 1323 G4cout << "cos(tlab)=" << cost << " std:: 1321 G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl; 1324 } 1322 } 1325 G4ThreeVector v1(sint*std::cos(phi),sint*st 1323 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost); 1326 v1 *= plab; 1324 v1 *= plab; 1327 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s 1325 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1)); 1328 1326 1329 nlv1.boost(-bst); 1327 nlv1.boost(-bst); 1330 1328 1331 G4ThreeVector np1 = nlv1.vect(); 1329 G4ThreeVector np1 = nlv1.vect(); 1332 1330 1333 1331 1334 G4double thetaCMS = np1.theta(); 1332 G4double thetaCMS = np1.theta(); 1335 1333 1336 return thetaCMS; 1334 return thetaCMS; 1337 } 1335 } 1338 1336 1339 ///////////////////////////////////////////// 1337 /////////////////////////////////////////////////////////////////////////////// 1340 // 1338 // 1341 // Test for given particle and element table 1339 // Test for given particle and element table of momentum, angle probability. 1342 // For the moment in lab system. 1340 // For the moment in lab system. 1343 1341 1344 void G4DiffuseElastic::TestAngleTable(const G 1342 void G4DiffuseElastic::TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom, 1345 G 1343 G4double Z, G4double A) 1346 { 1344 { 1347 fAtomicNumber = Z; // atomic number 1345 fAtomicNumber = Z; // atomic number 1348 fAtomicWeight = A; // number of nucleo 1346 fAtomicWeight = A; // number of nucleons 1349 fNuclearRadius = CalculateNuclearRad(fAtomi 1347 fNuclearRadius = CalculateNuclearRad(fAtomicWeight); 1350 1348 1351 1349 1352 1350 1353 G4cout<<"G4DiffuseElastic::TestAngleTable() 1351 G4cout<<"G4DiffuseElastic::TestAngleTable() init the element with Z = " 1354 <<Z<<"; and A = "<<A<<G4endl; 1352 <<Z<<"; and A = "<<A<<G4endl; 1355 1353 1356 fElementNumberVector.push_back(fAtomicNumbe 1354 fElementNumberVector.push_back(fAtomicNumber); 1357 1355 1358 1356 1359 1357 1360 1358 1361 G4int i=0, j; 1359 G4int i=0, j; 1362 G4double a = 0., z = theParticle->GetPDGCha 1360 G4double a = 0., z = theParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass(); 1363 G4double alpha1=0., alpha2=0., alphaMax=0., 1361 G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.; 1364 G4double deltaL10 = 0., deltaL96 = 0., delt 1362 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.; 1365 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0. 1363 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.; 1366 G4double epsilon = 0.001; 1364 G4double epsilon = 0.001; 1367 1365 1368 G4Integrator<G4DiffuseElastic,G4double(G4Di 1366 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral; 1369 1367 1370 fAngleTable = new G4PhysicsTable(fEnergyBin 1368 fAngleTable = new G4PhysicsTable(fEnergyBin); 1371 1369 1372 fWaveVector = partMom/hbarc; 1370 fWaveVector = partMom/hbarc; 1373 1371 1374 G4double kR = fWaveVector*fNuclearRadiu 1372 G4double kR = fWaveVector*fNuclearRadius; 1375 G4double kR2 = kR*kR; 1373 G4double kR2 = kR*kR; 1376 G4double kRmax = 10.6; // 10.6, 18, 10.174 1374 G4double kRmax = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25. 1377 G4double kRcoul = 1.2; // 1.4, 2.5; // on t 1375 G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1 1378 1376 1379 alphaMax = kRmax*kRmax/kR2; 1377 alphaMax = kRmax*kRmax/kR2; 1380 1378 1381 if (alphaMax > 4.) alphaMax = 4.; // vmg05 1379 if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2 1382 1380 1383 alphaCoulomb = kRcoul*kRcoul/kR2; 1381 alphaCoulomb = kRcoul*kRcoul/kR2; 1384 1382 1385 if( z ) 1383 if( z ) 1386 { 1384 { 1387 a = partMom/m1; // beta*gamma 1385 a = partMom/m1; // beta*gamma for m1 1388 fBeta = a/std::sqrt(1+a*a); 1386 fBeta = a/std::sqrt(1+a*a); 1389 fZommerfeld = CalculateZommerfeld( fBet 1387 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber); 1390 fAm = CalculateAm( partMom, fZo 1388 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber); 1391 } 1389 } 1392 G4PhysicsFreeVector* angleVector = new G4Ph 1390 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1); 1393 1391 1394 // G4PhysicsLogVector* angleBins = new G4P 1392 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin ); 1395 1393 1396 1394 1397 fAddCoulomb = false; 1395 fAddCoulomb = false; 1398 1396 1399 for(j = 1; j < fAngleBin; j++) 1397 for(j = 1; j < fAngleBin; j++) 1400 { 1398 { 1401 // alpha1 = angleBins->GetLowEdgeEnergy 1399 // alpha1 = angleBins->GetLowEdgeEnergy(j-1); 1402 // alpha2 = angleBins->GetLowEdgeEnergy 1400 // alpha2 = angleBins->GetLowEdgeEnergy(j); 1403 1401 1404 alpha1 = alphaMax*(j-1)/fAngleBin; 1402 alpha1 = alphaMax*(j-1)/fAngleBin; 1405 alpha2 = alphaMax*( j )/fAngleBin; 1403 alpha2 = alphaMax*( j )/fAngleBin; 1406 1404 1407 if( ( alpha2 > alphaCoulomb ) && z ) fAdd 1405 if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true; 1408 1406 1409 deltaL10 = integral.Legendre10(this, &G4D 1407 deltaL10 = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2); 1410 deltaL96 = integral.Legendre96(this, &G4D 1408 deltaL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2); 1411 deltaAG = integral.AdaptiveGauss(this, & 1409 deltaAG = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction, 1412 alpha1 1410 alpha1, alpha2,epsilon); 1413 1411 1414 // G4cout<<alpha1<<"\t"<<std::sqrt(alph 1412 // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t" 1415 // <<deltaL10<<"\t"<<deltaL96<<"\t" 1413 // <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl; 1416 1414 1417 sumL10 += deltaL10; 1415 sumL10 += deltaL10; 1418 sumL96 += deltaL96; 1416 sumL96 += deltaL96; 1419 sumAG += deltaAG; 1417 sumAG += deltaAG; 1420 1418 1421 G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/d 1419 G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t" 1422 <<sumL10<<"\t"<<sumL96<<"\t"<<sum 1420 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl; 1423 1421 1424 angleVector->PutValue( j-1 , alpha1, sumL 1422 angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2 1425 } 1423 } 1426 fAngleTable->insertAt(i,angleVector); 1424 fAngleTable->insertAt(i,angleVector); 1427 fAngleBank.push_back(fAngleTable); 1425 fAngleBank.push_back(fAngleTable); 1428 1426 1429 /* 1427 /* 1430 // Integral over all angle range - Bad accu 1428 // Integral over all angle range - Bad accuracy !!! 1431 1429 1432 sumL10 = integral.Legendre10(this, &G4Diffu 1430 sumL10 = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2); 1433 sumL96 = integral.Legendre96(this, &G4Diffu 1431 sumL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2); 1434 sumAG = integral.AdaptiveGauss(this, &G4Di 1432 sumAG = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction, 1435 0., al 1433 0., alpha2,epsilon); 1436 G4cout<<G4endl; 1434 G4cout<<G4endl; 1437 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/deg 1435 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t" 1438 <<sumL10<<"\t"<<sumL96<<"\t"<<sum 1436 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl; 1439 */ 1437 */ 1440 return; 1438 return; 1441 } 1439 } 1442 1440 1443 // 1441 // 1444 // 1442 // 1445 ///////////////////////////////////////////// 1443 ///////////////////////////////////////////////////////////////////////////////// 1446 1444