Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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