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

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

update processes
File size: 12.8 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.6 2008/11/18 18:26:30 ahoward Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
29//
30// Hadronic Process: Nuclear De-excitations
31// by V. Lara (Oct 1998)
32//
33
34#include "G4GEMChannel.hh"
35#include "G4PairingCorrection.hh"
36
37
38G4GEMChannel::G4GEMChannel(const G4int theA, const G4int theZ,
39                           G4GEMProbability * 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
52G4GEMChannel::G4GEMChannel(const G4int theA, const G4int theZ, const G4String & aName,
53                           G4GEMProbability * 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
67G4GEMChannel::G4GEMChannel(const G4int theA, const G4int theZ, const G4String * aName,
68                           G4GEMProbability * 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
83G4GEMChannel::~G4GEMChannel()
84{
85    if (MyOwnLevelDensity) delete theLevelDensityPtr;
86}
87
88
89
90
91G4GEMChannel::G4GEMChannel(const G4GEMChannel & ) : G4VEvaporationChannel()
92{
93    throw G4HadronicException(__FILE__, __LINE__, "G4GEMChannel::copy_costructor meant to not be accessable");
94}
95
96const G4GEMChannel & G4GEMChannel::operator=(const G4GEMChannel & )
97{
98    throw G4HadronicException(__FILE__, __LINE__, "G4GEMChannel::operator= meant to not be accessable");
99    return *this;
100}
101
102G4bool G4GEMChannel::operator==(const G4GEMChannel & right) const 
103{
104    return (this == (G4GEMChannel *) &right);
105    //  return false;
106}
107
108G4bool G4GEMChannel::operator!=(const G4GEMChannel & right) const 
109{
110    return (this != (G4GEMChannel *) &right);
111    //  return true;
112}
113
114
115
116void G4GEMChannel::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    G4double ExEnergy = fragment.GetExcitationEnergy() -
128      G4PairingCorrection::GetInstance()->GetPairingCorrection(AResidual,ZResidual);
129   
130    // We only take into account channels which are physically allowed
131    if (AResidual <= 0 || ZResidual <= 0 || AResidual < ZResidual ||
132        (AResidual == ZResidual && AResidual > 1) || ExEnergy <= 0.0) 
133      {
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
143//      std::cout << "\tfragment (" << A << ',' << Z << ") residual (" << AResidual << ',' << ZResidual << ')';
144//      std::cout << " U = " << fragment.GetExcitationEnergy();
145//      std::cout << " delta = " << G4PairingCorrection::GetInstance()->GetPairingCorrection(anA,aZ);
146//      std::cout << " U-delta = " << ExEnergy;
147//      std::cout << " V = " << CoulombBarrier;
148       
149        // Maximal Kinetic Energy
150        MaximalKineticEnergy = CalcMaximalKineticEnergy(G4ParticleTable::GetParticleTable()->
151                                                        GetIonTable()->GetNucleusMass(aZ,anA)+ExEnergy);
152
153//      std::cout << " Tmax-V = " << MaximalKineticEnergy << '\n';
154               
155        // Emission probability
156        if (MaximalKineticEnergy <= 0.0) 
157          {
158            EmissionProbability = 0.0;
159          }
160        else 
161          { 
162            // Total emission probability for this channel
163            EmissionProbability = theEvaporationProbabilityPtr->EmissionProbability(fragment,MaximalKineticEnergy);
164          }
165      }
166   
167    return;
168}
169
170
171G4FragmentVector * G4GEMChannel::BreakUp(const G4Fragment & theNucleus)
172{
173   
174    G4double EvaporatedKineticEnergy = CalcKineticEnergy(theNucleus);
175    G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A);
176    G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
177   
178    G4ThreeVector momentum(IsotropicVector(std::sqrt(EvaporatedKineticEnergy*
179                                                (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
180   
181    momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
182
183    G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
184    EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
185   
186    G4Fragment * EvaporatedFragment = new G4Fragment(A,Z,EvaporatedMomentum);
187#ifdef PRECOMPOUND_TEST
188    EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation"));
189#endif
190    // ** And now the residual nucleus **
191    G4double theExEnergy = theNucleus.GetExcitationEnergy();
192    G4double theMass = G4ParticleTable::GetParticleTable()->
193        GetIonTable()->GetNucleusMass(static_cast<G4int>(theNucleus.GetZ()),
194                                      static_cast<G4int>(theNucleus.GetA()));
195    G4double ResidualEnergy = theMass + (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
196       
197    G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
198    ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
199       
200    G4Fragment * ResidualFragment = new G4Fragment( AResidual, ZResidual, ResidualMomentum );
201
202#ifdef PRECOMPOUND_TEST
203    ResidualFragment->SetCreatorModel(G4String("ResidualNucleus"));
204#endif
205   
206   
207    G4FragmentVector * theResult = new G4FragmentVector;
208   
209#ifdef debug
210    G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
211    G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
212    if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 10.0*eV) {
213        G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4GEMChanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl;
214        G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV    Final :" 
215               <<Efinal/MeV << " MeV    Delta: " <<  (Efinal-theNucleus.GetMomentum().e())/MeV
216               << " MeV" << G4endl;
217    }
218    if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 10.0*eV ||
219        std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 10.0*eV ||
220        std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 10.0*eV ) {
221        G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4GEMChanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl;
222        G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV    Final :" 
223               <<Pfinal/MeV << " MeV    Delta: " <<  Pfinal-theNucleus.GetMomentum().vect()
224               << " MeV" << G4endl;   
225       
226    }
227#endif
228    theResult->push_back(EvaporatedFragment);
229    theResult->push_back(ResidualFragment);
230    return theResult; 
231} 
232
233
234
235G4double G4GEMChannel::CalcMaximalKineticEnergy(const G4double NucleusTotalE)
236    // Calculate maximal kinetic energy that can be carried by fragment.
237{
238    G4double ResidualMass = G4ParticleTable::GetParticleTable()->
239        GetIonTable()->GetNucleusMass( ZResidual, AResidual );
240    G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->
241        GetIonTable()->GetNucleusMass( Z, A );
242       
243    G4double T = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass - ResidualMass*ResidualMass)/
244        (2.0*NucleusTotalE) -
245        EvaporatedMass - CoulombBarrier;
246       
247    return T;
248}
249
250
251
252
253G4double G4GEMChannel::CalcKineticEnergy(const G4Fragment & fragment)
254    // Samples fragment kinetic energy.
255{
256    G4double U = fragment.GetExcitationEnergy();
257       
258    G4double NuclearMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(Z,A);
259
260    G4double a = theLevelDensityPtr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),
261                                                           static_cast<G4int>(fragment.GetZ()),U);
262    G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(fragment.GetA()),
263                                                                               static_cast<G4int>(fragment.GetZ()));
264   
265    G4double Alpha = theEvaporationProbabilityPtr->CalcAlphaParam(fragment);
266    G4double Beta = theEvaporationProbabilityPtr->CalcBetaParam(fragment);
267    G4double Normalization = theEvaporationProbabilityPtr->GetNormalization();
268   
269    G4double Ux = (2.5 + 150.0/AResidual)*MeV;
270    G4double Ex = Ux + delta0;
271    G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
272    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));
273
274    G4double InitialLevelDensity;
275    if ( U < Ex ) 
276      {
277        InitialLevelDensity = (pi/12.0)*std::exp((U-E0)/T)/T;
278      } 
279    else 
280      {
281        InitialLevelDensity = (pi/12.0)*std::exp(2*std::sqrt(a*(U-delta0)))/std::pow(a*std::pow(U-delta0,5.0),1.0/4.0);
282      }
283
284    const G4double Spin = theEvaporationProbabilityPtr->GetSpin();
285    G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
286    G4double RN = 0.0;
287    if (A > 4) 
288    {
289        G4double R1 = std::pow(G4double(fragment.GetA()-A),1.0/3.0);
290        G4double R2 = std::pow(G4double(A),1.0/3.0);
291        RN = 1.12*(R1 + R2) - 0.86*((R1+R2)/(R1*R2));
292        RN *= fermi;
293    }
294    else 
295    {
296        RN = 1.5*fermi;
297    }
298    G4double GeometricalXS = pi*RN*RN*std::pow(G4double(fragment.GetA()-A),2./3.); 
299   
300
301    G4double ConstantFactor = g*GeometricalXS*Alpha/InitialLevelDensity;
302    ConstantFactor *= pi/12.0;
303    G4double theEnergy = MaximalKineticEnergy + CoulombBarrier;
304    G4double KineticEnergy;
305    G4double Probability;
306
307//      std::cout << "\t\tEjectile (" << A << ',' << Z << ") V = " << CoulombBarrier
308//            << " Beta = " << Beta << " V+Beta = " << CoulombBarrier+Beta << '\n';
309
310    do 
311      {
312        KineticEnergy =  CoulombBarrier + G4UniformRand()*MaximalKineticEnergy;
313        Probability = ConstantFactor*(KineticEnergy + Beta);
314        if ( theEnergy-KineticEnergy < Ex) 
315          {
316            Probability *= std::exp((theEnergy-KineticEnergy-E0)/T)/T;
317          } 
318        else 
319          {
320            Probability *= std::exp(2*std::sqrt(a*(theEnergy-KineticEnergy-delta0)))/
321              std::pow(a*std::pow(theEnergy-KineticEnergy-delta0,5.0),1.0/4.0);
322          }
323        Probability /= Normalization;
324    }
325    while (G4UniformRand() > Probability);
326
327
328    //  ---------------------------------------------------------------------------------------------------
329   
330    return KineticEnergy;
331} 
332
333
334G4ThreeVector G4GEMChannel::IsotropicVector(const G4double Magnitude)
335    // Samples a isotropic random vectorwith a magnitud given by Magnitude.
336    // By default Magnitude = 1.0
337{
338    G4double CosTheta = 1.0 - 2.0*G4UniformRand();
339    G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
340    G4double Phi = twopi*G4UniformRand();
341    G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
342                         Magnitude*std::sin(Phi)*SinTheta,
343                         Magnitude*CosTheta);
344    return Vector;
345}
346
347
348
Note: See TracBrowser for help on using the repository browser.