source: trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/src/G4GEMChannel.cc @ 1348

Last change on this file since 1348 was 1347, checked in by garnier, 14 years ago

geant4 tag 9.4
File size: 12.9 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26//
27// $Id: G4GEMChannel.cc,v 1.12 2010/11/18 16:21:17 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// Hadronic Process: Nuclear De-excitations
31// by V. Lara (Oct 1998)
32//
33// J. M. Quesada (September 2009) bugs fixed in  probability distribution for kinetic
34//              energy sampling:
35//              -hbarc instead of hbar_Planck (BIG BUG)
36//              -quantities for initial nucleus and residual are calculated separately
37// V.Ivanchenko (September 2009) Added proper protection for unphysical final state
38// J. M. Quesada (October 2009) fixed bug in CoulombBarrier calculation
39
40#include "G4GEMChannel.hh"
41#include "G4PairingCorrection.hh"
42#include "G4Pow.hh"
43
44G4GEMChannel::G4GEMChannel(G4int theA, G4int theZ, const G4String & aName,
45                           G4GEMProbability * aEmissionStrategy,
46                           G4VCoulombBarrier * aCoulombBarrier) :
47  G4VEvaporationChannel(aName),
48  A(theA),
49  Z(theZ),
50  theEvaporationProbabilityPtr(aEmissionStrategy),
51  theCoulombBarrierPtr(aCoulombBarrier),
52  EmissionProbability(0.0),
53  MaximalKineticEnergy(-1000.0)
54{ 
55  theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
56  MyOwnLevelDensity = true;
57  EvaporatedMass = G4NucleiProperties::GetNuclearMass(A, Z);
58  ResidualMass = CoulombBarrier = 0.0;
59  fG4pow = G4Pow::GetInstance(); 
60}
61
62G4GEMChannel::G4GEMChannel() :
63  G4VEvaporationChannel("dummy"),
64  A(0),
65  Z(0),
66  theEvaporationProbabilityPtr(0),
67  theCoulombBarrierPtr(0),
68  EmissionProbability(0.0),
69  MaximalKineticEnergy(-1000.0)
70{ 
71  theLevelDensityPtr = 0;
72  MyOwnLevelDensity = true;
73  EvaporatedMass = 0.0;
74  ResidualMass = CoulombBarrier = 0.0;
75  fG4pow = G4Pow::GetInstance(); 
76}
77
78G4GEMChannel::~G4GEMChannel()
79{
80  if (MyOwnLevelDensity) { delete theLevelDensityPtr; }
81}
82
83void G4GEMChannel::Initialize(const G4Fragment & fragment)
84{
85  G4int anA = fragment.GetA_asInt();
86  G4int aZ  = fragment.GetZ_asInt();
87  ResidualA = anA - A;
88  ResidualZ = aZ - Z;
89  //G4cout << "G4GEMChannel::Initialize: Z= " << aZ << " A= " << anA
90  //       << " Zres= " << ResidualZ << " Ares= " << ResidualA << G4endl;
91
92  // We only take into account channels which are physically allowed
93  if (ResidualA <= 0 || ResidualZ <= 0 || ResidualA < ResidualZ ||
94      (ResidualA == ResidualZ && ResidualA > 1)) 
95    {
96      CoulombBarrier = 0.0;
97      MaximalKineticEnergy = -1000.0*MeV;
98      EmissionProbability = 0.0;
99    } 
100  else 
101    {
102      // Effective excitation energy
103      // JMQ 071009: pairing in ExEnergy should be the one of parent compound nucleus
104      // FIXED the bug causing reported crash by VI (negative Probabilities
105      // due to inconsistency in Coulomb barrier calculation (CoulombBarrier and -Beta
106      // param for protons must be the same)   
107      //    G4double ExEnergy = fragment.GetExcitationEnergy() -
108      //    G4PairingCorrection::GetInstance()->GetPairingCorrection(ResidualA,ResidualZ);
109      G4double ExEnergy = fragment.GetExcitationEnergy() -
110        G4PairingCorrection::GetInstance()->GetPairingCorrection(anA,aZ);
111
112      //G4cout << "Eexc(MeV)= " << ExEnergy/MeV << G4endl;
113
114      if( ExEnergy <= 0.0) {
115        CoulombBarrier = 0.0;
116        MaximalKineticEnergy = -1000.0*MeV;
117        EmissionProbability = 0.0;
118
119      } else {
120
121        ResidualMass = G4NucleiProperties::GetNuclearMass(ResidualA, ResidualZ);
122
123        // Coulomb Barrier calculation
124        CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(ResidualA,ResidualZ,ExEnergy);
125        //G4cout << "CBarrier(MeV)= " << CoulombBarrier/MeV << G4endl;
126
127        //Maximal kinetic energy (JMQ : at the Coulomb barrier)
128        MaximalKineticEnergy = 
129          CalcMaximalKineticEnergy(fragment.GetGroundStateMass()+ExEnergy);
130        //G4cout << "MaxE(MeV)= " << MaximalKineticEnergy/MeV << G4endl;
131               
132        // Emission probability
133        if (MaximalKineticEnergy <= 0.0) 
134          {
135            EmissionProbability = 0.0;
136          }
137        else 
138          { 
139            // Total emission probability for this channel
140            EmissionProbability = 
141              theEvaporationProbabilityPtr->EmissionProbability(fragment,MaximalKineticEnergy);
142          }
143      }
144    }   
145    //G4cout << "Prob= " << EmissionProbability << G4endl;
146    return;
147}
148
149G4FragmentVector * G4GEMChannel::BreakUp(const G4Fragment & theNucleus)
150{
151  G4double EvaporatedKineticEnergy = CalcKineticEnergy(theNucleus);
152  G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
153 
154  G4ThreeVector momentum(IsotropicVector(std::sqrt(EvaporatedKineticEnergy*
155                                                   (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
156   
157  momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
158
159  G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
160  EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
161  G4Fragment * EvaporatedFragment = new G4Fragment(A,Z,EvaporatedMomentum);
162  // ** And now the residual nucleus **
163  G4double theExEnergy = theNucleus.GetExcitationEnergy();
164  G4double theMass = theNucleus.GetGroundStateMass();
165  G4double ResidualEnergy = theMass + (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
166       
167  G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
168  ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
169       
170  G4Fragment * ResidualFragment = new G4Fragment( ResidualA, ResidualZ, ResidualMomentum );
171   
172  G4FragmentVector * theResult = new G4FragmentVector;
173   
174#ifdef debug
175  G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
176  G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
177  if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 10.0*eV) {
178    G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4GEMChanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl;
179    G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV    Final :" 
180           <<Efinal/MeV << " MeV    Delta: " <<  (Efinal-theNucleus.GetMomentum().e())/MeV
181           << " MeV" << G4endl;
182  }
183  if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 10.0*eV ||
184      std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 10.0*eV ||
185      std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 10.0*eV ) {
186    G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4GEMChanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl;
187    G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV    Final :" 
188           <<Pfinal/MeV << " MeV    Delta: " <<  Pfinal-theNucleus.GetMomentum().vect()
189           << " MeV" << G4endl;   
190       
191  }
192#endif
193  theResult->push_back(EvaporatedFragment);
194  theResult->push_back(ResidualFragment);
195  return theResult; 
196} 
197
198G4double G4GEMChannel::CalcMaximalKineticEnergy(const G4double NucleusTotalE)
199// Calculate maximal kinetic energy that can be carried by fragment.
200//JMQ this is not the assimptotic kinetic energy but the one at the Coulomb barrier
201{
202  G4double T = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass - ResidualMass*ResidualMass)/
203    (2.0*NucleusTotalE) - EvaporatedMass - CoulombBarrier;
204       
205  return T;
206}
207
208G4double G4GEMChannel::CalcKineticEnergy(const G4Fragment & fragment)
209// Samples fragment kinetic energy.
210{
211  G4double U = fragment.GetExcitationEnergy();
212 
213  G4double Alpha = theEvaporationProbabilityPtr->CalcAlphaParam(fragment);
214  G4double Beta = theEvaporationProbabilityPtr->CalcBetaParam(fragment);
215
216  G4double Normalization = theEvaporationProbabilityPtr->GetNormalization();
217
218  //                       ***RESIDUAL***
219  //JMQ (September 2009) the following quantities  refer to the RESIDUAL:
220  G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(ResidualA,ResidualZ);
221  G4double Ux = (2.5 + 150.0/ResidualA)*MeV;
222  G4double Ex = Ux + delta0;
223  G4double InitialLevelDensity;
224  //                    ***end RESIDUAL ***
225 
226  //                       ***PARENT***
227  //JMQ (September 2009) the following quantities   refer to the PARENT:
228 
229  G4double deltaCN = 
230    G4PairingCorrection::GetInstance()->GetPairingCorrection(fragment.GetA_asInt(),
231                                                             fragment.GetZ_asInt());
232  G4double aCN = theLevelDensityPtr->LevelDensityParameter(fragment.GetA_asInt(),
233                                                           fragment.GetZ_asInt(),U-deltaCN);   
234  G4double UxCN = (2.5 + 150.0/fragment.GetA())*MeV;
235  G4double ExCN = UxCN + deltaCN;
236  G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
237  //                       ***end PARENT***
238 
239  //JMQ quantities calculated for  CN in InitialLevelDensity
240  if ( U < ExCN ) 
241    {
242      G4double E0CN = ExCN - TCN*(std::log(TCN/MeV) - std::log(aCN*MeV)/4.0 
243                                  - 1.25*std::log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
244      InitialLevelDensity = (pi/12.0)*std::exp((U-E0CN)/TCN)/TCN;
245    } 
246  else 
247    {
248      G4double x  = U-deltaCN;
249      G4double x1 = std::sqrt(aCN*x);
250      InitialLevelDensity = (pi/12.0)*std::exp(2*x1)/(x*std::sqrt(x1));
251      //InitialLevelDensity =
252      //(pi/12.0)*std::exp(2*std::sqrt(aCN*(U-deltaCN)))/std::pow(aCN*std::pow(U-deltaCN,5.0),1.0/4.0);
253    }
254 
255  const G4double Spin = theEvaporationProbabilityPtr->GetSpin();
256  //JMQ  BIG BUG fixed: hbarc instead of hbar_Planck !!!!
257  //     it was fixed in total probability (for this channel) but remained still here!!
258  //    G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
259  G4double g = (2.0*Spin+1.0)*EvaporatedMass/(pi2* hbarc*hbarc);
260  //
261  //JMQ  fix on Rb and  geometrical cross sections according to Furihata's paper
262  //                      (JAERI-Data/Code 2001-105, p6)
263  G4double Rb = 0.0; 
264  if (A > 4) 
265    {
266      G4double Ad = fG4pow->Z13(ResidualA);
267      G4double Aj = fG4pow->Z13(A);
268      //        RN = 1.12*(R1 + R2) - 0.86*((R1+R2)/(R1*R2));
269      Rb = 1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85;
270      Rb *= fermi;
271    }
272  else if (A>1)
273    {
274      G4double Ad = fG4pow->Z13(ResidualA);
275      G4double Aj = fG4pow->Z13(A);
276      Rb=1.5*(Aj+Ad)*fermi;
277    }
278  else 
279    {
280      G4double Ad = fG4pow->Z13(ResidualA);
281      Rb = 1.5*Ad*fermi;
282    }
283  //   G4double GeometricalXS = pi*RN*RN*std::pow(G4double(fragment.GetA()-A),2./3.);
284  G4double GeometricalXS = pi*Rb*Rb; 
285   
286  G4double ConstantFactor = g*GeometricalXS*Alpha/InitialLevelDensity;
287  ConstantFactor *= pi/12.0;
288  //JMQ : this is the assimptotic maximal kinetic energy of the ejectile
289  //      (obtained by energy-momentum conservation).
290  //      In general smaller than U-theSeparationEnergy
291  G4double theEnergy = MaximalKineticEnergy + CoulombBarrier;
292  G4double KineticEnergy;
293  G4double Probability;
294
295  do 
296    {
297      KineticEnergy =  CoulombBarrier + G4UniformRand()*MaximalKineticEnergy;
298      Probability = ConstantFactor*(KineticEnergy + Beta);
299      G4double a = 
300        theLevelDensityPtr->LevelDensityParameter(ResidualA,ResidualZ,theEnergy-KineticEnergy-delta0);
301      G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
302      //JMQ fix in units
303       
304      if ( theEnergy-KineticEnergy < Ex) 
305        {
306          G4double E0 = Ex - T*(std::log(T/MeV) - std::log(a*MeV)/4.0 
307                                - 1.25*std::log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
308          Probability *= std::exp((theEnergy-KineticEnergy-E0)/T)/T;
309        } 
310      else 
311        {
312          Probability *= std::exp(2*std::sqrt(a*(theEnergy-KineticEnergy-delta0)))/
313            std::pow(a*fG4pow->powN(theEnergy-KineticEnergy-delta0,5), 0.25);
314        }
315    }
316  while (Normalization*G4UniformRand() > Probability);
317   
318  return KineticEnergy;
319} 
320
321
322G4ThreeVector G4GEMChannel::IsotropicVector(const G4double Magnitude)
323    // Samples a isotropic random vectorwith a magnitud given by Magnitude.
324    // By default Magnitude = 1.0
325{
326  G4double CosTheta = 1.0 - 2.0*G4UniformRand();
327  G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
328  G4double Phi = twopi*G4UniformRand();
329  G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
330                       Magnitude*std::sin(Phi)*SinTheta,
331                       Magnitude*CosTheta);
332  return Vector;
333}
334
335
336
Note: See TracBrowser for help on using the repository browser.