| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28 // Physics model class G4DiffuseElastic
29 //
30 //
31 // G4 Model: optical diffuse elastic scatterin
32 //
33 // 24-May-07 V. Grichine
34 //
35 // 21.10.15 V. Grichine
36 // Bug fixed in BuildAngleTable, i
37 // angle bins at high energies > 5
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
75 G4DiffuseElastic::G4DiffuseElastic()
76 : G4HadronElastic("DiffuseElastic"), fPartic
77 {
78 SetMinEnergy( 0.01*MeV );
79 SetMaxEnergy( G4HadronicParameters::Instance
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 ori
96 fAngleBin = 200;
97
98 fEnergyVector = new G4PhysicsLogVector( the
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
117 G4DiffuseElastic::~G4DiffuseElastic()
118 {
119 if ( fEnergyVector )
120 {
121 delete fEnergyVector;
122 fEnergyVector = 0;
123 }
124 for ( std::vector<G4PhysicsTable*>::iterator
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 ele
138
139 void G4DiffuseElastic::Initialise()
140 {
141
142 // fEnergyVector = new G4PhysicsLogVector( t
143
144 const G4ElementTable* theElementTable = G4El
145
146 std::size_t jEl, numOfEl = G4Element::GetNum
147
148 for( jEl = 0; jEl < numOfEl; ++jEl) // appli
149 {
150 fAtomicNumber = (*theElementTable)[jEl]->G
151 fAtomicWeight = G4NistManager::Instance()-
152 fNuclearRadius = CalculateNuclearRad(fAtom
153
154 if( verboseLevel > 0 )
155 {
156 G4cout<<"G4DiffuseElastic::Initialise()
157 <<(*theElementTable)[jEl]->GetName()<<G4
158 }
159 fElementNumberVector.push_back(fAtomicNumb
160 fElementNameVector.push_back((*theElementT
161
162 BuildAngleTable();
163 fAngleBank.push_back(fAngleTable);
164 }
165 return;
166 }
167
168 //////////////////////////////////////////////
169 //
170 // return differential elastic cross section d
171
172 G4double
173 G4DiffuseElastic::GetDiffuseElasticXsc( const
174 G4doub
175 G4double momentum,
176 G4doub
177 {
178 fParticle = particle;
179 fWaveVector = momentum/hbarc;
180 fAtomicWeight = A;
181 fAddCoulomb = false;
182 fNuclearRadius = CalculateNuclearRad(A);
183
184 G4double sigma = fNuclearRadius*fNuclearRadi
185
186 return sigma;
187 }
188
189 //////////////////////////////////////////////
190 //
191 // return invariant differential elastic cross
192
193 G4double
194 G4DiffuseElastic::GetInvElasticXsc( const G4Pa
195 G4doub
196 G4double plab,
197 G4doub
198 {
199 G4double m1 = particle->GetPDGMass();
200 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla
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 = thePro
207 else if (iZ == 1 && iA == 2) theDef = theDeu
208 else if (iZ == 1 && iA == 3) theDef = G4Trit
209 else if (iZ == 2 && iA == 3) theDef = G4He3:
210 else if (iZ == 2 && iA == 4) theDef = theAlp
211 else theDef = G4ParticleTable::GetParticleTa
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)/pto
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( parti
232
233 sigma *= pi/ptot2;
234
235 return sigma;
236 }
237
238 //////////////////////////////////////////////
239 //
240 // return differential elastic cross section d
241 // correction
242
243 G4double
244 G4DiffuseElastic::GetDiffuseElasticSumXsc( con
245 G4doub
246 G4double momentum,
247 G4doub
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*
259 G4double kRtC = 1.9;
260
261 if( z && (kRt > kRtC) )
262 {
263 fAddCoulomb = true;
264 fBeta = CalculateParticleBeta( parti
265 fZommerfeld = CalculateZommerfeld( fBeta,
266 fAm = CalculateAm( momentum, fZomm
267 }
268 G4double sigma = fNuclearRadius*fNuclearRadi
269
270 return sigma;
271 }
272
273 //////////////////////////////////////////////
274 //
275 // return invariant differential elastic cross
276 // correction
277
278 G4double
279 G4DiffuseElastic::GetInvElasticSumXsc( const G
280 G4doub
281 G4double plab,
282 G4doub
283 {
284 G4double m1 = particle->GetPDGMass();
285
286 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla
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 = thePro
294 else if (iZ == 1 && iA == 2) theDef = theDeu
295 else if (iZ == 1 && iA == 3) theDef = G4Trit
296 else if (iZ == 2 && iA == 3) theDef = G4He3:
297 else if (iZ == 2 && iA == 4) theDef = theAlp
298 else theDef = G4ParticleTable::GetParticleTa
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)/
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( pa
319
320 sigma *= pi/ptot2;
321
322 return sigma;
323 }
324
325 //////////////////////////////////////////////
326 //
327 // return invariant differential elastic cross
328 // correction
329
330 G4double
331 G4DiffuseElastic::GetInvCoulombElasticXsc( con
332 G4doub
333 G4double plab,
334 G4doub
335 {
336 G4double m1 = particle->GetPDGMass();
337 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla
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 = thePro
344 else if (iZ == 1 && iA == 2) theDef = theDeu
345 else if (iZ == 1 && iA == 3) theDef = G4Trit
346 else if (iZ == 2 && iA == 3) theDef = G4He3:
347 else if (iZ == 2 && iA == 4) theDef = theAlp
348 else theDef = G4ParticleTable::GetParticleTa
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)/pto
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( parti
369
370 sigma *= pi/ptot2;
371
372 return sigma;
373 }
374
375 //////////////////////////////////////////////
376 //
377 // return differential elastic probability d(p
378
379 G4double
380 G4DiffuseElastic::GetDiffElasticProb( // G4Par
381 G4doub
382 // G4double momentum,
383 // G4double A
384 )
385 {
386 G4double sigma, bzero, bzero2, bonebyarg, bo
387 G4double delta, diffuse, gamma;
388 G4double e1, e2, bone, bone2;
389
390 // G4double wavek = momentum/hbarc; // wave
391 // G4double r0 = 1.08*fermi;
392 // G4double rad = r0*G4Pow::GetInstance()-
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;
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;
435
436 G4double kgamma = lambda*(1.-G4Exp(-fWave
437 G4double kgamma2 = kgamma*kgamma;
438
439 // G4double dk2t = delta*fWaveVector*fWaveV
440 // G4double dk2t2 = dk2t*dk2t;
441 // G4double pikdt = pi*fWaveVector*diffuse*t
442
443 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa
444
445 damp = DampFactor(pikdt);
446 damp2 = damp*damp;
447
448 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector
449 G4double e2dk3t = -2.*e2*delta*fWaveVector*
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(p
465 // Coulomb correction
466
467 G4double
468 G4DiffuseElastic::GetDiffElasticSumProb( // G4
469 G4doub
470 // G4double momentum,
471 // G4double A
472 )
473 {
474 G4double sigma, bzero, bzero2, bonebyarg, bo
475 G4double delta, diffuse, gamma;
476 G4double e1, e2, bone, bone2;
477
478 // G4double wavek = momentum/hbarc; // wave
479 // G4double r0 = 1.08*fermi;
480 // G4double rad = r0*G4Pow::GetInstance()-
481
482 G4double kr = fWaveVector*fNuclearRadius;
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;
523 G4double kgamma = lambda*(1.-G4Exp(-fWave
524
525 // G4cout<<"kgamma = "<<kgamma<<G4endl;
526
527 if(fAddCoulomb) // add Coulomb correction
528 {
529 G4double sinHalfTheta = std::sin(0.5*thet
530 G4double sinHalfTheta2 = sinHalfTheta*sinH
531
532 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta
533 // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe
534 }
535
536 G4double kgamma2 = kgamma*kgamma;
537
538 // G4double dk2t = delta*fWaveVector*fWaveV
539 // G4cout<<"dk2t = "<<dk2t<<G4endl;
540 // G4double dk2t2 = dk2t*dk2t;
541 // G4double pikdt = pi*fWaveVector*diffuse*t
542
543 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa
544
545 // G4cout<<"pikdt = "<<pikdt<<G4endl;
546
547 damp = DampFactor(pikdt);
548 damp2 = damp*damp;
549
550 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector
551 G4double e2dk3t = -2.*e2*delta*fWaveVector*
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*fZommerf
560 sigma += kr2*bonebyarg2; // correction at J
561
562 sigma *= damp2; // *rad*rad;
563
564 return sigma;
565 }
566
567
568 //////////////////////////////////////////////
569 //
570 // return differential elastic probability d(p
571 // Coulomb correction. It is called from Build
572
573 G4double
574 G4DiffuseElastic::GetDiffElasticSumProbA( G4do
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, bo
583 G4double delta, diffuse, gamma;
584 G4double e1, e2, bone, bone2;
585
586 // G4double wavek = momentum/hbarc; // wave
587 // G4double r0 = 1.08*fermi;
588 // G4double rad = r0*G4Pow::GetInstance()-
589
590 G4double kr = fWaveVector*fNuclearRadius;
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;
631 G4double kgamma = lambda*(1.-G4Exp(-fWave
632
633 // G4cout<<"kgamma = "<<kgamma<<G4endl;
634
635 if( fAddCoulomb ) // add Coulomb correction
636 {
637 G4double sinHalfTheta = theta*0.5; // std
638 G4double sinHalfTheta2 = sinHalfTheta*sinH
639
640 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta
641 // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe
642 }
643 G4double kgamma2 = kgamma*kgamma;
644
645 // G4double dk2t = delta*fWaveVector*fWaveV
646 // G4cout<<"dk2t = "<<dk2t<<G4endl;
647 // G4double dk2t2 = dk2t*dk2t;
648 // G4double pikdt = pi*fWaveVector*diffuse*t
649
650 G4double pikdt = lambda*(1. - G4Exp( -pi*
651
652 // G4cout<<"pikdt = "<<pikdt<<G4endl;
653
654 damp = DampFactor( pikdt );
655 damp2 = damp*damp;
656
657 G4double mode2k2 = ( e1*e1 + e2*e2 )*fWaveVe
658 G4double e2dk3t = -2.*e2*delta*fWaveVector*
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*fZommerf
667 sigma += kr2*bonebyarg2; // correction at J
668
669 sigma *= damp2; // *rad*rad;
670
671 return sigma;
672 }
673
674
675 //////////////////////////////////////////////
676 //
677 // return differential elastic probability 2*p
678
679 G4double
680 G4DiffuseElastic::GetIntegrandFunction( G4doub
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(sig
694
695 G4double
696 G4DiffuseElastic::IntegralElasticProb( const
697 G4doub
698 G4double momentum,
699 G4doub
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(G4Dif
710
711 // result = integral.Legendre10(this,&G4Diff
712 result = integral.Legendre96(this,&G4Diffuse
713
714 return result;
715 }
716
717 //////////////////////////////////////////////
718 //
719 // Return inv momentum transfer -t > 0
720
721 G4double G4DiffuseElastic::SampleT( const G4Pa
722 {
723 G4double theta = SampleThetaCMS( aParticle,
724 G4double t = 2*p*p*( 1 - std::cos(theta)
725 return t;
726 }
727
728 //////////////////////////////////////////////
729 //
730 // Return scattering angle sampled in cms
731
732
733 G4double
734 G4DiffuseElastic::SampleThetaCMS(const G4Parti
735 G4doubl
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(G4Dif
752
753 // result = integral.Legendre10(this,&G4Diff
754 norm = integral.Legendre96(this,&G4DiffuseEl
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,&G4Diffu
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 r
784 //////////////////////////////////////////////
785 //
786 // Return inv momentum transfer -t > 0 from in
787
788 G4double G4DiffuseElastic::SampleInvariantT( c
789
790 {
791 fParticle = aParticle;
792 G4double m1 = fParticle->GetPDGMass(), t;
793 G4double totElab = std::sqrt(m1*m1+p*p);
794 G4double mass2 = G4NucleiProperties::GetNucl
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::
815 return t;
816 }
817 }
818
819 t = SampleTableT( aParticle, momentumCMS, G
820
821 return t;
822 }
823
824 //////////////////////////////////////////////
825
826 G4double G4DiffuseElastic::NeutronTuniform(G4i
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 in
839
840 G4double G4DiffuseElastic::SampleTableT( const
841
842 {
843 G4double alpha = SampleTableThetaCMS( aParti
844 G4double t = 2*p*p*( 1 - std::cos(std::s
845 // G4double t = p*p*alpha; /
846 return t;
847 }
848
849 //////////////////////////////////////////////
850 //
851 // Return scattering angle2 sampled in cms acc
852
853
854 G4double
855 G4DiffuseElastic::SampleTableThetaCMS(const G4
856 G4doubl
857 {
858 std::size_t iElement;
859 G4int iMomentum, iAngle;
860 G4double randAngle, position, theta1, theta2
861 G4double m1 = particle->GetPDGMass();
862
863 for(iElement = 0; iElement < fElementNumberV
864 {
865 if( std::fabs(Z - fElementNumberVector[iEl
866 }
867 if ( iElement == fElementNumberVector.size()
868 {
869 InitialiseOnFly(Z,A); // table preparation
870
871 // iElement--;
872
873 // G4cout << "G4DiffuseElastic: Element wi
874 // << " is not found, return zero angle" <
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
882
883 for( iMomentum = 0; iMomentum < fEnergyBin;
884 {
885 if( kinE < fEnergyVector->GetLowEdgeEnergy
886 }
887 if ( iMomentum >= fEnergyBin ) iMomentum = f
888 if ( iMomentum < 0 ) iMomentum = 0
889
890 // G4cout<<"iMomentum = "<<iMomentum<<G4endl
891
892 if (iMomentum == fEnergyBin -1 || iMomentum
893 {
894 position = (*(*fAngleTable)(iMomentum))(fA
895
896 // G4cout<<"position = "<<position<<G4endl
897
898 for(iAngle = 0; iAngle < fAngleBin-1; iAng
899 {
900 if( position > (*(*fAngleTable)(iMomentu
901 }
902 if (iAngle >= fAngleBin-1) iAngle = fAngle
903
904 // G4cout<<"iAngle = "<<iAngle<<G4endl;
905
906 randAngle = GetScatteringAngle(iMomentum,
907
908 // G4cout<<"randAngle = "<<randAngle<<G4en
909 }
910 else // kinE inside between energy table ed
911 {
912 // position = (*(*fAngleTable)(iMomentum))
913 position = (*(*fAngleTable)(iMomentum))(0)
914
915 // G4cout<<"position = "<<position<<G4endl
916
917 for(iAngle = 0; iAngle < fAngleBin-1; iAng
918 {
919 // if( position < (*(*fAngleTable)(iMome
920 if( position > (*(*fAngleTable)(iMomentu
921 }
922 if (iAngle >= fAngleBin-1) iAngle = fAngle
923
924 // G4cout<<"iAngle = "<<iAngle<<G4endl;
925
926 theta2 = GetScatteringAngle(iMomentum, iA
927
928 // G4cout<<"theta2 = "<<theta2<<G4endl;
929 E2 = fEnergyVector->GetLowEdgeEnergy(iMome
930
931 // G4cout<<"E2 = "<<E2<<G4endl;
932
933 iMomentum--;
934
935 // position = (*(*fAngleTable)(iMomentum))
936
937 // G4cout<<"position = "<<position<<G4endl
938
939 for(iAngle = 0; iAngle < fAngleBin-1; iAng
940 {
941 // if( position < (*(*fAngleTable)(iMome
942 if( position > (*(*fAngleTable)(iMomentu
943 }
944 if (iAngle >= fAngleBin-1) iAngle = fAngle
945
946 theta1 = GetScatteringAngle(iMomentum, iA
947
948 // G4cout<<"theta1 = "<<theta1<<G4endl;
949
950 E1 = fEnergyVector->GetLowEdgeEnergy(iMome
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<<G4en
962 }
963 // G4double angle = randAngle;
964 // if (randAngle > 0.) randAngle /= 2*pi*std
965
966 if(randAngle < 0.) randAngle = 0.;
967
968 return randAngle;
969 }
970
971 //////////////////////////////////////////////
972 //
973 // Initialisation for given particle on fly us
974
975 void G4DiffuseElastic::InitialiseOnFly(G4doubl
976 {
977 fAtomicNumber = Z; // atomic number
978 fAtomicWeight = G4NistManager::Instance()->
979
980 fNuclearRadius = CalculateNuclearRad(fAtomic
981
982 if( verboseLevel > 0 )
983 {
984 G4cout<<"G4DiffuseElastic::InitialiseOnFly
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
999 // For the moment in lab system.
1000
1001 void G4DiffuseElastic::BuildAngleTable()
1002 {
1003 G4int i, j;
1004 G4double partMom, kinE, a = 0., z = fPartic
1005 G4double alpha1, alpha2, alphaMax, alphaCou
1006
1007 G4Integrator<G4DiffuseElastic,G4double(G4Di
1008
1009 fAngleTable = new G4PhysicsTable( fEnergyBi
1010
1011 for( i = 0; i < fEnergyBin; i++)
1012 {
1013 kinE = fEnergyVector->GetLowEdgeEn
1014 partMom = std::sqrt( kinE*(kinE + 2*m
1015
1016 fWaveVector = partMom/hbarc;
1017
1018 G4double kR = fWaveVector*fNuclearRad
1019 G4double kR2 = kR*kR;
1020 G4double kRmax = 18.6; // 10.6; 10.6, 18
1021 G4double kRcoul = 1.9; // 1.2; 1.4, 2.5;
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.; //
1029 // if ( alphaMax > 4. || alphaMax < 1. )
1030
1031 // if ( alphaMax > 4. || alphaMax < 1. )
1032
1033 // G4cout<<"alphaMax = "<<alphaMax<<", ";
1034
1035 if ( alphaMax >= CLHEP::pi*CLHEP::pi ) al
1036
1037 alphaCoulomb = kRcoul*kRcoul/kR2;
1038
1039 if( z )
1040 {
1041 a = partMom/m1; // be
1042 fBeta = a/std::sqrt(1+a*a);
1043 fZommerfeld = CalculateZommerfeld( fBet
1044 fAm = CalculateAm( partMom, fZo
1045 }
1046 G4PhysicsFreeVector* angleVector = new G4
1047
1048 // G4PhysicsLogVector* angleBins = new G
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
1061 // alpha2 = angleBins->GetLowEdgeEnergy
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 )
1071 if( ( alpha1 < alphaCoulomb ) && z ) fA
1072
1073 delta = integral.Legendre10(this, &G4Di
1074 // delta = integral.Legendre96(this, &G
1075
1076 sum += delta;
1077
1078 angleVector->PutValue( j-1 , alpha1, su
1079 // G4cout<<"j-1 = "<<j-1<<"; alpha
1080 }
1081 fAngleTable->insertAt(i, angleVector);
1082
1083 // delete[] angleVector;
1084 // delete[] angleBins;
1085 }
1086 return;
1087 }
1088
1089 /////////////////////////////////////////////
1090 //
1091 //
1092
1093 G4double
1094 G4DiffuseElastic:: GetScatteringAngle( G4int
1095 {
1096 G4double x1, x2, y1, y2, randAngle;
1097
1098 if( iAngle == 0 )
1099 {
1100 randAngle = (*fAngleTable)(iMomentum)->Ge
1101 // iAngle++;
1102 }
1103 else
1104 {
1105 if ( iAngle >= G4int((*fAngleTable)(iMome
1106 {
1107 iAngle = G4int((*fAngleTable)(iMomentum
1108 }
1109 y1 = (*(*fAngleTable)(iMomentum))(iAngle-
1110 y2 = (*(*fAngleTable)(iMomentum))(iAngle)
1111
1112 x1 = (*fAngleTable)(iMomentum)->GetLowEdg
1113 x2 = (*fAngleTable)(iMomentum)->GetLowEdg
1114
1115 if ( x1 == x2 ) randAngle = x2;
1116 else
1117 {
1118 if ( y1 == y2 ) randAngle = x1 + ( x2 -
1119 else
1120 {
1121 randAngle = x1 + ( position - y1 )*(
1122 }
1123 }
1124 }
1125 return randAngle;
1126 }
1127
1128
1129
1130 /////////////////////////////////////////////
1131 //
1132 // Return scattering angle sampled in lab sys
1133
1134
1135
1136 G4double
1137 G4DiffuseElastic::SampleThetaLab( const G4Had
1138 G4dou
1139 {
1140 const G4ParticleDefinition* theParticle = a
1141 G4double m1 = theParticle->GetPDGMass();
1142 G4double plab = aParticle->GetTotalMomentum
1143 G4LorentzVector lv1 = aParticle->Get4Moment
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
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
1202 }
1203 G4ThreeVector v1(sint*std::cos(phi),sint*st
1204 v1 *= ptot;
1205 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s
1206
1207 nlv1.boost(bst);
1208
1209 G4ThreeVector np1 = nlv1.vect();
1210
1211 // G4double theta = std::acos( np1.z()/np
1212
1213 G4double theta = np1.theta();
1214
1215 return theta;
1216 }
1217
1218 /////////////////////////////////////////////
1219 //
1220 // Return scattering angle in lab system (tar
1221
1222
1223
1224 G4double
1225 G4DiffuseElastic::ThetaCMStoThetaLab( const G
1226 G4dou
1227 {
1228 const G4ParticleDefinition* theParticle = a
1229 G4double m1 = theParticle->GetPDGMass();
1230 // G4double plab = aParticle->GetTotalMomen
1231 G4LorentzVector lv1 = aParticle->Get4Moment
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::
1264 }
1265 G4ThreeVector v1(sint*std::cos(phi),sint*st
1266 v1 *= ptot;
1267 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s
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 (tar
1281
1282
1283
1284 G4double
1285 G4DiffuseElastic::ThetaLabToThetaCMS( const G
1286 G4dou
1287 {
1288 const G4ParticleDefinition* theParticle = a
1289 G4double m1 = theParticle->GetPDGMass();
1290 G4double plab = aParticle->GetTotalMomentum
1291 G4LorentzVector lv1 = aParticle->Get4Moment
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::
1324 }
1325 G4ThreeVector v1(sint*std::cos(phi),sint*st
1326 v1 *= plab;
1327 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s
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
1342 // For the moment in lab system.
1343
1344 void G4DiffuseElastic::TestAngleTable(const G
1345 G
1346 {
1347 fAtomicNumber = Z; // atomic number
1348 fAtomicWeight = A; // number of nucleo
1349 fNuclearRadius = CalculateNuclearRad(fAtomi
1350
1351
1352
1353 G4cout<<"G4DiffuseElastic::TestAngleTable()
1354 <<Z<<"; and A = "<<A<<G4endl;
1355
1356 fElementNumberVector.push_back(fAtomicNumbe
1357
1358
1359
1360
1361 G4int i=0, j;
1362 G4double a = 0., z = theParticle->GetPDGCha
1363 G4double alpha1=0., alpha2=0., alphaMax=0.,
1364 G4double deltaL10 = 0., deltaL96 = 0., delt
1365 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.
1366 G4double epsilon = 0.001;
1367
1368 G4Integrator<G4DiffuseElastic,G4double(G4Di
1369
1370 fAngleTable = new G4PhysicsTable(fEnergyBin
1371
1372 fWaveVector = partMom/hbarc;
1373
1374 G4double kR = fWaveVector*fNuclearRadiu
1375 G4double kR2 = kR*kR;
1376 G4double kRmax = 10.6; // 10.6, 18, 10.174
1377 G4double kRcoul = 1.2; // 1.4, 2.5; // on t
1378
1379 alphaMax = kRmax*kRmax/kR2;
1380
1381 if (alphaMax > 4.) alphaMax = 4.; // vmg05
1382
1383 alphaCoulomb = kRcoul*kRcoul/kR2;
1384
1385 if( z )
1386 {
1387 a = partMom/m1; // beta*gamma
1388 fBeta = a/std::sqrt(1+a*a);
1389 fZommerfeld = CalculateZommerfeld( fBet
1390 fAm = CalculateAm( partMom, fZo
1391 }
1392 G4PhysicsFreeVector* angleVector = new G4Ph
1393
1394 // G4PhysicsLogVector* angleBins = new G4P
1395
1396
1397 fAddCoulomb = false;
1398
1399 for(j = 1; j < fAngleBin; j++)
1400 {
1401 // alpha1 = angleBins->GetLowEdgeEnergy
1402 // alpha2 = angleBins->GetLowEdgeEnergy
1403
1404 alpha1 = alphaMax*(j-1)/fAngleBin;
1405 alpha2 = alphaMax*( j )/fAngleBin;
1406
1407 if( ( alpha2 > alphaCoulomb ) && z ) fAdd
1408
1409 deltaL10 = integral.Legendre10(this, &G4D
1410 deltaL96 = integral.Legendre96(this, &G4D
1411 deltaAG = integral.AdaptiveGauss(this, &
1412 alpha1
1413
1414 // G4cout<<alpha1<<"\t"<<std::sqrt(alph
1415 // <<deltaL10<<"\t"<<deltaL96<<"\t"
1416
1417 sumL10 += deltaL10;
1418 sumL96 += deltaL96;
1419 sumAG += deltaAG;
1420
1421 G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/d
1422 <<sumL10<<"\t"<<sumL96<<"\t"<<sum
1423
1424 angleVector->PutValue( j-1 , alpha1, sumL
1425 }
1426 fAngleTable->insertAt(i,angleVector);
1427 fAngleBank.push_back(fAngleTable);
1428
1429 /*
1430 // Integral over all angle range - Bad accu
1431
1432 sumL10 = integral.Legendre10(this, &G4Diffu
1433 sumL96 = integral.Legendre96(this, &G4Diffu
1434 sumAG = integral.AdaptiveGauss(this, &G4Di
1435 0., al
1436 G4cout<<G4endl;
1437 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/deg
1438 <<sumL10<<"\t"<<sumL96<<"\t"<<sum
1439 */
1440 return;
1441 }
1442
1443 //
1444 //
1445 /////////////////////////////////////////////
1446