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