Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // 28 // Physics model class G4DiffuseElasticV2 29 // 30 // 31 // G4 Model: optical diffuse elastic scatterin 32 // 33 // 24-May-07 V. Grichine 34 // 35 // 21.10.15 V. Grichine 36 // Bug fixed in BuildAngleTable, i 37 // angle bins at high energies > 5 38 // 39 // 24.11.17 W. Pokorski, code cleanup and perf 40 // 41 42 #include "G4DiffuseElasticV2.hh" 43 #include "G4ParticleTable.hh" 44 #include "G4ParticleDefinition.hh" 45 #include "G4IonTable.hh" 46 #include "G4NucleiProperties.hh" 47 48 #include "Randomize.hh" 49 #include "G4Integrator.hh" 50 #include "globals.hh" 51 #include "G4PhysicalConstants.hh" 52 #include "G4SystemOfUnits.hh" 53 54 #include "G4Proton.hh" 55 #include "G4Neutron.hh" 56 #include "G4Deuteron.hh" 57 #include "G4Alpha.hh" 58 #include "G4PionPlus.hh" 59 #include "G4PionMinus.hh" 60 61 #include "G4Element.hh" 62 #include "G4ElementTable.hh" 63 #include "G4NistManager.hh" 64 #include "G4PhysicsTable.hh" 65 #include "G4PhysicsLogVector.hh" 66 #include "G4PhysicsFreeVector.hh" 67 68 #include "G4Exp.hh" 69 70 #include "G4HadronicParameters.hh" 71 72 ////////////////////////////////////////////// 73 // 74 75 76 G4DiffuseElasticV2::G4DiffuseElasticV2() 77 : G4HadronElastic("DiffuseElasticV2"), fPart 78 { 79 SetMinEnergy( 0.01*MeV ); 80 SetMaxEnergy( G4HadronicParameters::Instance 81 82 verboseLevel = 0; 83 lowEnergyRecoilLimit = 100.*keV; 84 lowEnergyLimitQ = 0.0*GeV; 85 lowEnergyLimitHE = 0.0*GeV; 86 lowestEnergyLimit = 0.0*keV; 87 plabLowLimit = 20.0*MeV; 88 89 theProton = G4Proton::Proton(); 90 theNeutron = G4Neutron::Neutron(); 91 92 fEnergyBin = 300; // Increased from the ori 93 fAngleBin = 200; 94 95 fEnergyVector = new G4PhysicsLogVector( the 96 97 fEnergyAngleVector = 0; 98 fEnergySumVector = 0; 99 100 fParticle = 0; 101 fWaveVector = 0.; 102 fAtomicWeight = 0.; 103 fAtomicNumber = 0.; 104 fNuclearRadius = 0.; 105 fBeta = 0.; 106 fZommerfeld = 0.; 107 fAm = 0.; 108 fAddCoulomb = false; 109 } 110 111 ////////////////////////////////////////////// 112 // 113 // Destructor 114 115 G4DiffuseElasticV2::~G4DiffuseElasticV2() 116 { 117 if ( fEnergyVector ) 118 { 119 delete fEnergyVector; 120 fEnergyVector = 0; 121 } 122 } 123 124 ////////////////////////////////////////////// 125 // 126 // Initialisation for given particle using ele 127 128 void G4DiffuseElasticV2::Initialise() 129 { 130 131 const G4ElementTable* theElementTable = G4El 132 133 std::size_t jEl, numOfEl = G4Element::GetNum 134 135 for( jEl = 0; jEl < numOfEl; ++jEl) // appli 136 { 137 fAtomicNumber = (*theElementTable)[jEl]->G 138 fAtomicWeight = G4NistManager::Instance()- 139 fNuclearRadius = CalculateNuclearRad(fAtom 140 141 if( verboseLevel > 0 ) 142 { 143 G4cout<<"G4DiffuseElasticV2::Initialise( 144 <<(*theElementTable)[jEl]->GetName()<<G4 145 } 146 fElementNumberVector.push_back(fAtomicNumb 147 fElementNameVector.push_back((*theElementT 148 149 BuildAngleTable(); 150 151 fEnergyAngleVectorBank.push_back(fEnergyAn 152 fEnergySumVectorBank.push_back(fEnergySumV 153 154 } 155 return; 156 } 157 158 ////////////////////////////////////////////// 159 // 160 // return differential elastic probability d(p 161 // Coulomb correction. It is called from Build 162 163 G4double 164 G4DiffuseElasticV2::GetDiffElasticSumProbA( G4 165 { 166 167 G4double sigma, bzero, bzero2, bonebyarg, bo 168 G4double delta, diffuse, gamma; 169 G4double e1, e2, bone, bone2; 170 171 // G4double wavek = momentum/hbarc; // wave 172 // G4double r0 = 1.08*fermi; 173 // G4double rad = r0*G4Pow::GetInstance()- 174 175 G4double kr = fWaveVector*fNuclearRadius; 176 G4double kr2 = kr*kr; 177 G4double krt = kr*theta; 178 179 bzero = BesselJzero(krt); 180 bzero2 = bzero*bzero; 181 bone = BesselJone(krt); 182 bone2 = bone*bone; 183 bonebyarg = BesselOneByArg(krt); 184 bonebyarg2 = bonebyarg*bonebyarg; 185 186 if ( fParticle == theProton ) 187 { 188 diffuse = 0.63*fermi; 189 gamma = 0.3*fermi; 190 delta = 0.1*fermi*fermi; 191 e1 = 0.3*fermi; 192 e2 = 0.35*fermi; 193 } 194 else if ( fParticle == theNeutron ) 195 { 196 diffuse = 0.63*fermi; 197 gamma = 0.3*fermi; 198 delta = 0.1*fermi*fermi; 199 e1 = 0.3*fermi; 200 e2 = 0.35*fermi; 201 } 202 else // as proton, if were not defined 203 { 204 diffuse = 0.63*fermi; 205 gamma = 0.3*fermi; 206 delta = 0.1*fermi*fermi; 207 e1 = 0.3*fermi; 208 e2 = 0.35*fermi; 209 } 210 211 G4double lambda = 15; // 15 ok 212 // G4double kgamma = fWaveVector*gamma; 213 G4double kgamma = lambda*(1.-G4Exp(-fWave 214 215 if( fAddCoulomb ) // add Coulomb correction 216 { 217 G4double sinHalfTheta = std::sin(0.5*thet 218 G4double sinHalfTheta2 = sinHalfTheta*sinH 219 220 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta 221 } 222 223 G4double kgamma2 = kgamma*kgamma; 224 225 // G4double dk2t = delta*fWaveVector*fWaveV 226 // G4double dk2t2 = dk2t*dk2t; 227 228 // G4double pikdt = pi*fWaveVector*diffuse*t 229 G4double pikdt = lambda*(1. - G4Exp( -pi* 230 231 damp = DampFactor( pikdt ); 232 damp2 = damp*damp; 233 234 G4double mode2k2 = ( e1*e1 + e2*e2 )*fWaveVe 235 G4double e2dk3t = -2.*e2*delta*fWaveVector* 236 237 sigma = kgamma2; 238 // sigma += dk2t2; 239 sigma *= bzero2; 240 sigma += mode2k2*bone2; 241 sigma += e2dk3t*bzero*bone; 242 243 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerf 244 sigma += kr2*bonebyarg2; // correction at J 245 246 sigma *= damp2; // *rad*rad; 247 248 return sigma; 249 } 250 251 252 ////////////////////////////////////////////// 253 // 254 // return differential elastic probability 2*p 255 256 G4double 257 G4DiffuseElasticV2::GetIntegrandFunction( G4do 258 { 259 G4double result; 260 261 result = GetDiffElasticSumProbA(alpha) * 2 262 263 return result; 264 } 265 266 267 ////////////////////////////////////////////// 268 ///////////////////// Table preparation and r 269 ////////////////////////////////////////////// 270 // 271 // Return inv momentum transfer -t > 0 from in 272 273 G4double G4DiffuseElasticV2::SampleInvariantT( 274 275 { 276 fParticle = aParticle; 277 G4double m1 = fParticle->GetPDGMass(), t; 278 G4double totElab = std::sqrt(m1*m1+p*p); 279 G4double mass2 = G4NucleiProperties::GetNucl 280 G4LorentzVector lv1(p,0.0,0.0,totElab); 281 G4LorentzVector lv(0.0,0.0,0.0,mass2); 282 lv += lv1; 283 284 G4ThreeVector bst = lv.boostVector(); 285 lv1.boost(-bst); 286 287 G4ThreeVector p1 = lv1.vect(); 288 G4double momentumCMS = p1.mag(); 289 290 if( aParticle == theNeutron) 291 { 292 G4double Tmax = NeutronTuniform( Z ); 293 G4double pCMS2 = momentumCMS*momentumCMS; 294 G4double Tkin = std::sqrt(pCMS2+m1*m1)-m1; 295 296 if( Tkin <= Tmax ) 297 { 298 t = 4.*pCMS2*G4UniformRand(); 299 return t; 300 } 301 } 302 303 t = SampleTableT( aParticle, momentumCMS, G 304 305 return t; 306 } 307 308 ////////////////////////////////////////////// 309 310 G4double G4DiffuseElasticV2::NeutronTuniform(G 311 { 312 G4double elZ = G4double(Z); 313 elZ -= 1.; 314 G4double Tkin = 12.*G4Exp(-elZ/10.) + 1.; 315 316 return Tkin; 317 } 318 319 320 ////////////////////////////////////////////// 321 // 322 // Return inv momentum transfer -t > 0 from in 323 324 G4double G4DiffuseElasticV2::SampleTableT( con 325 326 { 327 G4double alpha = SampleTableThetaCMS( aParti 328 G4double t = 2*p*p*( 1 - std::cos(alpha) 329 330 return t; 331 } 332 333 ////////////////////////////////////////////// 334 // 335 // Return scattering angle2 sampled in cms acc 336 337 338 G4double 339 G4DiffuseElasticV2::SampleTableThetaCMS(const 340 G4doubl 341 { 342 std::size_t iElement; 343 G4int iMomentum; 344 unsigned long iAngle = 0; 345 G4double randAngle, position, theta1, theta2 346 G4double m1 = particle->GetPDGMass(); 347 348 for(iElement = 0; iElement < fElementNumberV 349 { 350 if( std::fabs(Z - fElementNumberVector[iEl 351 } 352 353 if ( iElement == fElementNumberVector.size() 354 { 355 InitialiseOnFly(Z,A); // table preparation 356 } 357 358 fEnergyAngleVector = fEnergyAngleVectorBank[ 359 fEnergySumVector = fEnergySumVectorBank[iEle 360 361 362 G4double kinE = std::sqrt(momentum*momentum 363 364 iMomentum = G4int(fEnergyVector->FindBin(kin 365 366 position = (*(*fEnergySumVector)[iMomentum]) 367 368 for(iAngle = 0; iAngle < fAngleBin; ++iAngle 369 { 370 if (position > (*(*fEnergySumVector)[iMo 371 } 372 373 374 if (iMomentum == fEnergyBin -1 || iMomentum 375 { 376 randAngle = GetScatteringAngle(iMomentum 377 } 378 else // kinE inside between energy table ed 379 { 380 theta2 = GetScatteringAngle(iMomentum, 381 382 E2 = fEnergyVector->Energy(iMomentum); 383 384 iMomentum--; 385 386 theta1 = GetScatteringAngle(iMomentum, 387 388 E1 = fEnergyVector->Energy(iMomentum); 389 390 W = 1.0/(E2 - E1); 391 W1 = (E2 - kinE)*W; 392 W2 = (kinE - E1)*W; 393 394 randAngle = W1*theta1 + W2*theta2; 395 } 396 397 398 399 if(randAngle < 0.) randAngle = 0.; 400 401 return randAngle; 402 } 403 404 ////////////////////////////////////////////// 405 // 406 // Initialisation for given particle on fly us 407 408 void G4DiffuseElasticV2::InitialiseOnFly(G4dou 409 { 410 fAtomicNumber = Z; // atomic number 411 fAtomicWeight = G4NistManager::Instance()-> 412 413 fNuclearRadius = CalculateNuclearRad(fAtomic 414 415 if( verboseLevel > 0 ) 416 { 417 G4cout<<"G4DiffuseElasticV2::InitialiseOnF 418 <<Z<<"; and A = "<<A<<G4endl; 419 } 420 fElementNumberVector.push_back(fAtomicNumber 421 422 BuildAngleTable(); 423 424 fEnergyAngleVectorBank.push_back(fEnergyAngl 425 fEnergySumVectorBank.push_back(fEnergySumVec 426 427 return; 428 } 429 430 ////////////////////////////////////////////// 431 // 432 // Build for given particle and element table 433 // For the moment in lab system. 434 435 void G4DiffuseElasticV2::BuildAngleTable() 436 { 437 G4double partMom, kinE, a = 0., z = fParticl 438 G4double alpha1, alpha2, alphaMax, alphaCoul 439 440 G4Integrator<G4DiffuseElasticV2,G4double(G4D 441 442 fEnergyAngleVector = new std::vector<std::ve 443 fEnergySumVector = new std::vector<std::vect 444 445 for( G4int i = 0; i < fEnergyBin; ++i) 446 { 447 kinE = fEnergyVector->Energy(i); 448 partMom = std::sqrt( kinE*(kinE + 2*m1 449 450 fWaveVector = partMom/hbarc; 451 452 G4double kR = fWaveVector*fNuclearRadi 453 G4double kRmax = 18.6; // 10.6; 10.6, 18, 454 G4double kRcoul = 1.9; // 1.2; 1.4, 2.5; / 455 456 alphaMax = kRmax/kR; 457 458 if ( alphaMax >= CLHEP::pi ) alphaMax = CL 459 460 alphaCoulomb = kRcoul/kR; 461 462 if( z ) 463 { 464 a = partMom/m1; // bet 465 fBeta = a/std::sqrt(1+a*a); 466 fZommerfeld = CalculateZommerfeld( fBeta 467 fAm = CalculateAm( partMom, fZom 468 fAddCoulomb = true; 469 } 470 471 std::vector<G4double>* angleVector = new s 472 std::vector<G4double>* sumVector = new std 473 474 475 G4double delth = alphaMax/fAngleBin; 476 477 sum = 0.; 478 479 for(G4int j = (G4int)fAngleBin-1; j >= 0; 480 { 481 alpha1 = delth*j; 482 alpha2 = alpha1 + delth; 483 484 if( fAddCoulomb && ( alpha2 < alphaCoulo 485 486 delta = integral.Legendre10(this, &G4Dif 487 488 sum += delta; 489 490 (*angleVector)[j] = alpha1; 491 (*sumVector)[j] = sum; 492 493 } 494 fEnergyAngleVector->push_back(angleVector) 495 fEnergySumVector->push_back(sumVector); 496 497 } 498 return; 499 } 500 501 ////////////////////////////////////////////// 502 // 503 // 504 505 G4double 506 G4DiffuseElasticV2::GetScatteringAngle( G4int 507 { 508 G4double x1, x2, y1, y2, randAngle = 0; 509 510 if( iAngle == 0 ) 511 { 512 randAngle = (*(*fEnergyAngleVector)[iMomen 513 } 514 else 515 { 516 if ( iAngle >= (*fEnergyAngleVector)[iMome 517 { 518 iAngle = (*fEnergyAngleVector)[iMomentum 519 } 520 521 y1 = (*(*fEnergySumVector)[iMomentum])[iAn 522 y2 = (*(*fEnergySumVector)[iMomentum])[iAn 523 524 x1 = (*(*fEnergyAngleVector)[iMomentum])[i 525 x2 = (*(*fEnergyAngleVector)[iMomentum])[i 526 527 if ( x1 == x2 ) randAngle = x2; 528 else 529 { 530 if ( y1 == y2 ) randAngle = x1 + ( x2 - 531 else 532 { 533 randAngle = x1 + ( position - y1 )*( x 534 } 535 } 536 } 537 538 return randAngle; 539 } 540 541 542 543 544 ////////////////////////////////////////////// 545 // 546 // Return scattering angle in lab system (targ 547 548 549 550 G4double 551 G4DiffuseElasticV2::ThetaCMStoThetaLab( const 552 G4doub 553 { 554 const G4ParticleDefinition* theParticle = aP 555 G4double m1 = theParticle->GetPDGMass(); 556 G4LorentzVector lv1 = aParticle->Get4Momentu 557 G4LorentzVector lv(0.0,0.0,0.0,tmass); 558 559 lv += lv1; 560 561 G4ThreeVector bst = lv.boostVector(); 562 563 lv1.boost(-bst); 564 565 G4ThreeVector p1 = lv1.vect(); 566 G4double ptot = p1.mag(); 567 568 G4double phi = G4UniformRand()*twopi; 569 G4double cost = std::cos(thetaCMS); 570 G4double sint; 571 572 if( cost >= 1.0 ) 573 { 574 cost = 1.0; 575 sint = 0.0; 576 } 577 else if( cost <= -1.0) 578 { 579 cost = -1.0; 580 sint = 0.0; 581 } 582 else 583 { 584 sint = std::sqrt((1.0-cost)*(1.0+cost)); 585 } 586 if (verboseLevel>1) 587 { 588 G4cout << "cos(tcms)=" << cost << " std::s 589 } 590 G4ThreeVector v1(sint*std::cos(phi),sint*std 591 v1 *= ptot; 592 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),st 593 594 nlv1.boost(bst); 595 596 G4ThreeVector np1 = nlv1.vect(); 597 598 G4double thetaLab = np1.theta(); 599 600 return thetaLab; 601 } 602 ////////////////////////////////////////////// 603 // 604 // Return scattering angle in CMS system (targ 605 606 607 608 G4double 609 G4DiffuseElasticV2::ThetaLabToThetaCMS( const 610 G4doub 611 { 612 const G4ParticleDefinition* theParticle = aP 613 G4double m1 = theParticle->GetPDGMass(); 614 G4double plab = aParticle->GetTotalMomentum( 615 G4LorentzVector lv1 = aParticle->Get4Momentu 616 G4LorentzVector lv(0.0,0.0,0.0,tmass); 617 618 lv += lv1; 619 620 G4ThreeVector bst = lv.boostVector(); 621 622 G4double phi = G4UniformRand()*twopi; 623 G4double cost = std::cos(thetaLab); 624 G4double sint; 625 626 if( cost >= 1.0 ) 627 { 628 cost = 1.0; 629 sint = 0.0; 630 } 631 else if( cost <= -1.0) 632 { 633 cost = -1.0; 634 sint = 0.0; 635 } 636 else 637 { 638 sint = std::sqrt((1.0-cost)*(1.0+cost)); 639 } 640 if (verboseLevel>1) 641 { 642 G4cout << "cos(tlab)=" << cost << " std::s 643 } 644 G4ThreeVector v1(sint*std::cos(phi),sint*std 645 v1 *= plab; 646 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),st 647 648 nlv1.boost(-bst); 649 650 G4ThreeVector np1 = nlv1.vect(); 651 G4double thetaCMS = np1.theta(); 652 653 return thetaCMS; 654 } 655 656