Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 // 28 // Physics model class G4DiffuseElasticV2 29 // 30 // 31 // G4 Model: optical diffuse elastic scattering with 4-momentum balance 32 // 33 // 24-May-07 V. Grichine 34 // 35 // 21.10.15 V. Grichine 36 // Bug fixed in BuildAngleTable, improving accuracy for 37 // angle bins at high energies > 50 GeV for pions. 38 // 39 // 24.11.17 W. Pokorski, code cleanup and performance improvements 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"), fParticle(0) 78 { 79 SetMinEnergy( 0.01*MeV ); 80 SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() ); 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 original 200 to have no wider log-energy-bins up to 10 PeV 93 fAngleBin = 200; 94 95 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin ); 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 element table of application 127 128 void G4DiffuseElasticV2::Initialise() 129 { 130 131 const G4ElementTable* theElementTable = G4Element::GetElementTable(); 132 133 std::size_t jEl, numOfEl = G4Element::GetNumberOfElements(); 134 135 for( jEl = 0; jEl < numOfEl; ++jEl) // application element loop 136 { 137 fAtomicNumber = (*theElementTable)[jEl]->GetZ(); // atomic number 138 fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( fAtomicNumber ) ); 139 fNuclearRadius = CalculateNuclearRad(fAtomicWeight); 140 141 if( verboseLevel > 0 ) 142 { 143 G4cout<<"G4DiffuseElasticV2::Initialise() the element: " 144 <<(*theElementTable)[jEl]->GetName()<<G4endl; 145 } 146 fElementNumberVector.push_back(fAtomicNumber); 147 fElementNameVector.push_back((*theElementTable)[jEl]->GetName()); 148 149 BuildAngleTable(); 150 151 fEnergyAngleVectorBank.push_back(fEnergyAngleVector); 152 fEnergySumVectorBank.push_back(fEnergySumVector); 153 154 } 155 return; 156 } 157 158 //////////////////////////////////////////////////////////////////////////// 159 // 160 // return differential elastic probability d(probability)/d(t) with 161 // Coulomb correction. It is called from BuildAngleTable() 162 163 G4double 164 G4DiffuseElasticV2::GetDiffElasticSumProbA( G4double theta ) 165 { 166 167 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2; 168 G4double delta, diffuse, gamma; 169 G4double e1, e2, bone, bone2; 170 171 // G4double wavek = momentum/hbarc; // wave vector 172 // G4double r0 = 1.08*fermi; 173 // G4double rad = r0*G4Pow::GetInstance()->A13(A); 174 175 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad; 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; // wavek*delta; 213 G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta; 214 215 if( fAddCoulomb ) // add Coulomb correction 216 { 217 G4double sinHalfTheta = std::sin(0.5*theta); 218 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta; 219 220 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0() 221 } 222 223 G4double kgamma2 = kgamma*kgamma; 224 225 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta; 226 // G4double dk2t2 = dk2t*dk2t; 227 228 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta; 229 G4double pikdt = lambda*(1. - G4Exp( -pi*fWaveVector*diffuse*theta/lambda ) ); // wavek*delta; 230 231 damp = DampFactor( pikdt ); 232 damp2 = damp*damp; 233 234 G4double mode2k2 = ( e1*e1 + e2*e2 )*fWaveVector*fWaveVector; 235 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta; 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*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/() 244 sigma += kr2*bonebyarg2; // correction at J1()/() 245 246 sigma *= damp2; // *rad*rad; 247 248 return sigma; 249 } 250 251 252 //////////////////////////////////////////////////////////////////////////// 253 // 254 // return differential elastic probability 2*pi*sin(theta)*d(probability)/d(omega) 255 256 G4double 257 G4DiffuseElasticV2::GetIntegrandFunction( G4double alpha ) 258 { 259 G4double result; 260 261 result = GetDiffElasticSumProbA(alpha) * 2 * CLHEP::pi * std::sin(alpha); 262 263 return result; 264 } 265 266 267 ///////////////////////////////////////////////////////////////////////////// 268 ///////////////////// Table preparation and reading //////////////////////// 269 //////////////////////////////////////////////////////////////////////////// 270 // 271 // Return inv momentum transfer -t > 0 from initialisation table 272 273 G4double G4DiffuseElasticV2::SampleInvariantT( const G4ParticleDefinition* aParticle, G4double p, 274 G4int Z, G4int A) 275 { 276 fParticle = aParticle; 277 G4double m1 = fParticle->GetPDGMass(), t; 278 G4double totElab = std::sqrt(m1*m1+p*p); 279 G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z); 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, G4double(Z), G4double(A) ); // sample theta in cms 304 305 return t; 306 } 307 308 /////////////////////////////////////////////////////// 309 310 G4double G4DiffuseElasticV2::NeutronTuniform(G4int Z) 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 initialisation table 323 324 G4double G4DiffuseElasticV2::SampleTableT( const G4ParticleDefinition* aParticle, G4double p, 325 G4double Z, G4double A) 326 { 327 G4double alpha = SampleTableThetaCMS( aParticle, p, Z, A); // sample theta in cms 328 G4double t = 2*p*p*( 1 - std::cos(alpha) ); // -t !!! 329 330 return t; 331 } 332 333 //////////////////////////////////////////////////////////////////////////// 334 // 335 // Return scattering angle2 sampled in cms according to precalculated table. 336 337 338 G4double 339 G4DiffuseElasticV2::SampleTableThetaCMS(const G4ParticleDefinition* particle, 340 G4double momentum, G4double Z, G4double A) 341 { 342 std::size_t iElement; 343 G4int iMomentum; 344 unsigned long iAngle = 0; 345 G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W; 346 G4double m1 = particle->GetPDGMass(); 347 348 for(iElement = 0; iElement < fElementNumberVector.size(); iElement++) 349 { 350 if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break; 351 } 352 353 if ( iElement == fElementNumberVector.size() ) 354 { 355 InitialiseOnFly(Z,A); // table preparation, if needed 356 } 357 358 fEnergyAngleVector = fEnergyAngleVectorBank[iElement]; 359 fEnergySumVector = fEnergySumVectorBank[iElement]; 360 361 362 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1; 363 364 iMomentum = G4int(fEnergyVector->FindBin(kinE,1000) + 1); 365 366 position = (*(*fEnergySumVector)[iMomentum])[0]*G4UniformRand(); 367 368 for(iAngle = 0; iAngle < fAngleBin; ++iAngle) 369 { 370 if (position > (*(*fEnergySumVector)[iMomentum])[iAngle]) break; 371 } 372 373 374 if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges 375 { 376 randAngle = GetScatteringAngle(iMomentum, iAngle, position); 377 } 378 else // kinE inside between energy table edges 379 { 380 theta2 = GetScatteringAngle(iMomentum, iAngle, position); 381 382 E2 = fEnergyVector->Energy(iMomentum); 383 384 iMomentum--; 385 386 theta1 = GetScatteringAngle(iMomentum, iAngle, position); 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 using new element number 407 408 void G4DiffuseElasticV2::InitialiseOnFly(G4double Z, G4double A) 409 { 410 fAtomicNumber = Z; // atomic number 411 fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( Z ) ); 412 413 fNuclearRadius = CalculateNuclearRad(fAtomicWeight); 414 415 if( verboseLevel > 0 ) 416 { 417 G4cout<<"G4DiffuseElasticV2::InitialiseOnFly() the element with Z = " 418 <<Z<<"; and A = "<<A<<G4endl; 419 } 420 fElementNumberVector.push_back(fAtomicNumber); 421 422 BuildAngleTable(); 423 424 fEnergyAngleVectorBank.push_back(fEnergyAngleVector); 425 fEnergySumVectorBank.push_back(fEnergySumVector); 426 427 return; 428 } 429 430 /////////////////////////////////////////////////////////////////////////////// 431 // 432 // Build for given particle and element table of momentum, angle probability. 433 // For the moment in lab system. 434 435 void G4DiffuseElasticV2::BuildAngleTable() 436 { 437 G4double partMom, kinE, a = 0., z = fParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass(); 438 G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.; 439 440 G4Integrator<G4DiffuseElasticV2,G4double(G4DiffuseElasticV2::*)(G4double)> integral; 441 442 fEnergyAngleVector = new std::vector<std::vector<G4double>*>; 443 fEnergySumVector = new std::vector<std::vector<G4double>*>; 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*fNuclearRadius; 453 G4double kRmax = 18.6; // 10.6; 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25. 454 G4double kRcoul = 1.9; // 1.2; 1.4, 2.5; // on the first slope of J1 455 456 alphaMax = kRmax/kR; 457 458 if ( alphaMax >= CLHEP::pi ) alphaMax = CLHEP::pi; // vmg21.10.15 459 460 alphaCoulomb = kRcoul/kR; 461 462 if( z ) 463 { 464 a = partMom/m1; // beta*gamma for m1 465 fBeta = a/std::sqrt(1+a*a); 466 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber); 467 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber); 468 fAddCoulomb = true; 469 } 470 471 std::vector<G4double>* angleVector = new std::vector<G4double>(fAngleBin); 472 std::vector<G4double>* sumVector = new std::vector<G4double>(fAngleBin); 473 474 475 G4double delth = alphaMax/fAngleBin; 476 477 sum = 0.; 478 479 for(G4int j = (G4int)fAngleBin-1; j >= 0; --j) 480 { 481 alpha1 = delth*j; 482 alpha2 = alpha1 + delth; 483 484 if( fAddCoulomb && ( alpha2 < alphaCoulomb)) fAddCoulomb = false; 485 486 delta = integral.Legendre10(this, &G4DiffuseElasticV2::GetIntegrandFunction, alpha1, alpha2); 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 iMomentum, unsigned long iAngle, G4double position ) 507 { 508 G4double x1, x2, y1, y2, randAngle = 0; 509 510 if( iAngle == 0 ) 511 { 512 randAngle = (*(*fEnergyAngleVector)[iMomentum])[iAngle]; 513 } 514 else 515 { 516 if ( iAngle >= (*fEnergyAngleVector)[iMomentum]->size() ) 517 { 518 iAngle = (*fEnergyAngleVector)[iMomentum]->size() - 1; 519 } 520 521 y1 = (*(*fEnergySumVector)[iMomentum])[iAngle-1]; 522 y2 = (*(*fEnergySumVector)[iMomentum])[iAngle]; 523 524 x1 = (*(*fEnergyAngleVector)[iMomentum])[iAngle-1]; 525 x2 = (*(*fEnergyAngleVector)[iMomentum])[iAngle]; 526 527 if ( x1 == x2 ) randAngle = x2; 528 else 529 { 530 if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand(); 531 else 532 { 533 randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 ); 534 } 535 } 536 } 537 538 return randAngle; 539 } 540 541 542 543 544 //////////////////////////////////////////////////////////////////////////// 545 // 546 // Return scattering angle in lab system (target at rest) knowing theta in CMS 547 548 549 550 G4double 551 G4DiffuseElasticV2::ThetaCMStoThetaLab( const G4DynamicParticle* aParticle, 552 G4double tmass, G4double thetaCMS) 553 { 554 const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); 555 G4double m1 = theParticle->GetPDGMass(); 556 G4LorentzVector lv1 = aParticle->Get4Momentum(); 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::sin(tcms)=" << sint << G4endl; 589 } 590 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost); 591 v1 *= ptot; 592 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1)); 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 (target at rest) knowing theta in Lab 605 606 607 608 G4double 609 G4DiffuseElasticV2::ThetaLabToThetaCMS( const G4DynamicParticle* aParticle, 610 G4double tmass, G4double thetaLab) 611 { 612 const G4ParticleDefinition* theParticle = aParticle->GetDefinition(); 613 G4double m1 = theParticle->GetPDGMass(); 614 G4double plab = aParticle->GetTotalMomentum(); 615 G4LorentzVector lv1 = aParticle->Get4Momentum(); 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::sin(tlab)=" << sint << G4endl; 643 } 644 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost); 645 v1 *= plab; 646 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1)); 647 648 nlv1.boost(-bst); 649 650 G4ThreeVector np1 = nlv1.vect(); 651 G4double thetaCMS = np1.theta(); 652 653 return thetaCMS; 654 } 655 656