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

Last change on this file since 1245 was 819, checked in by garnier, 16 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.