source: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundEmission.cc @ 1315

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 17.7 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: G4PreCompoundEmission.cc,v 1.28 2010/02/25 10:27:36 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
29//
30// -------------------------------------------------------------------
31//
32// GEANT4 Class file
33//
34//
35// File name:     G4PreCompoundEmission
36//
37// Author:         V.Lara
38//
39// Modified: 
40// 15.01.2010 J.M.Quesada  added protection against unphysical values of parameter an
41// 19.01.2010 V.Ivanchenko simplified computation of parameter an, sample cosTheta
42//                         instead of theta; protect all calls to sqrt
43//
44
45#include "G4PreCompoundEmission.hh"
46#include "G4PreCompoundParameters.hh"
47
48#include "G4PreCompoundEmissionFactory.hh"
49#include "G4HETCEmissionFactory.hh"
50
51const G4PreCompoundEmission & 
52G4PreCompoundEmission::operator=(const G4PreCompoundEmission &)
53{
54  throw G4HadronicException(__FILE__, __LINE__,
55                            "G4PreCompoundEmission::operator= meant to not be accessable");
56  return *this;
57}
58
59
60G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const
61{
62  return false;
63}
64
65G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const
66{
67  return true;
68}
69
70G4PreCompoundEmission::G4PreCompoundEmission()
71{
72  theFragmentsFactory = new G4PreCompoundEmissionFactory();
73  theFragmentsVector = new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
74}
75
76G4PreCompoundEmission::~G4PreCompoundEmission()
77{
78  if (theFragmentsFactory) delete theFragmentsFactory;
79  if (theFragmentsVector) delete theFragmentsVector;
80}
81
82void G4PreCompoundEmission::SetDefaultModel()
83{
84  if (theFragmentsFactory) delete theFragmentsFactory;
85  theFragmentsFactory = new G4PreCompoundEmissionFactory();
86  if (theFragmentsVector) 
87    {
88      theFragmentsVector->SetVector(theFragmentsFactory->GetFragmentVector());
89    }
90  else 
91    {
92      theFragmentsVector = new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
93    }
94  theFragmentsVector->ResetStage();
95  return;
96}
97
98void G4PreCompoundEmission::SetHETCModel()
99{
100  if (theFragmentsFactory) delete theFragmentsFactory;
101  theFragmentsFactory = new G4HETCEmissionFactory();
102  if (theFragmentsVector) 
103    {
104      theFragmentsVector->SetVector(theFragmentsFactory->GetFragmentVector());
105    }
106  else 
107    {
108      theFragmentsVector = new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
109    }
110  theFragmentsVector->ResetStage();
111  return;
112}
113
114G4ReactionProduct * G4PreCompoundEmission::PerformEmission(G4Fragment & aFragment)
115{
116#ifdef debug
117  G4Fragment InitialState(aFragment);
118#endif
119  // Choose a Fragment for emission
120  G4VPreCompoundFragment * theFragment = theFragmentsVector->ChooseFragment();
121  if (theFragment == 0)
122    {
123      G4cerr <<  "G4PreCompoundEmission::PerformEmission : I couldn't choose a fragment\n"
124             << "while trying to de-excite\n" 
125             << aFragment << '\n';
126      throw G4HadronicException(__FILE__, __LINE__, "");
127    }
128  // Kinetic Energy of emitted fragment
129  G4double KineticEnergyOfEmittedFragment = theFragment->GetKineticEnergy(aFragment);
130 
131  // Calculate the fragment momentum (three vector)
132  G4ThreeVector momentum = AngularDistribution(theFragment,aFragment,KineticEnergyOfEmittedFragment);
133 
134  // Mass of emittef fragment
135  G4double EmittedMass = theFragment->GetNuclearMass();
136 
137  // Now we can calculate the four momentum
138  // both options are valid and give the same result but 2nd one is faster
139  // G4LorentzVector EmittedMomentum(momentum,std::sqrt(momentum.mag2()+EmittedMass*EmittedMass));
140  G4LorentzVector EmittedMomentum(momentum,EmittedMass+KineticEnergyOfEmittedFragment);
141   
142  // Perform Lorentz boost
143  EmittedMomentum.boost(aFragment.GetMomentum().boostVector()); 
144
145  // Set emitted fragment momentum
146  theFragment->SetMomentum(EmittedMomentum);   
147
148
149  // NOW THE RESIDUAL NUCLEUS
150  // ------------------------
151   
152  // Now the residual nucleus.
153  // The energy conservation says that
154  G4double ResidualEcm = 
155    //    aFragment.GetGroundStateMass() + aFragment.GetExcitationEnergy() // initial energy in cm
156    aFragment.GetMomentum().m()
157    - (EmittedMass+KineticEnergyOfEmittedFragment); 
158
159  // Then the four momentum for residual is
160  G4LorentzVector RestMomentum(-momentum,ResidualEcm);
161  // This could save a Lorentz boost
162  // G4LorentzVector RestMomentum2(aFragment.GetMomentum()-EmittedMomentum);
163
164  // Just for test
165  // Excitation energy
166  //  G4double anU = ResidualEcm - theFragment->GetRestNuclearMass();
167  // This is equivalent
168  //  G4double anU = theFragment->GetMaximalKineticEnergy() - KineticEnergyOfEmittedFragment +
169  //    theFragment->GetCoulombBarrier();
170   
171  // check that Excitation energy is >= 0
172  G4double anU = RestMomentum.m()-theFragment->GetRestNuclearMass();
173  if (anU < 0.0) {
174    throw G4HadronicException(__FILE__, __LINE__, 
175                              "G4PreCompoundModel::DeExcite: Excitation energy less than 0!");
176  }
177     
178  // Update nucleus parameters:
179  // --------------------------
180
181  // Number of excitons
182  aFragment.SetNumberOfParticles(aFragment.GetNumberOfParticles()-
183                                 static_cast<G4int>(theFragment->GetA()));
184  // Number of charges
185  aFragment.SetNumberOfCharged(aFragment.GetNumberOfCharged()-
186                               static_cast<G4int>(theFragment->GetZ()));
187   
188  // Atomic number
189  aFragment.SetA(theFragment->GetRestA());
190   
191  // Charge
192  aFragment.SetZ(theFragment->GetRestZ());
193
194   
195  // Perform Lorentz boosts
196  RestMomentum.boost(aFragment.GetMomentum().boostVector());
197
198  // Update nucleus momentum
199  aFragment.SetMomentum(RestMomentum);
200       
201  // Create a G4ReactionProduct
202  G4ReactionProduct * MyRP = theFragment->GetReactionProduct();
203#ifdef PRECOMPOUND_TEST
204  MyRP->SetCreatorModel("G4PreCompoundModel");
205#endif
206#ifdef debug
207  CheckConservation(InitialState,aFragment,MyRP);
208#endif
209  return MyRP;
210}
211
212G4ThreeVector
213G4PreCompoundEmission::AngularDistribution(G4VPreCompoundFragment * theFragment,
214                                           const G4Fragment& aFragment,
215                                           const G4double kinEnergyOfEmittedFrag) const
216{
217  G4double p = aFragment.GetNumberOfParticles();
218  G4double h = aFragment.GetNumberOfHoles();
219  G4double U = aFragment.GetExcitationEnergy();
220
221  G4double ekin = std::max(0.0, kinEnergyOfEmittedFrag);
222       
223  // Emission particle separation energy
224  G4double Bemission = theFragment->GetBindingEnergy();
225
226  // Fermi energy
227  G4double Ef = G4PreCompoundParameters::GetAddress()->GetFermiEnergy();
228       
229  //
230  //  G4EvaporationLevelDensityParameter theLDP;
231  //  G4double g = (6.0/pi2)*aFragment.GetA()*
232
233  G4double g = (6.0/pi2)*aFragment.GetA()
234    *G4PreCompoundParameters::GetAddress()->GetLevelDensity();
235       
236  // Average exciton energy relative to bottom of nuclear well
237  G4double Eav = 2.0*p*(p+1.0)/((p+h)*g);
238       
239  // Excitation energy relative to the Fermi Level
240  G4double Uf = std::max(U - (p - h)*Ef , 0.0);
241  //  G4double Uf = U - KineticEnergyOfEmittedFragment - Bemission;
242
243  G4double w_num = rho(p+1,h,g,Uf,Ef);
244  G4double w_den = rho(p,h,g,Uf,Ef);
245  if (w_num > 0.0 && w_den > 0.0)
246    {
247      Eav *= (w_num/w_den);
248      Eav += - Uf/(p+h) + Ef;
249    }
250  else 
251    {
252      Eav = Ef;
253    }
254 
255
256  // VI + JMQ 19/01/2010 update computation of the parameter an
257  //
258  G4double an = 0.0;
259  G4double Eeff = ekin + Bemission + Ef;
260  if(ekin > DBL_MIN && Eeff > DBL_MIN) {
261
262    G4double zeta = std::max(1.0,9.3/std::sqrt(ekin/MeV));
263 
264    an = 3.0*std::sqrt((ProjEnergy+Ef)*Eeff)/(zeta*Eav);
265
266    G4int ne = aFragment.GetNumberOfExcitons() - 1;
267    if ( ne > 1 ) { an /= (G4double)ne; }
268                       
269    // protection of exponent
270    if ( an > 10. ) { an = 10.; }
271  }
272
273  // sample cosine of theta and not theta as in old versions 
274  G4double random = G4UniformRand();
275  G4double cost;
276 
277  if(an < 0.1) { cost = 1. - 2*random; }
278  else {
279    G4double exp2an = std::exp(-2*an);
280    cost = 1. + std::log(1-random*(1-exp2an))/an;
281    if(cost > 1.) { cost = 1.; }
282    else if(cost < -1.) {cost = -1.; }
283  } 
284
285  G4double phi = twopi*G4UniformRand();
286 
287  // Calculate the momentum magnitude of emitted fragment       
288  G4double pmag = std::sqrt(ekin*(ekin + 2.0*theFragment->GetNuclearMass()));
289 
290  G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
291
292  G4ThreeVector momentum(pmag*std::cos(phi)*sint,pmag*std::sin(phi)*sint,pmag*cost);
293  // theta is the angle wrt the incident direction
294  momentum.rotateUz(theIncidentDirection);
295
296  return momentum;
297}
298
299G4double G4PreCompoundEmission::rho(const G4double p, const G4double h, const G4double g, 
300                                    const G4double E, const G4double Ef) const
301{       
302  // 25.02.2010 V.Ivanchenko added more protections
303  G4double Aph   = (p*p + h*h + p - 3.0*h)/(4.0*g);
304  //  G4double alpha = (p*p + h*h)/(2.0*g);
305 
306  if ( E - Aph < 0.0) { return 0.0; }
307 
308  G4double logConst =  (p+h)*std::log(g) - logfactorial(p+h-1) - logfactorial(p) - logfactorial(h);
309
310  // initialise values using j=0
311
312  G4double t1=1;
313  G4double t2=1;
314  G4double logt3=(p+h-1) * std::log(E-Aph) + logConst;
315  const G4double logmax = 200.;
316  if(logt3 > logmax) { logt3 = logmax; }
317  G4double tot = std::exp( logt3 );
318
319  // and now sum rest of terms
320  // 25.02.2010 V.Ivanchenko change while to for loop and cleanup
321  G4double Eeff = E - Aph; 
322  for(G4int j=1; j<=h; ++j) 
323    {
324      Eeff -= Ef;
325      if(Eeff < 0.0) { break; }
326      t1 *= -1.;
327      t2 *= (G4double)(h+1-j)/(G4double)j;
328      logt3 = (p+h-1) * std::log( Eeff) + logConst;
329      if(logt3 > logmax) { logt3 = logmax; }
330      tot += t1*t2*std::exp(logt3);
331    }
332       
333  return tot;
334}
335
336G4double G4PreCompoundEmission::factorial(G4double a) const
337{
338  // Values of factorial function from 0 to 60
339  const G4int factablesize = 61;
340  static const G4double fact[factablesize] = 
341    {
342      1.0, // 0!
343      1.0, // 1!
344      2.0, // 2!
345      6.0, // 3!
346      24.0, // 4!
347      120.0, // 5!
348      720.0, // 6!
349      5040.0, // 7!
350      40320.0, // 8!
351      362880.0, // 9!
352      3628800.0, // 10!
353      39916800.0, // 11!
354      479001600.0, // 12!
355      6227020800.0, // 13!
356      87178291200.0, // 14!
357      1307674368000.0, // 15!
358      20922789888000.0, // 16!
359      355687428096000.0, // 17!
360      6402373705728000.0, // 18!
361      121645100408832000.0, // 19!
362      2432902008176640000.0, // 20!
363      51090942171709440000.0, // 21!
364      1124000727777607680000.0, // 22!
365      25852016738884976640000.0, // 23!
366      620448401733239439360000.0, // 24!
367      15511210043330985984000000.0, // 25!
368      403291461126605635584000000.0, // 26!
369      10888869450418352160768000000.0, // 27!
370      304888344611713860501504000000.0, // 28!
371      8841761993739701954543616000000.0, // 29!
372      265252859812191058636308480000000.0, // 30!
373      8222838654177922817725562880000000.0, // 31!
374      263130836933693530167218012160000000.0, // 32!
375      8683317618811886495518194401280000000.0, // 33!
376      295232799039604140847618609643520000000.0, // 34!
377      10333147966386144929666651337523200000000.0, // 35!
378      371993326789901217467999448150835200000000.0, // 36!
379      13763753091226345046315979581580902400000000.0, // 37!
380      523022617466601111760007224100074291200000000.0, // 38!
381      20397882081197443358640281739902897356800000000.0, // 39!
382      815915283247897734345611269596115894272000000000.0, // 40!
383      33452526613163807108170062053440751665152000000000.0, // 41!
384      1405006117752879898543142606244511569936384000000000.0, // 42!
385      60415263063373835637355132068513997507264512000000000.0, // 43!
386      2658271574788448768043625811014615890319638528000000000.0, // 44!
387      119622220865480194561963161495657715064383733760000000000.0, // 45!
388      5502622159812088949850305428800254892961651752960000000000.0, // 46!
389      258623241511168180642964355153611979969197632389120000000000.0, // 47!
390      12413915592536072670862289047373375038521486354677760000000000.0, // 48!
391      608281864034267560872252163321295376887552831379210240000000000.0, // 49!
392      30414093201713378043612608166064768844377641568960512000000000000.0, // 50!
393      1551118753287382280224243016469303211063259720016986112000000000000.0, // 51!
394      80658175170943878571660636856403766975289505440883277824000000000000.0, // 52!
395      4274883284060025564298013753389399649690343788366813724672000000000000.0, // 53!
396      230843697339241380472092742683027581083278564571807941132288000000000000.0, // 54!
397      12696403353658275925965100847566516959580321051449436762275840000000000000.0, // 55!
398      710998587804863451854045647463724949736497978881168458687447040000000000000.0, // 56!
399      40526919504877216755680601905432322134980384796226602145184481280000000000000.0, // 57!
400      2350561331282878571829474910515074683828862318181142924420699914240000000000000.0, // 58!
401      138683118545689835737939019720389406345902876772687432540821294940160000000000000.0, // 59!
402      8320987112741390144276341183223364380754172606361245952449277696409600000000000000.0  // 60!
403    };
404  //    fact[0] = 1;
405  //    for (G4int n = 1; n < 21; n++) {
406  //      fact[n] = fact[n-1]*static_cast<G4double>(n);
407  //    }
408  G4double result(0.0);
409  G4int ia = static_cast<G4int>(a);
410  if (ia < factablesize) 
411    {
412      result = fact[ia];
413    }
414  else
415    {
416      result = fact[factablesize-1];
417      for (G4int n = factablesize; n < ia+1; ++n)
418        {
419          result *= static_cast<G4double>(n);
420        }
421    }
422   
423    return result;
424}
425G4double G4PreCompoundEmission::logfactorial(G4double a) const
426{
427  // Values of logs of factorial function from 0 to 60
428
429  G4double result(0.0);
430  const G4int factablesize = 61;
431  const G4double halfLn2pi = 0.918938533;      // 0.5 log(2* pi)
432  static G4double logfact[factablesize];
433  static bool needinit=true;
434 
435  if (needinit)
436  {
437      needinit=false;
438      for ( G4int n=0; n < factablesize; ++n)
439      {
440         logfact[n]=std::log(factorial(n));
441      }
442  }
443
444  G4int ia = static_cast<G4int>(a);
445  if (ia < factablesize)
446  {
447      result = logfact[ia];
448  } else {
449      result = (ia+0.5)*std::log(G4double(ia)) - ia + halfLn2pi;
450  }
451   
452  return result;
453}
454
455#ifdef debug
456void G4PreCompoundEmission::CheckConservation(const G4Fragment & theInitialState,
457                                              const G4Fragment & theResidual,
458                                              G4ReactionProduct * theEmitted) const
459{
460  G4double ProductsEnergy = theEmitted->GetTotalEnergy() + theResidual.GetMomentum().e();
461  G4ThreeVector ProductsMomentum(theEmitted->GetMomentum()+theResidual.GetMomentum().vect());
462  G4int ProductsA = theEmitted->GetDefinition()->GetBaryonNumber() + theResidual.GetA();
463  G4int ProductsZ = theEmitted->GetDefinition()->GetPDGCharge() + theResidual.GetZ();
464
465  if (ProductsA != theInitialState.GetA()) {
466    G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl;
467    G4cout << "G4PreCompoundEmission.cc: Barionic Number Conservation"
468           << G4endl; 
469    G4cout << "Initial A = " << theInitialState.GetA() 
470           << "   Fragments A = " << ProductsA << "   Diference --> " 
471           << theInitialState.GetA() - ProductsA << G4endl;
472  }
473  if (ProductsZ != theInitialState.GetZ()) {
474    G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
475    G4cout << "G4PreCompoundEmission.cc: Charge Conservation test"
476           << G4endl; 
477    G4cout << "Initial Z = " << theInitialState.GetZ() 
478           << "   Fragments Z = " << ProductsZ << "   Diference --> " 
479           << theInitialState.GetZ() - ProductsZ << G4endl;
480  }
481  if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 10.0*eV) {
482    G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
483    G4cout << "G4PreCompoundEmission.cc: Energy Conservation test" 
484           << G4endl; 
485    G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
486           << "   Fragments E = " << ProductsEnergy/MeV  << " MeV   Diference --> " 
487           << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl;
488  } 
489  if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 10.0*eV || 
490      std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 10.0*eV ||
491      std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 10.0*eV) {
492    G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
493    G4cout << "G4PreCompoundEmission.cc: Momentum Conservation test" 
494           << G4endl; 
495    G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
496           << "   Fragments P = " << ProductsMomentum  << " MeV   Diference --> " 
497           << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl;
498  }
499  return;
500}
501
502#endif
Note: See TracBrowser for help on using the repository browser.