source: trunk/source/processes/hadronic/models/coherent_elastic/src/G4DiffuseElastic.cc@ 1200

Last change on this file since 1200 was 1196, checked in by garnier, 16 years ago

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