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