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

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

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