ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PreCompoundEmission.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PreCompoundEmission.cc
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"
53 #include "G4HETCEmissionFactory.hh"
54 #include "G4HadronicException.hh"
55 #include "G4NuclearLevelData.hh"
56 #include "G4DeexPrecoParameters.hh"
57 
59 {
66  fFermiEnergy = param->GetFermiEnergy();
68 }
69 
71 {
72  delete theFragmentsFactory;
73  delete theFragmentsVector;
74 }
75 
77 {
80  if (theFragmentsVector) {
82  } else {
85  }
86 }
87 
89 {
92  if (theFragmentsVector) {
94  } else {
97  }
98 }
99 
102 {
103  // Choose a Fragment for emission
104  G4VPreCompoundFragment * thePreFragment =
106  if (thePreFragment == nullptr)
107  {
108  G4cout << "G4PreCompoundEmission::PerformEmission : "
109  << "I couldn't choose a fragment\n"
110  << "while trying to de-excite\n"
111  << aFragment << G4endl;
112  throw G4HadronicException(__FILE__, __LINE__, "");
113  }
114 
115  //G4cout << "Chosen fragment: " << G4endl;
116  //G4cout << *thePreFragment << G4endl;
117 
118  // Kinetic Energy of emitted fragment
119  G4double kinEnergy = thePreFragment->SampleKineticEnergy(aFragment);
120  kinEnergy = std::max(kinEnergy, 0.0);
121 
122  // Calculate the fragment momentum (three vector)
124  AngularDistribution(thePreFragment,aFragment,kinEnergy);
125  } else {
126  G4double pmag =
127  std::sqrt(kinEnergy*(kinEnergy + 2.0*thePreFragment->GetNuclearMass()));
129  }
130 
131  // Mass of emittef fragment
132  G4double EmittedMass = thePreFragment->GetNuclearMass();
133  // Now we can calculate the four momentum
134  // both options are valid and give the same result but 2nd one is faster
135  G4LorentzVector Emitted4Momentum(theFinalMomentum,EmittedMass + kinEnergy);
136 
137  // Perform Lorentz boost
138  G4LorentzVector Rest4Momentum = aFragment.GetMomentum();
139  Emitted4Momentum.boost(Rest4Momentum.boostVector());
140 
141  // Set emitted fragment momentum
142  thePreFragment->SetMomentum(Emitted4Momentum);
143 
144  // NOW THE RESIDUAL NUCLEUS
145  // ------------------------
146 
147  Rest4Momentum -= Emitted4Momentum;
148 
149  // Update nucleus parameters:
150  // --------------------------
151 
152  // Z and A
153  aFragment.SetZandA_asInt(thePreFragment->GetRestZ(),
154  thePreFragment->GetRestA());
155 
156  // Number of excitons
157  aFragment.SetNumberOfParticles(aFragment.GetNumberOfParticles()-
158  thePreFragment->GetA());
159  // Number of charges
160  aFragment.SetNumberOfCharged(aFragment.GetNumberOfCharged()-
161  thePreFragment->GetZ());
162 
163  // Update nucleus momentum
164  // A check on consistence of Z, A, and mass will be performed
165  aFragment.SetMomentum(Rest4Momentum);
166 
167  // Create a G4ReactionProduct
168  G4ReactionProduct * MyRP = thePreFragment->GetReactionProduct();
169 
170  return MyRP;
171 }
172 
174  G4VPreCompoundFragment* thePreFragment,
175  const G4Fragment& aFragment,
176  G4double ekin)
177 {
178  G4int p = aFragment.GetNumberOfParticles();
179  G4int h = aFragment.GetNumberOfHoles();
180  G4double U = aFragment.GetExcitationEnergy();
181 
182  // Emission particle separation energy
183  G4double Bemission = thePreFragment->GetBindingEnergy();
184 
185  G4double gg = (6.0/pi2)*fNuclData->GetLevelDensity(aFragment.GetZ_asInt(),
186  aFragment.GetA_asInt(),U);
187 
188  // Average exciton energy relative to bottom of nuclear well
189  G4double Eav = 2*p*(p+1)/((p+h)*gg);
190 
191  // Excitation energy relative to the Fermi Level
192  G4double Uf = std::max(U - (p - h)*fFermiEnergy , 0.0);
193  // G4double Uf = U - KineticEnergyOfEmittedFragment - Bemission;
194 
195  G4double w_num = rho(p+1, h, gg, Uf, fFermiEnergy);
196  G4double w_den = rho(p, h, gg, Uf, fFermiEnergy);
197  if (w_num > 0.0 && w_den > 0.0)
198  {
199  Eav *= (w_num/w_den);
200  Eav += - Uf/(p+h) + fFermiEnergy;
201  }
202  else
203  {
204  Eav = fFermiEnergy;
205  }
206 
207  // VI + JMQ 19/01/2010 update computation of the parameter an
208  //
209  G4double an = 0.0;
210  G4double Eeff = ekin + Bemission + fFermiEnergy;
211  if(ekin > DBL_MIN && Eeff > DBL_MIN) {
212 
213  G4double zeta = std::max(1.0,9.3/std::sqrt(ekin/MeV));
214 
215  // This should be the projectile energy. If I would know which is
216  // the projectile (proton, neutron) I could remove the binding energy.
217  // But, what happens if INC precedes precompound? This approximation
218  // seems to work well enough
219  G4double ProjEnergy = aFragment.GetExcitationEnergy();
220 
221  an = 3*std::sqrt((ProjEnergy+fFermiEnergy)*Eeff)/(zeta*Eav);
222 
223  G4int ne = aFragment.GetNumberOfExcitons() - 1;
224  if ( ne > 1 ) { an /= (G4double)ne; }
225 
226  // protection of exponent
227  if ( an > 10. ) { an = 10.; }
228  }
229 
230  // sample cosine of theta and not theta as in old versions
231  G4double random = G4UniformRand();
232  G4double cost;
233 
234  if(an < 0.1) { cost = 1. - 2*random; }
235  else {
236  G4double exp2an = G4Exp(-2*an);
237  cost = 1. + G4Log(1-random*(1-exp2an))/an;
238  if(cost > 1.) { cost = 1.; }
239  else if(cost < -1.) {cost = -1.; }
240  }
241 
243 
244  // Calculate the momentum magnitude of emitted fragment
245  G4double pmag =
246  std::sqrt(ekin*(ekin + 2.0*thePreFragment->GetNuclearMass()));
247 
248  G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
249 
250  theFinalMomentum.set(pmag*std::cos(phi)*sint,pmag*std::sin(phi)*sint,
251  pmag*cost);
252 
253  // theta is the angle wrt the incident direction
254  G4ThreeVector theIncidentDirection = aFragment.GetMomentum().vect().unit();
255  theFinalMomentum.rotateUz(theIncidentDirection);
256 }
257 
259  G4double E, G4double Ef) const
260 {
261  // 25.02.2010 V.Ivanchenko added more protections
262  G4double Aph = (p*p + h*h + p - 3.0*h)/(4.0*gg);
263  // G4double alpha = (p*p + h*h)/(2.0*gg);
264 
265  if ( E - Aph < 0.0) { return 0.0; }
266 
267  G4double logConst = (p+h)*G4Log(gg)
268  - g4calc->logfactorial(p+h-1) - g4calc->logfactorial(p)
269  - g4calc->logfactorial(h);
270 
271  // initialise values using j=0
272 
273  G4double t1=1;
274  G4double t2=1;
275  G4double logt3 = (p+h-1) * G4Log(E-Aph) + logConst;
276  const G4double logmax = 200.;
277  if(logt3 > logmax) { logt3 = logmax; }
278  G4double tot = G4Exp( logt3 );
279 
280  // and now sum rest of terms
281  // 25.02.2010 V.Ivanchenko change while to for loop and cleanup
282  G4double Eeff = E - Aph;
283  for(G4int j=1; j<=h; ++j)
284  {
285  Eeff -= Ef;
286  if(Eeff < 0.0) { break; }
287  t1 *= -1.;
288  t2 *= (G4double)(h+1-j)/(G4double)j;
289  logt3 = (p+h-1) * G4Log( Eeff) + logConst;
290  if(logt3 > logmax) { logt3 = logmax; }
291  tot += t1*t2*G4Exp(logt3);
292  }
293 
294  return tot;
295 }