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