Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NuclNuclDiffuseElastic.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4NuclNuclDiffuseElastic.cc 98826 2016-08-12 12:36:07Z gcosmo $
27 //
28 //
29 // Physics model class G4NuclNuclDiffuseElastic
30 //
31 //
32 // G4 Model: optical diffuse elastic scattering with 4-momentum balance
33 //
34 // 24-May-07 V. Grichine
35 //
36 
38 #include "G4ParticleTable.hh"
39 #include "G4ParticleDefinition.hh"
40 #include "G4IonTable.hh"
41 #include "G4NucleiProperties.hh"
42 
43 #include "Randomize.hh"
44 #include "G4Integrator.hh"
45 #include "globals.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 
49 #include "G4Proton.hh"
50 #include "G4Neutron.hh"
51 #include "G4Deuteron.hh"
52 #include "G4Alpha.hh"
53 #include "G4PionPlus.hh"
54 #include "G4PionMinus.hh"
55 
56 #include "G4Element.hh"
57 #include "G4ElementTable.hh"
58 #include "G4NistManager.hh"
59 #include "G4PhysicsTable.hh"
60 #include "G4PhysicsLogVector.hh"
61 #include "G4PhysicsFreeVector.hh"
62 
64 //
65 // Test Constructor. Just to check xsc
66 
67 
69  : G4HadronElastic("NNDiffuseElastic"), fParticle(0)
70 {
71  SetMinEnergy( 50*MeV );
72  SetMaxEnergy( 1.*TeV );
73  verboseLevel = 0;
74  lowEnergyRecoilLimit = 100.*keV;
75  lowEnergyLimitQ = 0.0*GeV;
76  lowEnergyLimitHE = 0.0*GeV;
77  lowestEnergyLimit= 0.0*keV;
78  plabLowLimit = 20.0*MeV;
79 
80  theProton = G4Proton::Proton();
81  theNeutron = G4Neutron::Neutron();
82  theDeuteron = G4Deuteron::Deuteron();
83  theAlpha = G4Alpha::Alpha();
84  thePionPlus = G4PionPlus::PionPlus();
85  thePionMinus= G4PionMinus::PionMinus();
86 
87  fEnergyBin = 200;
88  fAngleBin = 200;
89 
90  fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
91  fAngleTable = 0;
92 
93  fParticle = 0;
94  fWaveVector = 0.;
95  fAtomicWeight = 0.;
96  fAtomicNumber = 0.;
97  fNuclearRadius = 0.;
98  fBeta = 0.;
99  fZommerfeld = 0.;
100  fAm = 0.;
101  fAddCoulomb = false;
102  // Ranges of angle table relative to current Rutherford (Coulomb grazing) angle
103 
104  // Empirical parameters
105 
106  fCofAlphaMax = 1.5;
107  fCofAlphaCoulomb = 0.5;
108 
109  fProfileDelta = 1.;
110  fProfileAlpha = 0.5;
111 
112  fCofLambda = 1.0;
113  fCofDelta = 0.04;
114  fCofAlpha = 0.095;
115 
116  fNuclearRadius1 = fNuclearRadius2 = fNuclearRadiusSquare
117  = fRutherfordRatio = fCoulombPhase0 = fHalfRutThetaTg = fHalfRutThetaTg2
118  = fRutherfordTheta = fProfileLambda = fCofPhase = fCofFar = fCofAlphaMax
119  = fCofAlphaCoulomb = fSumSigma = fEtaRatio = fReZ = 0.0;
120  fMaxL = 0;
121 
122  fNuclearRadiusCof = 1.0;
123 }
124 
126 //
127 // Destructor
128 
130 {
131  if ( fEnergyVector ) {
132  delete fEnergyVector;
133  fEnergyVector = 0;
134  }
135 
136  for ( std::vector<G4PhysicsTable*>::iterator it = fAngleBank.begin();
137  it != fAngleBank.end(); ++it ) {
138  if ( (*it) ) (*it)->clearAndDestroy();
139  delete *it;
140  *it = 0;
141  }
142  fAngleTable = 0;
143 }
144 
146 //
147 // Initialisation for given particle using element table of application
148 
150 {
151 
152  // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
153 
154  const G4ElementTable* theElementTable = G4Element::GetElementTable();
155  size_t jEl, numOfEl = G4Element::GetNumberOfElements();
156 
157  // projectile radius
158 
159  G4double A1 = G4double( fParticle->GetBaryonNumber() );
160  G4double R1 = CalculateNuclearRad(A1);
161 
162  for(jEl = 0 ; jEl < numOfEl; ++jEl) // application element loop
163  {
164  fAtomicNumber = (*theElementTable)[jEl]->GetZ(); // atomic number
165  fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( fAtomicNumber ) );
166 
167  fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
168  fNuclearRadius += R1;
169 
170  if(verboseLevel > 0)
171  {
172  G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element: "
173  <<(*theElementTable)[jEl]->GetName()<<G4endl;
174  }
175  fElementNumberVector.push_back(fAtomicNumber);
176  fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
177 
178  BuildAngleTable();
179  fAngleBank.push_back(fAngleTable);
180  }
181 }
182 
183 
185 //
186 // return differential elastic cross section d(sigma)/d(omega)
187 
188 G4double
190  G4double theta,
191  G4double momentum,
192  G4double A )
193 {
194  fParticle = particle;
195  fWaveVector = momentum/hbarc;
196  fAtomicWeight = A;
197  fAddCoulomb = false;
198  fNuclearRadius = CalculateNuclearRad(A);
199 
200  G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta);
201 
202  return sigma;
203 }
204 
206 //
207 // return invariant differential elastic cross section d(sigma)/d(tMand)
208 
209 G4double
211  G4double tMand,
212  G4double plab,
213  G4double A, G4double Z )
214 {
215  G4double m1 = particle->GetPDGMass();
216  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
217 
218  G4int iZ = static_cast<G4int>(Z+0.5);
219  G4int iA = static_cast<G4int>(A+0.5);
220  G4ParticleDefinition * theDef = 0;
221 
222  if (iZ == 1 && iA == 1) theDef = theProton;
223  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
224  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
225  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
226  else if (iZ == 2 && iA == 4) theDef = theAlpha;
227  else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
228 
229  G4double tmass = theDef->GetPDGMass();
230 
231  G4LorentzVector lv(0.0,0.0,0.0,tmass);
232  lv += lv1;
233 
234  G4ThreeVector bst = lv.boostVector();
235  lv1.boost(-bst);
236 
237  G4ThreeVector p1 = lv1.vect();
238  G4double ptot = p1.mag();
239  G4double ptot2 = ptot*ptot;
240  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
241 
242  if( cost >= 1.0 ) cost = 1.0;
243  else if( cost <= -1.0) cost = -1.0;
244 
245  G4double thetaCMS = std::acos(cost);
246 
247  G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A);
248 
249  sigma *= pi/ptot2;
250 
251  return sigma;
252 }
253 
255 //
256 // return differential elastic cross section d(sigma)/d(omega) with Coulomb
257 // correction
258 
259 G4double
261  G4double theta,
262  G4double momentum,
263  G4double A, G4double Z )
264 {
265  fParticle = particle;
266  fWaveVector = momentum/hbarc;
267  fAtomicWeight = A;
268  fAtomicNumber = Z;
269  fNuclearRadius = CalculateNuclearRad(A);
270  fAddCoulomb = false;
271 
272  G4double z = particle->GetPDGCharge();
273 
274  G4double kRt = fWaveVector*fNuclearRadius*theta;
275  G4double kRtC = 1.9;
276 
277  if( z && (kRt > kRtC) )
278  {
279  fAddCoulomb = true;
280  fBeta = CalculateParticleBeta( particle, momentum);
281  fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
282  fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
283  }
284  G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta);
285 
286  return sigma;
287 }
288 
290 //
291 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
292 // correction
293 
294 G4double
296  G4double tMand,
297  G4double plab,
298  G4double A, G4double Z )
299 {
300  G4double m1 = particle->GetPDGMass();
301 
302  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
303 
304  G4int iZ = static_cast<G4int>(Z+0.5);
305  G4int iA = static_cast<G4int>(A+0.5);
306 
307  G4ParticleDefinition* theDef = 0;
308 
309  if (iZ == 1 && iA == 1) theDef = theProton;
310  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
311  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
312  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
313  else if (iZ == 2 && iA == 4) theDef = theAlpha;
314  else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
315 
316  G4double tmass = theDef->GetPDGMass();
317 
318  G4LorentzVector lv(0.0,0.0,0.0,tmass);
319  lv += lv1;
320 
321  G4ThreeVector bst = lv.boostVector();
322  lv1.boost(-bst);
323 
324  G4ThreeVector p1 = lv1.vect();
325  G4double ptot = p1.mag();
326  G4double ptot2 = ptot*ptot;
327  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
328 
329  if( cost >= 1.0 ) cost = 1.0;
330  else if( cost <= -1.0) cost = -1.0;
331 
332  G4double thetaCMS = std::acos(cost);
333 
334  G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z );
335 
336  sigma *= pi/ptot2;
337 
338  return sigma;
339 }
340 
342 //
343 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
344 // correction
345 
346 G4double
348  G4double tMand,
349  G4double plab,
350  G4double A, G4double Z )
351 {
352  G4double m1 = particle->GetPDGMass();
353  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
354 
355  G4int iZ = static_cast<G4int>(Z+0.5);
356  G4int iA = static_cast<G4int>(A+0.5);
357  G4ParticleDefinition * theDef = 0;
358 
359  if (iZ == 1 && iA == 1) theDef = theProton;
360  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
361  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
362  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
363  else if (iZ == 2 && iA == 4) theDef = theAlpha;
364  else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
365 
366  G4double tmass = theDef->GetPDGMass();
367 
368  G4LorentzVector lv(0.0,0.0,0.0,tmass);
369  lv += lv1;
370 
371  G4ThreeVector bst = lv.boostVector();
372  lv1.boost(-bst);
373 
374  G4ThreeVector p1 = lv1.vect();
375  G4double ptot = p1.mag();
376  G4double ptot2 = ptot*ptot;
377  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
378 
379  if( cost >= 1.0 ) cost = 1.0;
380  else if( cost <= -1.0) cost = -1.0;
381 
382  G4double thetaCMS = std::acos(cost);
383 
384  G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z );
385 
386  sigma *= pi/ptot2;
387 
388  return sigma;
389 }
390 
392 //
393 // return differential elastic probability d(probability)/d(omega)
394 
395 G4double
396 G4NuclNuclDiffuseElastic::GetDiffElasticProb( // G4ParticleDefinition* particle,
397  G4double theta
398  // G4double momentum,
399  // G4double A
400  )
401 {
402  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
403  G4double delta, diffuse, gamma;
404  G4double e1, e2, bone, bone2;
405 
406  // G4double wavek = momentum/hbarc; // wave vector
407  // G4double r0 = 1.08*fermi;
408  // G4double rad = r0*G4Pow::GetInstance()->A13(A);
409 
410  G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
411  G4double kr2 = kr*kr;
412  G4double krt = kr*theta;
413 
414  bzero = BesselJzero(krt);
415  bzero2 = bzero*bzero;
416  bone = BesselJone(krt);
417  bone2 = bone*bone;
418  bonebyarg = BesselOneByArg(krt);
419  bonebyarg2 = bonebyarg*bonebyarg;
420 
421  // VI - Coverity complains
422  /*
423  if (fParticle == theProton)
424  {
425  diffuse = 0.63*fermi;
426  gamma = 0.3*fermi;
427  delta = 0.1*fermi*fermi;
428  e1 = 0.3*fermi;
429  e2 = 0.35*fermi;
430  }
431  else // as proton, if were not defined
432  {
433  */
434  diffuse = 0.63*fermi;
435  gamma = 0.3*fermi;
436  delta = 0.1*fermi*fermi;
437  e1 = 0.3*fermi;
438  e2 = 0.35*fermi;
439  //}
440  G4double lambda = 15.; // 15 ok
441 
442  // G4double kgamma = fWaveVector*gamma; // wavek*delta;
443 
444  G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
445  G4double kgamma2 = kgamma*kgamma;
446 
447  // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
448  // G4double dk2t2 = dk2t*dk2t;
449  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
450 
451  G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
452 
453  damp = DampFactor(pikdt);
454  damp2 = damp*damp;
455 
456  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
457  G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
458 
459 
460  sigma = kgamma2;
461  // sigma += dk2t2;
462  sigma *= bzero2;
463  sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
464  sigma += kr2*bonebyarg2;
465  sigma *= damp2; // *rad*rad;
466 
467  return sigma;
468 }
469 
471 //
472 // return differential elastic probability d(probability)/d(omega) with
473 // Coulomb correction
474 
475 G4double
476 G4NuclNuclDiffuseElastic::GetDiffElasticSumProb( // G4ParticleDefinition* particle,
477  G4double theta
478  // G4double momentum,
479  // G4double A
480  )
481 {
482  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
483  G4double delta, diffuse, gamma;
484  G4double e1, e2, bone, bone2;
485 
486  // G4double wavek = momentum/hbarc; // wave vector
487  // G4double r0 = 1.08*fermi;
488  // G4double rad = r0*G4Pow::GetInstance()->A13(A);
489 
490  G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
491  G4double kr2 = kr*kr;
492  G4double krt = kr*theta;
493 
494  bzero = BesselJzero(krt);
495  bzero2 = bzero*bzero;
496  bone = BesselJone(krt);
497  bone2 = bone*bone;
498  bonebyarg = BesselOneByArg(krt);
499  bonebyarg2 = bonebyarg*bonebyarg;
500 
501  if (fParticle == theProton)
502  {
503  diffuse = 0.63*fermi;
504  // diffuse = 0.6*fermi;
505  gamma = 0.3*fermi;
506  delta = 0.1*fermi*fermi;
507  e1 = 0.3*fermi;
508  e2 = 0.35*fermi;
509  }
510  else // as proton, if were not defined
511  {
512  diffuse = 0.63*fermi;
513  gamma = 0.3*fermi;
514  delta = 0.1*fermi*fermi;
515  e1 = 0.3*fermi;
516  e2 = 0.35*fermi;
517  }
518  G4double lambda = 15.; // 15 ok
519  // G4double kgamma = fWaveVector*gamma; // wavek*delta;
520  G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
521 
522  // G4cout<<"kgamma = "<<kgamma<<G4endl;
523 
524  if(fAddCoulomb) // add Coulomb correction
525  {
526  G4double sinHalfTheta = std::sin(0.5*theta);
527  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
528 
529  kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
530  // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
531  }
532 
533  G4double kgamma2 = kgamma*kgamma;
534 
535  // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
536  // G4cout<<"dk2t = "<<dk2t<<G4endl;
537  // G4double dk2t2 = dk2t*dk2t;
538  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
539 
540  G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
541 
542  // G4cout<<"pikdt = "<<pikdt<<G4endl;
543 
544  damp = DampFactor(pikdt);
545  damp2 = damp*damp;
546 
547  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
548  G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
549 
550  sigma = kgamma2;
551  // sigma += dk2t2;
552  sigma *= bzero2;
553  sigma += mode2k2*bone2;
554  sigma += e2dk3t*bzero*bone;
555 
556  // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
557  sigma += kr2*bonebyarg2; // correction at J1()/()
558 
559  sigma *= damp2; // *rad*rad;
560 
561  return sigma;
562 }
563 
564 
566 //
567 // return differential elastic probability d(probability)/d(t) with
568 // Coulomb correction
569 
570 G4double
572 {
573  G4double theta;
574 
575  theta = std::sqrt(alpha);
576 
577  // theta = std::acos( 1 - alpha/2. );
578 
579  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
580  G4double delta, diffuse, gamma;
581  G4double e1, e2, bone, bone2;
582 
583  // G4double wavek = momentum/hbarc; // wave vector
584  // G4double r0 = 1.08*fermi;
585  // G4double rad = r0*G4Pow::GetInstance()->A13(A);
586 
587  G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
588  G4double kr2 = kr*kr;
589  G4double krt = kr*theta;
590 
591  bzero = BesselJzero(krt);
592  bzero2 = bzero*bzero;
593  bone = BesselJone(krt);
594  bone2 = bone*bone;
595  bonebyarg = BesselOneByArg(krt);
596  bonebyarg2 = bonebyarg*bonebyarg;
597 
598  if (fParticle == theProton)
599  {
600  diffuse = 0.63*fermi;
601  // diffuse = 0.6*fermi;
602  gamma = 0.3*fermi;
603  delta = 0.1*fermi*fermi;
604  e1 = 0.3*fermi;
605  e2 = 0.35*fermi;
606  }
607  else // as proton, if were not defined
608  {
609  diffuse = 0.63*fermi;
610  gamma = 0.3*fermi;
611  delta = 0.1*fermi*fermi;
612  e1 = 0.3*fermi;
613  e2 = 0.35*fermi;
614  }
615  G4double lambda = 15.; // 15 ok
616  // G4double kgamma = fWaveVector*gamma; // wavek*delta;
617  G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
618 
619  // G4cout<<"kgamma = "<<kgamma<<G4endl;
620 
621  if(fAddCoulomb) // add Coulomb correction
622  {
623  G4double sinHalfTheta = theta*0.5; // std::sin(0.5*theta);
624  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
625 
626  kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
627  // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
628  }
629 
630  G4double kgamma2 = kgamma*kgamma;
631 
632  // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
633  // G4cout<<"dk2t = "<<dk2t<<G4endl;
634  // G4double dk2t2 = dk2t*dk2t;
635  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
636 
637  G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
638 
639  // G4cout<<"pikdt = "<<pikdt<<G4endl;
640 
641  damp = DampFactor(pikdt);
642  damp2 = damp*damp;
643 
644  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
645  G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
646 
647  sigma = kgamma2;
648  // sigma += dk2t2;
649  sigma *= bzero2;
650  sigma += mode2k2*bone2;
651  sigma += e2dk3t*bzero*bone;
652 
653  // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
654  sigma += kr2*bonebyarg2; // correction at J1()/()
655 
656  sigma *= damp2; // *rad*rad;
657 
658  return sigma;
659 }
660 
661 
663 //
664 // return differential elastic probability 2*pi*sin(theta)*d(probability)/d(omega)
665 
666 G4double
668 {
670 
671  result = GetDiffElasticSumProbA(alpha);
672 
673  // result *= 2*pi*std::sin(theta);
674 
675  return result;
676 }
677 
679 //
680 // return integral elastic cross section d(sigma)/d(omega) integrated 0 - theta
681 
682 G4double
684  G4double theta,
685  G4double momentum,
686  G4double A )
687 {
689  fParticle = particle;
690  fWaveVector = momentum/hbarc;
691  fAtomicWeight = A;
692 
693  fNuclearRadius = CalculateNuclearRad(A);
694 
695 
697 
698  // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
699  result = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
700 
701  return result;
702 }
703 
705 //
706 // Return inv momentum transfer -t > 0
707 
709  G4double p, G4double A)
710 {
711  G4double theta = SampleThetaCMS( aParticle, p, A); // sample theta in cms
712  G4double t = 2*p*p*( 1 - std::cos(theta) ); // -t !!!
713  return t;
714 }
715 
717 //
718 // Return scattering angle sampled in cms
719 
720 
721 G4double
723  G4double momentum, G4double A)
724 {
725  G4int i, iMax = 100;
726  G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
727 
728  fParticle = particle;
729  fWaveVector = momentum/hbarc;
730  fAtomicWeight = A;
731 
732  fNuclearRadius = CalculateNuclearRad(A);
733 
734  thetaMax = 10.174/fWaveVector/fNuclearRadius;
735 
736  if (thetaMax > pi) thetaMax = pi;
737 
739 
740  // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
741  norm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., thetaMax );
742 
743  norm *= G4UniformRand();
744 
745  for(i = 1; i <= iMax; i++)
746  {
747  theta1 = (i-1)*thetaMax/iMax;
748  theta2 = i*thetaMax/iMax;
749  sum += integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, theta1, theta2);
750 
751  if ( sum >= norm )
752  {
753  result = 0.5*(theta1 + theta2);
754  break;
755  }
756  }
757  if (i > iMax ) result = 0.5*(theta1 + theta2);
758 
759  G4double sigma = pi*thetaMax/iMax;
760 
761  result += G4RandGauss::shoot(0.,sigma);
762 
763  if(result < 0.) result = 0.;
764  if(result > thetaMax) result = thetaMax;
765 
766  return result;
767 }
768 
771 
773 //
774 // Return inv momentum transfer -t > 0 from initialisation table
775 
777  G4int Z, G4int A)
778 {
779  fParticle = aParticle;
780  G4double m1 = fParticle->GetPDGMass();
781  G4double totElab = std::sqrt(m1*m1+p*p);
783  G4LorentzVector lv1(p,0.0,0.0,totElab);
784  G4LorentzVector lv(0.0,0.0,0.0,mass2);
785  lv += lv1;
786 
787  G4ThreeVector bst = lv.boostVector();
788  lv1.boost(-bst);
789 
790  G4ThreeVector p1 = lv1.vect();
791  G4double momentumCMS = p1.mag();
792 
793  G4double t = SampleTableT( aParticle, momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms
794  return t;
795 }
796 
798 //
799 // Return inv momentum transfer -t > 0 from initialisation table
800 
802  G4double Z, G4double A)
803 {
804  G4double alpha = SampleTableThetaCMS( aParticle, p, Z, A); // sample theta2 in cms
805  // G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) ); // -t !!!
806  G4double t = p*p*alpha; // -t !!!
807  return t;
808 }
809 
811 //
812 // Return scattering angle2 sampled in cms according to precalculated table.
813 
814 
815 G4double
817  G4double momentum, G4double Z, G4double A)
818 {
819  size_t iElement;
820  G4int iMomentum, iAngle;
821  G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W;
822  G4double m1 = particle->GetPDGMass();
823 
824  for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
825  {
826  if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
827  }
828  if ( iElement == fElementNumberVector.size() )
829  {
830  InitialiseOnFly(Z,A); // table preparation, if needed
831 
832  // iElement--;
833 
834  // G4cout << "G4NuclNuclDiffuseElastic: Element with atomic number " << Z
835  // << " is not found, return zero angle" << G4endl;
836  // return 0.; // no table for this element
837  }
838  // G4cout<<"iElement = "<<iElement<<G4endl;
839 
840  fAngleTable = fAngleBank[iElement];
841 
842  G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
843 
844  for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
845  {
846  // G4cout<<iMomentum<<", kinE = "<<kinE/MeV<<", vectorE = "<<fEnergyVector->GetLowEdgeEnergy(iMomentum)/MeV<<G4endl;
847 
848  if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
849  }
850  // G4cout<<"iMomentum = "<<iMomentum<<", fEnergyBin -1 = "<<fEnergyBin -1<<G4endl;
851 
852 
853  if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1; // kinE is more then theMaxEnergy
854  if ( iMomentum < 0 ) iMomentum = 0; // against negative index, kinE < theMinEnergy
855 
856 
857  if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges
858  {
859  position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
860 
861  // G4cout<<"position = "<<position<<G4endl;
862 
863  for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
864  {
865  if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
866  }
867  if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
868 
869  // G4cout<<"iAngle = "<<iAngle<<G4endl;
870 
871  randAngle = GetScatteringAngle(iMomentum, iAngle, position);
872 
873  // G4cout<<"randAngle = "<<randAngle<<G4endl;
874  }
875  else // kinE inside between energy table edges
876  {
877  // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
878  position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand();
879 
880  // G4cout<<"position = "<<position<<G4endl;
881 
882  for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
883  {
884  // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
885  if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
886  }
887  if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
888 
889  // G4cout<<"iAngle = "<<iAngle<<G4endl;
890 
891  theta2 = GetScatteringAngle(iMomentum, iAngle, position);
892 
893  // G4cout<<"theta2 = "<<theta2<<G4endl;
894 
895  E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
896 
897  // G4cout<<"E2 = "<<E2<<G4endl;
898 
899  iMomentum--;
900 
901  // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
902 
903  // G4cout<<"position = "<<position<<G4endl;
904 
905  for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
906  {
907  // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
908  if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
909  }
910  if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
911 
912  theta1 = GetScatteringAngle(iMomentum, iAngle, position);
913 
914  // G4cout<<"theta1 = "<<theta1<<G4endl;
915 
916  E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
917 
918  // G4cout<<"E1 = "<<E1<<G4endl;
919 
920  W = 1.0/(E2 - E1);
921  W1 = (E2 - kinE)*W;
922  W2 = (kinE - E1)*W;
923 
924  randAngle = W1*theta1 + W2*theta2;
925 
926  // randAngle = theta2;
927  }
928  // G4double angle = randAngle;
929  // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
930  // G4cout<<"randAngle = "<<randAngle/degree<<G4endl;
931 
932  return randAngle;
933 }
934 
936 //
937 // Initialisation for given particle on fly using new element number
938 
940 {
941  fAtomicNumber = Z; // atomic number
942  fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( Z ) );
943 
944  G4double A1 = G4double( fParticle->GetBaryonNumber() );
945  G4double R1 = CalculateNuclearRad(A1);
946 
947  fNuclearRadius = CalculateNuclearRad(fAtomicWeight) + R1;
948 
949  if( verboseLevel > 0 )
950  {
951  G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element with Z = "
952  <<Z<<"; and A = "<<A<<G4endl;
953  }
954  fElementNumberVector.push_back(fAtomicNumber);
955 
956  BuildAngleTable();
957 
958  fAngleBank.push_back(fAngleTable);
959 
960  return;
961 }
962 
964 //
965 // Build for given particle and element table of momentum, angle probability.
966 // For the moment in lab system.
967 
969 {
970  G4int i, j;
971  G4double partMom, kinE, m1 = fParticle->GetPDGMass();
972  G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
973 
974  // G4cout<<"particle z = "<<z<<"; particle m1 = "<<m1/GeV<<" GeV"<<G4endl;
975 
977 
978  fAngleTable = new G4PhysicsTable(fEnergyBin);
979 
980  for( i = 0; i < fEnergyBin; i++)
981  {
982  kinE = fEnergyVector->GetLowEdgeEnergy(i);
983 
984  // G4cout<<G4endl;
985  // G4cout<<"kinE = "<<kinE/MeV<<" MeV"<<G4endl;
986 
987  partMom = std::sqrt( kinE*(kinE + 2*m1) );
988 
989  InitDynParameters(fParticle, partMom);
990 
991  alphaMax = fRutherfordTheta*fCofAlphaMax;
992 
993  if(alphaMax > pi) alphaMax = pi;
994 
995  // VI: Coverity complain
996  //alphaMax = pi2;
997 
998  alphaCoulomb = fRutherfordTheta*fCofAlphaCoulomb;
999 
1000  // G4cout<<"alphaCoulomb = "<<alphaCoulomb/degree<<"; alphaMax = "<<alphaMax/degree<<G4endl;
1001 
1002  G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1003 
1004  // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1005 
1006  G4double delth = (alphaMax-alphaCoulomb)/fAngleBin;
1007 
1008  sum = 0.;
1009 
1010  // fAddCoulomb = false;
1011  fAddCoulomb = true;
1012 
1013  // for(j = 1; j < fAngleBin; j++)
1014  for(j = fAngleBin-1; j >= 1; j--)
1015  {
1016  // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1017  // alpha2 = angleBins->GetLowEdgeEnergy(j);
1018 
1019  // alpha1 = alphaMax*(j-1)/fAngleBin;
1020  // alpha2 = alphaMax*( j )/fAngleBin;
1021 
1022  alpha1 = alphaCoulomb + delth*(j-1);
1023  // if(alpha1 < kRlim2) alpha1 = kRlim2;
1024  alpha2 = alpha1 + delth;
1025 
1026  delta = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc, alpha1, alpha2);
1027  // delta = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1028 
1029  sum += delta;
1030 
1031  angleVector->PutValue( j-1 , alpha1, sum ); // alpha2
1032  // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2/degree<<"; sum = "<<sum<<G4endl;
1033  }
1034  fAngleTable->insertAt(i,angleVector);
1035 
1036  // delete[] angleVector;
1037  // delete[] angleBins;
1038  }
1039  return;
1040 }
1041 
1043 //
1044 //
1045 
1046 G4double
1048 {
1049  G4double x1, x2, y1, y2, randAngle;
1050 
1051  if( iAngle == 0 )
1052  {
1053  randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1054  // iAngle++;
1055  }
1056  else
1057  {
1058  if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
1059  {
1060  iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1061  }
1062  y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1063  y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1064 
1065  x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1066  x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1067 
1068  if ( x1 == x2 ) randAngle = x2;
1069  else
1070  {
1071  if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand();
1072  else
1073  {
1074  randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
1075  }
1076  }
1077  }
1078  return randAngle;
1079 }
1080 
1081 
1082 
1084 //
1085 // Return scattering angle sampled in lab system (target at rest)
1086 
1087 
1088 
1089 G4double
1091  G4double tmass, G4double A)
1092 {
1093  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1094  G4double m1 = theParticle->GetPDGMass();
1095  G4double plab = aParticle->GetTotalMomentum();
1096  G4LorentzVector lv1 = aParticle->Get4Momentum();
1097  G4LorentzVector lv(0.0,0.0,0.0,tmass);
1098  lv += lv1;
1099 
1100  G4ThreeVector bst = lv.boostVector();
1101  lv1.boost(-bst);
1102 
1103  G4ThreeVector p1 = lv1.vect();
1104  G4double ptot = p1.mag();
1105  G4double tmax = 4.0*ptot*ptot;
1106  G4double t = 0.0;
1107 
1108 
1109  //
1110  // Sample t
1111  //
1112 
1113  t = SampleT( theParticle, ptot, A);
1114 
1115  // NaN finder
1116  if(!(t < 0.0 || t >= 0.0))
1117  {
1118  if (verboseLevel > 0)
1119  {
1120  G4cout << "G4NuclNuclDiffuseElastic:WARNING: A = " << A
1121  << " mom(GeV)= " << plab/GeV
1122  << " S-wave will be sampled"
1123  << G4endl;
1124  }
1125  t = G4UniformRand()*tmax;
1126  }
1127  if(verboseLevel>1)
1128  {
1129  G4cout <<" t= " << t << " tmax= " << tmax
1130  << " ptot= " << ptot << G4endl;
1131  }
1132  // Sampling of angles in CM system
1133 
1134  G4double phi = G4UniformRand()*twopi;
1135  G4double cost = 1. - 2.0*t/tmax;
1136  G4double sint;
1137 
1138  if( cost >= 1.0 )
1139  {
1140  cost = 1.0;
1141  sint = 0.0;
1142  }
1143  else if( cost <= -1.0)
1144  {
1145  cost = -1.0;
1146  sint = 0.0;
1147  }
1148  else
1149  {
1150  sint = std::sqrt((1.0-cost)*(1.0+cost));
1151  }
1152  if (verboseLevel>1)
1153  {
1154  G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
1155  }
1156  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1157  v1 *= ptot;
1158  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1159 
1160  nlv1.boost(bst);
1161 
1162  G4ThreeVector np1 = nlv1.vect();
1163 
1164  // G4double theta = std::acos( np1.z()/np1.mag() ); // degree;
1165 
1166  G4double theta = np1.theta();
1167 
1168  return theta;
1169 }
1170 
1172 //
1173 // Return scattering angle in lab system (target at rest) knowing theta in CMS
1174 
1175 
1176 
1177 G4double
1179  G4double tmass, G4double thetaCMS)
1180 {
1181  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1182  G4double m1 = theParticle->GetPDGMass();
1183  // G4double plab = aParticle->GetTotalMomentum();
1184  G4LorentzVector lv1 = aParticle->Get4Momentum();
1185  G4LorentzVector lv(0.0,0.0,0.0,tmass);
1186 
1187  lv += lv1;
1188 
1189  G4ThreeVector bst = lv.boostVector();
1190 
1191  lv1.boost(-bst);
1192 
1193  G4ThreeVector p1 = lv1.vect();
1194  G4double ptot = p1.mag();
1195 
1196  G4double phi = G4UniformRand()*twopi;
1197  G4double cost = std::cos(thetaCMS);
1198  G4double sint;
1199 
1200  if( cost >= 1.0 )
1201  {
1202  cost = 1.0;
1203  sint = 0.0;
1204  }
1205  else if( cost <= -1.0)
1206  {
1207  cost = -1.0;
1208  sint = 0.0;
1209  }
1210  else
1211  {
1212  sint = std::sqrt((1.0-cost)*(1.0+cost));
1213  }
1214  if (verboseLevel>1)
1215  {
1216  G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
1217  }
1218  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1219  v1 *= ptot;
1220  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1221 
1222  nlv1.boost(bst);
1223 
1224  G4ThreeVector np1 = nlv1.vect();
1225 
1226 
1227  G4double thetaLab = np1.theta();
1228 
1229  return thetaLab;
1230 }
1231 
1233 //
1234 // Return scattering angle in CMS system (target at rest) knowing theta in Lab
1235 
1236 
1237 
1238 G4double
1240  G4double tmass, G4double thetaLab)
1241 {
1242  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1243  G4double m1 = theParticle->GetPDGMass();
1244  G4double plab = aParticle->GetTotalMomentum();
1245  G4LorentzVector lv1 = aParticle->Get4Momentum();
1246  G4LorentzVector lv(0.0,0.0,0.0,tmass);
1247 
1248  lv += lv1;
1249 
1250  G4ThreeVector bst = lv.boostVector();
1251 
1252  // lv1.boost(-bst);
1253 
1254  // G4ThreeVector p1 = lv1.vect();
1255  // G4double ptot = p1.mag();
1256 
1257  G4double phi = G4UniformRand()*twopi;
1258  G4double cost = std::cos(thetaLab);
1259  G4double sint;
1260 
1261  if( cost >= 1.0 )
1262  {
1263  cost = 1.0;
1264  sint = 0.0;
1265  }
1266  else if( cost <= -1.0)
1267  {
1268  cost = -1.0;
1269  sint = 0.0;
1270  }
1271  else
1272  {
1273  sint = std::sqrt((1.0-cost)*(1.0+cost));
1274  }
1275  if (verboseLevel>1)
1276  {
1277  G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
1278  }
1279  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1280  v1 *= plab;
1281  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
1282 
1283  nlv1.boost(-bst);
1284 
1285  G4ThreeVector np1 = nlv1.vect();
1286 
1287 
1288  G4double thetaCMS = np1.theta();
1289 
1290  return thetaCMS;
1291 }
1292 
1294 //
1295 // Test for given particle and element table of momentum, angle probability.
1296 // For the moment in lab system.
1297 
1299  G4double Z, G4double A)
1300 {
1301  fAtomicNumber = Z; // atomic number
1302  fAtomicWeight = A; // number of nucleons
1303  fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
1304 
1305 
1306 
1307  G4cout<<"G4NuclNuclDiffuseElastic::TestAngleTable() init the element with Z = "
1308  <<Z<<"; and A = "<<A<<G4endl;
1309 
1310  fElementNumberVector.push_back(fAtomicNumber);
1311 
1312 
1313 
1314 
1315  G4int i=0, j;
1316  G4double a = 0., z = theParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
1317  G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
1318  G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1319  G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1320  G4double epsilon = 0.001;
1321 
1323 
1324  fAngleTable = new G4PhysicsTable(fEnergyBin);
1325 
1326  fWaveVector = partMom/hbarc;
1327 
1328  G4double kR = fWaveVector*fNuclearRadius;
1329  G4double kR2 = kR*kR;
1330  G4double kRmax = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
1331  G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1
1332 
1333  alphaMax = kRmax*kRmax/kR2;
1334 
1335  if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2
1336 
1337  alphaCoulomb = kRcoul*kRcoul/kR2;
1338 
1339  if( z )
1340  {
1341  a = partMom/m1; // beta*gamma for m1
1342  fBeta = a/std::sqrt(1+a*a);
1343  fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1344  fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1345  }
1346  G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1347 
1348  // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1349 
1350 
1351  fAddCoulomb = false;
1352 
1353  for(j = 1; j < fAngleBin; j++)
1354  {
1355  // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1356  // alpha2 = angleBins->GetLowEdgeEnergy(j);
1357 
1358  alpha1 = alphaMax*(j-1)/fAngleBin;
1359  alpha2 = alphaMax*( j )/fAngleBin;
1360 
1361  if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1362 
1363  deltaL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1364  deltaL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1366  alpha1, alpha2,epsilon);
1367 
1368  // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1369  // <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl;
1370 
1371  sumL10 += deltaL10;
1372  sumL96 += deltaL96;
1373  sumAG += deltaAG;
1374 
1375  G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1376  <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1377 
1378  angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2
1379  }
1380  fAngleTable->insertAt(i,angleVector);
1381  fAngleBank.push_back(fAngleTable);
1382 
1383  /*
1384  // Integral over all angle range - Bad accuracy !!!
1385 
1386  sumL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1387  sumL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1388  sumAG = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction,
1389  0., alpha2,epsilon);
1390  G4cout<<G4endl;
1391  G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
1392  <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1393  */
1394  return;
1395 }
1396 
1398 //
1399 //
1400 
1402 {
1403  G4double legPol, epsilon = 1.e-6;
1404  G4double x = std::cos(theta);
1405 
1406  if ( n < 0 ) legPol = 0.;
1407  else if( n == 0 ) legPol = 1.;
1408  else if( n == 1 ) legPol = x;
1409  else if( n == 2 ) legPol = (3.*x*x-1.)/2.;
1410  else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/2.;
1411  else if( n == 4 ) legPol = (35.*x*x*x*x-30.*x*x+3.)/8.;
1412  else if( n == 5 ) legPol = (63.*x*x*x*x*x-70.*x*x*x+15.*x)/8.;
1413  else if( n == 6 ) legPol = (231.*x*x*x*x*x*x-315.*x*x*x*x+105.*x*x-5.)/16.;
1414  else
1415  {
1416  // legPol = ( (2*n-1)*x*GetLegendrePol(n-1,x) - (n-1)*GetLegendrePol(n-2,x) )/n;
1417 
1418  legPol = std::sqrt( 2./(n*CLHEP::pi*std::sin(theta+epsilon)) )*std::sin( (n+0.5)*theta+0.25*CLHEP::pi );
1419  }
1420  return legPol;
1421 }
1422 
1424 //
1425 //
1426 
1428 {
1429  G4int n;
1430  G4double n2, cofn, shny, chny, fn, gn;
1431 
1432  G4double x = z.real();
1433  G4double y = z.imag();
1434 
1435  G4double outRe = 0., outIm = 0.;
1436 
1437  G4double twox = 2.*x;
1438  G4double twoxy = twox*y;
1439  G4double twox2 = twox*twox;
1440 
1441  G4double cof1 = G4Exp(-x*x)/CLHEP::pi;
1442 
1443  G4double cos2xy = std::cos(twoxy);
1444  G4double sin2xy = std::sin(twoxy);
1445 
1446  G4double twoxcos2xy = twox*cos2xy;
1447  G4double twoxsin2xy = twox*sin2xy;
1448 
1449  for( n = 1; n <= nMax; n++)
1450  {
1451  n2 = n*n;
1452 
1453  cofn = G4Exp(-0.5*n2)/(n2+twox2); // /(n2+0.5*twox2);
1454 
1455  chny = std::cosh(n*y);
1456  shny = std::sinh(n*y);
1457 
1458  fn = twox - twoxcos2xy*chny + n*sin2xy*shny;
1459  gn = twoxsin2xy*chny + n*cos2xy*shny;
1460 
1461  fn *= cofn;
1462  gn *= cofn;
1463 
1464  outRe += fn;
1465  outIm += gn;
1466  }
1467  outRe *= 2*cof1;
1468  outIm *= 2*cof1;
1469 
1470  if(std::abs(x) < 0.0001)
1471  {
1472  outRe += GetErf(x);
1473  outIm += cof1*y;
1474  }
1475  else
1476  {
1477  outRe += GetErf(x) + cof1*(1-cos2xy)/twox;
1478  outIm += cof1*sin2xy/twox;
1479  }
1480  return G4complex(outRe, outIm);
1481 }
1482 
1483 
1485 //
1486 //
1487 
1489 {
1490  G4double outRe, outIm;
1491 
1492  G4double x = z.real();
1493  G4double y = z.imag();
1494  fReZ = x;
1495 
1497 
1498  outRe = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpSin, 0., y );
1499  outIm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpCos, 0., y );
1500 
1501  outRe *= 2./std::sqrt(CLHEP::pi);
1502  outIm *= 2./std::sqrt(CLHEP::pi);
1503 
1504  outRe += GetErf(x);
1505 
1506  return G4complex(outRe, outIm);
1507 }
1508 
1510 //
1511 //
1512 
1513 
1515 {
1516  G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1517  G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1518 
1519  G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1520  G4double kappa = u/std::sqrt(CLHEP::pi);
1521  G4double dTheta = theta - fRutherfordTheta;
1522  u *= dTheta;
1523  G4double u2 = u*u;
1524  G4double u2m2p3 = u2*2./3.;
1525 
1526  G4complex im = G4complex(0.,1.);
1527  G4complex order = G4complex(u,u);
1528  order /= std::sqrt(2.);
1529 
1530  G4complex gamma = CLHEP::pi*kappa*GetErfcInt(-order)*std::exp(im*(u*u+0.25*CLHEP::pi));
1531  G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1532  G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1533  G4complex out = gamma*(1. - a1*dTheta) - a0;
1534 
1535  return out;
1536 }
1537 
1539 //
1540 //
1541 
1543 {
1544  G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1545  G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1546 
1547  G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1548  G4double kappa = u/std::sqrt(CLHEP::pi);
1549  G4double dTheta = theta - fRutherfordTheta;
1550  u *= dTheta;
1551  G4double u2 = u*u;
1552  G4double u2m2p3 = u2*2./3.;
1553 
1554  G4complex im = G4complex(0.,1.);
1555  G4complex order = G4complex(u,u);
1556  order /= std::sqrt(2.);
1557  G4complex gamma = CLHEP::pi*kappa*GetErfcInt(order)*std::exp(im*(u*u+0.25*CLHEP::pi));
1558  G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1559  G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1560  G4complex out = -gamma*(1. - a1*dTheta) - a0;
1561 
1562  return out;
1563 }
1564 
1566 //
1567 //
1568 
1570 {
1571  G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
1572  G4complex out = G4complex(kappa/fWaveVector,0.);
1573 
1574  out *= PhaseNear(theta);
1575 
1576  if( theta <= fRutherfordTheta )
1577  {
1578  out *= GammaLess(theta) + ProfileNear(theta);
1579  // out *= GammaMore(theta) + ProfileNear(theta);
1580  out += CoulombAmplitude(theta);
1581  }
1582  else
1583  {
1584  out *= GammaMore(theta) + ProfileNear(theta);
1585  // out *= GammaLess(theta) + ProfileNear(theta);
1586  }
1587  return out;
1588 }
1589 
1591 //
1592 //
1593 
1595 {
1596  G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1597  G4double dTheta = 0.5*(theta - fRutherfordTheta);
1598  G4double sindTheta = std::sin(dTheta);
1599  G4double persqrt2 = std::sqrt(0.5);
1600 
1601  G4complex order = G4complex(persqrt2,persqrt2);
1602  order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*sindTheta;
1603  // order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*dTheta;
1604 
1605  G4complex out;
1606 
1607  if ( theta <= fRutherfordTheta )
1608  {
1609  out = 1. - 0.5*GetErfcInt(-order)*ProfileNear(theta);
1610  }
1611  else
1612  {
1613  out = 0.5*GetErfcInt(order)*ProfileNear(theta);
1614  }
1615 
1616  out *= CoulombAmplitude(theta);
1617 
1618  return out;
1619 }
1620 
1622 //
1623 //
1624 
1626 {
1627  G4int n;
1628  G4double T12b, b, b2; // cosTheta = std::cos(theta);
1629  G4complex out = G4complex(0.,0.), shiftC, shiftN;
1630  G4complex im = G4complex(0.,1.);
1631 
1632  for( n = 0; n < fMaxL; n++)
1633  {
1634  shiftC = std::exp( im*2.*CalculateCoulombPhase(n) );
1635  // b = ( fZommerfeld + std::sqrt( fZommerfeld*fZommerfeld + n*(n+1) ) )/fWaveVector;
1636  b = ( std::sqrt( G4double(n*(n+1)) ) )/fWaveVector;
1637  b2 = b*b;
1638  T12b = fSumSigma*G4Exp(-b2/fNuclearRadiusSquare)/CLHEP::pi/fNuclearRadiusSquare;
1639  shiftN = std::exp( -0.5*(1.-im*fEtaRatio)*T12b ) - 1.;
1640  out += (2.*n+1.)*shiftC*shiftN*GetLegendrePol(n, theta);
1641  }
1642  out /= 2.*im*fWaveVector;
1643  out += CoulombAmplitude(theta);
1644  return out;
1645 }
1646 
1647 
1649 //
1650 //
1651 
1653 {
1654  G4int n;
1655  G4double T12b, a, aTemp, b2, sinThetaH = std::sin(0.5*theta);
1656  G4double sinThetaH2 = sinThetaH*sinThetaH;
1657  G4complex out = G4complex(0.,0.);
1658  G4complex im = G4complex(0.,1.);
1659 
1660  a = -fSumSigma/CLHEP::twopi/fNuclearRadiusSquare;
1661  b2 = fWaveVector*fWaveVector*fNuclearRadiusSquare*sinThetaH2;
1662 
1663  aTemp = a;
1664 
1665  for( n = 1; n < fMaxL; n++)
1666  {
1667  T12b = aTemp*G4Exp(-b2/n)/n;
1668  aTemp *= a;
1669  out += T12b;
1670  G4cout<<"out = "<<out<<G4endl;
1671  }
1672  out *= -4.*im*fWaveVector/CLHEP::pi;
1673  out += CoulombAmplitude(theta);
1674  return out;
1675 }
1676 
1677 
1679 //
1680 // Test for given particle and element table of momentum, angle probability.
1681 // For the partMom in CMS.
1682 
1684  G4double partMom, G4double Z, G4double A)
1685 {
1686  fAtomicNumber = Z; // atomic number
1687  fAtomicWeight = A; // number of nucleons
1688 
1689  fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight);
1690  G4double A1 = G4double( theParticle->GetBaryonNumber() );
1691  fNuclearRadius1 = CalculateNuclearRad(A1);
1692  // fNuclearRadius = std::sqrt(fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2);
1693  fNuclearRadius = fNuclearRadius1 + fNuclearRadius2;
1694 
1695  G4double a = 0.;
1696  G4double z = theParticle->GetPDGCharge();
1697  G4double m1 = theParticle->GetPDGMass();
1698 
1699  fWaveVector = partMom/CLHEP::hbarc;
1700 
1701  G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
1702  G4cout<<"kR = "<<lambda<<G4endl;
1703 
1704  if( z )
1705  {
1706  a = partMom/m1; // beta*gamma for m1
1707  fBeta = a/std::sqrt(1+a*a);
1708  fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1709  fRutherfordRatio = fZommerfeld/fWaveVector;
1710  fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1711  }
1712  G4cout<<"fZommerfeld = "<<fZommerfeld<<G4endl;
1713  fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
1714  G4cout<<"fProfileLambda = "<<fProfileLambda<<G4endl;
1715  fProfileDelta = fCofDelta*fProfileLambda;
1716  fProfileAlpha = fCofAlpha*fProfileLambda;
1717 
1720 
1721  return;
1722 }
1724 //
1725 // Test for given particle and element table of momentum, angle probability.
1726 // For the partMom in CMS.
1727 
1729  G4double partMom)
1730 {
1731  G4double a = 0.;
1732  G4double z = theParticle->GetPDGCharge();
1733  G4double m1 = theParticle->GetPDGMass();
1734 
1735  fWaveVector = partMom/CLHEP::hbarc;
1736 
1737  G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
1738 
1739  if( z )
1740  {
1741  a = partMom/m1; // beta*gamma for m1
1742  fBeta = a/std::sqrt(1+a*a);
1743  fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1744  fRutherfordRatio = fZommerfeld/fWaveVector;
1745  fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1746  }
1747  fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
1748  fProfileDelta = fCofDelta*fProfileLambda;
1749  fProfileAlpha = fCofAlpha*fProfileLambda;
1750 
1753 
1754  return;
1755 }
1756 
1757 
1759 //
1760 // Test for given particle and element table of momentum, angle probability.
1761 // For the partMom in CMS.
1762 
1764  G4double partMom, G4double Z, G4double A)
1765 {
1766  fAtomicNumber = Z; // target atomic number
1767  fAtomicWeight = A; // target number of nucleons
1768 
1769  fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight); // target nucleus radius
1770  G4double A1 = G4double( aParticle->GetDefinition()->GetBaryonNumber() );
1771  fNuclearRadius1 = CalculateNuclearRad(A1); // projectile nucleus radius
1772  fNuclearRadiusSquare = fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2;
1773 
1774 
1775  G4double a = 0., kR12;
1776  G4double z = aParticle->GetDefinition()->GetPDGCharge();
1777  G4double m1 = aParticle->GetDefinition()->GetPDGMass();
1778 
1779  fWaveVector = partMom/CLHEP::hbarc;
1780 
1781  G4double pN = A1 - z;
1782  if( pN < 0. ) pN = 0.;
1783 
1784  G4double tN = A - Z;
1785  if( tN < 0. ) tN = 0.;
1786 
1787  G4double pTkin = aParticle->GetKineticEnergy();
1788  pTkin /= A1;
1789 
1790 
1791  fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
1792  (z*tN+pN*Z)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
1793 
1794  G4cout<<"fSumSigma = "<<fSumSigma/CLHEP::millibarn<<" mb"<<G4endl;
1795  G4cout<<"pi*R2 = "<<CLHEP::pi*fNuclearRadiusSquare/CLHEP::millibarn<<" mb"<<G4endl;
1796  kR12 = fWaveVector*std::sqrt(fNuclearRadiusSquare);
1797  G4cout<<"k*sqrt(R2) = "<<kR12<<" "<<G4endl;
1798  fMaxL = (G4int(kR12)+1)*4;
1799  G4cout<<"fMaxL = "<<fMaxL<<" "<<G4endl;
1800 
1801  if( z )
1802  {
1803  a = partMom/m1; // beta*gamma for m1
1804  fBeta = a/std::sqrt(1+a*a);
1805  fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1806  fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1807  }
1808 
1810 
1811 
1812  return;
1813 }
1814 
1815 
1817 //
1818 // Returns nucleon-nucleon cross-section based on N. Starkov parametrisation of
1819 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
1820 // projectile nucleon is pParticle with pTkin shooting target nucleon tParticle
1821 
1822 G4double
1824  G4double pTkin,
1825  G4ParticleDefinition* tParticle)
1826 {
1827  G4double xsection(0), /*Delta,*/ A0, B0;
1828  G4double hpXsc(0);
1829  G4double hnXsc(0);
1830 
1831 
1832  G4double targ_mass = tParticle->GetPDGMass();
1833  G4double proj_mass = pParticle->GetPDGMass();
1834 
1835  G4double proj_energy = proj_mass + pTkin;
1836  G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
1837 
1838  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
1839 
1840  sMand /= CLHEP::GeV*CLHEP::GeV; // in GeV for parametrisation
1841  proj_momentum /= CLHEP::GeV;
1842  proj_energy /= CLHEP::GeV;
1843  proj_mass /= CLHEP::GeV;
1844  G4double logS = G4Log(sMand);
1845 
1846  // General PDG fit constants
1847 
1848 
1849  // fEtaRatio=Re[f(0)]/Im[f(0)]
1850 
1851  if( proj_momentum >= 1.2 )
1852  {
1853  fEtaRatio = 0.13*(logS - 5.8579332)*G4Pow::GetInstance()->powA(sMand,-0.18);
1854  }
1855  else if( proj_momentum >= 0.6 )
1856  {
1857  fEtaRatio = -75.5*(G4Pow::GetInstance()->powA(proj_momentum,0.25)-0.95)/
1858  (G4Pow::GetInstance()->powA(3*proj_momentum,2.2)+1);
1859  }
1860  else
1861  {
1862  fEtaRatio = 15.5*proj_momentum/(27*proj_momentum*proj_momentum*proj_momentum+2);
1863  }
1864  G4cout<<"fEtaRatio = "<<fEtaRatio<<G4endl;
1865 
1866  // xsc
1867 
1868  if( proj_momentum >= 10. ) // high energy: pp = nn = np
1869  // if( proj_momentum >= 2.)
1870  {
1871  //Delta = 1.;
1872 
1873  //if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
1874 
1875  //AR-12Aug2016 if( proj_momentum >= 10.)
1876  {
1877  B0 = 7.5;
1878  A0 = 100. - B0*G4Log(3.0e7);
1879 
1880  xsection = A0 + B0*G4Log(proj_energy) - 11
1881  + 103*G4Pow::GetInstance()->powA(2*0.93827*proj_energy + proj_mass*proj_mass+
1882  0.93827*0.93827,-0.165); // mb
1883  }
1884  }
1885  else // low energy pp = nn != np
1886  {
1887  if(pParticle == tParticle) // pp or nn // nn to be pp
1888  {
1889  if( proj_momentum < 0.73 )
1890  {
1891  hnXsc = 23 + 50*( G4Pow::GetInstance()->powA( G4Log(0.73/proj_momentum), 3.5 ) );
1892  }
1893  else if( proj_momentum < 1.05 )
1894  {
1895  hnXsc = 23 + 40*(G4Log(proj_momentum/0.73))*
1896  (G4Log(proj_momentum/0.73));
1897  }
1898  else // if( proj_momentum < 10. )
1899  {
1900  hnXsc = 39.0 +
1901  75*(proj_momentum - 1.2)/(G4Pow::GetInstance()->powA(proj_momentum,3.0) + 0.15);
1902  }
1903  xsection = hnXsc;
1904  }
1905  else // pn to be np
1906  {
1907  if( proj_momentum < 0.8 )
1908  {
1909  hpXsc = 33+30*G4Pow::GetInstance()->powA(G4Log(proj_momentum/1.3),4.0);
1910  }
1911  else if( proj_momentum < 1.4 )
1912  {
1913  hpXsc = 33+30*G4Pow::GetInstance()->powA(G4Log(proj_momentum/0.95),2.0);
1914  }
1915  else // if( proj_momentum < 10. )
1916  {
1917  hpXsc = 33.3+
1918  20.8*(G4Pow::GetInstance()->powA(proj_momentum,2.0)-1.35)/
1919  (G4Pow::GetInstance()->powA(proj_momentum,2.50)+0.95);
1920  }
1921  xsection = hpXsc;
1922  }
1923  }
1924  xsection *= CLHEP::millibarn; // parametrised in mb
1925  G4cout<<"xsection = "<<xsection/CLHEP::millibarn<<" mb"<<G4endl;
1926  return xsection;
1927 }
1928 
1930 //
1931 // The ratio el/ruth for Fresnel smooth nucleus profile
1932 
1934 {
1935  G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1936  G4double dTheta = 0.5*(theta - fRutherfordTheta);
1937  G4double sindTheta = std::sin(dTheta);
1938 
1939  G4double prof = Profile(theta);
1940  G4double prof2 = prof*prof;
1941  // G4double profmod = std::abs(prof);
1942  G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
1943 
1944  order = std::abs(order); // since sin changes sign!
1945  // G4cout<<"order = "<<order<<G4endl;
1946 
1947  G4double cosFresnel = GetCint(order);
1948  G4double sinFresnel = GetSint(order);
1949 
1950  G4double out;
1951 
1952  if ( theta <= fRutherfordTheta )
1953  {
1954  out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1955  out += ( cosFresnel + sinFresnel - 1. )*prof;
1956  }
1957  else
1958  {
1959  out = 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1960  }
1961 
1962  return out;
1963 }
1964 
1966 //
1967 // For the calculation of arg Gamma(z) one needs complex extension
1968 // of ln(Gamma(z))
1969 
1971 {
1972  const G4double cof[6] = { 76.18009172947146, -86.50532032941677,
1973  24.01409824083091, -1.231739572450155,
1974  0.1208650973866179e-2, -0.5395239384953e-5 } ;
1975  G4int j;
1976  G4complex z = zz - 1.0;
1977  G4complex tmp = z + 5.5;
1978  tmp -= (z + 0.5) * std::log(tmp);
1979  G4complex ser = G4complex(1.000000000190015,0.);
1980 
1981  for ( j = 0; j <= 5; j++ )
1982  {
1983  z += 1.0;
1984  ser += cof[j]/z;
1985  }
1986  return -tmp + std::log(2.5066282746310005*ser);
1987 }
1988 
1990 //
1991 // Bessel J0 function based on rational approximation from
1992 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
1993 
1995 {
1996  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
1997 
1998  modvalue = std::fabs(value);
1999 
2000  if ( value < 8.0 && value > -8.0 )
2001  {
2002  value2 = value*value;
2003 
2004  fact1 = 57568490574.0 + value2*(-13362590354.0
2005  + value2*( 651619640.7
2006  + value2*(-11214424.18
2007  + value2*( 77392.33017
2008  + value2*(-184.9052456 ) ) ) ) );
2009 
2010  fact2 = 57568490411.0 + value2*( 1029532985.0
2011  + value2*( 9494680.718
2012  + value2*(59272.64853
2013  + value2*(267.8532712
2014  + value2*1.0 ) ) ) );
2015 
2016  bessel = fact1/fact2;
2017  }
2018  else
2019  {
2020  arg = 8.0/modvalue;
2021 
2022  value2 = arg*arg;
2023 
2024  shift = modvalue-0.785398164;
2025 
2026  fact1 = 1.0 + value2*(-0.1098628627e-2
2027  + value2*(0.2734510407e-4
2028  + value2*(-0.2073370639e-5
2029  + value2*0.2093887211e-6 ) ) );
2030 
2031  fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
2032  + value2*(-0.6911147651e-5
2033  + value2*(0.7621095161e-6
2034  - value2*0.934945152e-7 ) ) );
2035 
2036  bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
2037  }
2038  return bessel;
2039 }
2040 
2042 //
2043 // Bessel J1 function based on rational approximation from
2044 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
2045 
2047 {
2048  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
2049 
2050  modvalue = std::fabs(value);
2051 
2052  if ( modvalue < 8.0 )
2053  {
2054  value2 = value*value;
2055 
2056  fact1 = value*(72362614232.0 + value2*(-7895059235.0
2057  + value2*( 242396853.1
2058  + value2*(-2972611.439
2059  + value2*( 15704.48260
2060  + value2*(-30.16036606 ) ) ) ) ) );
2061 
2062  fact2 = 144725228442.0 + value2*(2300535178.0
2063  + value2*(18583304.74
2064  + value2*(99447.43394
2065  + value2*(376.9991397
2066  + value2*1.0 ) ) ) );
2067  bessel = fact1/fact2;
2068  }
2069  else
2070  {
2071  arg = 8.0/modvalue;
2072 
2073  value2 = arg*arg;
2074 
2075  shift = modvalue - 2.356194491;
2076 
2077  fact1 = 1.0 + value2*( 0.183105e-2
2078  + value2*(-0.3516396496e-4
2079  + value2*(0.2457520174e-5
2080  + value2*(-0.240337019e-6 ) ) ) );
2081 
2082  fact2 = 0.04687499995 + value2*(-0.2002690873e-3
2083  + value2*( 0.8449199096e-5
2084  + value2*(-0.88228987e-6
2085  + value2*0.105787412e-6 ) ) );
2086 
2087  bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
2088 
2089  if (value < 0.0) bessel = -bessel;
2090  }
2091  return bessel;
2092 }
2093 
2094 //
2095 //
G4double G4ParticleHPJENDLHEData::G4double result
const G4double a0
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
ThreeVector shoot(const G4int Ap, const G4int Af)
Hep3Vector boostVector() const
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double CalculateNuclearRad(G4double A)
G4double Legendre96(T &typeT, F f, G4double a, G4double b)
G4double GetKineticEnergy() const
double x() const
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
static constexpr double millibarn
Definition: SystemOfUnits.h:86
static constexpr double hbarc
const char * p
Definition: xmltok.h:285
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:503
G4double SampleTableThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
void PutValue(size_t index, G4double energy, G4double dataValue)
G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position)
G4double GetDiffElasticProb(G4double theta)
G4complex AmplitudeNear(G4double theta)
G4double GetInvCoulombElasticXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
G4complex PhaseNear(G4double theta)
G4complex CoulombAmplitude(G4double theta)
G4ParticleDefinition * GetDefinition() const
tuple x
Definition: test.py:50
G4complex AmplitudeGla(G4double theta)
G4double GetDiffuseElasticSumXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4double GetCoulombElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
G4double GetLowEdgeEnergy(size_t binNumber) const
G4double GetDiffuseElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
double z() const
G4double SampleThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double A)
static constexpr double twopi
Definition: G4SIunits.hh:76
G4double AdaptiveGauss(T &typeT, F f, G4double a, G4double b, G4double e)
G4double Profile(G4double theta)
G4complex AmplitudeGG(G4double theta)
G4double GetTotalMomentum() const
static constexpr double TeV
Definition: G4SIunits.hh:218
tuple b
Definition: test.py:12
void SetMinEnergy(G4double anEnergy)
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4IonTable * GetIonTable() const
Hep3Vector vect() const
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
G4double GetHadronNucleonXscNS(G4ParticleDefinition *pParticle, G4double pTkin, G4ParticleDefinition *tParticle)
static size_t GetNumberOfElements()
Definition: G4Element.cc:405
const XML_Char int const XML_Char * value
Definition: expat.h:331
const G4ParticleDefinition * GetDefinition() const
static constexpr double degree
Definition: G4SIunits.hh:144
void TestAngleTable(const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
void InitParametersGla(const G4DynamicParticle *aParticle, G4double partMom, G4double Z, G4double A)
HepLorentzVector & boost(double, double, double)
G4double GetIntegrandFunction(G4double theta)
static G4Triton * Triton()
Definition: G4Triton.cc:95
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
G4double GetInvElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4double SampleT(const G4ParticleDefinition *aParticle, G4double p, G4double A)
G4double GetDiffElasticSumProb(G4double theta)
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
const G4int n
G4complex GammaMore(G4double theta)
const G4LorentzVector & Get4Momentum() const
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4LorentzVector Get4Momentum() const
double theta() const
G4double ThetaCMStoThetaLab(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaCMS)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static constexpr double GeV
void InitialiseOnFly(G4double Z, G4double A)
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
void InitDynParameters(const G4ParticleDefinition *theParticle, G4double partMom)
G4double GetRatioGen(G4double theta)
G4complex GammaLess(G4double theta)
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
G4double ProfileNear(G4double theta)
G4double GetFresnelIntegrandXsc(G4double alpha)
G4double SampleTableT(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
double y() const
G4double GetAtomicMassAmu(const G4String &symb) const
G4complex GetErfcInt(G4complex z)
static constexpr double GeV
Definition: G4SIunits.hh:217
tuple z
Definition: test.py:28
void insertAt(size_t, G4PhysicsVector *)
void SetMaxEnergy(const G4double anEnergy)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double CalcMandelstamS(const G4double mp, const G4double mt, const G4double Plab)
static constexpr double pi
Definition: G4SIunits.hh:75
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4double GetLegendrePol(G4int n, G4double x)
std::vector< G4Element * > G4ElementTable
double G4double
Definition: G4Types.hh:76
static constexpr double fermi
Definition: G4SIunits.hh:103
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:398
G4double GetPDGCharge() const
static const G4double alpha
G4complex AmplitudeSim(G4double theta)
static G4He3 * He3()
Definition: G4He3.cc:94
G4double GetInvElasticSumXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
double mag() const
static constexpr double keV
Definition: G4SIunits.hh:216
G4complex GetErfComp(G4complex z, G4int nMax)
G4double ThetaLabToThetaCMS(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaLab)
G4double GetDiffElasticSumProbA(G4double alpha)
G4double IntegralElasticProb(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
static constexpr double twopi
Definition: SystemOfUnits.h:55
double epsilon(double density, double temperature)
static constexpr double pi
Definition: SystemOfUnits.h:54
G4double SampleThetaLab(const G4HadProjectile *aParticle, G4double tmass, G4double A)
G4double GetTotalMomentum() const
G4complex GammaLogarithm(G4complex xx)
void InitParameters(const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)