Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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, result, theta1, theta2, thetaMax, sum = 0.;
728
729 fParticle = particle;
730 fWaveVector = momentum/hbarc;
731 fAtomicWeight = A;
732
733 fNuclearRadius = CalculateNuclearRad(A);
734
735 thetaMax = 10.174/fWaveVector/fNuclearRadius;
736
737 if (thetaMax > pi) thetaMax = pi;
738
739 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
740
741 // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
742 norm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., thetaMax );
743
744 norm *= G4UniformRand();
745
746 for(i = 1; i <= iMax; i++)
747 {
748 theta1 = (i-1)*thetaMax/iMax;
749 theta2 = i*thetaMax/iMax;
750 sum += integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, theta1, theta2);
751
752 if ( sum >= norm )
753 {
754 result = 0.5*(theta1 + theta2);
755 break;
756 }
757 }
758 if (i > iMax ) result = 0.5*(theta1 + theta2);
759
760 G4double sigma = pi*thetaMax/iMax;
761
762 result += G4RandGauss::shoot(0.,sigma);
763
764 if(result < 0.) result = 0.;
765 if(result > thetaMax) result = thetaMax;
766
767 return result;
768}
769
770/////////////////////////////////////////////////////////////////////////////
771///////////////////// Table preparation and reading ////////////////////////
772
773////////////////////////////////////////////////////////////////////////////
774//
775// Interface function. Return inv momentum transfer -t > 0 from initialisation table
776
777G4double G4NuclNuclDiffuseElastic::SampleInvariantT( const G4ParticleDefinition* aParticle, G4double p,
778 G4int Z, G4int A)
779{
780 fParticle = aParticle;
781 fAtomicNumber = G4double(Z);
782 fAtomicWeight = G4double(A);
783
784 G4double t(0.), m1 = fParticle->GetPDGMass();
785 G4double totElab = std::sqrt(m1*m1+p*p);
786 G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z);
787 G4LorentzVector lv1(p,0.0,0.0,totElab);
788 G4LorentzVector lv(0.0,0.0,0.0,mass2);
789 lv += lv1;
790
791 G4ThreeVector bst = lv.boostVector();
792 lv1.boost(-bst);
793
794 G4ThreeVector p1 = lv1.vect();
795 G4double momentumCMS = p1.mag();
796
797 // t = SampleTableT( aParticle, momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms
798
799 t = SampleCoulombMuCMS( aParticle, momentumCMS); // sample theta2 in cms
800
801
802 return t;
803}
804
805////////////////////////////////////////////////////////////////////////////
806//
807// Return inv momentum transfer -t > 0 as Coulomb scattering <= thetaC
808
809G4double G4NuclNuclDiffuseElastic::SampleCoulombMuCMS( const G4ParticleDefinition* aParticle, G4double p)
810{
811 G4double t(0.), rand(0.), mu(0.);
812
813 G4double A1 = G4double( aParticle->GetBaryonNumber() );
814 G4double R1 = CalculateNuclearRad(A1);
815
816 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
817 fNuclearRadius += R1;
818
819 InitDynParameters(fParticle, p);
820
821 fCoulombMuC = fHalfRutThetaTg2/(1.+fHalfRutThetaTg2);
822
823 rand = G4UniformRand();
824
825 // sample (1-cosTheta) in 0 < Theta < ThetaC as Coulomb scattering
826
827 mu = fCoulombMuC*rand*fAm;
828 mu /= fAm + fCoulombMuC*( 1. - rand );
829
830 t = 4.*p*p*mu;
831
832 return t;
833}
834
835
836////////////////////////////////////////////////////////////////////////////
837//
838// Return inv momentum transfer -t > 0 from initialisation table
839
840G4double G4NuclNuclDiffuseElastic::SampleTableT( const G4ParticleDefinition* aParticle, G4double p,
841 G4double Z, G4double A)
842{
843 G4double alpha = SampleTableThetaCMS( aParticle, p, Z, A); // sample theta2 in cms
844 // G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) ); // -t !!!
845 G4double t = p*p*alpha; // -t !!!
846 return t;
847}
848
849////////////////////////////////////////////////////////////////////////////
850//
851// Return scattering angle2 sampled in cms according to precalculated table.
852
853
854G4double
855G4NuclNuclDiffuseElastic::SampleTableThetaCMS(const G4ParticleDefinition* particle,
856 G4double momentum, G4double Z, G4double A)
857{
858 std::size_t iElement;
859 G4int iMomentum, iAngle;
860 G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W;
861 G4double m1 = particle->GetPDGMass();
862
863 for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
864 {
865 if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
866 }
867 if ( iElement == fElementNumberVector.size() )
868 {
869 InitialiseOnFly(Z,A); // table preparation, if needed
870
871 // iElement--;
872
873 // G4cout << "G4NuclNuclDiffuseElastic: Element with atomic number " << Z
874 // << " is not found, return zero angle" << G4endl;
875 // return 0.; // no table for this element
876 }
877 // G4cout<<"iElement = "<<iElement<<G4endl;
878
879 fAngleTable = fAngleBank[iElement];
880
881 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
882
883 for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
884 {
885 // G4cout<<iMomentum<<", kinE = "<<kinE/MeV<<", vectorE = "<<fEnergyVector->GetLowEdgeEnergy(iMomentum)/MeV<<G4endl;
886
887 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
888 }
889 // G4cout<<"iMomentum = "<<iMomentum<<", fEnergyBin -1 = "<<fEnergyBin -1<<G4endl;
890
891
892 if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1; // kinE is more then theMaxEnergy
893 if ( iMomentum < 0 ) iMomentum = 0; // against negative index, kinE < theMinEnergy
894
895
896 if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges
897 {
898 position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
899
900 // G4cout<<"position = "<<position<<G4endl;
901
902 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
903 {
904 if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
905 }
906 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
907
908 // G4cout<<"iAngle = "<<iAngle<<G4endl;
909
910 randAngle = GetScatteringAngle(iMomentum, iAngle, position);
911
912 // G4cout<<"randAngle = "<<randAngle<<G4endl;
913 }
914 else // kinE inside between energy table edges
915 {
916 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
917 position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand();
918
919 // G4cout<<"position = "<<position<<G4endl;
920
921 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
922 {
923 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
924 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
925 }
926 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
927
928 // G4cout<<"iAngle = "<<iAngle<<G4endl;
929
930 theta2 = GetScatteringAngle(iMomentum, iAngle, position);
931
932 // G4cout<<"theta2 = "<<theta2<<G4endl;
933
934 E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
935
936 // G4cout<<"E2 = "<<E2<<G4endl;
937
938 iMomentum--;
939
940 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
941
942 // G4cout<<"position = "<<position<<G4endl;
943
944 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
945 {
946 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
947 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
948 }
949 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
950
951 theta1 = GetScatteringAngle(iMomentum, iAngle, position);
952
953 // G4cout<<"theta1 = "<<theta1<<G4endl;
954
955 E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
956
957 // G4cout<<"E1 = "<<E1<<G4endl;
958
959 W = 1.0/(E2 - E1);
960 W1 = (E2 - kinE)*W;
961 W2 = (kinE - E1)*W;
962
963 randAngle = W1*theta1 + W2*theta2;
964
965 // randAngle = theta2;
966 }
967 // G4double angle = randAngle;
968 // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
969 // G4cout<<"randAngle = "<<randAngle/degree<<G4endl;
970
971 return randAngle;
972}
973
974//////////////////////////////////////////////////////////////////////////////
975//
976// Initialisation for given particle on fly using new element number
977
978void G4NuclNuclDiffuseElastic::InitialiseOnFly(G4double Z, G4double A)
979{
980 fAtomicNumber = Z; // atomic number
981 fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( Z ) );
982
983 G4double A1 = G4double( fParticle->GetBaryonNumber() );
984 G4double R1 = CalculateNuclearRad(A1);
985
986 fNuclearRadius = CalculateNuclearRad(fAtomicWeight) + R1;
987
988 if( verboseLevel > 0 )
989 {
990 G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element with Z = "
991 <<Z<<"; and A = "<<A<<G4endl;
992 }
993 fElementNumberVector.push_back(fAtomicNumber);
994
995 BuildAngleTable();
996
997 fAngleBank.push_back(fAngleTable);
998
999 return;
1000}
1001
1002///////////////////////////////////////////////////////////////////////////////
1003//
1004// Build for given particle and element table of momentum, angle probability.
1005// For the moment in lab system.
1006
1007void G4NuclNuclDiffuseElastic::BuildAngleTable()
1008{
1009 G4int i, j;
1010 G4double partMom, kinE, m1 = fParticle->GetPDGMass();
1011 G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
1012
1013 // G4cout<<"particle z = "<<z<<"; particle m1 = "<<m1/GeV<<" GeV"<<G4endl;
1014
1015 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1016
1017 fAngleTable = new G4PhysicsTable(fEnergyBin);
1018
1019 for( i = 0; i < fEnergyBin; i++)
1020 {
1021 kinE = fEnergyVector->GetLowEdgeEnergy(i);
1022
1023 // G4cout<<G4endl;
1024 // G4cout<<"kinE = "<<kinE/MeV<<" MeV"<<G4endl;
1025
1026 partMom = std::sqrt( kinE*(kinE + 2*m1) );
1027
1028 InitDynParameters(fParticle, partMom);
1029
1030 alphaMax = fRutherfordTheta*fCofAlphaMax;
1031
1032 if(alphaMax > pi) alphaMax = pi;
1033
1034 // VI: Coverity complain
1035 //alphaMax = pi2;
1036
1037 alphaCoulomb = fRutherfordTheta*fCofAlphaCoulomb;
1038
1039 // G4cout<<"alphaCoulomb = "<<alphaCoulomb/degree<<"; alphaMax = "<<alphaMax/degree<<G4endl;
1040
1041 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1042
1043 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1044
1045 G4double delth = (alphaMax-alphaCoulomb)/fAngleBin;
1046
1047 sum = 0.;
1048
1049 // fAddCoulomb = false;
1050 fAddCoulomb = true;
1051
1052 // for(j = 1; j < fAngleBin; j++)
1053 for(j = fAngleBin-1; j >= 1; j--)
1054 {
1055 // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1056 // alpha2 = angleBins->GetLowEdgeEnergy(j);
1057
1058 // alpha1 = alphaMax*(j-1)/fAngleBin;
1059 // alpha2 = alphaMax*( j )/fAngleBin;
1060
1061 alpha1 = alphaCoulomb + delth*(j-1);
1062 // if(alpha1 < kRlim2) alpha1 = kRlim2;
1063 alpha2 = alpha1 + delth;
1064
1065 delta = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc, alpha1, alpha2);
1066 // delta = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1067
1068 sum += delta;
1069
1070 angleVector->PutValue( j-1 , alpha1, sum ); // alpha2
1071 // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2/degree<<"; sum = "<<sum<<G4endl;
1072 }
1073 fAngleTable->insertAt(i,angleVector);
1074
1075 // delete[] angleVector;
1076 // delete[] angleBins;
1077 }
1078 return;
1079}
1080
1081/////////////////////////////////////////////////////////////////////////////////
1082//
1083//
1084
1085G4double
1086G4NuclNuclDiffuseElastic:: GetScatteringAngle( G4int iMomentum, G4int iAngle, G4double position )
1087{
1088 G4double x1, x2, y1, y2, randAngle;
1089
1090 if( iAngle == 0 )
1091 {
1092 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1093 // iAngle++;
1094 }
1095 else
1096 {
1097 if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
1098 {
1099 iAngle = G4int((*fAngleTable)(iMomentum)->GetVectorLength() - 1);
1100 }
1101 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1102 y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1103
1104 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1105 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1106
1107 if ( x1 == x2 ) randAngle = x2;
1108 else
1109 {
1110 if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand();
1111 else
1112 {
1113 randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
1114 }
1115 }
1116 }
1117 return randAngle;
1118}
1119
1120
1121
1122////////////////////////////////////////////////////////////////////////////
1123//
1124// Return scattering angle sampled in lab system (target at rest)
1125
1126
1127
1128G4double
1129G4NuclNuclDiffuseElastic::SampleThetaLab( const G4HadProjectile* aParticle,
1130 G4double tmass, G4double A)
1131{
1132 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1133 G4double m1 = theParticle->GetPDGMass();
1134 G4double plab = aParticle->GetTotalMomentum();
1135 G4LorentzVector lv1 = aParticle->Get4Momentum();
1136 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1137 lv += lv1;
1138
1139 G4ThreeVector bst = lv.boostVector();
1140 lv1.boost(-bst);
1141
1142 G4ThreeVector p1 = lv1.vect();
1143 G4double ptot = p1.mag();
1144 G4double tmax = 4.0*ptot*ptot;
1145 G4double t = 0.0;
1146
1147
1148 //
1149 // Sample t
1150 //
1151
1152 t = SampleT( theParticle, ptot, A);
1153
1154 // NaN finder
1155 if(!(t < 0.0 || t >= 0.0))
1156 {
1157 if (verboseLevel > 0)
1158 {
1159 G4cout << "G4NuclNuclDiffuseElastic:WARNING: A = " << A
1160 << " mom(GeV)= " << plab/GeV
1161 << " S-wave will be sampled"
1162 << G4endl;
1163 }
1164 t = G4UniformRand()*tmax;
1165 }
1166 if(verboseLevel>1)
1167 {
1168 G4cout <<" t= " << t << " tmax= " << tmax
1169 << " ptot= " << ptot << G4endl;
1170 }
1171 // Sampling of angles in CM system
1172
1173 G4double phi = G4UniformRand()*twopi;
1174 G4double cost = 1. - 2.0*t/tmax;
1175 G4double sint;
1176
1177 if( cost >= 1.0 )
1178 {
1179 cost = 1.0;
1180 sint = 0.0;
1181 }
1182 else if( cost <= -1.0)
1183 {
1184 cost = -1.0;
1185 sint = 0.0;
1186 }
1187 else
1188 {
1189 sint = std::sqrt((1.0-cost)*(1.0+cost));
1190 }
1191 if (verboseLevel>1)
1192 {
1193 G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
1194 }
1195 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1196 v1 *= ptot;
1197 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1198
1199 nlv1.boost(bst);
1200
1201 G4ThreeVector np1 = nlv1.vect();
1202
1203 // G4double theta = std::acos( np1.z()/np1.mag() ); // degree;
1204
1205 G4double theta = np1.theta();
1206
1207 return theta;
1208}
1209
1210////////////////////////////////////////////////////////////////////////////
1211//
1212// Return scattering angle in lab system (target at rest) knowing theta in CMS
1213
1214
1215
1216G4double
1217G4NuclNuclDiffuseElastic::ThetaCMStoThetaLab( const G4DynamicParticle* aParticle,
1218 G4double tmass, G4double thetaCMS)
1219{
1220 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1221 G4double m1 = theParticle->GetPDGMass();
1222 // G4double plab = aParticle->GetTotalMomentum();
1223 G4LorentzVector lv1 = aParticle->Get4Momentum();
1224 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1225
1226 lv += lv1;
1227
1228 G4ThreeVector bst = lv.boostVector();
1229
1230 lv1.boost(-bst);
1231
1232 G4ThreeVector p1 = lv1.vect();
1233 G4double ptot = p1.mag();
1234
1235 G4double phi = G4UniformRand()*twopi;
1236 G4double cost = std::cos(thetaCMS);
1237 G4double sint;
1238
1239 if( cost >= 1.0 )
1240 {
1241 cost = 1.0;
1242 sint = 0.0;
1243 }
1244 else if( cost <= -1.0)
1245 {
1246 cost = -1.0;
1247 sint = 0.0;
1248 }
1249 else
1250 {
1251 sint = std::sqrt((1.0-cost)*(1.0+cost));
1252 }
1253 if (verboseLevel>1)
1254 {
1255 G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
1256 }
1257 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1258 v1 *= ptot;
1259 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1260
1261 nlv1.boost(bst);
1262
1263 G4ThreeVector np1 = nlv1.vect();
1264
1265
1266 G4double thetaLab = np1.theta();
1267
1268 return thetaLab;
1269}
1270
1271////////////////////////////////////////////////////////////////////////////
1272//
1273// Return scattering angle in CMS system (target at rest) knowing theta in Lab
1274
1275
1276
1277G4double
1278G4NuclNuclDiffuseElastic::ThetaLabToThetaCMS( const G4DynamicParticle* aParticle,
1279 G4double tmass, G4double thetaLab)
1280{
1281 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1282 G4double m1 = theParticle->GetPDGMass();
1283 G4double plab = aParticle->GetTotalMomentum();
1284 G4LorentzVector lv1 = aParticle->Get4Momentum();
1285 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1286
1287 lv += lv1;
1288
1289 G4ThreeVector bst = lv.boostVector();
1290
1291 // lv1.boost(-bst);
1292
1293 // G4ThreeVector p1 = lv1.vect();
1294 // G4double ptot = p1.mag();
1295
1296 G4double phi = G4UniformRand()*twopi;
1297 G4double cost = std::cos(thetaLab);
1298 G4double sint;
1299
1300 if( cost >= 1.0 )
1301 {
1302 cost = 1.0;
1303 sint = 0.0;
1304 }
1305 else if( cost <= -1.0)
1306 {
1307 cost = -1.0;
1308 sint = 0.0;
1309 }
1310 else
1311 {
1312 sint = std::sqrt((1.0-cost)*(1.0+cost));
1313 }
1314 if (verboseLevel>1)
1315 {
1316 G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
1317 }
1318 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1319 v1 *= plab;
1320 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
1321
1322 nlv1.boost(-bst);
1323
1324 G4ThreeVector np1 = nlv1.vect();
1325
1326
1327 G4double thetaCMS = np1.theta();
1328
1329 return thetaCMS;
1330}
1331
1332///////////////////////////////////////////////////////////////////////////////
1333//
1334// Test for given particle and element table of momentum, angle probability.
1335// For the moment in lab system.
1336
1337void G4NuclNuclDiffuseElastic::TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
1338 G4double Z, G4double A)
1339{
1340 fAtomicNumber = Z; // atomic number
1341 fAtomicWeight = A; // number of nucleons
1342 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
1343
1344
1345
1346 G4cout<<"G4NuclNuclDiffuseElastic::TestAngleTable() init the element with Z = "
1347 <<Z<<"; and A = "<<A<<G4endl;
1348
1349 fElementNumberVector.push_back(fAtomicNumber);
1350
1351
1352
1353
1354 G4int i=0, j;
1355 G4double a = 0., z = theParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
1356 G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
1357 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1358 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1359 G4double epsilon = 0.001;
1360
1361 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1362
1363 fAngleTable = new G4PhysicsTable(fEnergyBin);
1364
1365 fWaveVector = partMom/hbarc;
1366
1367 G4double kR = fWaveVector*fNuclearRadius;
1368 G4double kR2 = kR*kR;
1369 G4double kRmax = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
1370 G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1
1371
1372 alphaMax = kRmax*kRmax/kR2;
1373
1374 if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2
1375
1376 alphaCoulomb = kRcoul*kRcoul/kR2;
1377
1378 if( z )
1379 {
1380 a = partMom/m1; // beta*gamma for m1
1381 fBeta = a/std::sqrt(1+a*a);
1382 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1383 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1384 }
1385 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1386
1387 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1388
1389
1390 fAddCoulomb = false;
1391
1392 for(j = 1; j < fAngleBin; j++)
1393 {
1394 // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1395 // alpha2 = angleBins->GetLowEdgeEnergy(j);
1396
1397 alpha1 = alphaMax*(j-1)/fAngleBin;
1398 alpha2 = alphaMax*( j )/fAngleBin;
1399
1400 if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1401
1402 deltaL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1403 deltaL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1404 deltaAG = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction,
1405 alpha1, alpha2,epsilon);
1406
1407 // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1408 // <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl;
1409
1410 sumL10 += deltaL10;
1411 sumL96 += deltaL96;
1412 sumAG += deltaAG;
1413
1414 G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1415 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1416
1417 angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2
1418 }
1419 fAngleTable->insertAt(i,angleVector);
1420 fAngleBank.push_back(fAngleTable);
1421
1422 /*
1423 // Integral over all angle range - Bad accuracy !!!
1424
1425 sumL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1426 sumL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1427 sumAG = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction,
1428 0., alpha2,epsilon);
1429 G4cout<<G4endl;
1430 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
1431 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1432 */
1433 return;
1434}
1435
1436/////////////////////////////////////////////////////////////////
1437//
1438//
1439
1440 G4double G4NuclNuclDiffuseElastic::GetLegendrePol(G4int n, G4double theta)
1441{
1442 G4double legPol, epsilon = 1.e-6;
1443 G4double x = std::cos(theta);
1444
1445 if ( n < 0 ) legPol = 0.;
1446 else if( n == 0 ) legPol = 1.;
1447 else if( n == 1 ) legPol = x;
1448 else if( n == 2 ) legPol = (3.*x*x-1.)/2.;
1449 else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/2.;
1450 else if( n == 4 ) legPol = (35.*x*x*x*x-30.*x*x+3.)/8.;
1451 else if( n == 5 ) legPol = (63.*x*x*x*x*x-70.*x*x*x+15.*x)/8.;
1452 else if( n == 6 ) legPol = (231.*x*x*x*x*x*x-315.*x*x*x*x+105.*x*x-5.)/16.;
1453 else
1454 {
1455 // legPol = ( (2*n-1)*x*GetLegendrePol(n-1,x) - (n-1)*GetLegendrePol(n-2,x) )/n;
1456
1457 legPol = std::sqrt( 2./(n*CLHEP::pi*std::sin(theta+epsilon)) )*std::sin( (n+0.5)*theta+0.25*CLHEP::pi );
1458 }
1459 return legPol;
1460}
1461
1462/////////////////////////////////////////////////////////////////
1463//
1464//
1465
1466G4complex G4NuclNuclDiffuseElastic::GetErfComp(G4complex z, G4int nMax)
1467{
1468 G4int n;
1469 G4double n2, cofn, shny, chny, fn, gn;
1470
1471 G4double x = z.real();
1472 G4double y = z.imag();
1473
1474 G4double outRe = 0., outIm = 0.;
1475
1476 G4double twox = 2.*x;
1477 G4double twoxy = twox*y;
1478 G4double twox2 = twox*twox;
1479
1480 G4double cof1 = G4Exp(-x*x)/CLHEP::pi;
1481
1482 G4double cos2xy = std::cos(twoxy);
1483 G4double sin2xy = std::sin(twoxy);
1484
1485 G4double twoxcos2xy = twox*cos2xy;
1486 G4double twoxsin2xy = twox*sin2xy;
1487
1488 for( n = 1; n <= nMax; n++)
1489 {
1490 n2 = n*n;
1491
1492 cofn = G4Exp(-0.5*n2)/(n2+twox2); // /(n2+0.5*twox2);
1493
1494 chny = std::cosh(n*y);
1495 shny = std::sinh(n*y);
1496
1497 fn = twox - twoxcos2xy*chny + n*sin2xy*shny;
1498 gn = twoxsin2xy*chny + n*cos2xy*shny;
1499
1500 fn *= cofn;
1501 gn *= cofn;
1502
1503 outRe += fn;
1504 outIm += gn;
1505 }
1506 outRe *= 2*cof1;
1507 outIm *= 2*cof1;
1508
1509 if(std::abs(x) < 0.0001)
1510 {
1511 outRe += GetErf(x);
1512 outIm += cof1*y;
1513 }
1514 else
1515 {
1516 outRe += GetErf(x) + cof1*(1-cos2xy)/twox;
1517 outIm += cof1*sin2xy/twox;
1518 }
1519 return G4complex(outRe, outIm);
1520}
1521
1522
1523/////////////////////////////////////////////////////////////////
1524//
1525//
1526
1527G4complex G4NuclNuclDiffuseElastic::GetErfInt(G4complex z) // , G4int nMax)
1528{
1529 G4double outRe, outIm;
1530
1531 G4double x = z.real();
1532 G4double y = z.imag();
1533 fReZ = x;
1534
1535 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1536
1537 outRe = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpSin, 0., y );
1538 outIm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpCos, 0., y );
1539
1540 outRe *= 2./std::sqrt(CLHEP::pi);
1541 outIm *= 2./std::sqrt(CLHEP::pi);
1542
1543 outRe += GetErf(x);
1544
1545 return G4complex(outRe, outIm);
1546}
1547
1548/////////////////////////////////////////////////////////////////
1549//
1550//
1551
1552
1553G4complex G4NuclNuclDiffuseElastic::GammaLess(G4double theta)
1554{
1555 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1556 G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1557
1558 G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1559 G4double kappa = u/std::sqrt(CLHEP::pi);
1560 G4double dTheta = theta - fRutherfordTheta;
1561 u *= dTheta;
1562 G4double u2 = u*u;
1563 G4double u2m2p3 = u2*2./3.;
1564
1565 G4complex im = G4complex(0.,1.);
1566 G4complex order = G4complex(u,u);
1567 order /= std::sqrt(2.);
1568
1569 G4complex gamma = CLHEP::pi*kappa*GetErfcInt(-order)*std::exp(im*(u*u+0.25*CLHEP::pi));
1570 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1571 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1572 G4complex out = gamma*(1. - a1*dTheta) - a0;
1573
1574 return out;
1575}
1576
1577/////////////////////////////////////////////////////////////////
1578//
1579//
1580
1581G4complex G4NuclNuclDiffuseElastic::GammaMore(G4double theta)
1582{
1583 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1584 G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1585
1586 G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1587 G4double kappa = u/std::sqrt(CLHEP::pi);
1588 G4double dTheta = theta - fRutherfordTheta;
1589 u *= dTheta;
1590 G4double u2 = u*u;
1591 G4double u2m2p3 = u2*2./3.;
1592
1593 G4complex im = G4complex(0.,1.);
1594 G4complex order = G4complex(u,u);
1595 order /= std::sqrt(2.);
1596 G4complex gamma = CLHEP::pi*kappa*GetErfcInt(order)*std::exp(im*(u*u+0.25*CLHEP::pi));
1597 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1598 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1599 G4complex out = -gamma*(1. - a1*dTheta) - a0;
1600
1601 return out;
1602}
1603
1604/////////////////////////////////////////////////////////////////
1605//
1606//
1607
1608 G4complex G4NuclNuclDiffuseElastic::AmplitudeNear(G4double theta)
1609{
1610 G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
1611 G4complex out = G4complex(kappa/fWaveVector,0.);
1612
1613 out *= PhaseNear(theta);
1614
1615 if( theta <= fRutherfordTheta )
1616 {
1617 out *= GammaLess(theta) + ProfileNear(theta);
1618 // out *= GammaMore(theta) + ProfileNear(theta);
1619 out += CoulombAmplitude(theta);
1620 }
1621 else
1622 {
1623 out *= GammaMore(theta) + ProfileNear(theta);
1624 // out *= GammaLess(theta) + ProfileNear(theta);
1625 }
1626 return out;
1627}
1628
1629/////////////////////////////////////////////////////////////////
1630//
1631//
1632
1633G4complex G4NuclNuclDiffuseElastic::AmplitudeSim(G4double theta)
1634{
1635 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1636 G4double dTheta = 0.5*(theta - fRutherfordTheta);
1637 G4double sindTheta = std::sin(dTheta);
1638 G4double persqrt2 = std::sqrt(0.5);
1639
1640 G4complex order = G4complex(persqrt2,persqrt2);
1641 order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*sindTheta;
1642 // order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*dTheta;
1643
1644 G4complex out;
1645
1646 if ( theta <= fRutherfordTheta )
1647 {
1648 out = 1. - 0.5*GetErfcInt(-order)*ProfileNear(theta);
1649 }
1650 else
1651 {
1652 out = 0.5*GetErfcInt(order)*ProfileNear(theta);
1653 }
1654
1655 out *= CoulombAmplitude(theta);
1656
1657 return out;
1658}
1659
1660/////////////////////////////////////////////////////////////////
1661//
1662//
1663
1664 G4complex G4NuclNuclDiffuseElastic::AmplitudeGla(G4double theta)
1665{
1666 G4int n;
1667 G4double T12b, b, b2; // cosTheta = std::cos(theta);
1668 G4complex out = G4complex(0.,0.), shiftC, shiftN;
1669 G4complex im = G4complex(0.,1.);
1670
1671 for( n = 0; n < fMaxL; n++)
1672 {
1673 shiftC = std::exp( im*2.*CalculateCoulombPhase(n) );
1674 // b = ( fZommerfeld + std::sqrt( fZommerfeld*fZommerfeld + n*(n+1) ) )/fWaveVector;
1675 b = ( std::sqrt( G4double(n*(n+1)) ) )/fWaveVector;
1676 b2 = b*b;
1677 T12b = fSumSigma*G4Exp(-b2/fNuclearRadiusSquare)/CLHEP::pi/fNuclearRadiusSquare;
1678 shiftN = std::exp( -0.5*(1.-im*fEtaRatio)*T12b ) - 1.;
1679 out += (2.*n+1.)*shiftC*shiftN*GetLegendrePol(n, theta);
1680 }
1681 out /= 2.*im*fWaveVector;
1682 out += CoulombAmplitude(theta);
1683 return out;
1684}
1685
1686
1687/////////////////////////////////////////////////////////////////
1688//
1689//
1690
1691 G4complex G4NuclNuclDiffuseElastic::AmplitudeGG(G4double theta)
1692{
1693 G4int n;
1694 G4double T12b, a, aTemp, b2, sinThetaH = std::sin(0.5*theta);
1695 G4double sinThetaH2 = sinThetaH*sinThetaH;
1696 G4complex out = G4complex(0.,0.);
1697 G4complex im = G4complex(0.,1.);
1698
1699 a = -fSumSigma/CLHEP::twopi/fNuclearRadiusSquare;
1700 b2 = fWaveVector*fWaveVector*fNuclearRadiusSquare*sinThetaH2;
1701
1702 aTemp = a;
1703
1704 for( n = 1; n < fMaxL; n++)
1705 {
1706 T12b = aTemp*G4Exp(-b2/n)/n;
1707 aTemp *= a;
1708 out += T12b;
1709 G4cout<<"out = "<<out<<G4endl;
1710 }
1711 out *= -4.*im*fWaveVector/CLHEP::pi;
1712 out += CoulombAmplitude(theta);
1713 return out;
1714}
1715
1716
1717///////////////////////////////////////////////////////////////////////////////
1718//
1719// Test for given particle and element table of momentum, angle probability.
1720// For the partMom in CMS.
1721
1722void G4NuclNuclDiffuseElastic::InitParameters(const G4ParticleDefinition* theParticle,
1723 G4double partMom, G4double Z, G4double A)
1724{
1725 fAtomicNumber = Z; // atomic number
1726 fAtomicWeight = A; // number of nucleons
1727
1728 fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight);
1729 G4double A1 = G4double( theParticle->GetBaryonNumber() );
1730 fNuclearRadius1 = CalculateNuclearRad(A1);
1731 // fNuclearRadius = std::sqrt(fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2);
1732 fNuclearRadius = fNuclearRadius1 + fNuclearRadius2;
1733
1734 G4double a = 0.;
1735 G4double z = theParticle->GetPDGCharge();
1736 G4double m1 = theParticle->GetPDGMass();
1737
1738 fWaveVector = partMom/CLHEP::hbarc;
1739
1740 G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
1741 G4cout<<"kR = "<<lambda<<G4endl;
1742
1743 if( z )
1744 {
1745 a = partMom/m1; // beta*gamma for m1
1746 fBeta = a/std::sqrt(1+a*a);
1747 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1748 fRutherfordRatio = fZommerfeld/fWaveVector;
1749 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1750 }
1751 G4cout<<"fZommerfeld = "<<fZommerfeld<<G4endl;
1752 fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
1753 G4cout<<"fProfileLambda = "<<fProfileLambda<<G4endl;
1754 fProfileDelta = fCofDelta*fProfileLambda;
1755 fProfileAlpha = fCofAlpha*fProfileLambda;
1756
1757 CalculateCoulombPhaseZero();
1758 CalculateRutherfordAnglePar();
1759
1760 return;
1761}
1762///////////////////////////////////////////////////////////////////////////////
1763//
1764// Test for given particle and element table of momentum, angle probability.
1765// For the partMom in CMS.
1766
1767void G4NuclNuclDiffuseElastic::InitDynParameters(const G4ParticleDefinition* theParticle,
1768 G4double partMom)
1769{
1770 G4double a = 0.;
1771 G4double z = theParticle->GetPDGCharge();
1772 G4double m1 = theParticle->GetPDGMass();
1773
1774 fWaveVector = partMom/CLHEP::hbarc;
1775
1776 G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
1777
1778 if( z )
1779 {
1780 a = partMom/m1; // beta*gamma for m1
1781 fBeta = a/std::sqrt(1+a*a);
1782 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1783 fRutherfordRatio = fZommerfeld/fWaveVector;
1784 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1785 }
1786 fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
1787 fProfileDelta = fCofDelta*fProfileLambda;
1788 fProfileAlpha = fCofAlpha*fProfileLambda;
1789
1790 CalculateCoulombPhaseZero();
1791 CalculateRutherfordAnglePar();
1792
1793 return;
1794}
1795
1796
1797///////////////////////////////////////////////////////////////////////////////
1798//
1799// Test for given particle and element table of momentum, angle probability.
1800// For the partMom in CMS.
1801
1802void G4NuclNuclDiffuseElastic::InitParametersGla(const G4DynamicParticle* aParticle,
1803 G4double partMom, G4double Z, G4double A)
1804{
1805 fAtomicNumber = Z; // target atomic number
1806 fAtomicWeight = A; // target number of nucleons
1807
1808 fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight); // target nucleus radius
1809 G4double A1 = G4double( aParticle->GetDefinition()->GetBaryonNumber() );
1810 fNuclearRadius1 = CalculateNuclearRad(A1); // projectile nucleus radius
1811 fNuclearRadiusSquare = fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2;
1812
1813
1814 G4double a = 0., kR12;
1815 G4double z = aParticle->GetDefinition()->GetPDGCharge();
1816 G4double m1 = aParticle->GetDefinition()->GetPDGMass();
1817
1818 fWaveVector = partMom/CLHEP::hbarc;
1819
1820 G4double pN = A1 - z;
1821 if( pN < 0. ) pN = 0.;
1822
1823 G4double tN = A - Z;
1824 if( tN < 0. ) tN = 0.;
1825
1826 G4double pTkin = aParticle->GetKineticEnergy();
1827 pTkin /= A1;
1828
1829
1830 fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
1831 (z*tN+pN*Z)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
1832
1833 G4cout<<"fSumSigma = "<<fSumSigma/CLHEP::millibarn<<" mb"<<G4endl;
1834 G4cout<<"pi*R2 = "<<CLHEP::pi*fNuclearRadiusSquare/CLHEP::millibarn<<" mb"<<G4endl;
1835 kR12 = fWaveVector*std::sqrt(fNuclearRadiusSquare);
1836 G4cout<<"k*sqrt(R2) = "<<kR12<<" "<<G4endl;
1837 fMaxL = (G4int(kR12)+1)*4;
1838 G4cout<<"fMaxL = "<<fMaxL<<" "<<G4endl;
1839
1840 if( z )
1841 {
1842 a = partMom/m1; // beta*gamma for m1
1843 fBeta = a/std::sqrt(1+a*a);
1844 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1845 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1846 }
1847
1848 CalculateCoulombPhaseZero();
1849
1850
1851 return;
1852}
1853
1854
1855/////////////////////////////////////////////////////////////////////////////////////
1856//
1857// Returns nucleon-nucleon cross-section based on N. Starkov parametrisation of
1858// data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
1859// projectile nucleon is pParticle with pTkin shooting target nucleon tParticle
1860
1861G4double
1862G4NuclNuclDiffuseElastic::GetHadronNucleonXscNS( G4ParticleDefinition* pParticle,
1863 G4double pTkin,
1864 G4ParticleDefinition* tParticle)
1865{
1866 G4double xsection(0), /*Delta,*/ A0, B0;
1867 G4double hpXsc(0);
1868 G4double hnXsc(0);
1869
1870
1871 G4double targ_mass = tParticle->GetPDGMass();
1872 G4double proj_mass = pParticle->GetPDGMass();
1873
1874 G4double proj_energy = proj_mass + pTkin;
1875 G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
1876
1877 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
1878
1879 sMand /= CLHEP::GeV*CLHEP::GeV; // in GeV for parametrisation
1880 proj_momentum /= CLHEP::GeV;
1881 proj_energy /= CLHEP::GeV;
1882 proj_mass /= CLHEP::GeV;
1883 G4double logS = G4Log(sMand);
1884
1885 // General PDG fit constants
1886
1887
1888 // fEtaRatio=Re[f(0)]/Im[f(0)]
1889
1890 if( proj_momentum >= 1.2 )
1891 {
1892 fEtaRatio = 0.13*(logS - 5.8579332)*G4Pow::GetInstance()->powA(sMand,-0.18);
1893 }
1894 else if( proj_momentum >= 0.6 )
1895 {
1896 fEtaRatio = -75.5*(G4Pow::GetInstance()->powA(proj_momentum,0.25)-0.95)/
1897 (G4Pow::GetInstance()->powA(3*proj_momentum,2.2)+1);
1898 }
1899 else
1900 {
1901 fEtaRatio = 15.5*proj_momentum/(27*proj_momentum*proj_momentum*proj_momentum+2);
1902 }
1903 G4cout<<"fEtaRatio = "<<fEtaRatio<<G4endl;
1904
1905 // xsc
1906
1907 if( proj_momentum >= 10. ) // high energy: pp = nn = np
1908 // if( proj_momentum >= 2.)
1909 {
1910 //Delta = 1.;
1911
1912 //if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
1913
1914 //AR-12Aug2016 if( proj_momentum >= 10.)
1915 {
1916 B0 = 7.5;
1917 A0 = 100. - B0*G4Log(3.0e7);
1918
1919 xsection = A0 + B0*G4Log(proj_energy) - 11
1920 + 103*G4Pow::GetInstance()->powA(2*0.93827*proj_energy + proj_mass*proj_mass+
1921 0.93827*0.93827,-0.165); // mb
1922 }
1923 }
1924 else // low energy pp = nn != np
1925 {
1926 if(pParticle == tParticle) // pp or nn // nn to be pp
1927 {
1928 if( proj_momentum < 0.73 )
1929 {
1930 hnXsc = 23 + 50*( G4Pow::GetInstance()->powA( G4Log(0.73/proj_momentum), 3.5 ) );
1931 }
1932 else if( proj_momentum < 1.05 )
1933 {
1934 hnXsc = 23 + 40*(G4Log(proj_momentum/0.73))*
1935 (G4Log(proj_momentum/0.73));
1936 }
1937 else // if( proj_momentum < 10. )
1938 {
1939 hnXsc = 39.0 +
1940 75*(proj_momentum - 1.2)/(G4Pow::GetInstance()->powA(proj_momentum,3.0) + 0.15);
1941 }
1942 xsection = hnXsc;
1943 }
1944 else // pn to be np
1945 {
1946 if( proj_momentum < 0.8 )
1947 {
1948 hpXsc = 33+30*G4Pow::GetInstance()->powA(G4Log(proj_momentum/1.3),4.0);
1949 }
1950 else if( proj_momentum < 1.4 )
1951 {
1952 hpXsc = 33+30*G4Pow::GetInstance()->powA(G4Log(proj_momentum/0.95),2.0);
1953 }
1954 else // if( proj_momentum < 10. )
1955 {
1956 hpXsc = 33.3+
1957 20.8*(G4Pow::GetInstance()->powA(proj_momentum,2.0)-1.35)/
1958 (G4Pow::GetInstance()->powA(proj_momentum,2.50)+0.95);
1959 }
1960 xsection = hpXsc;
1961 }
1962 }
1963 xsection *= CLHEP::millibarn; // parametrised in mb
1964 G4cout<<"xsection = "<<xsection/CLHEP::millibarn<<" mb"<<G4endl;
1965 return xsection;
1966}
1967
1968/////////////////////////////////////////////////////////////////
1969//
1970// The ratio el/ruth for Fresnel smooth nucleus profile
1971
1972G4double G4NuclNuclDiffuseElastic::GetRatioGen(G4double theta)
1973{
1974 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1975 G4double dTheta = 0.5*(theta - fRutherfordTheta);
1976 G4double sindTheta = std::sin(dTheta);
1977
1978 G4double prof = Profile(theta);
1979 G4double prof2 = prof*prof;
1980 // G4double profmod = std::abs(prof);
1981 G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
1982
1983 order = std::abs(order); // since sin changes sign!
1984 // G4cout<<"order = "<<order<<G4endl;
1985
1986 G4double cosFresnel = GetCint(order);
1987 G4double sinFresnel = GetSint(order);
1988
1989 G4double out;
1990
1991 if ( theta <= fRutherfordTheta )
1992 {
1993 out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1994 out += ( cosFresnel + sinFresnel - 1. )*prof;
1995 }
1996 else
1997 {
1998 out = 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1999 }
2000
2001 return out;
2002}
2003
2004///////////////////////////////////////////////////////////////////
2005//
2006// For the calculation of arg Gamma(z) one needs complex extension
2007// of ln(Gamma(z))
2008
2009G4complex G4NuclNuclDiffuseElastic::GammaLogarithm(G4complex zz)
2010{
2011 const G4double cof[6] = { 76.18009172947146, -86.50532032941677,
2012 24.01409824083091, -1.231739572450155,
2013 0.1208650973866179e-2, -0.5395239384953e-5 } ;
2014 G4int j;
2015 G4complex z = zz - 1.0;
2016 G4complex tmp = z + 5.5;
2017 tmp -= (z + 0.5) * std::log(tmp);
2018 G4complex ser = G4complex(1.000000000190015,0.);
2019
2020 for ( j = 0; j <= 5; j++ )
2021 {
2022 z += 1.0;
2023 ser += cof[j]/z;
2024 }
2025 return -tmp + std::log(2.5066282746310005*ser);
2026}
2027
2028/////////////////////////////////////////////////////////////
2029//
2030// Bessel J0 function based on rational approximation from
2031// J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
2032
2033G4double G4NuclNuclDiffuseElastic::BesselJzero(G4double value)
2034{
2035 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
2036
2037 modvalue = std::fabs(value);
2038
2039 if ( value < 8.0 && value > -8.0 )
2040 {
2041 value2 = value*value;
2042
2043 fact1 = 57568490574.0 + value2*(-13362590354.0
2044 + value2*( 651619640.7
2045 + value2*(-11214424.18
2046 + value2*( 77392.33017
2047 + value2*(-184.9052456 ) ) ) ) );
2048
2049 fact2 = 57568490411.0 + value2*( 1029532985.0
2050 + value2*( 9494680.718
2051 + value2*(59272.64853
2052 + value2*(267.8532712
2053 + value2*1.0 ) ) ) );
2054
2055 bessel = fact1/fact2;
2056 }
2057 else
2058 {
2059 arg = 8.0/modvalue;
2060
2061 value2 = arg*arg;
2062
2063 shift = modvalue-0.785398164;
2064
2065 fact1 = 1.0 + value2*(-0.1098628627e-2
2066 + value2*(0.2734510407e-4
2067 + value2*(-0.2073370639e-5
2068 + value2*0.2093887211e-6 ) ) );
2069
2070 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
2071 + value2*(-0.6911147651e-5
2072 + value2*(0.7621095161e-6
2073 - value2*0.934945152e-7 ) ) );
2074
2075 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
2076 }
2077 return bessel;
2078}
2079
2080/////////////////////////////////////////////////////////////
2081//
2082// Bessel J1 function based on rational approximation from
2083// J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
2084
2085G4double G4NuclNuclDiffuseElastic::BesselJone(G4double value)
2086{
2087 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
2088
2089 modvalue = std::fabs(value);
2090
2091 if ( modvalue < 8.0 )
2092 {
2093 value2 = value*value;
2094
2095 fact1 = value*(72362614232.0 + value2*(-7895059235.0
2096 + value2*( 242396853.1
2097 + value2*(-2972611.439
2098 + value2*( 15704.48260
2099 + value2*(-30.16036606 ) ) ) ) ) );
2100
2101 fact2 = 144725228442.0 + value2*(2300535178.0
2102 + value2*(18583304.74
2103 + value2*(99447.43394
2104 + value2*(376.9991397
2105 + value2*1.0 ) ) ) );
2106 bessel = fact1/fact2;
2107 }
2108 else
2109 {
2110 arg = 8.0/modvalue;
2111
2112 value2 = arg*arg;
2113
2114 shift = modvalue - 2.356194491;
2115
2116 fact1 = 1.0 + value2*( 0.183105e-2
2117 + value2*(-0.3516396496e-4
2118 + value2*(0.2457520174e-5
2119 + value2*(-0.240337019e-6 ) ) ) );
2120
2121 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
2122 + value2*( 0.8449199096e-5
2123 + value2*(-0.88228987e-6
2124 + value2*0.105787412e-6 ) ) );
2125
2126 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
2127
2128 if (value < 0.0) bessel = -bessel;
2129 }
2130 return bessel;
2131}
2132
2133//
2134//
2135/////////////////////////////////////////////////////////////////////////////////
epsilon
G4double epsilon(G4double density, G4double temperature)
Definition: G4DNAElectronHoleRecombination.cc:95
G4ElementTable
std::vector< G4Element * > G4ElementTable
Definition: G4ElementTable.hh:44
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
a0
const G4double a0
Definition: G4IonCoulombCrossSection.cc:68
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:227
G4double
double G4double
Definition: G4Types.hh:83
G4complex
std::complex< G4double > G4complex
Definition: G4Types.hh:88
G4int
int G4int
Definition: G4Types.hh:85
Z
const G4int Z[17]
Definition: G4WaterStopping.cc:51
A
const G4double A[17]
Definition: G4WaterStopping.cc:53
alpha2
const G4double alpha2
Definition: G4WentzelOKandVIxSection.cc:70
G4endl
#define G4endl
Definition: G4ios.hh:57
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:88
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:93
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::Get4Momentum
G4LorentzVector Get4Momentum() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4DynamicParticle::GetTotalMomentum
G4double GetTotalMomentum() const
G4Element::GetElementTable
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:403
G4Element::GetNumberOfElements
static size_t GetNumberOfElements()
Definition: G4Element.cc:410
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
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:54
G4He3::He3
static G4He3 * He3()
Definition: G4He3.cc:93
G4IonTable::GetIon
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:522
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:103
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:62
G4NistManager::GetAtomicMassAmu
G4double GetAtomicMassAmu(const G4String &symb) const
Definition: G4NistManager.hh:371
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:1278
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:1129
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:1337
G4NuclNuclDiffuseElastic::GetRatioGen
G4double GetRatioGen(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:1972
G4NuclNuclDiffuseElastic::Profile
G4double Profile(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:908
G4NuclNuclDiffuseElastic::AmplitudeSim
G4complex AmplitudeSim(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:1633
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:978
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:1664
G4NuclNuclDiffuseElastic::GetErfComp
G4complex GetErfComp(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.cc:1466
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:1440
G4NuclNuclDiffuseElastic::GetIntegrandFunction
G4double GetIntegrandFunction(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:668
G4NuclNuclDiffuseElastic::AmplitudeNear
G4complex AmplitudeNear(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:1608
G4NuclNuclDiffuseElastic::AmplitudeGG
G4complex AmplitudeGG(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:1691
G4NuclNuclDiffuseElastic::GetHadronNucleonXscNS
G4double GetHadronNucleonXscNS(G4ParticleDefinition *pParticle, G4double pTkin, G4ParticleDefinition *tParticle)
Definition: G4NuclNuclDiffuseElastic.cc:1862
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:855
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:1802
G4NuclNuclDiffuseElastic::PhaseNear
G4complex PhaseNear(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:927
G4NuclNuclDiffuseElastic::GammaLogarithm
G4complex GammaLogarithm(G4complex xx)
Definition: G4NuclNuclDiffuseElastic.cc:2009
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:1086
G4NuclNuclDiffuseElastic::GetDiffElasticSumProbA
G4double GetDiffElasticSumProbA(G4double alpha)
Definition: G4NuclNuclDiffuseElastic.cc:572
G4NuclNuclDiffuseElastic::ThetaCMStoThetaLab
G4double ThetaCMStoThetaLab(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaCMS)
Definition: G4NuclNuclDiffuseElastic.cc:1217
G4NuclNuclDiffuseElastic::SampleTableT
G4double SampleTableT(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:840
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:1767
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:1722
G4NuclNuclDiffuseElastic::GetErfcInt
G4complex GetErfcInt(G4complex z)
Definition: G4NuclNuclDiffuseElastic.hh:672
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:777
G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc
G4double GetFresnelIntegrandXsc(G4double alpha)
Definition: G4NuclNuclDiffuseElastic.hh:1035
G4NuclNuclDiffuseElastic::SampleCoulombMuCMS
G4double SampleCoulombMuCMS(const G4ParticleDefinition *aParticle, G4double p)
Definition: G4NuclNuclDiffuseElastic.cc:809
G4NuclNuclDiffuseElastic::GammaMore
G4complex GammaMore(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:1581
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:1553
G4NucleiProperties::GetNuclearMass
static G4double GetNuclearMass(const G4double A, const G4double Z)
Definition: G4NucleiProperties.cc:50
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:113
G4ParticleTable::GetIonTable
G4IonTable * GetIonTable() const
Definition: G4ParticleTable.cc:698
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:87
G4PhysicsFreeVector::PutValue
void PutValue(const std::size_t index, const G4double e, const G4double value)
Definition: G4PhysicsFreeVector.hh:95
G4PhysicsVector::GetLowEdgeEnergy
G4double GetLowEdgeEnergy(const std::size_t index) const
G4PionMinus::PionMinus
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:97
G4PionPlus::PionPlus
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:97
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:92
G4Triton::Triton
static G4Triton * Triton()
Definition: G4Triton.cc:93
W
#define W
Definition: crc32.c:84