source: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationChannel.cc @ 1316

Last change on this file since 1316 was 962, checked in by garnier, 15 years ago

update processes
File size: 13.3 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//J.M. Quesada (August2008). Based on:
28//
29// Hadronic Process: Nuclear De-excitations
30// by V. Lara (Oct 1998)
31//
32// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
33// cross section option
34// JMQ (06 September 2008) Also external choices have been added for
35// superimposed Coulomb barrier (if useSICB is set true, by default is false)
36
37
38#include "G4EvaporationChannel.hh"
39#include "G4PairingCorrection.hh"
40
41
42
43G4EvaporationChannel::G4EvaporationChannel(const G4int anA, const G4int aZ, const G4String & aName,
44                                           G4VEmissionProbability * aEmissionStrategy,
45                                           G4VCoulombBarrier * aCoulombBarrier):
46    G4VEvaporationChannel(aName),
47    theA(anA),
48    theZ(aZ),
49    theEvaporationProbabilityPtr(aEmissionStrategy),
50    theCoulombBarrierPtr(aCoulombBarrier),
51    EmissionProbability(0.0),
52    MaximalKineticEnergy(-1000.0)
53{ 
54    theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
55//    MyOwnLevelDensity = true; 
56}
57
58G4EvaporationChannel::~G4EvaporationChannel()
59{
60  delete theLevelDensityPtr;
61}
62
63G4EvaporationChannel::G4EvaporationChannel(const G4EvaporationChannel & ) : G4VEvaporationChannel()
64{
65    throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationChannel::copy_costructor meant to not be accessable");
66}
67
68const G4EvaporationChannel & G4EvaporationChannel::operator=(const G4EvaporationChannel & )
69{
70    throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationChannel::operator= meant to not be accessable");
71    return *this;
72}
73
74G4bool G4EvaporationChannel::operator==(const G4EvaporationChannel & right) const 
75{
76    return (this == (G4EvaporationChannel *) &right);
77    //  return false;
78}
79
80G4bool G4EvaporationChannel::operator!=(const G4EvaporationChannel & right) const 
81{
82    return (this != (G4EvaporationChannel *) &right);
83    //  return true;
84}
85
86//void G4EvaporationChannel::SetLevelDensityParameter(G4VLevelDensityParameter * aLevelDensity)
87//  {
88//    if (MyOwnLevelDensity) delete theLevelDensityPtr;
89//    theLevelDensityPtr = aLevelDensity;
90//    MyOwnLevelDensity = false;
91//  }
92
93void G4EvaporationChannel::Initialize(const G4Fragment & fragment)
94{
95  //for inverse cross section choice
96  theEvaporationProbabilityPtr->SetOPTxs(OPTxs);
97  // for superimposed Coulomb Barrier for inverse cross sections
98  theEvaporationProbabilityPtr->UseSICB(useSICB);
99 
100 
101  G4int FragmentA = static_cast<G4int>(fragment.GetA());
102  G4int FragmentZ = static_cast<G4int>(fragment.GetZ());
103  ResidualA = FragmentA - theA;
104  ResidualZ = FragmentZ - theZ;
105 
106  //Effective excitation energy
107  G4double ExEnergy = fragment.GetExcitationEnergy() - 
108    G4PairingCorrection::GetInstance()->GetPairingCorrection(FragmentA,FragmentZ);
109 
110 
111  // Only channels which are physically allowed are taken into account
112  if (ResidualA <= 0 || ResidualZ <= 0 || ResidualA < ResidualZ ||
113      (ResidualA == ResidualZ && ResidualA > 1) || ExEnergy <= 0.0) {
114    CoulombBarrier=0.0;
115    MaximalKineticEnergy = -1000.0*MeV;
116    EmissionProbability = 0.0;
117  } else {
118    CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(ResidualA,ResidualZ,ExEnergy);
119    // Maximal Kinetic Energy
120    MaximalKineticEnergy = CalcMaximalKineticEnergy
121      (G4ParticleTable::GetParticleTable()->
122       GetIonTable()->GetNucleusMass(FragmentZ,FragmentA)+ExEnergy);
123   
124    // Emission probability
125    // Protection for the case Tmax<V. If not set in this way we could end up in an
126    // infinite loop in  the method GetKineticEnergy if OPTxs!=0 && useSICB=true.
127    // Of course for OPTxs=0 we have the Coulomb barrier
128   
129    G4double limit;
130    if (OPTxs==0 || (OPTxs!=0 && useSICB)) 
131      limit= CoulombBarrier;
132    else limit=0.;
133 
134    // The threshold for charged particle emission must be  set to 0 if Coulomb
135    //cutoff  is included in the cross sections
136    //if (MaximalKineticEnergy <= 0.0) EmissionProbability = 0.0; 
137    if (MaximalKineticEnergy <= limit) EmissionProbability = 0.0;
138    else { 
139      // Total emission probability for this channel
140      EmissionProbability = theEvaporationProbabilityPtr->
141        EmissionProbability(fragment, MaximalKineticEnergy);
142    }
143  }
144 
145  return;
146}
147
148G4FragmentVector * G4EvaporationChannel::BreakUp(const G4Fragment & theNucleus)
149{
150  G4double EvaporatedKineticEnergy=GetKineticEnergy(theNucleus);
151 
152  G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
153    GetIonMass(theZ,theA);
154  G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
155 
156  G4ThreeVector momentum(IsotropicVector
157                         (std::sqrt(EvaporatedKineticEnergy*
158                                    (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
159 
160  momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
161 
162  G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
163  EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
164 
165  G4Fragment * EvaporatedFragment = new G4Fragment(theA,theZ,EvaporatedMomentum);
166#ifdef PRECOMPOUND_TEST
167  EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation"));
168#endif
169  // ** And now the residual nucleus **
170  G4double theExEnergy = theNucleus.GetExcitationEnergy();
171  G4double theMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
172    GetNucleusMass(static_cast<G4int>(theNucleus.GetZ()),
173                   static_cast<G4int>(theNucleus.GetA()));
174  G4double ResidualEnergy = theMass + 
175    (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
176 
177  G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
178  ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
179 
180  G4Fragment * ResidualFragment = new G4Fragment(ResidualA, ResidualZ, ResidualMomentum );
181 
182#ifdef PRECOMPOUND_TEST
183  ResidualFragment->SetCreatorModel(G4String("ResidualNucleus"));
184#endif
185  G4FragmentVector * theResult = new G4FragmentVector;
186 
187#ifdef debug
188  G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
189  G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
190  if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 1.0*keV) {
191    G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl;
192    G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV    Final :" 
193           <<Efinal/MeV << " MeV    Delta: " <<  (Efinal-theNucleus.GetMomentum().e())/MeV
194           << " MeV" << G4endl;
195  }
196  if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 1.0*keV ||
197      std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 1.0*keV ||
198      std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 1.0*keV ) {
199    G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl;
200    G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV    Final :" 
201           <<Pfinal/MeV << " MeV    Delta: " <<  Pfinal-theNucleus.GetMomentum().vect()
202           << " MeV" << G4endl;   
203   
204  }
205#endif
206  theResult->push_back(EvaporatedFragment);
207  theResult->push_back(ResidualFragment);
208  return theResult; 
209} 
210
211/////////////////////////////////////////
212// Calculates the maximal kinetic energy that can be carried by fragment.
213G4double G4EvaporationChannel::CalcMaximalKineticEnergy(const G4double NucleusTotalE)
214{
215  G4double ResidualMass = G4ParticleTable::GetParticleTable()->
216    GetIonTable()->GetNucleusMass( ResidualZ, ResidualA );
217  G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->
218    GetIonTable()->GetNucleusMass( theZ, theA );
219
220  // This is the "true" assimptotic kinetic energy (from energy conservation)   
221  G4double Tmax = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass -               
222                   ResidualMass*ResidualMass)/(2.0*NucleusTotalE) - EvaporatedMass;
223 
224  //JMQ (13-09-08) bug fixed: in the original version the Tmax is calculated
225  //at the Coulomb barrier
226  //IMPORTANT: meaning of Tmax differs in OPTxs=0 and OPTxs!=0
227  //When OPTxs!=0 Tmax is the TRUE (assimptotic) maximal kinetic energy
228 
229  if(OPTxs==0) 
230    Tmax=Tmax- CoulombBarrier;
231 
232  return Tmax;
233}
234
235///////////////////////////////////////////
236//JMQ: New method for MC sampling of kinetic energy. Substitutes old CalcKineticEnergy
237G4double G4EvaporationChannel::GetKineticEnergy(const G4Fragment & aFragment)
238{
239 
240  if (OPTxs==0) {
241    // It uses Dostrovsky's approximation for the inverse reaction cross
242    // in the probability for fragment emission
243    // MaximalKineticEnergy energy in the original version (V.Lara) was calculated at
244    //the Coulomb barrier.
245   
246   
247    if (MaximalKineticEnergy < 0.0)   
248      throw G4HadronicException(__FILE__, __LINE__, 
249                                "G4EvaporationChannel::CalcKineticEnergy: maximal kinetic at the Coulomb barrier is less than 0");
250   
251    G4double Rb = 4.0*theLevelDensityPtr->
252      LevelDensityParameter(ResidualA+theA,ResidualZ+theZ,MaximalKineticEnergy)*
253      MaximalKineticEnergy;
254    G4double RbSqrt = std::sqrt(Rb);
255    G4double PEX1 = 0.0;
256    if (RbSqrt < 160.0) PEX1 = std::exp(-RbSqrt);
257    G4double Rk = 0.0;
258    G4double FRk = 0.0;
259    do {
260      G4double RandNumber = G4UniformRand();
261      Rk = 1.0 + (1./RbSqrt)*std::log(RandNumber + (1.0-RandNumber)*PEX1);
262      G4double Q1 = 1.0;
263      G4double Q2 = 1.0;
264      if (theZ == 0) { // for emitted neutron
265        G4double Beta = (2.12/std::pow(G4double(ResidualA),2./3.) - 0.05)*MeV/
266          (0.76 + 2.2/std::pow(G4double(ResidualA),1./3.));
267        Q1 = 1.0 + Beta/(MaximalKineticEnergy);
268        Q2 = Q1*std::sqrt(Q1);
269      } 
270     
271      FRk = (3.0*std::sqrt(3.0)/2.0)/Q2 * Rk * (Q1 - Rk*Rk);
272     
273    } while (FRk < G4UniformRand());
274   
275    G4double result =  MaximalKineticEnergy * (1.0-Rk*Rk) + CoulombBarrier;
276    return result;
277  } else if (OPTxs==1 || OPTxs==2 || OPTxs==3 || OPTxs==4) {
278   
279    G4double V;
280    if(useSICB) V= CoulombBarrier;
281    else V=0.;
282    //If Coulomb barrier is just included  in the cross sections
283    //  G4double V=0.;
284
285    G4double Tmax=MaximalKineticEnergy;
286    G4double T(0.0);
287    G4double NormalizedProbability(1.0);
288   
289    // A pointer is created in order to access the distribution function.
290    G4EvaporationProbability * G4EPtemp;
291   
292    if (theA==1 && theZ==0) G4EPtemp=new G4NeutronEvaporationProbability();
293    else if (theA==1 && theZ==1) G4EPtemp=new G4ProtonEvaporationProbability();
294    else if (theA==2 && theZ==1 ) G4EPtemp=new G4DeuteronEvaporationProbability();
295    else if (theA==3 && theZ==1 ) G4EPtemp=new G4TritonEvaporationProbability();
296    else if (theA==3 && theZ==2 ) G4EPtemp=new G4He3EvaporationProbability();
297    else if (theA==4 && theZ==2) G4EPtemp=new G4AlphaEvaporationProbability(); 
298    else {
299      std::ostringstream errOs;
300      errOs << "ejected particle out of range in G4EvaporationChannel"  << G4endl;
301      throw G4HadronicException(__FILE__, __LINE__, errOs.str());
302    }
303
304      //for cross section selection and superimposed Coulom Barrier for xs
305      G4EPtemp->SetOPTxs(OPTxs);
306      G4EPtemp->UseSICB(useSICB);
307
308    do
309      { 
310        T=V+G4UniformRand()*(Tmax-V);
311        NormalizedProbability=G4EPtemp->ProbabilityDistributionFunction(aFragment,T)/
312          (this->GetEmissionProbability());
313       
314      }
315    while (G4UniformRand() > NormalizedProbability);
316   delete G4EPtemp;
317   return T;
318  } else{
319    std::ostringstream errOs;
320    errOs << "Bad option for energy sampling in evaporation"  <<G4endl;
321    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
322  }
323}
324
325G4ThreeVector G4EvaporationChannel::IsotropicVector(const G4double Magnitude)
326    // Samples a isotropic random vectorwith a magnitud given by Magnitude.
327    // By default Magnitude = 1.0
328{
329    G4double CosTheta = 1.0 - 2.0*G4UniformRand();
330    G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
331    G4double Phi = twopi*G4UniformRand();
332    G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
333                         Magnitude*std::sin(Phi)*SinTheta,
334                         Magnitude*CosTheta);
335    return Vector;
336            }
337
338
339
340
341
342   
343
344
345 
346
Note: See TracBrowser for help on using the repository browser.