source: trunk/source/processes/hadronic/models/low_energy/src/G4LElastic.cc@ 1201

Last change on this file since 1201 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 16.3 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//
27//
28// Physics model class G4LElastic
29//
30//
31// G4 Model: Low-energy Elastic scattering
32// F.W. Jones, TRIUMF, 04-JUN-96
33//
34// use -scheme for elastic scattering: HPW, 20th June 1997
35// most of the code comes from the old Low-energy Elastic class
36//
37// 25-JUN-98 FWJ: replaced missing Initialize for ParticleChange.
38// 14-DEC-05 V.Ivanchenko: restore 1.19 version (7.0)
39// 23-JAN-07 V.Ivanchenko: add protection inside sqrt
40//
41
42#include "globals.hh"
43#include "G4LElastic.hh"
44#include "Randomize.hh"
45#include "G4ParticleTable.hh"
46#include "G4IonTable.hh"
47
48
49G4HadFinalState*
50G4LElastic::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
51{
52 if(getenv("debug_LElastic")) verboseLevel = 5;
53 theParticleChange.Clear();
54 const G4HadProjectile* aParticle = &aTrack;
55 G4double atno2 = targetNucleus.GetN();
56 G4double zTarget = targetNucleus.GetZ();
57 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
58 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
59
60// Elastic scattering off Hydrogen
61
62 G4DynamicParticle* aSecondary = 0;
63 if (atno2 < 1.5) {
64 const G4ParticleDefinition* aParticleType = aParticle->GetDefinition();
65 if (aParticleType == G4PionPlus::PionPlus())
66 aSecondary = LightMedia.PionPlusExchange(aParticle, targetNucleus);
67 else if (aParticleType == G4PionMinus::PionMinus())
68 aSecondary = LightMedia.PionMinusExchange(aParticle, targetNucleus);
69 else if (aParticleType == G4KaonPlus::KaonPlus())
70 aSecondary = LightMedia.KaonPlusExchange(aParticle, targetNucleus);
71 else if (aParticleType == G4KaonZeroShort::KaonZeroShort())
72 aSecondary = LightMedia.KaonZeroShortExchange(aParticle,targetNucleus);
73 else if (aParticleType == G4KaonZeroLong::KaonZeroLong())
74 aSecondary = LightMedia.KaonZeroLongExchange(aParticle, targetNucleus);
75 else if (aParticleType == G4KaonMinus::KaonMinus())
76 aSecondary = LightMedia.KaonMinusExchange(aParticle, targetNucleus);
77 else if (aParticleType == G4Proton::Proton())
78 aSecondary = LightMedia.ProtonExchange(aParticle, targetNucleus);
79 else if (aParticleType == G4AntiProton::AntiProton())
80 aSecondary = LightMedia.AntiProtonExchange(aParticle, targetNucleus);
81 else if (aParticleType == G4Neutron::Neutron())
82 aSecondary = LightMedia.NeutronExchange(aParticle, targetNucleus);
83 else if (aParticleType == G4AntiNeutron::AntiNeutron())
84 aSecondary = LightMedia.AntiNeutronExchange(aParticle, targetNucleus);
85 else if (aParticleType == G4Lambda::Lambda())
86 aSecondary = LightMedia.LambdaExchange(aParticle, targetNucleus);
87 else if (aParticleType == G4AntiLambda::AntiLambda())
88 aSecondary = LightMedia.AntiLambdaExchange(aParticle, targetNucleus);
89 else if (aParticleType == G4SigmaPlus::SigmaPlus())
90 aSecondary = LightMedia.SigmaPlusExchange(aParticle, targetNucleus);
91 else if (aParticleType == G4SigmaMinus::SigmaMinus())
92 aSecondary = LightMedia.SigmaMinusExchange(aParticle, targetNucleus);
93 else if (aParticleType == G4AntiSigmaPlus::AntiSigmaPlus())
94 aSecondary = LightMedia.AntiSigmaPlusExchange(aParticle,targetNucleus);
95 else if (aParticleType == G4AntiSigmaMinus::AntiSigmaMinus())
96 aSecondary= LightMedia.AntiSigmaMinusExchange(aParticle,targetNucleus);
97 else if (aParticleType == G4XiZero::XiZero())
98 aSecondary = LightMedia.XiZeroExchange(aParticle, targetNucleus);
99 else if (aParticleType == G4XiMinus::XiMinus())
100 aSecondary = LightMedia.XiMinusExchange(aParticle, targetNucleus);
101 else if (aParticleType == G4AntiXiZero::AntiXiZero())
102 aSecondary = LightMedia.AntiXiZeroExchange(aParticle, targetNucleus);
103 else if (aParticleType == G4AntiXiMinus::AntiXiMinus())
104 aSecondary = LightMedia.AntiXiMinusExchange(aParticle, targetNucleus);
105 else if (aParticleType == G4OmegaMinus::OmegaMinus())
106 aSecondary = LightMedia.OmegaMinusExchange(aParticle, targetNucleus);
107 else if (aParticleType == G4AntiOmegaMinus::AntiOmegaMinus())
108 aSecondary= LightMedia.AntiOmegaMinusExchange(aParticle,targetNucleus);
109 else if (aParticleType == G4KaonPlus::KaonPlus())
110 aSecondary = LightMedia.KaonPlusExchange(aParticle, targetNucleus);
111 }
112
113// Has a charge or strangeness exchange occurred?
114 if (aSecondary) {
115 aSecondary->SetMomentum(aParticle->Get4Momentum().vect());
116 theParticleChange.SetStatusChange(stopAndKill);
117 theParticleChange.AddSecondary(aSecondary);
118 }
119 // G4cout << "Entering elastic scattering 1"<<G4endl;
120
121 G4double p = aParticle->GetTotalMomentum()/GeV;
122 if (verboseLevel > 1)
123 G4cout << "G4LElastic::DoIt: Incident particle p=" << p << " GeV" << G4endl;
124
125 if (p < 0.01) return &theParticleChange;
126
127 // G4cout << "Entering elastic scattering 2"<<G4endl;
128// Compute the direction of elastic scattering.
129// It is planned to replace this code with a method based on
130// parameterized functions and a Monte Carlo method to invert the CDF.
131
132 G4double ran = G4UniformRand();
133 G4double aa, bb, cc, dd, rr;
134 if (atno2 <= 62.) {
135 aa = std::pow(atno2, 1.63);
136 bb = 14.5*std::pow(atno2, 0.66);
137 cc = 1.4*std::pow(atno2, 0.33);
138 dd = 10.;
139 }
140 else {
141 aa = std::pow(atno2, 1.33);
142 bb = 60.*std::pow(atno2, 0.33);
143 cc = 0.4*std::pow(atno2, 0.40);
144 dd = 10.;
145 }
146 aa = aa/bb;
147 cc = cc/dd;
148 rr = (aa + cc)*ran;
149 if (verboseLevel > 1) {
150 G4cout << "DoIt: aa,bb,cc,dd,rr" << G4endl;
151 G4cout << aa << " " << bb << " " << cc << " " << dd << " " << rr << G4endl;
152 }
153 G4double t1 = -std::log(ran)/bb;
154 G4double t2 = -std::log(ran)/dd;
155 if (verboseLevel > 1) {
156 G4cout << "t1,Fctcos " << t1 << " " << Fctcos(t1, aa, bb, cc, dd, rr) <<
157 G4endl;
158 G4cout << "t2,Fctcos " << t2 << " " << Fctcos(t2, aa, bb, cc, dd, rr) <<
159 G4endl;
160 }
161 G4double eps = 0.001;
162 G4int ind1 = 10;
163 G4double t;
164 G4int ier1;
165 ier1 = Rtmi(&t, t1, t2, eps, ind1,
166 aa, bb, cc, dd, rr);
167 if (verboseLevel > 1) {
168 G4cout << "From Rtmi, ier1=" << ier1 << G4endl;
169 G4cout << "t, Fctcos " << t << " " << Fctcos(t, aa, bb, cc, dd, rr) <<
170 G4endl;
171 }
172 if (ier1 != 0) t = 0.25*(3.*t1 + t2);
173 if (verboseLevel > 1) {
174 G4cout << "t, Fctcos " << t << " " << Fctcos(t, aa, bb, cc, dd, rr) <<
175 G4endl;
176 }
177 G4double phi = G4UniformRand()*twopi;
178 rr = 0.5*t/(p*p);
179 if (rr > 1.) rr = 0.;
180 if (verboseLevel > 1)
181 G4cout << "rr=" << rr << G4endl;
182 G4double cost = 1. - rr;
183 G4double sint = std::sqrt(std::max(rr*(2. - rr), 0.));
184 if (sint == 0.) return &theParticleChange;
185 // G4cout << "Entering elastic scattering 3"<<G4endl;
186 if (verboseLevel > 1) G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
187// Scattered particle referred to axis of incident particle
188 G4double m1=aParticle->GetDefinition()->GetPDGMass();
189 G4int Z=static_cast<G4int>(zTarget+.5);
190 G4int A=static_cast<G4int>(atno2);
191 if(G4UniformRand()<atno2-A) A++;
192 //G4cout << " ion info "<<atno2 << " "<<A<<" "<<Z<<" "<<zTarget<<G4endl;
193 G4double m2=G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z)->GetPDGMass();
194// non relativistic approximation
195 G4double a=1+m2/m1;
196 G4double b=-2.*p*cost;
197 G4double c=p*p*(1-m2/m1);
198 G4double p1 = (-b+std::sqrt(std::max(0.0,b*b-4.*a*c)))/(2.*a);
199 G4double px = p1*sint*std::sin(phi);
200 G4double py = p1*sint*std::cos(phi);
201 G4double pz = p1*cost;
202
203// relativistic calculation
204 G4double etot = std::sqrt(m1*m1+p*p)+m2;
205 a = etot*etot-p*p*cost*cost;
206 b = 2*p*p*(m1*cost*cost-etot);
207 c = p*p*p*p*sint*sint;
208
209 G4double de = (-b-std::sqrt(std::max(0.0,b*b-4.*a*c)))/(2.*a);
210 G4double e1 = std::sqrt(p*p+m1*m1)-de;
211 G4double p12=e1*e1-m1*m1;
212 p1 = std::sqrt(std::max(1.*eV*eV,p12));
213 px = p1*sint*std::sin(phi);
214 py = p1*sint*std::cos(phi);
215 pz = p1*cost;
216
217 if (verboseLevel > 1)
218 {
219 G4cout << "Relevant test "<<p<<" "<<p1<<" "<<cost<<" "<<de<<G4endl;
220 G4cout << "p1/p = "<<p1/p<<" "<<m1<<" "<<m2<<" "<<a<<" "<<b<<" "<<c<<G4endl;
221 G4cout << "rest = "<< b*b<<" "<<4.*a*c<<" "<<G4endl;
222 G4cout << "make p1 = "<< p12<<" "<<e1*e1<<" "<<m1*m1<<" "<<G4endl;
223 }
224// Incident particle
225 G4double pxinc = p*aParticle->Get4Momentum().vect().unit().x();
226 G4double pyinc = p*aParticle->Get4Momentum().vect().unit().y();
227 G4double pzinc = p*aParticle->Get4Momentum().vect().unit().z();
228 if (verboseLevel > 1) {
229 G4cout << "NOM SCAT " << px << " " << py << " " << pz << G4endl;
230 G4cout << "INCIDENT " << pxinc << " " << pyinc << " " << pzinc << G4endl;
231 }
232
233// Transform scattered particle to reflect direction of incident particle
234 G4double pxnew, pynew, pznew;
235 Defs1(p, px, py, pz, pxinc, pyinc, pzinc, &pxnew, &pynew, &pznew);
236// Normalize:
237 G4double pxre=pxinc-pxnew;
238 G4double pyre=pyinc-pynew;
239 G4double pzre=pzinc-pznew;
240 G4ThreeVector it0(pxnew*GeV, pynew*GeV, pznew*GeV);
241 if(p1>0)
242 {
243 pxnew = pxnew/p1;
244 pynew = pynew/p1;
245 pznew = pznew/p1;
246 }
247 else
248 {
249 //G4double pphi = 2*pi*G4UniformRand();
250 //G4double ccth = 2*G4UniformRand()-1;
251 pxnew = 0;//std::sin(std::acos(ccth))*std::sin(pphi);
252 pynew = 0;//std::sin(std::acos(ccth))*std::cos(phi);
253 pznew = 1;//ccth;
254 }
255 if (verboseLevel > 1) {
256 G4cout << "DoIt: returning new momentum vector" << G4endl;
257 G4cout << "DoIt: "<<pxinc << " " << pyinc << " " << pzinc <<" "<<p<< G4endl;
258 G4cout << "DoIt: "<<pxnew << " " << pynew << " " << pznew <<" "<<p<< G4endl;
259 }
260
261 if (aSecondary)
262 {
263 aSecondary->SetMomentumDirection(pxnew, pynew, pznew);
264 }
265 else
266 {
267 try
268 {
269 theParticleChange.SetMomentumChange(pxnew, pynew, pznew);
270 theParticleChange.SetEnergyChange(std::sqrt(m1*m1+it0.mag2())-m1);
271 }
272 catch(G4HadronicException)
273 {
274 std::cerr << "GHADException originating from components of G4LElastic"<<std::cout;
275 throw;
276 }
277 G4ParticleDefinition * theDef = G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
278 G4ThreeVector it(pxre*GeV, pyre*GeV, pzre*GeV);
279 G4DynamicParticle * aSec =
280 new G4DynamicParticle(theDef, it.unit(), it.mag2()/(2.*m2));
281 theParticleChange.AddSecondary(aSec);
282 // G4cout << "Final check ###### "<<p<<" "<<it.mag()<<" "<<p1<<G4endl;
283 }
284 return &theParticleChange;
285}
286
287
288// The following is a "translation" of a root-finding routine
289// from GEANT3.21/GHEISHA. Some of the labelled block structure has
290// been retained for clarity. This routine will not be needed after
291// the planned revisions to DoIt().
292
293G4int
294G4LElastic::Rtmi(G4double* x, G4double xli, G4double xri, G4double eps,
295 G4int iend,
296 G4double aa, G4double bb, G4double cc, G4double dd,
297 G4double rr)
298{
299 G4int ier = 0;
300 G4double xl = xli;
301 G4double xr = xri;
302 *x = xl;
303 G4double tol = *x;
304 G4double f = Fctcos(tol, aa, bb, cc, dd, rr);
305 if (f == 0.) return ier;
306 G4double fl, fr;
307 fl = f;
308 *x = xr;
309 tol = *x;
310 f = Fctcos(tol, aa, bb, cc, dd, rr);
311 if (f == 0.) return ier;
312 fr = f;
313
314// Error return in case of wrong input data
315 if (fl*fr >= 0.) {
316 ier = 2;
317 return ier;
318 }
319
320// Basic assumption fl*fr less than 0 is satisfied.
321// Generate tolerance for function values.
322 G4int i = 0;
323 G4double tolf = 100.*eps;
324
325// Start iteration loop
326label4:
327 i++;
328
329// Start bisection loop
330 for (G4int k = 1; k <= iend; k++) {
331 *x = 0.5*(xl + xr);
332 tol = *x;
333 f = Fctcos(tol, aa, bb, cc, dd, rr);
334 if (f == 0.) return 0;
335 if (f*fr < 0.) { // Interchange xl and xr in order to get the
336 tol = xl; // same Sign in f and fr
337 xl = xr;
338 xr = tol;
339 tol = fl;
340 fl = fr;
341 fr = tol;
342 }
343 tol = f - fl;
344 G4double a = f*tol;
345 a = a + a;
346 if (a < fr*(fr - fl) && i <= iend) goto label17;
347 xr = *x;
348 fr = f;
349
350// Test on satisfactory accuracy in bisection loop
351 tol = eps;
352 a = std::abs(xr);
353 if (a > 1.) tol = tol*a;
354 if (std::abs(xr - xl) <= tol && std::abs(fr - fl) <= tolf) goto label14;
355 }
356// End of bisection loop
357
358// No convergence after iend iteration steps followed by iend
359// successive steps of bisection or steadily increasing function
360// values at right bounds. Error return.
361 ier = 1;
362
363label14:
364 if (std::abs(fr) > std::abs(fl)) {
365 *x = xl;
366 f = fl;
367 }
368 return ier;
369
370// Computation of iterated x-value by inverse parabolic interp
371label17:
372 G4double a = fr - f;
373 G4double dx = (*x - xl)*fl*(1. + f*(a - tol)/(a*(fr - fl)))/tol;
374 G4double xm = *x;
375 G4double fm = f;
376 *x = xl - dx;
377 tol = *x;
378 f = Fctcos(tol, aa, bb, cc, dd, rr);
379 if (f == 0.) return ier;
380
381// Test on satisfactory accuracy in iteration loop
382 tol = eps;
383 a = std::abs(*x);
384 if (a > 1) tol = tol*a;
385 if (std::abs(dx) <= tol && std::abs(f) <= tolf) return ier;
386
387// Preparation of next bisection loop
388 if (f*fl < 0.) {
389 xr = *x;
390 fr = f;
391 }
392 else {
393 xl = *x;
394 fl = f;
395 xr = xm;
396 fr = fm;
397 }
398 goto label4;
399}
400
401
402// Test function for root-finder
403
404G4double
405G4LElastic::Fctcos(G4double t,
406 G4double aa, G4double bb, G4double cc, G4double dd,
407 G4double rr)
408{
409 const G4double expxl = -82.;
410 const G4double expxu = 82.;
411
412 G4double test1 = -bb*t;
413 if (test1 > expxu) test1 = expxu;
414 if (test1 < expxl) test1 = expxl;
415
416 G4double test2 = -dd*t;
417 if (test2 > expxu) test2 = expxu;
418 if (test2 < expxl) test2 = expxl;
419
420 return aa*std::exp(test1) + cc*std::exp(test2) - rr;
421}
422
423
424void
425G4LElastic::Defs1(G4double p, G4double px, G4double py, G4double pz,
426 G4double pxinc, G4double pyinc, G4double pzinc,
427 G4double* pxnew, G4double* pynew, G4double* pznew)
428{
429// Transform scattered particle to reflect direction of incident particle
430 G4double pt2 = pxinc*pxinc + pyinc*pyinc;
431 if (pt2 > 0.) {
432 G4double cost = pzinc/p;
433 G4double sint1 = std::sqrt(std::abs((1. - cost )*(1.+cost)));
434 G4double sint2 = std::sqrt(pt2)/p;
435 G4double sint = 0.5*(sint1 + sint2);
436 G4double ph = pi*0.5;
437 if (pyinc < 0.) ph = pi*1.5;
438 if (std::abs(pxinc) > 1.e-6) ph = std::atan2(pyinc, pxinc);
439 G4double cosp = std::cos(ph);
440 G4double sinp = std::sin(ph);
441 if (verboseLevel > 1) {
442 G4cout << "cost sint " << cost << " " << sint << G4endl;
443 G4cout << "cosp sinp " << cosp << " " << sinp << G4endl;
444 }
445 *pxnew = cost*cosp*px - sinp*py + sint*cosp*pz;
446 *pynew = cost*sinp*px + cosp*py + sint*sinp*pz;
447 *pznew = -sint*px +cost*pz;
448 }
449 else {
450 G4double cost=pzinc/p;
451 *pxnew = cost*px;
452 *pynew = py;
453 *pznew = cost*pz;
454 }
455}
Note: See TracBrowser for help on using the repository browser.