ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DiffuseElasticV2.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DiffuseElasticV2.cc
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 G4DiffuseElasticV2
29 //
30 //
31 // G4 Model: optical diffuse elastic scattering with 4-momentum balance
32 //
33 // 24-May-07 V. Grichine
34 //
35 // 21.10.15 V. Grichine
36 // Bug fixed in BuildAngleTable, improving accuracy for
37 // angle bins at high energies > 50 GeV for pions.
38 //
39 // 24.11.17 W. Pokorski, code cleanup and performance improvements
40 //
41 
42 #include "G4DiffuseElasticV2.hh"
43 #include "G4ParticleTable.hh"
44 #include "G4ParticleDefinition.hh"
45 #include "G4IonTable.hh"
46 #include "G4NucleiProperties.hh"
47 
48 #include "Randomize.hh"
49 #include "G4Integrator.hh"
50 #include "globals.hh"
51 #include "G4PhysicalConstants.hh"
52 #include "G4SystemOfUnits.hh"
53 
54 #include "G4Proton.hh"
55 #include "G4Neutron.hh"
56 #include "G4Deuteron.hh"
57 #include "G4Alpha.hh"
58 #include "G4PionPlus.hh"
59 #include "G4PionMinus.hh"
60 
61 #include "G4Element.hh"
62 #include "G4ElementTable.hh"
63 #include "G4NistManager.hh"
64 #include "G4PhysicsTable.hh"
65 #include "G4PhysicsLogVector.hh"
66 #include "G4PhysicsFreeVector.hh"
67 
68 #include "G4Exp.hh"
69 
70 #include "G4HadronicParameters.hh"
71 
73 //
74 
75 
77  : G4HadronElastic("DiffuseElasticV2"), fParticle(0)
78 {
79  SetMinEnergy( 0.01*MeV );
81 
82  verboseLevel = 0;
83  lowEnergyRecoilLimit = 100.*keV;
84  lowEnergyLimitQ = 0.0*GeV;
85  lowEnergyLimitHE = 0.0*GeV;
86  lowestEnergyLimit = 0.0*keV;
87  plabLowLimit = 20.0*MeV;
88 
91 
92  fEnergyBin = 300; // Increased from the original 200 to have no wider log-energy-bins up to 10 PeV
93  fAngleBin = 200;
94 
96 
98  fEnergySumVector = 0;
99 
100  fParticle = 0;
101  fWaveVector = 0.;
102  fAtomicWeight = 0.;
103  fAtomicNumber = 0.;
104  fNuclearRadius = 0.;
105  fBeta = 0.;
106  fZommerfeld = 0.;
107  fAm = 0.;
108  fAddCoulomb = false;
109 }
110 
112 //
113 // Destructor
114 
116 {
117  if ( fEnergyVector )
118  {
119  delete fEnergyVector;
120  fEnergyVector = 0;
121  }
122 }
123 
125 //
126 // Initialisation for given particle using element table of application
127 
129 {
130 
131  const G4ElementTable* theElementTable = G4Element::GetElementTable();
132 
133  size_t jEl, numOfEl = G4Element::GetNumberOfElements();
134 
135  for( jEl = 0; jEl < numOfEl; ++jEl) // application element loop
136  {
137  fAtomicNumber = (*theElementTable)[jEl]->GetZ(); // atomic number
140 
141  if( verboseLevel > 0 )
142  {
143  G4cout<<"G4DiffuseElasticV2::Initialise() the element: "
144  <<(*theElementTable)[jEl]->GetName()<<G4endl;
145  }
147  fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
148 
149  BuildAngleTable();
150 
153 
154  }
155  return;
156 }
157 
159 //
160 // return differential elastic probability d(probability)/d(t) with
161 // Coulomb correction. It is called from BuildAngleTable()
162 
163 G4double
165 {
166 
167  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
168  G4double delta, diffuse, gamma;
169  G4double e1, e2, bone, bone2;
170 
171  // G4double wavek = momentum/hbarc; // wave vector
172  // G4double r0 = 1.08*fermi;
173  // G4double rad = r0*G4Pow::GetInstance()->A13(A);
174 
175  G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
176  G4double kr2 = kr*kr;
177  G4double krt = kr*theta;
178 
179  bzero = BesselJzero(krt);
180  bzero2 = bzero*bzero;
181  bone = BesselJone(krt);
182  bone2 = bone*bone;
183  bonebyarg = BesselOneByArg(krt);
184  bonebyarg2 = bonebyarg*bonebyarg;
185 
186  if ( fParticle == theProton )
187  {
188  diffuse = 0.63*fermi;
189  gamma = 0.3*fermi;
190  delta = 0.1*fermi*fermi;
191  e1 = 0.3*fermi;
192  e2 = 0.35*fermi;
193  }
194  else if ( fParticle == theNeutron )
195  {
196  diffuse = 0.63*fermi;
197  gamma = 0.3*fermi;
198  delta = 0.1*fermi*fermi;
199  e1 = 0.3*fermi;
200  e2 = 0.35*fermi;
201  }
202  else // as proton, if were not defined
203  {
204  diffuse = 0.63*fermi;
205  gamma = 0.3*fermi;
206  delta = 0.1*fermi*fermi;
207  e1 = 0.3*fermi;
208  e2 = 0.35*fermi;
209  }
210 
211  G4double lambda = 15; // 15 ok
212  // G4double kgamma = fWaveVector*gamma; // wavek*delta;
213  G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
214 
215  if( fAddCoulomb ) // add Coulomb correction
216  {
217  G4double sinHalfTheta = std::sin(0.5*theta);
218  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
219 
220  kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
221  }
222 
223  G4double kgamma2 = kgamma*kgamma;
224 
225  // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
226  // G4double dk2t2 = dk2t*dk2t;
227 
228  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
229  G4double pikdt = lambda*(1. - G4Exp( -pi*fWaveVector*diffuse*theta/lambda ) ); // wavek*delta;
230 
231  damp = DampFactor( pikdt );
232  damp2 = damp*damp;
233 
234  G4double mode2k2 = ( e1*e1 + e2*e2 )*fWaveVector*fWaveVector;
235  G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
236 
237  sigma = kgamma2;
238  // sigma += dk2t2;
239  sigma *= bzero2;
240  sigma += mode2k2*bone2;
241  sigma += e2dk3t*bzero*bone;
242 
243  // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
244  sigma += kr2*bonebyarg2; // correction at J1()/()
245 
246  sigma *= damp2; // *rad*rad;
247 
248  return sigma;
249 }
250 
251 
253 //
254 // return differential elastic probability 2*pi*sin(theta)*d(probability)/d(omega)
255 
256 G4double
258 {
259  G4double result;
260 
261  result = GetDiffElasticSumProbA(alpha) * 2 * CLHEP::pi * std::sin(alpha);
262 
263  return result;
264 }
265 
266 
270 //
271 // Return inv momentum transfer -t > 0 from initialisation table
272 
274  G4int Z, G4int A)
275 {
276  fParticle = aParticle;
277  G4double m1 = fParticle->GetPDGMass(), t;
278  G4double totElab = std::sqrt(m1*m1+p*p);
280  G4LorentzVector lv1(p,0.0,0.0,totElab);
281  G4LorentzVector lv(0.0,0.0,0.0,mass2);
282  lv += lv1;
283 
284  G4ThreeVector bst = lv.boostVector();
285  lv1.boost(-bst);
286 
287  G4ThreeVector p1 = lv1.vect();
288  G4double momentumCMS = p1.mag();
289 
290  if( aParticle == theNeutron)
291  {
292  G4double Tmax = NeutronTuniform( Z );
293  G4double pCMS2 = momentumCMS*momentumCMS;
294  G4double Tkin = std::sqrt(pCMS2+m1*m1)-m1;
295 
296  if( Tkin <= Tmax )
297  {
298  t = 4.*pCMS2*G4UniformRand();
299  return t;
300  }
301  }
302 
303  t = SampleTableT( aParticle, momentumCMS, G4double(Z), G4double(A) ); // sample theta in cms
304 
305  return t;
306 }
307 
309 
311 {
312  G4double elZ = G4double(Z);
313  elZ -= 1.;
314  G4double Tkin = 12.*G4Exp(-elZ/10.) + 1.;
315 
316  return Tkin;
317 }
318 
319 
321 //
322 // Return inv momentum transfer -t > 0 from initialisation table
323 
325  G4double Z, G4double A)
326 {
327  G4double alpha = SampleTableThetaCMS( aParticle, p, Z, A); // sample theta in cms
328  G4double t = 2*p*p*( 1 - std::cos(alpha) ); // -t !!!
329 
330  return t;
331 }
332 
334 //
335 // Return scattering angle2 sampled in cms according to precalculated table.
336 
337 
338 G4double
341 {
342  size_t iElement;
343  G4int iMomentum;
344  unsigned long iAngle = 0;
345  G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W;
346  G4double m1 = particle->GetPDGMass();
347 
348  for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
349  {
350  if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
351  }
352 
353  if ( iElement == fElementNumberVector.size() )
354  {
355  InitialiseOnFly(Z,A); // table preparation, if needed
356  }
357 
360 
361 
362  G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
363 
364  iMomentum = fEnergyVector->FindBin(kinE,1000) + 1;
365 
366  position = (*(*fEnergySumVector)[iMomentum])[0]*G4UniformRand();
367 
368  for(iAngle = 0; iAngle < fAngleBin; iAngle++)
369  {
370  if (position > (*(*fEnergySumVector)[iMomentum])[iAngle]) break;
371  }
372 
373 
374  if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges
375  {
376  randAngle = GetScatteringAngle(iMomentum, iAngle, position);
377  }
378  else // kinE inside between energy table edges
379  {
380  theta2 = GetScatteringAngle(iMomentum, iAngle, position);
381 
382  E2 = fEnergyVector->Energy(iMomentum);
383 
384  iMomentum--;
385 
386  theta1 = GetScatteringAngle(iMomentum, iAngle, position);
387 
388  E1 = fEnergyVector->Energy(iMomentum);
389 
390  W = 1.0/(E2 - E1);
391  W1 = (E2 - kinE)*W;
392  W2 = (kinE - E1)*W;
393 
394  randAngle = W1*theta1 + W2*theta2;
395  }
396 
397 
398 
399  if(randAngle < 0.) randAngle = 0.;
400 
401  return randAngle;
402 }
403 
405 //
406 // Initialisation for given particle on fly using new element number
407 
409 {
410  fAtomicNumber = Z; // atomic number
411  fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( Z ) );
412 
414 
415  if( verboseLevel > 0 )
416  {
417  G4cout<<"G4DiffuseElasticV2::InitialiseOnFly() the element with Z = "
418  <<Z<<"; and A = "<<A<<G4endl;
419  }
421 
422  BuildAngleTable();
423 
426 
427  return;
428 }
429 
431 //
432 // Build for given particle and element table of momentum, angle probability.
433 // For the moment in lab system.
434 
436 {
437  G4int i, j;
438  G4double partMom, kinE, a = 0., z = fParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
439  G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
440 
442 
443  fEnergyAngleVector = new std::vector<std::vector<double>*>;
444  fEnergySumVector = new std::vector<std::vector<double>*>;
445 
446  for( i = 0; i < fEnergyBin; i++)
447  {
448  kinE = fEnergyVector->Energy(i);
449  partMom = std::sqrt( kinE*(kinE + 2*m1) );
450 
451  fWaveVector = partMom/hbarc;
452 
454  G4double kRmax = 18.6; // 10.6; 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
455  G4double kRcoul = 1.9; // 1.2; 1.4, 2.5; // on the first slope of J1
456 
457  alphaMax = kRmax/kR;
458 
459  if ( alphaMax >= CLHEP::pi ) alphaMax = CLHEP::pi; // vmg21.10.15
460 
461  alphaCoulomb = kRcoul/kR;
462 
463  if( z )
464  {
465  a = partMom/m1; // beta*gamma for m1
466  fBeta = a/std::sqrt(1+a*a);
469  fAddCoulomb = true;
470  }
471 
472  std::vector<double>* angleVector = new std::vector<double>(fAngleBin);
473  std::vector<double>* sumVector = new std::vector<double>(fAngleBin);
474 
475 
476  G4double delth = alphaMax/fAngleBin;
477 
478  sum = 0.;
479 
480  for(j = fAngleBin-1; j >= 0; j--)
481  {
482  alpha1 = delth*j;
483  alpha2 = alpha1 + delth;
484 
485  if( fAddCoulomb && ( alpha2 < alphaCoulomb)) fAddCoulomb = false;
486 
487  delta = integral.Legendre10(this, &G4DiffuseElasticV2::GetIntegrandFunction, alpha1, alpha2);
488 
489  sum += delta;
490 
491  (*angleVector)[j] = alpha1;
492  (*sumVector)[j] = sum;
493 
494  }
495  fEnergyAngleVector->push_back(angleVector);
496  fEnergySumVector->push_back(sumVector);
497 
498  }
499  return;
500 }
501 
503 //
504 //
505 
506 G4double
508 {
509  G4double x1, x2, y1, y2, randAngle = 0;
510 
511  if( iAngle == 0 )
512  {
513  randAngle = (*(*fEnergyAngleVector)[iMomentum])[iAngle];
514  }
515  else
516  {
517  if ( iAngle >= (*fEnergyAngleVector)[iMomentum]->size() )
518  {
519  iAngle = (*fEnergyAngleVector)[iMomentum]->size() - 1;
520  }
521 
522  y1 = (*(*fEnergySumVector)[iMomentum])[iAngle-1];
523  y2 = (*(*fEnergySumVector)[iMomentum])[iAngle];
524 
525  x1 = (*(*fEnergyAngleVector)[iMomentum])[iAngle-1];
526  x2 = (*(*fEnergyAngleVector)[iMomentum])[iAngle];
527 
528  if ( x1 == x2 ) randAngle = x2;
529  else
530  {
531  if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand();
532  else
533  {
534  randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
535  }
536  }
537  }
538 
539  return randAngle;
540 }
541 
542 
543 
544 
546 //
547 // Return scattering angle in lab system (target at rest) knowing theta in CMS
548 
549 
550 
551 G4double
553  G4double tmass, G4double thetaCMS)
554 {
555  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
556  G4double m1 = theParticle->GetPDGMass();
557  G4LorentzVector lv1 = aParticle->Get4Momentum();
558  G4LorentzVector lv(0.0,0.0,0.0,tmass);
559 
560  lv += lv1;
561 
562  G4ThreeVector bst = lv.boostVector();
563 
564  lv1.boost(-bst);
565 
566  G4ThreeVector p1 = lv1.vect();
567  G4double ptot = p1.mag();
568 
570  G4double cost = std::cos(thetaCMS);
571  G4double sint;
572 
573  if( cost >= 1.0 )
574  {
575  cost = 1.0;
576  sint = 0.0;
577  }
578  else if( cost <= -1.0)
579  {
580  cost = -1.0;
581  sint = 0.0;
582  }
583  else
584  {
585  sint = std::sqrt((1.0-cost)*(1.0+cost));
586  }
587  if (verboseLevel>1)
588  {
589  G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
590  }
591  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
592  v1 *= ptot;
593  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
594 
595  nlv1.boost(bst);
596 
597  G4ThreeVector np1 = nlv1.vect();
598 
599  G4double thetaLab = np1.theta();
600 
601  return thetaLab;
602 }
604 //
605 // Return scattering angle in CMS system (target at rest) knowing theta in Lab
606 
607 
608 
609 G4double
611  G4double tmass, G4double thetaLab)
612 {
613  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
614  G4double m1 = theParticle->GetPDGMass();
615  G4double plab = aParticle->GetTotalMomentum();
616  G4LorentzVector lv1 = aParticle->Get4Momentum();
617  G4LorentzVector lv(0.0,0.0,0.0,tmass);
618 
619  lv += lv1;
620 
621  G4ThreeVector bst = lv.boostVector();
622 
624  G4double cost = std::cos(thetaLab);
625  G4double sint;
626 
627  if( cost >= 1.0 )
628  {
629  cost = 1.0;
630  sint = 0.0;
631  }
632  else if( cost <= -1.0)
633  {
634  cost = -1.0;
635  sint = 0.0;
636  }
637  else
638  {
639  sint = std::sqrt((1.0-cost)*(1.0+cost));
640  }
641  if (verboseLevel>1)
642  {
643  G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
644  }
645  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
646  v1 *= plab;
647  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
648 
649  nlv1.boost(-bst);
650 
651  G4ThreeVector np1 = nlv1.vect();
652  G4double thetaCMS = np1.theta();
653 
654  return thetaCMS;
655 }
656