ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RKFieldIntegrator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RKFieldIntegrator.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 // G4RKFieldIntegrator
27 #include "G4RKFieldIntegrator.hh"
28 #include "G4PhysicalConstants.hh"
29 #include "G4SystemOfUnits.hh"
30 #include "G4NucleiProperties.hh"
31 #include "G4FermiMomentum.hh"
32 #include "G4NuclearFermiDensity.hh"
34 #include "G4Nucleon.hh"
35 #include "G4Exp.hh"
36 #include "G4Log.hh"
37 #include "G4Pow.hh"
38 
39 // Class G4RKFieldIntegrator
40 //*************************************************************************************************************************************
41 
42 // only theActive are propagated, nothing else
43 // only theSpectators define the field, nothing else
44 
45 void G4RKFieldIntegrator::Transport(G4KineticTrackVector &theActive, const G4KineticTrackVector &theSpectators, G4double theTimeStep)
46 {
47  (void)theActive;
48  (void)theSpectators;
49  (void)theTimeStep;
50 }
51 
52 
54 {
55  const G4double Alpha = 0.25/fermi/fermi;
56  const G4double t1 = -7264.04*fermi*fermi*fermi;
57  const G4double tGamma = 87.65*fermi*fermi*fermi*fermi*fermi*fermi;
58 // const G4double Gamma = 1.676;
59  const G4double Vo = -0.498*fermi;
60  const G4double GammaY = 1.4*fermi;
61 
62  G4double Etot = 0;
63  G4int nBarion = Barions.size();
64  for(G4int c1 = 0; c1 < nBarion; c1++)
65  {
66  G4KineticTrack* p1 = Barions.operator[](c1);
67  // Ekin
68  Etot += p1->Get4Momentum().e();
69  for(G4int c2 = c1 + 1; c2 < nBarion; c2++)
70  {
71  G4KineticTrack* p2 = Barions.operator[](c2);
72  G4double r12 = (p1->GetPosition() - p2->GetPosition()).mag()*fermi;
73 
74  // Esk2
75  Etot += t1*G4Pow::GetInstance()->A23(Alpha/pi)*G4Exp(-Alpha*r12*r12);
76 
77  // Eyuk
78  Etot += Vo*0.5/r12*G4Exp(1/(4*Alpha*GammaY*GammaY))*
79  (G4Exp(-r12/GammaY)*(1 - Erf(0.5/GammaY/std::sqrt(Alpha) - std::sqrt(Alpha)*r12)) -
80  G4Exp( r12/GammaY)*(1 - Erf(0.5/GammaY/std::sqrt(Alpha) + std::sqrt(Alpha)*r12)));
81 
82  // Ecoul
83  Etot += 1.44*p1->GetDefinition()->GetPDGCharge()*p2->GetDefinition()->GetPDGCharge()/r12*Erf(std::sqrt(Alpha)*r12);
84 
85  // Epaul
86  Etot = 0;
87 
88  for(G4int c3 = c2 + 1; c3 < nBarion; c3++)
89  {
90  G4KineticTrack* p3 = Barions.operator[](c3);
91  G4double r13 = (p1->GetPosition() - p3->GetPosition()).mag()*fermi;
92 
93  // Esk3
94  Etot = tGamma*G4Pow::GetInstance()->powA(4*Alpha*Alpha/3/pi/pi, 1.5)*G4Exp(-Alpha*(r12*r12 + r13*r13));
95  }
96  }
97  }
98  return Etot;
99 }
100 
101 //************************************************************************************************
102 // originated from the Numerical recipes error function
104 {
105  const G4double Z1 = 1;
106  const G4double HF = Z1/2;
107  const G4double C1 = 0.56418958;
108 
109  const G4double P10 = +3.6767877;
110  const G4double Q10 = +3.2584593;
111  const G4double P11 = -9.7970465E-2;
112 
113 // static G4ThreadLocal G4double P2[5] = { 7.3738883, 6.8650185, 3.0317993, 0.56316962, 4.3187787e-5 };
114 // static G4ThreadLocal G4double Q2[5] = { 7.3739609, 15.184908, 12.79553, 5.3542168, 1. };
115  const G4double P2[5] = { 7.3738883, 6.8650185, 3.0317993, 0.56316962, 4.3187787e-5 };
116  const G4double Q2[5] = { 7.3739609, 15.184908, 12.79553, 5.3542168, 1. };
117 
118  const G4double P30 = -1.2436854E-1;
119  const G4double Q30 = +4.4091706E-1;
120  const G4double P31 = -9.6821036E-2;
121 
122  G4double V = std::abs(X);
123  G4double H;
124  G4double Y;
125  G4int c1;
126 
127  if(V < HF)
128  {
129  Y = V*V;
130  H = X*(P10 + P11*Y)/(Q10+Y);
131  }
132  else
133  {
134  if(V < 4)
135  {
136  G4double AP = P2[4];
137  G4double AQ = Q2[4];
138  for(c1 = 3; c1 >= 0; c1--)
139  {
140  AP = P2[c1] + V*AP;
141  AQ = Q2[c1] + V*AQ;
142  }
143  H = 1 - G4Exp(-V*V)*AP/AQ;
144  }
145  else
146  {
147  Y = 1./V*V;
148  H = 1 - G4Exp(-V*V)*(C1+Y*(P30 + P31*Y)/(Q30 + Y))/V;
149  }
150  if (X < 0)
151  H = -H;
152  }
153  return H;
154 }
155 
156 //************************************************************************************************
157 //This is a QMD version to calculate excitation energy of a fragment,
158 //which consists from G4KTV &the Particles
159 /*
160 G4double G4RKFieldIntegrator::GetExcitationEnergy(const G4KineticTrackVector &theParticles)
161 {
162  // Excitation energy of a fragment consisting from A nucleons and Z protons
163  // is Etot - Z*Mp - (A - Z)*Mn - B(A, Z), where B(A,Z) is the binding energy of fragment
164  // and Mp, Mn are proton and neutron mass, respectively.
165  G4int NZ = 0;
166  G4int NA = 0;
167  G4double Etot = CalculateTotalEnergy(theParticles);
168  for(G4int cParticle = 0; cParticle < theParticles.length(); cParticle++)
169  {
170  G4KineticTrack* pKineticTrack = theParticles.at(cParticle);
171  G4int Encoding = std::abs(pKineticTrack->GetDefinition()->GetPDGEncoding());
172  if (Encoding == 2212)
173  NZ++, NA++;
174  if (Encoding == 2112)
175  NA++;
176  Etot -= pKineticTrack->GetDefinition()->GetPDGMass();
177  }
178  return Etot - G4NucleiProperties::GetBindingEnergy(NZ, NA);
179 }
180 */
181 
182 //*************************************************************************************************************************************
183 //This is a simplified method to get excitation energy of a residual
184 // nucleus with nHitNucleons.
186 {
187  const G4double MeanE = 50;
188  G4double Sum = 0;
189  for(G4int c1 = 0; c1 < nHitNucleons; c1++)
190  {
191  Sum += -MeanE*G4Log(G4UniformRand());
192  }
193  return Sum;
194 }
195 //*************************************************************************************************************************************
196 
197 /*
198 //This is free propagation of particles for CASCADE mode. Target nucleons should be frozen
199 void G4RKFieldIntegrator::Integrate(G4KineticTrackVector& theParticles)
200  {
201  for(G4int cParticle = 0; cParticle < theParticles.length(); cParticle++)
202  {
203  G4KineticTrack* pKineticTrack = theParticles.at(cParticle);
204  pKineticTrack->SetPosition(pKineticTrack->GetPosition() + theTimeStep*pKineticTrack->Get4Momentum().boostVector());
205  }
206  }
207 */
208 //*************************************************************************************************************************************
209 
211 {
212  for(size_t cParticle = 0; cParticle < theBarions.size(); cParticle++)
213  {
214  G4KineticTrack* pKineticTrack = theBarions[cParticle];
215  pKineticTrack->SetPosition(pKineticTrack->GetPosition() + theTimeStep*pKineticTrack->Get4Momentum().boostVector());
216  }
217 }
218 
219 //*************************************************************************************************************************************
220 
221 // constant to calculate theCoulomb barrier
222 const G4double G4RKFieldIntegrator::coulomb = 1.44 / 1.14 * MeV;
223 
224 // kaon's potential constant (real part only)
225 // 0.35 + i0.82 or 0.63 + i0.89 fermi
227 
228 // pion's potential constant (real part only)
230 // 0.35 + i0.82 or 0.63 + i0.89 fermi
232 
233 // antiproton's potential constant (real part only)
234 // 1.53 + i2.50 fermi
236 
237 // methods for calculating potentials for different types of particles
238 // aPosition is relative to the nucleus center
240 {
241  /*
242  const G4double Mn = 939.56563 * MeV; // mass of nuetron
243 
244  G4VNuclearDensity *theDencity;
245  if(theA < 17) theDencity = new G4NuclearShellModelDensity(theA, theZ);
246  else theDencity = new G4NuclearFermiDensity(theA, theZ);
247 
248  // GetDencity() accepts only G4ThreeVector so build it:
249  G4ThreeVector aPosition(0.0, 0.0, radius);
250  G4double density = theDencity->GetDensity(aPosition);
251  delete theDencity;
252 
253  G4FermiMomentum *fm = new G4FermiMomentum();
254  fm->Init(theA, theZ);
255  G4double fermiMomentum = fm->GetFermiMomentum(density);
256  delete fm;
257 
258  return sqr(fermiMomentum)/(2 * Mn)
259  + G4CreateNucleus::GetBindingEnergy(theZ, theA)/theA;
260  //+ G4NucleiProperties::GetBindingEnergy(theZ, theA)/theA;
261  */
262 
263  return 0.0;
264 }
265 
267 {
268  /*
269  // calculate Coulomb barrier value
270  G4double theCoulombBarrier = coulomb * theZ/(1. + G4Pow::GetInstance()->Z13(theA));
271  const G4double Mp = 938.27231 * MeV; // mass of proton
272 
273  G4VNuclearDensity *theDencity;
274  if(theA < 17) theDencity = new G4NuclearShellModelDensity(theA, theZ);
275  else theDencity = new G4NuclearFermiDensity(theA, theZ);
276 
277  // GetDencity() accepts only G4ThreeVector so build it:
278  G4ThreeVector aPosition(0.0, 0.0, radius);
279  G4double density = theDencity->GetDensity(aPosition);
280  delete theDencity;
281 
282  G4FermiMomentum *fm = new G4FermiMomentum();
283  fm->Init(theA, theZ);
284  G4double fermiMomentum = fm->GetFermiMomentum(density);
285  delete fm;
286 
287  return sqr(fermiMomentum)/ (2 * Mp)
288  + G4CreateNucleus::GetBindingEnergy(theZ, theA)/theA;
289  //+ G4NucleiProperties::GetBindingEnergy(theZ, theA)/theA
290  + theCoulombBarrier;
291  */
292 
293  return 0.0;
294 }
295 
297 {
298  /*
299  //G4double theM = G4NucleiProperties::GetAtomicMass(theA, theZ);
300  G4double theM = theZ * G4Proton::Proton()->GetPDGMass()
301  + (theA - theZ) * G4Neutron::Neutron()->GetPDGMass()
302  + G4CreateNucleus::GetBindingEnergy(theZ, theA);
303 
304  const G4double Mp = 938.27231 * MeV; // mass of proton
305  G4double mu = (theM * Mp)/(theM + Mp);
306 
307  // antiproton's potential coefficient
308  // V = coeff_antiproton * nucleus_density
309  G4double coeff_antiproton = -2.*pi/mu * (1. + Mp) * a_antiproton;
310 
311  G4VNuclearDensity *theDencity;
312  if(theA < 17) theDencity = new G4NuclearShellModelDensity(theA, theZ);
313  else theDencity = new G4NuclearFermiDensity(theA, theZ);
314 
315  // GetDencity() accepts only G4ThreeVector so build it:
316  G4ThreeVector aPosition(0.0, 0.0, radius);
317  G4double density = theDencity->GetDensity(aPosition);
318  delete theDencity;
319 
320  return coeff_antiproton * density;
321  */
322 
323  return 0.0;
324 }
325 
327 {
328  /*
329  //G4double theM = G4NucleiProperties::GetAtomicMass(theA, theZ);
330  G4double theM = theZ * G4Proton::Proton()->GetPDGMass()
331  + (theA - theZ) * G4Neutron::Neutron()->GetPDGMass()
332  + G4CreateNucleus::GetBindingEnergy(theZ, theA);
333 
334  const G4double Mk = 496. * MeV; // mass of "kaon"
335  G4double mu = (theM * Mk)/(theM + Mk);
336 
337  // kaon's potential coefficient
338  // V = coeff_kaon * nucleus_density
339  G4double coeff_kaon = -2.*pi/mu * (1. + Mk/theM) * a_kaon;
340 
341  G4VNuclearDensity *theDencity;
342  if(theA < 17) theDencity = new G4NuclearShellModelDensity(theA, theZ);
343  else theDencity = new G4NuclearFermiDensity(theA, theZ);
344 
345  // GetDencity() accepts only G4ThreeVector so build it:
346  G4ThreeVector aPosition(0.0, 0.0, radius);
347  G4double density = theDencity->GetDensity(aPosition);
348  delete theDencity;
349 
350  return coeff_kaon * density;
351  */
352 
353  return 0.0;
354 }
355 
357 {
358  /*
359  //G4double theM = G4NucleiProperties::GetAtomicMass(theA, theZ);
360  G4double theM = theZ * G4Proton::Proton()->GetPDGMass()
361  + (theA - theZ) * G4Neutron::Neutron()->GetPDGMass()
362  + G4CreateNucleus::GetBindingEnergy(theZ, theA);
363 
364  const G4double Mpi = 139. * MeV; // mass of "pion"
365  G4double mu = (theM * Mpi)/(theM + Mpi);
366 
367  // pion's potential coefficient
368  // V = coeff_pion * nucleus_density
369  G4double coeff_pion = -2.*pi/mu * (1. + Mpi) * a_pion;
370 
371  G4VNuclearDensity *theDencity;
372  if(theA < 17) theDencity = new G4NuclearShellModelDensity(theA, theZ);
373  else theDencity = new G4NuclearFermiDensity(theA, theZ);
374 
375  // GetDencity() accepts only G4ThreeVector so build it:
376  G4ThreeVector aPosition(0.0, 0.0, radius);
377  G4double density = theDencity->GetDensity(aPosition);
378  delete theDencity;
379 
380  return coeff_pion * density;
381  */
382 
383  return 0.0;
384 }