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

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

import all except CVS
File size: 12.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.5 2006/06/29 20:22:07 gunter Exp $
28// GEANT4 tag $Name:  $
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   
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      {
132        CoulombBarrier = 0.0;
133        MaximalKineticEnergy = -1000.0*MeV;
134        EmissionProbability = 0.0;
135      } 
136    else 
137      {
138        // Coulomb Barrier calculation
139        CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(AResidual,ZResidual,ExEnergy);
140
141//      std::cout << "\tfragment (" << A << ',' << Z << ") residual (" << AResidual << ',' << ZResidual << ')';
142//      std::cout << " U = " << fragment.GetExcitationEnergy();
143//      std::cout << " delta = " << G4PairingCorrection::GetInstance()->GetPairingCorrection(anA,aZ);
144//      std::cout << " U-delta = " << ExEnergy;
145//      std::cout << " V = " << CoulombBarrier;
146       
147        // Maximal Kinetic Energy
148        MaximalKineticEnergy = CalcMaximalKineticEnergy(G4ParticleTable::GetParticleTable()->
149                                                        GetIonTable()->GetNucleusMass(aZ,anA)+ExEnergy);
150
151//      std::cout << " Tmax-V = " << MaximalKineticEnergy << '\n';
152               
153        // Emission probability
154        if (MaximalKineticEnergy <= 0.0) 
155          {
156            EmissionProbability = 0.0;
157          }
158        else 
159          { 
160            // Total emission probability for this channel
161            EmissionProbability = theEvaporationProbabilityPtr->EmissionProbability(fragment,MaximalKineticEnergy);
162          }
163      }
164   
165    return;
166}
167
168
169G4FragmentVector * G4GEMChannel::BreakUp(const G4Fragment & theNucleus)
170{
171   
172    G4double EvaporatedKineticEnergy = CalcKineticEnergy(theNucleus);
173    G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A);
174    G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
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   
184    G4Fragment * EvaporatedFragment = new G4Fragment(A,Z,EvaporatedMomentum);
185#ifdef PRECOMPOUND_TEST
186    EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation"));
187#endif
188    // ** And now the residual nucleus **
189    G4double theExEnergy = theNucleus.GetExcitationEnergy();
190    G4double theMass = G4ParticleTable::GetParticleTable()->
191        GetIonTable()->GetNucleusMass(static_cast<G4int>(theNucleus.GetZ()),
192                                      static_cast<G4int>(theNucleus.GetA()));
193    G4double ResidualEnergy = theMass + (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
194       
195    G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
196    ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
197       
198    G4Fragment * ResidualFragment = new G4Fragment( AResidual, ZResidual, ResidualMomentum );
199
200#ifdef PRECOMPOUND_TEST
201    ResidualFragment->SetCreatorModel(G4String("ResidualNucleus"));
202#endif
203   
204   
205    G4FragmentVector * theResult = new G4FragmentVector;
206   
207#ifdef debug
208    G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
209    G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
210    if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 10.0*eV) {
211        G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4GEMChanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl;
212        G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV    Final :" 
213               <<Efinal/MeV << " MeV    Delta: " <<  (Efinal-theNucleus.GetMomentum().e())/MeV
214               << " MeV" << G4endl;
215    }
216    if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 10.0*eV ||
217        std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 10.0*eV ||
218        std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 10.0*eV ) {
219        G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4GEMChanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl;
220        G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV    Final :" 
221               <<Pfinal/MeV << " MeV    Delta: " <<  Pfinal-theNucleus.GetMomentum().vect()
222               << " MeV" << G4endl;   
223       
224    }
225#endif
226    theResult->push_back(EvaporatedFragment);
227    theResult->push_back(ResidualFragment);
228    return theResult; 
229} 
230
231
232
233G4double G4GEMChannel::CalcMaximalKineticEnergy(const G4double NucleusTotalE)
234    // Calculate maximal kinetic energy that can be carried by fragment.
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       
256    G4double NuclearMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(Z,A);
257
258    G4double a = theLevelDensityPtr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),
259                                                           static_cast<G4int>(fragment.GetZ()),U);
260    G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(fragment.GetA()),
261                                                                               static_cast<G4int>(fragment.GetZ()));
262   
263    G4double Alpha = theEvaporationProbabilityPtr->CalcAlphaParam(fragment);
264    G4double Beta = theEvaporationProbabilityPtr->CalcBetaParam(fragment);
265    G4double Normalization = theEvaporationProbabilityPtr->GetNormalization();
266   
267    G4double Ux = (2.5 + 150.0/AResidual)*MeV;
268    G4double Ex = Ux + delta0;
269    G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
270    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));
271
272    G4double InitialLevelDensity;
273    if ( U < Ex ) 
274      {
275        InitialLevelDensity = (pi/12.0)*std::exp((U-E0)/T)/T;
276      } 
277    else 
278      {
279        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);
280      }
281
282    const G4double Spin = theEvaporationProbabilityPtr->GetSpin();
283    G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
284    G4double RN = 0.0;
285    if (A > 4) 
286    {
287        G4double R1 = std::pow(G4double(fragment.GetA()-A),1.0/3.0);
288        G4double R2 = std::pow(G4double(A),1.0/3.0);
289        RN = 1.12*(R1 + R2) - 0.86*((R1+R2)/(R1*R2));
290        RN *= fermi;
291    }
292    else 
293    {
294        RN = 1.5*fermi;
295    }
296    G4double GeometricalXS = pi*RN*RN*std::pow(G4double(fragment.GetA()-A),2./3.); 
297   
298
299    G4double ConstantFactor = g*GeometricalXS*Alpha/InitialLevelDensity;
300    ConstantFactor *= pi/12.0;
301    G4double theEnergy = MaximalKineticEnergy + CoulombBarrier;
302    G4double KineticEnergy;
303    G4double Probability;
304
305//      std::cout << "\t\tEjectile (" << A << ',' << Z << ") V = " << CoulombBarrier
306//            << " Beta = " << Beta << " V+Beta = " << CoulombBarrier+Beta << '\n';
307
308    do 
309      {
310        KineticEnergy =  CoulombBarrier + G4UniformRand()*MaximalKineticEnergy;
311        Probability = ConstantFactor*(KineticEnergy + Beta);
312        if ( theEnergy-KineticEnergy < Ex) 
313          {
314            Probability *= std::exp((theEnergy-KineticEnergy-E0)/T)/T;
315          } 
316        else 
317          {
318            Probability *= std::exp(2*std::sqrt(a*(theEnergy-KineticEnergy-delta0)))/
319              std::pow(a*std::pow(theEnergy-KineticEnergy-delta0,5.0),1.0/4.0);
320          }
321        Probability /= Normalization;
322    }
323    while (G4UniformRand() > Probability);
324
325
326    //  ---------------------------------------------------------------------------------------------------
327   
328    return KineticEnergy;
329} 
330
331
332G4ThreeVector G4GEMChannel::IsotropicVector(const G4double Magnitude)
333    // Samples a isotropic random vectorwith a magnitud given by Magnitude.
334    // By default Magnitude = 1.0
335{
336    G4double CosTheta = 1.0 - 2.0*G4UniformRand();
337    G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
338    G4double Phi = twopi*G4UniformRand();
339    G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
340                         Magnitude*std::sin(Phi)*SinTheta,
341                         Magnitude*CosTheta);
342    return Vector;
343}
344
345
346
Note: See TracBrowser for help on using the repository browser.