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

Last change on this file since 1253 was 1228, checked in by garnier, 15 years ago

update geant4.9.3 tag
File size: 14.6 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.10 2009/10/08 07:55:39 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-03 $
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
43
44G4GEMChannel::G4GEMChannel(const G4int theA, const G4int theZ,
45                           G4GEMProbability * aEmissionStrategy,
46                           G4VCoulombBarrier * aCoulombBarrier):
47    A(theA),
48    Z(theZ),
49    theEvaporationProbabilityPtr(aEmissionStrategy),
50    theCoulombBarrierPtr(aCoulombBarrier),
51    EmissionProbability(0.0),
52    MaximalKineticEnergy(-1000.0)
53{ 
54    theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
55    MyOwnLevelDensity = true;
56}
57
58G4GEMChannel::G4GEMChannel(const G4int theA, const G4int theZ, const G4String & aName,
59                           G4GEMProbability * aEmissionStrategy,
60                           G4VCoulombBarrier * aCoulombBarrier) :
61  G4VEvaporationChannel(aName),
62  A(theA),
63  Z(theZ),
64  theEvaporationProbabilityPtr(aEmissionStrategy),
65  theCoulombBarrierPtr(aCoulombBarrier),
66  EmissionProbability(0.0),
67  MaximalKineticEnergy(-1000.0)
68{ 
69  theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
70  MyOwnLevelDensity = true;
71}
72
73G4GEMChannel::~G4GEMChannel()
74{
75    if (MyOwnLevelDensity) delete theLevelDensityPtr;
76}
77
78G4GEMChannel::G4GEMChannel(const G4GEMChannel & ) : G4VEvaporationChannel()
79{
80    throw G4HadronicException(__FILE__, __LINE__, "G4GEMChannel::copy_costructor meant to not be accessable");
81}
82
83const G4GEMChannel & G4GEMChannel::operator=(const G4GEMChannel & )
84{
85    throw G4HadronicException(__FILE__, __LINE__, "G4GEMChannel::operator= meant to not be accessable");
86    return *this;
87}
88
89G4bool G4GEMChannel::operator==(const G4GEMChannel & right) const 
90{
91    return (this == (G4GEMChannel *) &right);
92    //  return false;
93}
94
95G4bool G4GEMChannel::operator!=(const G4GEMChannel & right) const 
96{
97    return (this != (G4GEMChannel *) &right);
98    //  return true;
99}
100
101void G4GEMChannel::Initialize(const G4Fragment & fragment)
102{
103  G4int anA = static_cast<G4int>(fragment.GetA());
104  G4int aZ = static_cast<G4int>(fragment.GetZ());
105  AResidual = anA - A;
106  ZResidual = aZ - Z;
107  //G4cout << "G4GEMChannel::Initialize: Z= " << aZ << " A= " << anA
108  //       << " Zres= " << ZResidual << " Ares= " << AResidual << G4endl;
109
110  // We only take into account channels which are physically allowed
111  if (AResidual <= 0 || ZResidual <= 0 || AResidual < ZResidual ||
112      (AResidual == ZResidual && AResidual > 1)) 
113    {
114      CoulombBarrier = 0.0;
115      MaximalKineticEnergy = -1000.0*MeV;
116      EmissionProbability = 0.0;
117    } 
118  else 
119    {
120
121      // Effective excitation energy
122      // JMQ 071009: pairing in ExEnergy should be the one of parent compound nucleus
123      // FIXED the bug causing reported crash by VI (negative Probabilities
124      // due to inconsistency in Coulomb barrier calculation (CoulombBarrier and -Beta
125      // param for protons must be the same)   
126      //    G4double ExEnergy = fragment.GetExcitationEnergy() -
127      //    G4PairingCorrection::GetInstance()->GetPairingCorrection(AResidual,ZResidual);
128      G4double ExEnergy = fragment.GetExcitationEnergy() -
129        G4PairingCorrection::GetInstance()->GetPairingCorrection(anA,aZ);
130
131      //G4cout << "Eexc(MeV)= " << ExEnergy/MeV << G4endl;
132
133      if( ExEnergy <= 0.0) {
134        CoulombBarrier = 0.0;
135        MaximalKineticEnergy = -1000.0*MeV;
136        EmissionProbability = 0.0;
137
138      } else {
139
140        // Coulomb Barrier calculation
141        CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(AResidual,ZResidual,ExEnergy);
142        //G4cout << "CBarrier(MeV)= " << CoulombBarrier/MeV << G4endl;
143
144        //Maximal kinetic energy (JMQ : at the Coulomb barrier)
145        MaximalKineticEnergy = 
146          CalcMaximalKineticEnergy(G4ParticleTable::GetParticleTable()->
147                                   GetIonTable()->GetNucleusMass(aZ,anA)+ExEnergy);
148        //G4cout << "MaxE(MeV)= " << MaximalKineticEnergy/MeV << G4endl;
149
150               
151        // Emission probability
152        if (MaximalKineticEnergy <= 0.0) 
153          {
154            EmissionProbability = 0.0;
155          }
156        else 
157          { 
158            // Total emission probability for this channel
159            EmissionProbability = 
160              theEvaporationProbabilityPtr->EmissionProbability(fragment,MaximalKineticEnergy);
161          }
162      }
163    }
164   
165    //G4cout << "Prob= " << EmissionProbability << G4endl;
166    return;
167}
168
169G4FragmentVector * G4GEMChannel::BreakUp(const G4Fragment & theNucleus)
170{
171    G4double EvaporatedKineticEnergy = CalcKineticEnergy(theNucleus);
172    G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A);
173    G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
174
175   
176    G4ThreeVector momentum(IsotropicVector(std::sqrt(EvaporatedKineticEnergy*
177                                                (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
178   
179    momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
180
181    G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
182    EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
183    G4Fragment * EvaporatedFragment = new G4Fragment(A,Z,EvaporatedMomentum);
184#ifdef PRECOMPOUND_TEST
185    EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation"));
186#endif
187    // ** And now the residual nucleus **
188    G4double theExEnergy = theNucleus.GetExcitationEnergy();
189    G4double theMass = G4ParticleTable::GetParticleTable()->
190        GetIonTable()->GetNucleusMass(static_cast<G4int>(theNucleus.GetZ()),
191                                      static_cast<G4int>(theNucleus.GetA()));
192    G4double ResidualEnergy = theMass + (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
193       
194    G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
195    ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
196       
197    G4Fragment * ResidualFragment = new G4Fragment( AResidual, ZResidual, ResidualMomentum );
198
199#ifdef PRECOMPOUND_TEST
200    ResidualFragment->SetCreatorModel(G4String("ResidualNucleus"));
201#endif
202   
203   
204    G4FragmentVector * theResult = new G4FragmentVector;
205   
206#ifdef debug
207    G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
208    G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
209    if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 10.0*eV) {
210        G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4GEMChanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl;
211        G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV    Final :" 
212               <<Efinal/MeV << " MeV    Delta: " <<  (Efinal-theNucleus.GetMomentum().e())/MeV
213               << " MeV" << G4endl;
214    }
215    if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 10.0*eV ||
216        std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 10.0*eV ||
217        std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 10.0*eV ) {
218        G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4GEMChanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl;
219        G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV    Final :" 
220               <<Pfinal/MeV << " MeV    Delta: " <<  Pfinal-theNucleus.GetMomentum().vect()
221               << " MeV" << G4endl;   
222       
223    }
224#endif
225    theResult->push_back(EvaporatedFragment);
226    theResult->push_back(ResidualFragment);
227    return theResult; 
228} 
229
230
231
232G4double G4GEMChannel::CalcMaximalKineticEnergy(const G4double NucleusTotalE)
233// Calculate maximal kinetic energy that can be carried by fragment.
234//JMQ this is not the assimptotic kinetic energy but the one at the Coulomb barrier
235{
236  G4double ResidualMass = G4ParticleTable::GetParticleTable()->
237    GetIonTable()->GetNucleusMass( ZResidual, AResidual );
238  G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->
239    GetIonTable()->GetNucleusMass( Z, A );
240       
241  G4double T = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass - ResidualMass*ResidualMass)/
242    (2.0*NucleusTotalE) -
243    EvaporatedMass - CoulombBarrier;
244       
245  return T;
246}
247
248
249
250
251G4double G4GEMChannel::CalcKineticEnergy(const G4Fragment & fragment)
252// Samples fragment kinetic energy.
253{
254  G4double U = fragment.GetExcitationEnergy();
255  G4double NuclearMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(Z,A);
256 
257  G4double Alpha = theEvaporationProbabilityPtr->CalcAlphaParam(fragment);
258  G4double Beta = theEvaporationProbabilityPtr->CalcBetaParam(fragment);
259
260  G4double Normalization = theEvaporationProbabilityPtr->GetNormalization();
261
262//                       ***RESIDUAL***
263  //JMQ (September 2009) the following quantities  refer to the RESIDUAL:
264  G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(AResidual,ZResidual);
265  G4double Ux = (2.5 + 150.0/AResidual)*MeV;
266  G4double Ex = Ux + delta0;
267  G4double InitialLevelDensity;
268  //                    ***end RESIDUAL ***
269 
270  //                       ***PARENT***
271  //JMQ (September 2009) the following quantities   refer to the PARENT:
272 
273  G4double deltaCN = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(fragment.GetA()),
274                                                                              static_cast<G4int>(fragment.GetZ()));
275  G4double aCN = theLevelDensityPtr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),
276                                                           static_cast<G4int>(fragment.GetZ()),U-deltaCN);   
277  G4double UxCN = (2.5 + 150.0/fragment.GetA())*MeV;
278  G4double ExCN = UxCN + deltaCN;
279  G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
280  G4double E0CN = ExCN - TCN*(std::log(TCN/MeV) - std::log(aCN*MeV)/4.0 - 1.25*std::log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
281  //                       ***end PARENT***
282 
283  //JMQ quantities calculated for  CN in InitialLevelDensity
284  if ( U < ExCN ) 
285    {
286      InitialLevelDensity = (pi/12.0)*std::exp((U-E0CN)/TCN)/TCN;
287    } 
288  else 
289    {
290      InitialLevelDensity = (pi/12.0)*std::exp(2*std::sqrt(aCN*(U-deltaCN)))/std::pow(aCN*std::pow(U-deltaCN,5.0),1.0/4.0);
291    }
292 
293  const G4double Spin = theEvaporationProbabilityPtr->GetSpin();
294//JMQ  BIG BUG fixed: hbarc instead of hbar_Planck !!!!
295//     it was fixed in total probability (for this channel) but remained still here!!
296//    G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
297    G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbarc*hbarc);
298//
299//JMQ  fix on Rb and  geometrical cross sections according to Furihata's paper
300//                      (JAERI-Data/Code 2001-105, p6)
301  G4double Rb = 0.0; 
302    if (A > 4) 
303    {
304      G4double Ad = std::pow(G4double(fragment.GetA()-A),1.0/3.0);
305      G4double Aj = std::pow(G4double(A),1.0/3.0);
306      Rb = 1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85;
307      Rb *= fermi;
308    }
309    else if (A>1)
310      {
311        //      G4double R1 = std::pow(G4double(fragment.GetA()-A),1.0/3.0);
312        //      G4double R2 = std::pow(G4double(A),1.0/3.0);
313      G4double Ad = std::pow(G4double(fragment.GetA()-A),1.0/3.0);
314      G4double Aj = std::pow(G4double(A),1.0/3.0);
315      //        RN = 1.12*(R1 + R2) - 0.86*((R1+R2)/(R1*R2));
316      Rb=1.5*(Aj+Ad)*fermi;
317    }
318    else 
319    {
320      G4double Ad = std::pow(G4double(fragment.GetA()-A),1.0/3.0);
321      //        RN = 1.5*fermi;
322      Rb = 1.5*Ad*fermi;
323    }
324    //   G4double GeometricalXS = pi*RN*RN*std::pow(G4double(fragment.GetA()-A),2./3.);
325    G4double GeometricalXS = pi*Rb*Rb; 
326   
327
328    G4double ConstantFactor = g*GeometricalXS*Alpha/InitialLevelDensity;
329    ConstantFactor *= pi/12.0;
330    //JMQ : this is the assimptotic maximal kinetic energy of the ejectile
331    //      (obtained by energy-momentum conservation).
332    //      In general smaller than U-theSeparationEnergy
333    G4double theEnergy = MaximalKineticEnergy + CoulombBarrier;
334    G4double KineticEnergy;
335    G4double Probability;
336
337    do 
338      {
339        KineticEnergy =  CoulombBarrier + G4UniformRand()*MaximalKineticEnergy;
340        Probability = ConstantFactor*(KineticEnergy + Beta);
341        G4double a = theLevelDensityPtr->LevelDensityParameter(AResidual,ZResidual,theEnergy-KineticEnergy-delta0);
342        G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
343        //JMQ fix in units
344        G4double E0 = Ex - T*(std::log(T/MeV) - std::log(a*MeV)/4.0 - 1.25*std::log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
345       
346        if ( theEnergy-KineticEnergy < Ex) 
347          {
348            Probability *= std::exp((theEnergy-KineticEnergy-E0)/T)/T;
349          } 
350        else 
351          {
352            Probability *= std::exp(2*std::sqrt(a*(theEnergy-KineticEnergy-delta0)))/
353              std::pow(a*std::pow(theEnergy-KineticEnergy-delta0,5.0),1.0/4.0);
354          }
355        Probability /= Normalization;
356      }
357    while (G4UniformRand() > Probability);
358   
359    return KineticEnergy;
360} 
361
362
363G4ThreeVector G4GEMChannel::IsotropicVector(const G4double Magnitude)
364    // Samples a isotropic random vectorwith a magnitud given by Magnitude.
365    // By default Magnitude = 1.0
366{
367    G4double CosTheta = 1.0 - 2.0*G4UniformRand();
368    G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
369    G4double Phi = twopi*G4UniformRand();
370    G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
371                         Magnitude*std::sin(Phi)*SinTheta,
372                         Magnitude*CosTheta);
373    return Vector;
374}
375
376
377
Note: See TracBrowser for help on using the repository browser.