source: trunk/source/processes/hadronic/models/binary_cascade/src/G4RKFieldIntegrator.cc @ 1347

Last change on this file since 1347 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 13.5 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26// G4RKFieldIntegrator
27#include "G4RKFieldIntegrator.hh"
28#include "G4NucleiProperties.hh"
29#include "G4FermiMomentum.hh"
30#include "G4NuclearFermiDensity.hh"
31#include "G4NuclearShellModelDensity.hh"
32#include "G4Nucleon.hh"
33
34// Class G4RKFieldIntegrator
35//*************************************************************************************************************************************
36
37// only theActive are propagated, nothing else
38// only theSpectators define the field, nothing else
39
40void G4RKFieldIntegrator::Transport(G4KineticTrackVector &theActive, const G4KineticTrackVector &theSpectators, G4double theTimeStep)
41{
42   (void)theActive;
43   (void)theSpectators;
44   (void)theTimeStep;
45}
46
47
48G4double G4RKFieldIntegrator::CalculateTotalEnergy(const G4KineticTrackVector& Barions)
49{
50   const G4double Alpha  =  0.25/fermi/fermi;
51   const G4double t1     = -7264.04*fermi*fermi*fermi;
52   const G4double tGamma =  87.65*fermi*fermi*fermi*fermi*fermi*fermi;
53//   const G4double Gamma  =  1.676;
54   const G4double Vo     = -0.498*fermi;
55   const G4double GammaY =  1.4*fermi;
56
57   G4double Etot = 0;
58   G4int nBarion = Barions.size();
59   for(G4int c1 = 0; c1 < nBarion; c1++)
60      {
61      G4KineticTrack* p1 = Barions.operator[](c1);
62   // Ekin
63      Etot += p1->Get4Momentum().e();
64      for(G4int c2 = c1 + 1; c2 < nBarion; c2++)
65         {
66         G4KineticTrack* p2 = Barions.operator[](c2);
67         G4ThreeVector   rv = p1->GetPosition() - p2->GetPosition();
68         G4double r12 = std::sqrt(rv*rv)*fermi;
69
70         //  Esk2
71         Etot += t1*std::pow(Alpha/pi, 3/2)*std::exp(-Alpha*r12*r12); 
72
73         // Eyuk
74         Etot += Vo*0.5/r12*std::exp(1/(4*Alpha*GammaY*GammaY))*
75            (std::exp(-r12/GammaY)*(1 - Erf(0.5/GammaY/std::sqrt(Alpha) - std::sqrt(Alpha)*r12)) - 
76             std::exp( r12/GammaY)*(1 - Erf(0.5/GammaY/std::sqrt(Alpha) + std::sqrt(Alpha)*r12)));
77
78         // Ecoul       
79         Etot += 1.44*p1->GetDefinition()->GetPDGCharge()*p2->GetDefinition()->GetPDGCharge()/r12*Erf(std::sqrt(Alpha)*r12);
80
81         // Epaul
82         Etot = 0;
83
84         for(G4int c3 = c2 + 1; c3 < nBarion; c3++)
85            {
86            G4KineticTrack* p3 = Barions.operator[](c3);
87            G4ThreeVector rv = p1->GetPosition() - p3->GetPosition();
88            G4double r13 = std::sqrt(rv*rv)*fermi;
89
90            // Esk3
91            Etot  = tGamma*std::pow(4*Alpha*Alpha/3/pi/pi, 1.5)*std::exp(-Alpha*(r12*r12 + r13*r13));
92            }
93         }
94      }     
95   return Etot;
96} 
97
98//************************************************************************************************       
99// originated from the Numerical recipes error function
100G4double G4RKFieldIntegrator::Erf(G4double X)
101{
102   const G4double Z1 = 1;
103   const G4double HF = Z1/2;
104   const G4double C1 = 0.56418958;
105
106   const G4double P10 = +3.6767877;
107   const G4double Q10 = +3.2584593;
108   const G4double P11 = -9.7970465E-2;
109
110   static G4double P2[5] = { 7.3738883, 6.8650185,  3.0317993, 0.56316962, 4.3187787e-5 };
111   static G4double Q2[5] = { 7.3739609, 15.184908, 12.79553,   5.3542168,  1. };
112   
113   const G4double P30 = -1.2436854E-1;
114   const G4double Q30 = +4.4091706E-1;
115   const G4double P31 = -9.6821036E-2;
116
117   G4double V = std::abs(X);
118   G4double H;   
119   G4double Y;
120   G4int c1;
121     
122   if(V < HF)
123      {
124      Y = V*V;
125      H = X*(P10 + P11*Y)/(Q10+Y);
126      }
127   else
128      {
129      if(V < 4) 
130         {
131         G4double AP = P2[4];
132         G4double AQ = Q2[4];
133         for(c1 = 3; c1 >= 0; c1--)
134            {
135            AP = P2[c1] + V*AP;
136            AQ = Q2[c1] + V*AQ;
137            }
138         H = 1 - std::exp(-V*V)*AP/AQ;
139         }
140      else
141        {
142        Y = 1./V*V;
143        H = 1 - std::exp(-V*V)*(C1+Y*(P30 + P31*Y)/(Q30 + Y))/V;
144        }
145     if (X < 0) 
146        H =- H;
147     }
148   return H;
149}
150   
151//************************************************************************************************       
152//This is a QMD version to calculate excitation energy of a fragment,
153//which consists from G4KTV &the Particles
154/*
155G4double G4RKFieldIntegrator::GetExcitationEnergy(const G4KineticTrackVector &theParticles)
156{
157   // Excitation energy of a fragment consisting from A nucleons and Z protons
158   // is Etot - Z*Mp - (A - Z)*Mn - B(A, Z), where B(A,Z) is the binding energy of fragment
159   //  and Mp, Mn are proton and neutron mass, respectively.
160   G4int NZ = 0;
161   G4int NA = 0;
162   G4double Etot = CalculateTotalEnergy(theParticles);
163   for(G4int cParticle = 0; cParticle < theParticles.length(); cParticle++)
164      {
165      G4KineticTrack* pKineticTrack = theParticles.at(cParticle);
166      G4int Encoding =  std::abs(pKineticTrack->GetDefinition()->GetPDGEncoding());
167      if (Encoding == 2212)
168          NZ++, NA++;
169      if (Encoding == 2112)
170          NA++;
171      Etot -= pKineticTrack->GetDefinition()->GetPDGMass();
172      }
173   return Etot - G4NucleiProperties::GetBindingEnergy(NZ, NA);
174}
175*/
176
177//*************************************************************************************************************************************
178//This is a simplified method to get excitation energy of a residual
179// nucleus with nHitNucleons.
180G4double G4RKFieldIntegrator::GetExcitationEnergy(G4int nHitNucleons, const G4KineticTrackVector &)
181{
182   const G4double MeanE = 50;
183   G4double Sum = 0;
184   for(G4int c1 = 0; c1 < nHitNucleons; c1++)
185       {
186       Sum += -MeanE*std::log(G4UniformRand());
187       }
188   return Sum;
189}
190//*************************************************************************************************************************************
191
192/*
193//This is free propagation of particles for CASCADE mode. Target nucleons should be frozen
194void G4RKFieldIntegrator::Integrate(G4KineticTrackVector& theParticles)
195   {
196   for(G4int cParticle = 0; cParticle < theParticles.length(); cParticle++)
197      {
198      G4KineticTrack* pKineticTrack = theParticles.at(cParticle);
199      pKineticTrack->SetPosition(pKineticTrack->GetPosition() + theTimeStep*pKineticTrack->Get4Momentum().boostVector());
200      }
201   }
202*/   
203//*************************************************************************************************************************************
204
205void G4RKFieldIntegrator::Integrate(const G4KineticTrackVector& theBarions, G4double theTimeStep)
206{
207   for(size_t cParticle = 0; cParticle < theBarions.size(); cParticle++)
208      {
209      G4KineticTrack* pKineticTrack = theBarions[cParticle];
210      pKineticTrack->SetPosition(pKineticTrack->GetPosition() + theTimeStep*pKineticTrack->Get4Momentum().boostVector());
211      } 
212}
213   
214//*************************************************************************************************************************************
215
216// constant to calculate theCoulomb barrier
217const G4double G4RKFieldIntegrator::coulomb = 1.44 / 1.14 * MeV;
218
219// kaon's potential constant (real part only)
220// 0.35 + i0.82 or 0.63 + i0.89 fermi
221const G4double G4RKFieldIntegrator::a_kaon = 0.35;
222
223// pion's potential constant (real part only)
224//!! for pions it has todiffer from kaons
225// 0.35 + i0.82 or 0.63 + i0.89 fermi
226const G4double G4RKFieldIntegrator::a_pion = 0.35;
227
228// antiproton's potential constant (real part only)
229// 1.53 + i2.50 fermi
230const G4double G4RKFieldIntegrator::a_antiproton = 1.53;
231
232// methods for calculating potentials for different types of particles
233// aPosition is relative to the nucleus center
234G4double G4RKFieldIntegrator::GetNeutronPotential(G4double )
235{
236   /*
237   const G4double Mn  = 939.56563 * MeV; // mass of nuetron
238
239   G4VNuclearDensity *theDencity;
240   if(theA < 17) theDencity = new G4NuclearShellModelDensity(theA, theZ);
241   else          theDencity = new G4NuclearFermiDensity(theA, theZ);
242   
243   // GetDencity() accepts only G4ThreeVector so build it:   
244   G4ThreeVector aPosition(0.0, 0.0, radius);
245   G4double density = theDencity->GetDensity(aPosition);
246   delete theDencity;
247   
248   G4FermiMomentum *fm = new G4FermiMomentum();
249   fm->Init(theA, theZ);
250   G4double fermiMomentum = fm->GetFermiMomentum(density);
251   delete fm;
252   
253   return sqr(fermiMomentum)/(2 * Mn)
254      + G4CreateNucleus::GetBindingEnergy(theZ, theA)/theA;
255      //+ G4NucleiProperties::GetBindingEnergy(theZ, theA)/theA;
256   */
257   
258   return 0.0;
259}
260
261G4double G4RKFieldIntegrator::GetProtonPotential(G4double )
262{
263   /*
264   // calculate Coulomb barrier value
265   G4double theCoulombBarrier = coulomb * theZ/(1. + std::pow(theA, 1./3.));
266   const G4double Mp  = 938.27231 * MeV; // mass of proton
267
268   G4VNuclearDensity *theDencity;
269   if(theA < 17) theDencity = new G4NuclearShellModelDensity(theA, theZ);
270   else          theDencity = new G4NuclearFermiDensity(theA, theZ);
271   
272   // GetDencity() accepts only G4ThreeVector so build it:   
273   G4ThreeVector aPosition(0.0, 0.0, radius);
274   G4double density = theDencity->GetDensity(aPosition);
275   delete theDencity;
276   
277   G4FermiMomentum *fm = new G4FermiMomentum();
278   fm->Init(theA, theZ);
279   G4double fermiMomentum = fm->GetFermiMomentum(density);
280   delete fm;
281
282   return sqr(fermiMomentum)/ (2 * Mp)
283      + G4CreateNucleus::GetBindingEnergy(theZ, theA)/theA;
284      //+ G4NucleiProperties::GetBindingEnergy(theZ, theA)/theA
285      + theCoulombBarrier;
286   */
287
288   return 0.0;
289}
290
291G4double G4RKFieldIntegrator::GetAntiprotonPotential(G4double )
292{
293   /*
294   //G4double theM = G4NucleiProperties::GetAtomicMass(theA, theZ);
295   G4double theM = theZ * G4Proton::Proton()->GetPDGMass()
296      + (theA - theZ) * G4Neutron::Neutron()->GetPDGMass()
297      + G4CreateNucleus::GetBindingEnergy(theZ, theA);
298     
299   const G4double Mp  = 938.27231 * MeV; // mass of proton
300   G4double mu = (theM * Mp)/(theM + Mp);
301   
302   // antiproton's potential coefficient
303   //   V = coeff_antiproton * nucleus_density
304   G4double coeff_antiproton = -2.*pi/mu * (1. + Mp) * a_antiproton;
305   
306   G4VNuclearDensity *theDencity;
307   if(theA < 17) theDencity = new G4NuclearShellModelDensity(theA, theZ);
308   else          theDencity = new G4NuclearFermiDensity(theA, theZ);
309   
310   // GetDencity() accepts only G4ThreeVector so build it:   
311   G4ThreeVector aPosition(0.0, 0.0, radius);
312   G4double density = theDencity->GetDensity(aPosition);
313   delete theDencity;
314   
315   return coeff_antiproton * density;
316   */
317
318   return 0.0;
319}
320
321G4double G4RKFieldIntegrator::GetKaonPotential(G4double )
322{
323   /*
324   //G4double theM = G4NucleiProperties::GetAtomicMass(theA, theZ);
325   G4double theM = theZ * G4Proton::Proton()->GetPDGMass()
326      + (theA - theZ) * G4Neutron::Neutron()->GetPDGMass()
327      + G4CreateNucleus::GetBindingEnergy(theZ, theA);
328     
329   const G4double Mk  = 496. * MeV;      // mass of "kaon"
330   G4double mu = (theM * Mk)/(theM + Mk);
331   
332   // kaon's potential coefficient
333   //   V = coeff_kaon * nucleus_density
334   G4double coeff_kaon = -2.*pi/mu * (1. + Mk/theM) * a_kaon;
335   
336   G4VNuclearDensity *theDencity;
337   if(theA < 17) theDencity = new G4NuclearShellModelDensity(theA, theZ);
338   else          theDencity = new G4NuclearFermiDensity(theA, theZ);
339   
340   // GetDencity() accepts only G4ThreeVector so build it:   
341   G4ThreeVector aPosition(0.0, 0.0, radius);
342   G4double density = theDencity->GetDensity(aPosition);
343   delete theDencity;
344   
345   return coeff_kaon * density;
346   */
347
348   return 0.0;
349}
350
351G4double G4RKFieldIntegrator::GetPionPotential(G4double )
352{
353   /*
354   //G4double theM = G4NucleiProperties::GetAtomicMass(theA, theZ);
355   G4double theM = theZ * G4Proton::Proton()->GetPDGMass()
356      + (theA - theZ) * G4Neutron::Neutron()->GetPDGMass()
357      + G4CreateNucleus::GetBindingEnergy(theZ, theA);
358     
359   const G4double Mpi = 139. * MeV;      // mass of "pion"
360   G4double mu = (theM * Mpi)/(theM + Mpi);
361
362   // pion's potential coefficient
363   //   V = coeff_pion * nucleus_density
364   G4double coeff_pion = -2.*pi/mu * (1. + Mpi) * a_pion;
365   
366   G4VNuclearDensity *theDencity;
367   if(theA < 17) theDencity = new G4NuclearShellModelDensity(theA, theZ);
368   else          theDencity = new G4NuclearFermiDensity(theA, theZ);
369   
370   // GetDencity() accepts only G4ThreeVector so build it:   
371   G4ThreeVector aPosition(0.0, 0.0, radius);
372   G4double density = theDencity->GetDensity(aPosition);
373   delete theDencity;
374   
375   return coeff_pion * density;
376   */
377
378   return 0.0;
379}
Note: See TracBrowser for help on using the repository browser.