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