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

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

import all except CVS
File size: 11.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// $Id: G4EvaporationChannel.cc,v 1.5 2006/06/29 20:10:27 gunter Exp $
28// GEANT4 tag $Name:  $
29//
30// Hadronic Process: Nuclear De-excitations
31// by V. Lara (Oct 1998)
32//
33
34#include "G4EvaporationChannel.hh"
35#include "G4PairingCorrection.hh"
36
37
38G4EvaporationChannel::G4EvaporationChannel(const G4int theA, const G4int theZ,
39                                           G4VEmissionProbability * aEmissionStrategy,
40                                           G4VCoulombBarrier * aCoulombBarrier):
41    A(theA),
42    Z(theZ),
43    theEvaporationProbabilityPtr(aEmissionStrategy),
44    theCoulombBarrierPtr(aCoulombBarrier),
45    EmissionProbability(0.0),
46    MaximalKineticEnergy(-1000.0)
47{ 
48    theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
49    MyOwnLevelDensity = true;
50}
51
52G4EvaporationChannel::G4EvaporationChannel(const G4int theA, const G4int theZ, const G4String & aName,
53                                           G4VEmissionProbability * aEmissionStrategy,
54                                           G4VCoulombBarrier * aCoulombBarrier):
55    G4VEvaporationChannel(aName),
56    A(theA),
57    Z(theZ),
58    theEvaporationProbabilityPtr(aEmissionStrategy),
59    theCoulombBarrierPtr(aCoulombBarrier),
60    EmissionProbability(0.0),
61    MaximalKineticEnergy(-1000.0)
62{ 
63    theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
64    MyOwnLevelDensity = true;
65}
66
67G4EvaporationChannel::G4EvaporationChannel(const G4int theA, const G4int theZ, const G4String * aName,
68                                           G4VEmissionProbability * aEmissionStrategy,
69                                           G4VCoulombBarrier * aCoulombBarrier):
70    G4VEvaporationChannel(aName),
71    A(theA),
72    Z(theZ),
73    theEvaporationProbabilityPtr(aEmissionStrategy),
74    theCoulombBarrierPtr(aCoulombBarrier),
75    EmissionProbability(0.0),
76    MaximalKineticEnergy(-1000.0)
77{ 
78    theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
79    MyOwnLevelDensity = true;
80}
81
82
83G4EvaporationChannel::~G4EvaporationChannel()
84{
85    if (MyOwnLevelDensity) delete theLevelDensityPtr;
86}
87
88
89
90
91G4EvaporationChannel::G4EvaporationChannel(const G4EvaporationChannel & ) : G4VEvaporationChannel()
92{
93    throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationChannel::copy_costructor meant to not be accessable");
94}
95
96const G4EvaporationChannel & G4EvaporationChannel::operator=(const G4EvaporationChannel & )
97{
98    throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationChannel::operator= meant to not be accessable");
99    return *this;
100}
101
102G4bool G4EvaporationChannel::operator==(const G4EvaporationChannel & right) const 
103{
104    return (this == (G4EvaporationChannel *) &right);
105    //  return false;
106}
107
108G4bool G4EvaporationChannel::operator!=(const G4EvaporationChannel & right) const 
109{
110    return (this != (G4EvaporationChannel *) &right);
111    //  return true;
112}
113
114
115
116void G4EvaporationChannel::Initialize(const G4Fragment & fragment)
117{
118
119    G4int anA = static_cast<G4int>(fragment.GetA());
120    G4int aZ = static_cast<G4int>(fragment.GetZ());
121    AResidual = anA - A;
122    ZResidual = aZ - Z;
123
124    // Effective excitation energy
125    G4double ExEnergy = fragment.GetExcitationEnergy() - 
126      G4PairingCorrection::GetInstance()->GetPairingCorrection(anA,aZ);
127
128    // We only take into account channels which are physically allowed
129    if (AResidual <= 0 || ZResidual <= 0 || AResidual < ZResidual ||
130        (AResidual == ZResidual && AResidual > 1) || ExEnergy <= 0.0) {
131        //              LevelDensityParameter = 0.0;
132        CoulombBarrier = 0.0;
133        //              BindingEnergy = 0.0;
134        MaximalKineticEnergy = -1000.0*MeV;
135        EmissionProbability = 0.0;
136    } else {
137        // // Get Level Density
138        // LevelDensityParameter = theLevelDensityPtr->LevelDensityParameter(anA,aZ,ExEnergy);
139
140        // Coulomb Barrier calculation
141        CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(AResidual,ZResidual,ExEnergy);
142       
143        // // Binding Enegy (for separate fragment from nucleus)
144        // BindingEnergy = CalcBindingEnergy(anA,aZ);
145
146        // Maximal Kinetic Energy
147        MaximalKineticEnergy = CalcMaximalKineticEnergy(G4ParticleTable::GetParticleTable()->
148                                                        GetIonTable()->GetNucleusMass(aZ,anA)+ExEnergy);
149               
150        // Emission probability
151        if (MaximalKineticEnergy <= 0.0) EmissionProbability = 0.0;
152        else { 
153            // Total emission probability for this channel
154            EmissionProbability = theEvaporationProbabilityPtr->EmissionProbability(fragment,MaximalKineticEnergy);
155
156        }
157    }
158       
159    return;
160}
161
162
163G4FragmentVector * G4EvaporationChannel::BreakUp(const G4Fragment & theNucleus)
164{
165
166    G4double EvaporatedKineticEnergy = CalcKineticEnergy();
167    G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A);
168    G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
169
170    G4ThreeVector momentum(IsotropicVector(std::sqrt(EvaporatedKineticEnergy*
171                                                (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
172 
173    momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
174
175    G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
176    EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
177
178    G4Fragment * EvaporatedFragment = new G4Fragment(A,Z,EvaporatedMomentum);
179#ifdef PRECOMPOUND_TEST
180    EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation"));
181#endif
182    // ** And now the residual nucleus **
183    G4double theExEnergy = theNucleus.GetExcitationEnergy();
184    G4double theMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
185        GetNucleusMass(static_cast<G4int>(theNucleus.GetZ()),
186                       static_cast<G4int>(theNucleus.GetA()));
187    G4double ResidualEnergy = theMass + (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
188       
189    G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
190    ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
191       
192    G4Fragment * ResidualFragment = new G4Fragment( AResidual, ZResidual, ResidualMomentum );
193#ifdef PRECOMPOUND_TEST
194    ResidualFragment->SetCreatorModel(G4String("ResidualNucleus"));
195#endif
196    G4FragmentVector * theResult = new G4FragmentVector;
197
198#ifdef debug
199    G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
200    G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
201    if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 1.0*keV) {
202        G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl;
203        G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV    Final :" 
204               <<Efinal/MeV << " MeV    Delta: " <<  (Efinal-theNucleus.GetMomentum().e())/MeV
205               << " MeV" << G4endl;
206    }
207    if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 1.0*keV ||
208        std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 1.0*keV ||
209        std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 1.0*keV ) {
210        G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4Evaporation Chanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl;
211        G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV    Final :" 
212               <<Pfinal/MeV << " MeV    Delta: " <<  Pfinal-theNucleus.GetMomentum().vect()
213               << " MeV" << G4endl;   
214
215    }
216#endif
217    theResult->push_back(EvaporatedFragment);
218    theResult->push_back(ResidualFragment);
219    return theResult; 
220} 
221
222
223
224// G4double G4EvaporationChannel::CalcBindingEnergy(const G4int anA, const G4int aZ)
225// // Calculate Binding Energy for separate fragment from nucleus
226// {
227//      // Mass Excess for residual nucleus
228//      G4double ResNucMassExcess = G4NucleiProperties::GetNuclearMass(AResidual,ZResidual);
229//      // Mass Excess for fragment
230//      G4double FragmentMassExcess = G4NucleiProperties::GetNuclearMass(A,Z);
231//      // Mass Excess for Compound Nucleus
232//      G4double NucleusMassExcess = G4NucleiProperties::GetNuclearMass(anA,aZ);
233//
234//      return ResNucMassExcess + FragmentMassExcess - NucleusMassExcess;
235// }
236
237
238G4double G4EvaporationChannel::CalcMaximalKineticEnergy(const G4double NucleusTotalE)
239    // Calculate maximal kinetic energy that can be carried by fragment.
240{
241    G4double ResidualMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass( ZResidual, AResidual );
242    G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass( Z, A );
243       
244    G4double T = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass - ResidualMass*ResidualMass)/
245        (2.0*NucleusTotalE) -
246        EvaporatedMass - CoulombBarrier;
247       
248    return T;
249}
250
251
252
253
254G4double G4EvaporationChannel::CalcKineticEnergy(void)
255    // Samples fragment kinetic energy.
256    // It uses Dostrovsky's approximation for the inverse reaction cross
257    // in the probability for fragment emisson
258{
259    if (MaximalKineticEnergy < 0.0) 
260        throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationChannel::CalcKineticEnergy: maximal kinetic energy is less than 0");
261
262    G4double Rb = 4.0*theLevelDensityPtr->LevelDensityParameter(AResidual+A,ZResidual+Z,MaximalKineticEnergy)*
263        MaximalKineticEnergy;
264    G4double RbSqrt = std::sqrt(Rb);
265    G4double PEX1 = 0.0;
266    if (RbSqrt < 160.0) PEX1 = std::exp(-RbSqrt);
267    G4double Rk = 0.0;
268    G4double FRk = 0.0;
269    do {
270        G4double RandNumber = G4UniformRand();
271        Rk = 1.0 + (1./RbSqrt)*std::log(RandNumber + (1.0-RandNumber)*PEX1);
272        G4double Q1 = 1.0;
273        G4double Q2 = 1.0;
274        if (Z == 0) { // for emitted neutron
275            G4double Beta = (2.12/std::pow(G4double(AResidual),2./3.) - 0.05)*MeV/
276                (0.76 + 2.2/std::pow(G4double(AResidual),1./3.));
277            Q1 = 1.0 + Beta/(MaximalKineticEnergy);
278            Q2 = Q1*std::sqrt(Q1);
279        } 
280   
281        FRk = (3.0*std::sqrt(3.0)/2.0)/Q2 * Rk * (Q1 - Rk*Rk);
282   
283    } while (FRk < G4UniformRand());
284
285    G4double result =  MaximalKineticEnergy * (1.0-Rk*Rk) + CoulombBarrier;
286
287    return result;
288} 
289 
290
291G4ThreeVector G4EvaporationChannel::IsotropicVector(const G4double Magnitude)
292    // Samples a isotropic random vectorwith a magnitud given by Magnitude.
293    // By default Magnitude = 1.0
294{
295    G4double CosTheta = 1.0 - 2.0*G4UniformRand();
296    G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
297    G4double Phi = twopi*G4UniformRand();
298    G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
299                         Magnitude*std::sin(Phi)*SinTheta,
300                         Magnitude*CosTheta);
301    return Vector;
302}
303
304
305
Note: See TracBrowser for help on using the repository browser.