| Geant4 Cross Reference |
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 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4PreCompoundEmission
33 //
34 // Author: V.Lara
35 //
36 // Modified:
37 // 15.01.2010 J.M.Quesada added protection against unphysical values of parameter an
38 // 19.01.2010 V.Ivanchenko simplified computation of parameter an, sample cosTheta
39 // instead of theta; protect all calls to sqrt
40 // 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers
41 // use int Z and A and cleanup
42 //
43
44 #include "G4PreCompoundEmission.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47 #include "G4Pow.hh"
48 #include "G4Exp.hh"
49 #include "G4Log.hh"
50 #include "Randomize.hh"
51 #include "G4RandomDirection.hh"
52 #include "G4PreCompoundEmissionFactory.hh"
53 #include "G4HETCEmissionFactory.hh"
54 #include "G4HadronicException.hh"
55 #include "G4NuclearLevelData.hh"
56 #include "G4DeexPrecoParameters.hh"
57 #include "G4PhysicsModelCatalog.hh"
58
59
60 G4PreCompoundEmission::G4PreCompoundEmission()
61 {
62 theFragmentsFactory = new G4PreCompoundEmissionFactory();
63 theFragmentsVector =
64 new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
65 g4calc = G4Pow::GetInstance();
66 fNuclData = G4NuclearLevelData::GetInstance();
67 G4DeexPrecoParameters* param = fNuclData->GetParameters();
68 fFermiEnergy = param->GetFermiEnergy();
69 fUseAngularGenerator = param->UseAngularGen();
70 fModelID = G4PhysicsModelCatalog::GetModelID("model_PRECO");
71 }
72
73 G4PreCompoundEmission::~G4PreCompoundEmission()
74 {
75 delete theFragmentsFactory;
76 delete theFragmentsVector;
77 }
78
79 void G4PreCompoundEmission::SetDefaultModel()
80 {
81 if (theFragmentsFactory) { delete theFragmentsFactory; }
82 theFragmentsFactory = new G4PreCompoundEmissionFactory();
83 if (theFragmentsVector) {
84 theFragmentsVector->SetVector(theFragmentsFactory->GetFragmentVector());
85 } else {
86 theFragmentsVector =
87 new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
88 }
89 }
90
91 void G4PreCompoundEmission::SetHETCModel()
92 {
93 if (theFragmentsFactory) delete theFragmentsFactory;
94 theFragmentsFactory = new G4HETCEmissionFactory();
95 if (theFragmentsVector) {
96 theFragmentsVector->SetVector(theFragmentsFactory->GetFragmentVector());
97 } else {
98 theFragmentsVector =
99 new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
100 }
101 }
102
103 G4ReactionProduct*
104 G4PreCompoundEmission::PerformEmission(G4Fragment & aFragment)
105 {
106 // Choose a Fragment for emission
107 G4VPreCompoundFragment * thePreFragment =
108 theFragmentsVector->ChooseFragment();
109 if (thePreFragment == nullptr)
110 {
111 G4cout << "G4PreCompoundEmission::PerformEmission : "
112 << "I couldn't choose a fragment\n"
113 << "while trying to de-excite\n"
114 << aFragment << G4endl;
115 throw G4HadronicException(__FILE__, __LINE__, "");
116 }
117
118 //G4cout << "Chosen fragment: " << G4endl;
119 //G4cout << *thePreFragment << G4endl;
120
121 // Kinetic Energy of emitted fragment
122 G4double kinEnergy = thePreFragment->SampleKineticEnergy(aFragment);
123 kinEnergy = std::max(kinEnergy, 0.0);
124
125 // Calculate the fragment momentum (three vector)
126 if(fUseAngularGenerator) {
127 AngularDistribution(thePreFragment,aFragment,kinEnergy);
128 } else {
129 G4double pmag =
130 std::sqrt(kinEnergy*(kinEnergy + 2.0*thePreFragment->GetNuclearMass()));
131 theFinalMomentum = pmag*G4RandomDirection();
132 }
133
134 // Mass of emittef fragment
135 G4double EmittedMass = thePreFragment->GetNuclearMass();
136 // Now we can calculate the four momentum
137 // both options are valid and give the same result but 2nd one is faster
138 G4LorentzVector Emitted4Momentum(theFinalMomentum,EmittedMass + kinEnergy);
139
140 // Perform Lorentz boost
141 G4LorentzVector Rest4Momentum = aFragment.GetMomentum();
142 Emitted4Momentum.boost(Rest4Momentum.boostVector());
143
144 // Set emitted fragment momentum
145 thePreFragment->SetMomentum(Emitted4Momentum);
146
147 // NOW THE RESIDUAL NUCLEUS
148 // ------------------------
149
150 Rest4Momentum -= Emitted4Momentum;
151
152 // Update nucleus parameters:
153 // --------------------------
154
155 // Z and A
156 aFragment.SetZandA_asInt(thePreFragment->GetRestZ(),
157 thePreFragment->GetRestA());
158
159 // Number of excitons
160 aFragment.SetNumberOfParticles(aFragment.GetNumberOfParticles()-
161 thePreFragment->GetA());
162 // Number of charges
163 aFragment.SetNumberOfCharged(aFragment.GetNumberOfCharged()-
164 thePreFragment->GetZ());
165
166 // Update nucleus momentum
167 // A check on consistence of Z, A, and mass will be performed
168 aFragment.SetMomentum(Rest4Momentum);
169
170 // Create a G4ReactionProduct
171 G4ReactionProduct * MyRP = thePreFragment->GetReactionProduct();
172
173 // Set the creator model ID
174 aFragment.SetCreatorModelID(fModelID);
175 if (MyRP != nullptr) MyRP->SetCreatorModelID(fModelID);
176
177 return MyRP;
178 }
179
180 void G4PreCompoundEmission::AngularDistribution(
181 G4VPreCompoundFragment* thePreFragment,
182 const G4Fragment& aFragment,
183 G4double ekin)
184 {
185 G4int p = aFragment.GetNumberOfParticles();
186 G4int h = aFragment.GetNumberOfHoles();
187 G4double U = aFragment.GetExcitationEnergy();
188
189 // Emission particle separation energy
190 G4double Bemission = thePreFragment->GetBindingEnergy();
191
192 G4double gg = (6.0/pi2)*fNuclData->GetLevelDensity(aFragment.GetZ_asInt(),
193 aFragment.GetA_asInt(),U);
194
195 // Average exciton energy relative to bottom of nuclear well
196 G4double Eav = 2*p*(p+1)/((p+h)*gg);
197
198 // Excitation energy relative to the Fermi Level
199 G4double Uf = std::max(U - (p - h)*fFermiEnergy , 0.0);
200 // G4double Uf = U - KineticEnergyOfEmittedFragment - Bemission;
201
202 G4double w_num = rho(p+1, h, gg, Uf, fFermiEnergy);
203 G4double w_den = rho(p, h, gg, Uf, fFermiEnergy);
204 if (w_num > 0.0 && w_den > 0.0)
205 {
206 Eav *= (w_num/w_den);
207 Eav += - Uf/(p+h) + fFermiEnergy;
208 }
209 else
210 {
211 Eav = fFermiEnergy;
212 }
213
214 // VI + JMQ 19/01/2010 update computation of the parameter an
215 //
216 G4double an = 0.0;
217 G4double Eeff = ekin + Bemission + fFermiEnergy;
218 if(ekin > DBL_MIN && Eeff > DBL_MIN) {
219
220 G4double zeta = std::max(1.0,9.3/std::sqrt(ekin/CLHEP::MeV));
221
222 // This should be the projectile energy. If I would know which is
223 // the projectile (proton, neutron) I could remove the binding energy.
224 // But, what happens if INC precedes precompound? This approximation
225 // seems to work well enough
226 G4double ProjEnergy = aFragment.GetExcitationEnergy();
227
228 an = 3*std::sqrt((ProjEnergy+fFermiEnergy)*Eeff)/(zeta*Eav);
229
230 G4int ne = aFragment.GetNumberOfExcitons() - 1;
231 if ( ne > 1 ) { an /= static_cast<G4double>(ne); }
232
233 // protection of exponent
234 if ( an > 10. ) { an = 10.; }
235 }
236
237 // sample cosine of theta and not theta as in old versions
238 G4double random = G4UniformRand();
239 G4double cost;
240
241 if(an < 0.1) { cost = 1. - 2*random; }
242 else {
243 G4double exp2an = G4Exp(-2*an);
244 cost = 1. + G4Log(1-random*(1-exp2an))/an;
245 if(cost > 1.) { cost = 1.; }
246 else if(cost < -1.) {cost = -1.; }
247 }
248
249 G4double phi = CLHEP::twopi*G4UniformRand();
250
251 // Calculate the momentum magnitude of emitted fragment
252 G4double pmag =
253 std::sqrt(ekin*(ekin + 2.0*thePreFragment->GetNuclearMass()));
254
255 G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
256
257 theFinalMomentum.set(pmag*std::cos(phi)*sint,pmag*std::sin(phi)*sint,
258 pmag*cost);
259
260 // theta is the angle wrt the incident direction
261 G4ThreeVector theIncidentDirection = aFragment.GetMomentum().vect().unit();
262 theFinalMomentum.rotateUz(theIncidentDirection);
263 }
264
265 G4double G4PreCompoundEmission::rho(G4int p, G4int h, G4double gg,
266 G4double E, G4double Ef) const
267 {
268 // 25.02.2010 V.Ivanchenko added more protections
269 G4double Aph = (p*p + h*h + p - 3.0*h)/(4.0*gg);
270
271 if ( E - Aph < 0.0) { return 0.0; }
272
273 G4double logConst = (p+h)*G4Log(gg)
274 - g4calc->logfactorial(p+h-1) - g4calc->logfactorial(p)
275 - g4calc->logfactorial(h);
276
277 // initialise values using j=0
278
279 G4double t1=1;
280 G4double t2=1;
281 G4double logt3 = (p+h-1) * G4Log(E-Aph) + logConst;
282 const G4double logmax = 200.;
283 if(logt3 > logmax) { logt3 = logmax; }
284 G4double tot = G4Exp( logt3 );
285
286 // and now sum rest of terms
287 // 25.02.2010 V.Ivanchenko change while to for loop and cleanup
288 G4double Eeff = E - Aph;
289 for(G4int j=1; j<=h; ++j)
290 {
291 Eeff -= Ef;
292 if(Eeff < 0.0) { break; }
293 t1 *= -1.;
294 t2 *= static_cast<G4double>(h+1-j)/static_cast<G4double>(j);
295 logt3 = (p+h-1) * G4Log( Eeff) + logConst;
296 if(logt3 > logmax) { logt3 = logmax; }
297 tot += t1*t2*G4Exp(logt3);
298 }
299
300 return tot;
301 }
302