source: trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundModel.cc @ 1358

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

update ti head
File size: 15.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// $Id: G4PreCompoundModel.cc,v 1.25 2010/10/11 13:54:59 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-ref-09 $
28//
29// by V. Lara
30//
31// Modified:
32// 01.04.2008 J.M.Quesada Several changes. Soft cut-off switched off.
33// 01.05.2008 J.M.Quesada Protection against non-physical preeq.
34//                        transitional regime has been set
35// 03.09.2008 J.M.Quesada for external choice of inverse cross section option
36// 06.09.2008 J.M.Quesada Also external choices have been added for:
37//                      - superimposed Coulomb barrier (useSICB=true)
38//                      - "never go back"  hipothesis (useNGB=true)
39//                      - soft cutoff from preeq. to equlibrium (useSCO=true)
40//                      - CEM transition probabilities (useCEMtr=true) 
41// 20.08.2010 V.Ivanchenko Cleanup of the code:
42//                      - integer Z and A;
43//                      - emission and transition classes created at initialisation
44//                      - options are set at initialisation
45//                      - do not use copy-constructors for G4Fragment 
46
47#include "G4PreCompoundModel.hh"
48#include "G4PreCompoundEmission.hh"
49#include "G4PreCompoundTransitions.hh"
50#include "G4GNASHTransitions.hh"
51#include "G4ParticleDefinition.hh"
52#include "G4Proton.hh"
53#include "G4Neutron.hh"
54
55#include "G4NucleiProperties.hh"
56#include "G4PreCompoundParameters.hh"
57#include "Randomize.hh"
58#include "G4DynamicParticle.hh"
59#include "G4ParticleTypes.hh"
60#include "G4ParticleTable.hh"
61#include "G4LorentzVector.hh"
62
63#ifdef PRECOMPOUND_TEST
64G4Fragment G4PreCompoundModel::theInitialFragmentForTest;
65std::vector<G4String*> G4PreCompoundModel::theCreatorModels;
66#endif
67
68G4PreCompoundModel::G4PreCompoundModel(G4ExcitationHandler * const value) 
69  : G4VPreCompoundModel(value), useHETCEmission(false), useGNASHTransition(false), 
70    OPTxs(3), useSICB(false), useNGB(false), useSCO(false), useCEMtr(true) 
71{
72  theParameters = G4PreCompoundParameters::GetAddress();
73
74  theEmission = new G4PreCompoundEmission();
75  if(useHETCEmission) { theEmission->SetHETCModel(); }
76  else { theEmission->SetDefaultModel(); }
77  theEmission->SetOPTxs(OPTxs);
78  theEmission->UseSICB(useSICB);
79
80  if(useGNASHTransition) { theTransition = new G4GNASHTransitions; }
81  else { theTransition = new G4PreCompoundTransitions(); }
82  theTransition->UseNGB(useNGB);
83  theTransition->UseCEMtr(useCEMtr);
84
85  proton = G4Proton::Proton();
86  neutron = G4Neutron::Neutron();
87}
88
89G4PreCompoundModel::~G4PreCompoundModel() 
90{
91  delete theEmission;
92  delete theTransition;
93}
94
95/////////////////////////////////////////////////////////////////////////////////////////
96
97G4HadFinalState* G4PreCompoundModel::ApplyYourself(const G4HadProjectile & thePrimary,
98                                                   G4Nucleus & theNucleus)
99{ 
100  const G4ParticleDefinition* primary = thePrimary.GetDefinition();
101  if(primary != neutron && primary != proton) {
102    std::ostringstream errOs;
103    errOs << "BAD primary type in G4PreCompoundModel: " 
104          << primary->GetParticleName() <<G4endl;
105    throw G4HadronicException(__FILE__, __LINE__, errOs.str());
106  }
107  G4int Zp = 0;
108  G4int Ap = 1;
109  if(primary == proton) { Zp = 1; }
110
111  G4int A = theNucleus.GetA_asInt();
112  G4int Z = theNucleus.GetZ_asInt();
113
114  // 4-Momentum
115  G4LorentzVector p = thePrimary.Get4Momentum();
116  G4double mass = G4NucleiProperties::GetNuclearMass(A, Z);
117  p += G4LorentzVector(0.0,0.0,0.0,mass);
118
119  // prepare fragment
120  G4Fragment anInitialState(A + Ap, Z + Zp, p);
121  anInitialState.SetNumberOfParticles(2);
122  anInitialState.SetNumberOfHoles(1);
123  anInitialState.SetNumberOfCharged(1);
124  anInitialState.SetCreationTime(thePrimary.GetGlobalTime());
125  /*
126  // Number of Excited Particles
127  anInitialState.SetNumberOfParticles(1+thePrimary.GetDefinition()->GetBaryonNumber());
128 
129  // Number of Charged Excited Particles
130  // JMQ/AH modify number of charged particles with probability of the Z/A ratio of the nucleus:
131  //  if(G4UniformRand() <= aZ/anA) BUG! - integer arithmetic
132  if(G4UniformRand() <= (static_cast<G4double>(aZ))/(static_cast<G4double>(anA)))
133      anInitialState.SetNumberOfCharged(static_cast<G4int>(thePrimary.GetDefinition()->GetPDGCharge()+.01) + 1);
134  else
135      anInitialState.SetNumberOfCharged(static_cast<G4int>(thePrimary.GetDefinition()->GetPDGCharge()+.01));
136   
137//AH     anInitialState.SetNumberOfCharged(static_cast<G4int>(thePrimary.GetDefinition()->GetPDGCharge()+.01) +
138//AH                                static_cast<G4int>(0.5+G4UniformRand()));
139
140  // Number of Holes
141  anInitialState.SetNumberOfHoles(1);
142  */
143 
144  // call excitation handler
145  //  const G4Fragment aFragment(anInitialState);
146  G4ReactionProductVector * result = DeExcite(anInitialState);
147
148#ifdef PRECOMPOUND_TEST
149  for (std::vector<G4String*>::iterator icm = theCreatorModels.begin(); 
150       icm != theCreatorModels.end(); ++icm )
151    {
152      delete (*icm);
153    }
154  theCreatorModels.clear();
155#endif
156  // fill particle change
157  theResult.Clear();
158  theResult.SetStatusChange(stopAndKill);
159  for(G4ReactionProductVector::iterator i= result->begin(); i != result->end(); ++i)
160    {
161      G4DynamicParticle * aNew = 
162        new G4DynamicParticle((*i)->GetDefinition(),
163                              (*i)->GetTotalEnergy(),
164                              (*i)->GetMomentum());
165#ifdef PRECOMPOUND_TEST
166      theCreatorModels.push_back(new G4String((*i)->GetCreatorModel()));
167#endif
168      delete (*i);
169      theResult.AddSecondary(aNew);
170    }
171  delete result;
172 
173  //return the filled particle change
174  return &theResult;
175}
176
177/////////////////////////////////////////////////////////////////////////////////////////
178
179G4ReactionProductVector* G4PreCompoundModel::DeExcite(G4Fragment& aFragment)
180{
181  G4ReactionProductVector * Result = new G4ReactionProductVector;
182  G4double Eex = aFragment.GetExcitationEnergy();
183  G4int A = aFragment.GetA_asInt(); 
184
185  //G4cout << "### G4PreCompoundModel::DeExcite" << G4endl;
186  //G4cout << aFragment << G4endl;
187 
188  // Perform Equilibrium Emission
189  if (A < 5 || Eex < keV /*|| Eex > 3.*MeV*A*/) {
190    PerformEquilibriumEmission(aFragment, Result);
191    return Result;
192  }
193 
194  // main loop 
195  for (;;) {
196   
197    //fragment++;
198    //G4cout<<"-------------------------------------------------------------------"<<G4endl;
199    //G4cout<<"Fragment number .. "<<fragment<<G4endl;
200   
201    // Initialize fragment according with the nucleus parameters
202    theEmission->Initialize(aFragment);
203   
204    G4double g = (6.0/pi2)*aFragment.GetA_asInt()*theParameters->GetLevelDensity();
205   
206    G4int EquilibriumExcitonNumber = 
207      static_cast<G4int>(std::sqrt(2.0*g*aFragment.GetExcitationEnergy())+ 0.5);
208    //   
209    //    G4cout<<"Neq="<<EquilibriumExcitonNumber<<G4endl;
210    //
211    // J. M. Quesada (Jan. 08)  equilibrium hole number could be used as preeq.
212    // evap. delimiter (IAEA report)
213   
214    // Loop for transitions, it is performed while there are preequilibrium transitions.
215    G4bool ThereIsTransition = false;
216   
217    //        G4cout<<"----------------------------------------"<<G4endl;
218    //        G4double NP=aFragment.GetNumberOfParticles();
219    //        G4double NH=aFragment.GetNumberOfHoles();
220    //        G4double NE=aFragment.GetNumberOfExcitons();
221    //        G4cout<<" Ex. Energy="<<aFragment.GetExcitationEnergy()<<G4endl;
222    //        G4cout<<"N. excitons="<<NE<<"  N. Part="<<NP<<"N. Holes ="<<NH<<G4endl;
223    //G4int transition=0;
224    do {
225      //transition++;
226      //G4cout<<"transition number .."<<transition<<G4endl;
227      //G4cout<<" n ="<<aFragment.GetNumberOfExcitons()<<G4endl;
228      G4bool go_ahead = false;
229      // soft cutoff criterium as an "ad-hoc" solution to force increase in  evaporation 
230      //       G4double test = static_cast<G4double>(aFragment.GetNumberOfHoles());
231      G4int test = aFragment.GetNumberOfExcitons();
232      if (test < EquilibriumExcitonNumber) { go_ahead=true; }
233
234      //J. M. Quesada (Apr. 08): soft-cutoff switched off by default
235      if (useSCO) {
236        if (test < EquilibriumExcitonNumber)
237          {
238            G4double x = G4double(test)/G4double(EquilibriumExcitonNumber) - 1; 
239            if( G4UniformRand() < 1.0 -  std::exp(-x*x/0.32) ) { go_ahead = true; }
240            /*
241              test = test*test;
242              test /= 0.32;
243              test = 1.0 - std::exp(-test);
244              go_ahead = (G4UniformRand() < test);
245            */
246          }
247      } 
248       
249      // JMQ: WARNING:  CalculateProbability MUST be called prior to Get methods !!
250      // (O values would be returned otherwise)
251      G4double TotalTransitionProbability = 
252        theTransition->CalculateProbability(aFragment);
253      G4double P1 = theTransition->GetTransitionProb1();
254      G4double P2 = theTransition->GetTransitionProb2();
255      G4double P3 = theTransition->GetTransitionProb3();
256      //G4cout<<"P1="<<P1<<" P2="<<P2<<"  P3="<<P3<<G4endl;
257     
258      //J.M. Quesada (May. 08). Physical criterium (lamdas)  PREVAILS over
259      //                        approximation (critical exciton number)
260      if(P1 <= P2+P3) { go_ahead = false; }
261       
262      if (go_ahead &&  aFragment.GetA_asInt() > 4) 
263        {
264                               
265          G4double TotalEmissionProbability = 
266            theEmission->GetTotalProbability(aFragment);
267          //
268          //  G4cout<<"TotalEmissionProbability="<<TotalEmissionProbability<<G4endl;
269          //
270          // Check if number of excitons is greater than 0
271          // else perform equilibrium emission
272          if (aFragment.GetNumberOfExcitons() <= 0) 
273            {
274              PerformEquilibriumEmission(aFragment,Result);
275              return Result;
276            }
277           
278          //J.M.Quesada (May 08) this has already been done in order to decide 
279          //                     what to do (preeq-eq)
280          // Sum of all probabilities
281          G4double TotalProbability = TotalEmissionProbability
282            + TotalTransitionProbability;
283           
284          // Select subprocess
285          if (TotalProbability*G4UniformRand() > TotalEmissionProbability) 
286            {
287              // It will be transition to state with a new number of excitons
288              ThereIsTransition = true;         
289              // Perform the transition
290              theTransition->PerformTransition(aFragment);
291            } 
292          else 
293            {
294              // It will be fragment emission
295              ThereIsTransition = false;
296              Result->push_back(theEmission->PerformEmission(aFragment));
297            }
298        } 
299      else 
300        {
301          PerformEquilibriumEmission(aFragment,Result);
302          return Result;
303        }
304    } while (ThereIsTransition);   // end of do loop
305  } // end of for (;;) loop
306  return Result;
307}
308
309/////////////////////////////////////////////////////////////////////////////////////////
310//       Initialisation
311/////////////////////////////////////////////////////////////////////////////////////////
312
313void G4PreCompoundModel::UseHETCEmission() 
314{ 
315  useHETCEmission = true; 
316  theEmission->SetHETCModel();
317}
318
319void G4PreCompoundModel::UseDefaultEmission() 
320{ 
321  useHETCEmission = false; 
322  theEmission->SetDefaultModel();
323}
324
325void G4PreCompoundModel::UseGNASHTransition() { 
326  useGNASHTransition = true; 
327  delete theTransition;
328  theTransition = new G4GNASHTransitions;
329  theTransition->UseNGB(useNGB);
330  theTransition->UseCEMtr(useCEMtr);
331}
332
333void G4PreCompoundModel::UseDefaultTransition() { 
334  useGNASHTransition = false; 
335  delete theTransition;
336  theTransition = new G4PreCompoundTransitions();
337  theTransition->UseNGB(useNGB);
338  theTransition->UseCEMtr(useCEMtr);
339}
340
341void G4PreCompoundModel::SetOPTxs(G4int opt) 
342{ 
343  OPTxs = opt; 
344  theEmission->SetOPTxs(OPTxs);
345}
346
347void G4PreCompoundModel::UseSICB() 
348{ 
349  useSICB = true; 
350  theEmission->UseSICB(useSICB);
351}
352
353void G4PreCompoundModel::UseNGB() 
354{ 
355  useNGB = true; 
356}
357
358void G4PreCompoundModel::UseSCO() 
359{ 
360  useSCO = true; 
361}
362
363void G4PreCompoundModel::UseCEMtr() 
364{ 
365  useCEMtr = true; 
366}
367
368/////////////////////////////////////////////////////////////////////////////////////////
369
370#ifdef debug
371void G4PreCompoundModel::CheckConservation(const G4Fragment & theInitialState,
372                                           const G4Fragment & aFragment,
373                                           G4ReactionProductVector * Result) const
374{
375  G4double ProductsEnergy = aFragment.GetMomentum().e();
376  G4ThreeVector ProductsMomentum = aFragment.GetMomentum();
377  G4int ProductsA = static_cast<G4int>(aFragment.GetA());
378  G4int ProductsZ = static_cast<G4int>(aFragment.GetZ());
379  for (G4ReactionProductVector::iterator h = Result->begin(); 
380       h != Result->end(); ++h) 
381    {
382      ProductsEnergy += (*h)->GetTotalEnergy();
383      ProductsMomentum += (*h)->GetMomentum();
384      ProductsA += static_cast<G4int>((*h)->GetDefinition()->GetBaryonNumber());
385      ProductsZ += static_cast<G4int>((*h)->GetDefinition()->GetPDGCharge());
386    }
387
388  if (ProductsA != theInitialState.GetA()) 
389    {
390      G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!\n"
391             << "G4PreCompoundModel.cc: Barionic Number Conservation test for just preequilibrium fragments\n" 
392             << "Initial A = " << theInitialState.GetA() 
393             << "   Fragments A = " << ProductsA << "   Diference --> " 
394             << theInitialState.GetA() - ProductsA << '\n';
395    }
396  if (ProductsZ != theInitialState.GetZ()) 
397    {
398      G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!\n"
399             << "G4PreCompoundModel.cc: Charge Conservation test for just preequilibrium fragments\n" 
400             << "Initial Z = " << theInitialState.GetZ() 
401             << "   Fragments Z = " << ProductsZ << "   Diference --> " 
402             << theInitialState.GetZ() - ProductsZ << '\n';
403    }
404  if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) 
405    {
406      G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!\n" 
407             << "G4PreCompoundModel.cc: Energy Conservation test for just preequilibrium fragments\n" 
408             << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
409             << "   Fragments E = " << ProductsEnergy/MeV  << " MeV   Diference --> " 
410             << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV\n";
411    } 
412  if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV || 
413      std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV ||
414      std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) 
415    {
416      G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!\n"
417             << "G4PreCompoundModel.cc: Momentum Conservation test for just preequilibrium fragments\n" 
418             << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
419             << "   Fragments P = " << ProductsMomentum  << " MeV   Diference --> " 
420             << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV\n";
421    }
422  return;
423}
424
425#endif
426
427
428
Note: See TracBrowser for help on using the repository browser.