source: trunk/source/processes/hadronic/models/de_excitation/handler/src/G4ExcitationHandler.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

geant4 tag 9.4
File size: 19.9 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: G4ExcitationHandler.cc,v 1.40 2010/11/17 16:20:38 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
28//
29// Hadronic Process: Nuclear De-excitations
30// by V. Lara (May 1998)
31//
32//
33// Modified:
34// (30 June 1998) by V. Lara:
35//      -Modified the Transform method for use G4ParticleTable and
36//       therefore G4IonTable. It makes possible to convert all kind
37//       of fragments (G4Fragment) produced in deexcitation to
38//       G4DynamicParticle
39//      -It uses default algorithms for:
40//              Evaporation: G4Evaporation
41//              MultiFragmentation: G4StatMF
42//              Fermi Breakup model: G4FermiBreakUp
43// (24 Jul 2008) by M. A. Cortes Giraldo:
44//      -Max Z,A for Fermi Break-Up turns to 9,17 by default
45//      -BreakItUp() reorganised and bug in Evaporation loop fixed
46//      -Transform() optimised
47// (September 2008) by J. M. Quesada. External choices have been added for :
48//      -inverse cross section option (default OPTxs=3)
49//      -superimposed Coulomb barrier (if useSICB is set true, by default it is false)
50// (September 2009) by J. M. Quesada:
51//      -according to Igor Pshenichnov, SMM will be applied (just in case) only once.
52// (27 Nov 2009) by V.Ivanchenko:
53//      -cleanup the logic, reduce number internal vectors, fixed memory leak.
54// (11 May 2010) by V.Ivanchenko:
55//      -FermiBreakUp activated, used integer Z and A, used BreakUpFragment method for
56//       final photon deexcitation; used check on adundance of a fragment, decay
57//       unstable fragments with A <5
58//
59
60#include "G4ExcitationHandler.hh"
61#include "globals.hh"
62#include "G4LorentzVector.hh"
63#include "G4NistManager.hh"
64#include "G4ParticleTable.hh"
65#include "G4IonTable.hh"
66#include "G4IonConstructor.hh"
67#include "G4ParticleTypes.hh"
68
69#include <list>
70
71G4ExcitationHandler::G4ExcitationHandler():
72  // JMQ 160909 Fermi BreakUp & MultiFrag are on by default
73  // This is needed for activation of such models when G4BinaryLightIonReaction is used
74  // since no interface (for external activation via macro input file) is still available.
75  maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(4.0*GeV),
76  //  maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(3.0*MeV),
77  minExcitation(CLHEP::keV),
78  MyOwnEvaporationClass(true), MyOwnMultiFragmentationClass(true),MyOwnFermiBreakUpClass(true),
79  MyOwnPhotonEvaporationClass(true),OPTxs(3),useSICB(false)
80{                                                                         
81  theTableOfIons = G4ParticleTable::GetParticleTable()->GetIonTable();
82 
83  theEvaporation = new G4Evaporation;
84  theMultiFragmentation = new G4StatMF;
85  theFermiModel = new G4FermiBreakUp;
86  thePhotonEvaporation = new G4PhotonEvaporation;
87  SetParameters();
88}
89
90G4ExcitationHandler::~G4ExcitationHandler()
91{
92  if (MyOwnEvaporationClass) delete theEvaporation;
93  if (MyOwnMultiFragmentationClass) delete theMultiFragmentation;
94  if (MyOwnFermiBreakUpClass) delete theFermiModel;
95  if (MyOwnPhotonEvaporationClass) delete thePhotonEvaporation;
96}
97
98void G4ExcitationHandler::SetParameters()
99{
100  //for inverse cross section choice
101  theEvaporation->SetOPTxs(OPTxs);
102  //for the choice of superimposed Coulomb Barrier for inverse cross sections
103  theEvaporation->UseSICB(useSICB);
104  theEvaporation->Initialise();
105}
106
107G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) const
108{       
109 
110  // Variables existing until end of method
111  G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
112
113  G4FragmentVector * theTempResult = 0;      // pointer which receives temporal results
114  std::list<G4Fragment*> theEvapList;        // list to apply Evaporation or Fermi Break-Up
115  std::list<G4Fragment*> theEvapStableList;  // list to apply PhotonEvaporation
116  std::list<G4Fragment*> theResults;         // list to store final result
117  std::list<G4Fragment*>::iterator iList;
118  //
119  //G4cout << "@@@@@@@@@@ Start G4Excitation Handler @@@@@@@@@@@@@" << G4endl; 
120  //G4cout << theInitialState << G4endl; 
121 
122  // Variables to describe the excited configuration
123  G4double exEnergy = theInitialState.GetExcitationEnergy();
124  G4int A = theInitialState.GetA_asInt();
125  G4int Z = theInitialState.GetZ_asInt();
126
127  G4NistManager* nist = G4NistManager::Instance();
128 
129  // JMQ 150909:  first step in de-excitation chain (SMM will be used only here)
130 
131  // In case A <= 1 the fragment will not perform any nucleon emission
132  if (A <= 1)
133    {
134      theResults.push_back( theInitialStatePtr );
135    }
136  // check if a fragment is stable
137  else if(exEnergy < minExcitation && nist->GetIsotopeAbundance(Z, A) > 0.0)
138    {
139      theResults.push_back( theInitialStatePtr );
140    } 
141  else  // In all cases apply once theFermiModel, theMultiFragmentation or theEvaporation
142    {     
143      // JMQ 150909: first step in de-excitation is treated separately
144      // Fragments after the first step are stored in theEvapList
145      // Statistical Multifragmentation will take place only once
146      //
147      // Fermi Break-Up always first
148      if((A < 5) || (A<GetMaxA() && Z<GetMaxZ())) 
149        {
150          theTempResult = theFermiModel->BreakItUp(theInitialState);
151          // fragment was not decaied try to evaporate
152          if(1 == theTempResult->size()) {
153            if((*theTempResult)[0] !=  theInitialStatePtr) { 
154              delete (*theTempResult)[0]; 
155            }
156            delete theTempResult;
157            theTempResult = theEvaporation->BreakItUp(theInitialState);
158          }
159        }
160      else   if (exEnergy>GetMinE()*A) 
161        {
162          theTempResult = theMultiFragmentation->BreakItUp(theInitialState);
163          // fragment was not decaied try to evaporate
164          if(1 == theTempResult->size()) {
165            if((*theTempResult)[0] !=  theInitialStatePtr) { 
166              delete (*theTempResult)[0]; 
167            }
168            delete theTempResult;
169            theTempResult = theEvaporation->BreakItUp(theInitialState);
170          }
171        }
172      else 
173        {
174          theTempResult = theEvaporation->BreakItUp(theInitialState);
175        }
176     
177      G4bool deletePrimary = true;
178      G4int nsec=theTempResult->size();
179      if(nsec > 0) 
180        {     
181          // Sort out secondary fragments
182          G4FragmentVector::iterator j;
183          for (j = theTempResult->begin(); j != theTempResult->end(); ++j)
184            { 
185              if((*j) == theInitialStatePtr) { deletePrimary = false; }
186              A = (*j)->GetA_asInt(); 
187
188              if(A <= 1) { theResults.push_back(*j); }   // gamma, p, n
189              else if(1 == nsec) { theEvapStableList.push_back(*j); } // evaporation is not possible 
190              else  // Analyse fragment
191                {
192                  G4double exEnergy = (*j)->GetExcitationEnergy();
193                  if(exEnergy < minExcitation) {
194                    Z = (*j)->GetZ_asInt(); 
195                    if(nist->GetIsotopeAbundance(Z, A) > 0.0) { 
196                      theResults.push_back(*j); // stable fragment
197                    } else {   
198                      theEvapList.push_back(*j); // unstable cold fragment
199                    }
200                  } else { theEvapList.push_back(*j); } // hot fragment
201                }
202            }
203        }
204      if( deletePrimary ) { delete theInitialStatePtr; }
205      delete theTempResult;
206    }
207  //
208  // JMQ 150909: Further steps in de-excitation chain follow ..
209     
210  //G4cout << "## After first step " << theEvapList.size() << " for evap;  "
211  //     << theEvapStableList.size() << " for photo-evap; "
212  //     << theResults.size() << " results. " << G4endl;
213
214  // ------------------------------
215  // De-excitation loop
216  // ------------------------------
217     
218  for (iList = theEvapList.begin(); iList != theEvapList.end(); ++iList)
219    {
220      A = (*iList)->GetA_asInt(); 
221      Z = (*iList)->GetZ_asInt();
222         
223      // Fermi Break-Up
224      if ((A < 5) || (A < GetMaxA() && Z < GetMaxZ())) 
225        {
226          theTempResult = theFermiModel->BreakItUp(*(*iList));
227          // fragment was not decaied try to evaporate
228          if(1 == theTempResult->size()) {
229            if((*theTempResult)[0] !=  theInitialStatePtr) { delete (*theTempResult)[0]; }
230            delete theTempResult;
231            theTempResult = theEvaporation->BreakItUp(*(*iList));
232          }
233        }
234      else // apply Evaporation in another case
235        {
236          theTempResult = theEvaporation->BreakItUp(*(*iList));
237        }
238     
239      G4bool deletePrimary = true;
240      G4int nsec = theTempResult->size();
241                 
242      // Sort out secondary fragments
243      if ( nsec > 0 )
244        {
245          G4FragmentVector::iterator j;
246          for (j = theTempResult->begin(); j != theTempResult->end(); ++j)
247            {
248              if((*j) == (*iList)) { deletePrimary = false; }
249              A = (*j)->GetA_asInt();
250
251              if(A <= 1) { theResults.push_back(*j); }                // gamma, p, n
252              else if(1 == nsec) { theEvapStableList.push_back(*j); } // no evaporation
253              else  // Analyse fragment
254                {
255                  G4double exEnergy = (*j)->GetExcitationEnergy();
256                  if(exEnergy < minExcitation) {
257                    Z = (*j)->GetZ_asInt();
258                    if(nist->GetIsotopeAbundance(Z, A) > 0.0) { 
259                      theResults.push_back(*j); // stable fragment
260                    } else {   
261                      theEvapList.push_back(*j); // unstable cold fragment
262                    }
263                  } else { theEvapList.push_back(*j); } // hot fragment
264                }
265            }
266        }
267      if( deletePrimary ) { delete (*iList); }
268      delete theTempResult;
269    } // end of the loop over theEvapList
270
271  //G4cout << "## After 2nd step " << theEvapList.size() << " was evap;  "
272  //     << theEvapStableList.size() << " for photo-evap; "
273  //     << theResults.size() << " results. " << G4endl;
274     
275  // -----------------------
276  // Photon-Evaporation loop
277  // -----------------------
278 
279  // normally should not reach this point - it is kind of work around         
280  for (iList = theEvapStableList.begin(); iList != theEvapStableList.end(); ++iList)
281    {
282      // photon-evaporation is applied
283      theTempResult = thePhotonEvaporation->BreakUpFragment(*iList);     
284      G4int nsec = theTempResult->size();
285         
286      // if there is a gamma emission then
287      if (nsec > 0)
288        {
289          G4FragmentVector::iterator j;
290          for (j = theTempResult->begin(); j != theTempResult->end(); ++j)
291            {
292              theResults.push_back(*j); 
293            }
294        }
295      delete theTempResult;
296
297      // priamry fragment is kept
298      theResults.push_back(*iList); 
299
300    } // end of photon-evaporation loop
301
302  //G4cout << "## After 3d step " << theEvapList.size() << " was evap;  "
303  //     << theEvapStableList.size() << " was photo-evap; "
304  //     << theResults.size() << " results. " << G4endl;
305   
306#ifdef debug
307  CheckConservation(theInitialState,*theResults);
308#endif
309
310  G4ReactionProductVector * theReactionProductVector = new G4ReactionProductVector;
311
312  // MAC (24/07/08)
313  // To optimise the storing speed, we reserve space in memory for the vector
314  theReactionProductVector->reserve( theResults.size() );
315
316  G4int theFragmentA, theFragmentZ;
317  G4LorentzVector theFragmentMomentum;
318
319  std::list<G4Fragment*>::iterator i;
320  for (i = theResults.begin(); i != theResults.end(); ++i) 
321    {
322      theFragmentA = static_cast<G4int>((*i)->GetA());
323      theFragmentZ = static_cast<G4int>((*i)->GetZ());
324      theFragmentMomentum = (*i)->GetMomentum();
325      G4ParticleDefinition* theKindOfFragment = 0;
326      if (theFragmentA == 0) {       // photon or e-
327        theKindOfFragment = (*i)->GetParticleDefinition();   
328      } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
329        theKindOfFragment = G4Neutron::NeutronDefinition();
330      } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
331        theKindOfFragment = G4Proton::ProtonDefinition();
332      } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
333        theKindOfFragment = G4Deuteron::DeuteronDefinition();
334      } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
335        theKindOfFragment = G4Triton::TritonDefinition();
336      } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
337        theKindOfFragment = G4He3::He3Definition();
338      } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha
339        theKindOfFragment = G4Alpha::AlphaDefinition();;
340      } else {
341        theKindOfFragment = 
342          theTableOfIons->GetIon(theFragmentZ,theFragmentA,0.0);
343      }
344      if (theKindOfFragment != 0) 
345        {
346          G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
347          theNew->SetMomentum(theFragmentMomentum.vect());
348          theNew->SetTotalEnergy(theFragmentMomentum.e());
349          theNew->SetFormationTime((*i)->GetCreationTime());
350          theReactionProductVector->push_back(theNew);
351        }
352      delete (*i);
353    }
354
355  return theReactionProductVector;
356}
357
358G4ReactionProductVector * 
359G4ExcitationHandler::Transform(G4FragmentVector * theFragmentVector) const
360{
361  if (theFragmentVector == 0) return 0;
362 
363  // Conversion from G4FragmentVector to G4ReactionProductVector
364  G4ParticleDefinition *theGamma = G4Gamma::GammaDefinition();
365  G4ParticleDefinition *theNeutron = G4Neutron::NeutronDefinition();
366  G4ParticleDefinition *theProton = G4Proton::ProtonDefinition();   
367  G4ParticleDefinition *theDeuteron = G4Deuteron::DeuteronDefinition();
368  G4ParticleDefinition *theTriton = G4Triton::TritonDefinition();
369  G4ParticleDefinition *theHelium3 = G4He3::He3Definition();
370  G4ParticleDefinition *theAlpha = G4Alpha::AlphaDefinition();
371  G4ParticleDefinition *theKindOfFragment = 0;
372  theNeutron->SetVerboseLevel(2);
373  G4ReactionProductVector * theReactionProductVector = new G4ReactionProductVector;
374
375  // MAC (24/07/08)
376  // To optimise the storing speed, we reserve space in memory for the vector
377  theReactionProductVector->reserve( theFragmentVector->size() * sizeof(G4ReactionProduct*) );
378
379  G4int theFragmentA, theFragmentZ;
380  G4LorentzVector theFragmentMomentum;
381
382  G4FragmentVector::iterator i;
383  for (i = theFragmentVector->begin(); i != theFragmentVector->end(); i++) {
384    //    std::cout << (*i) <<'\n';
385    theFragmentA = static_cast<G4int>((*i)->GetA());
386    theFragmentZ = static_cast<G4int>((*i)->GetZ());
387    theFragmentMomentum = (*i)->GetMomentum();
388    theKindOfFragment = 0;
389    if (theFragmentA == 0 && theFragmentZ == 0) {       // photon
390      theKindOfFragment = theGamma;     
391    } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
392      theKindOfFragment = theNeutron;
393    } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
394      theKindOfFragment = theProton;
395    } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
396      theKindOfFragment = theDeuteron;
397    } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
398      theKindOfFragment = theTriton;
399    } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
400      theKindOfFragment = theHelium3;
401    } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha
402      theKindOfFragment = theAlpha;
403    } else {
404      theKindOfFragment = 
405        theTableOfIons->GetIon(theFragmentZ,theFragmentA,0.0);
406    }
407    if (theKindOfFragment != 0) 
408      {
409        G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
410        theNew->SetMomentum(theFragmentMomentum.vect());
411        theNew->SetTotalEnergy(theFragmentMomentum.e());
412        theNew->SetFormationTime((*i)->GetCreationTime());
413        theReactionProductVector->push_back(theNew);
414      }
415  }
416  if (theFragmentVector != 0)
417    { 
418      std::for_each(theFragmentVector->begin(), theFragmentVector->end(), DeleteFragment());
419      delete theFragmentVector;
420    }
421  G4ReactionProductVector::iterator debugit;
422  for(debugit=theReactionProductVector->begin(); 
423      debugit!=theReactionProductVector->end(); debugit++)
424    {
425    if((*debugit)->GetTotalEnergy()<1.*eV)
426      {
427        if(getenv("G4DebugPhotonevaporationData"))
428          {
429            G4cerr << "G4ExcitationHandler: Warning: Photonevaporation data not exact."<<G4endl;
430            G4cerr << "G4ExcitationHandler: Warning: Found gamma with energy = "
431                   << (*debugit)->GetTotalEnergy()/MeV << "MeV"
432                   << G4endl;
433          }
434        delete (*debugit);
435        *debugit = 0;
436      }
437  }
438  G4ReactionProduct* tmpPtr=0;
439  theReactionProductVector->erase(std::remove_if(theReactionProductVector->begin(),
440                                                 theReactionProductVector->end(),
441                                                 std::bind2nd(std::equal_to<G4ReactionProduct*>(),
442                                                              tmpPtr)),
443                                  theReactionProductVector->end());
444  return theReactionProductVector;
445}
446
447
448#ifdef debug
449void G4ExcitationHandler::CheckConservation(const G4Fragment & theInitialState,
450                                            G4FragmentVector * Result) const
451{
452  G4double ProductsEnergy =0;
453  G4ThreeVector ProductsMomentum;
454  G4int ProductsA = 0;
455  G4int ProductsZ = 0;
456  G4FragmentVector::iterator h;
457  for (h = Result->begin(); h != Result->end(); h++) {
458    G4LorentzVector tmp = (*h)->GetMomentum();
459    ProductsEnergy += tmp.e();
460    ProductsMomentum += tmp.vect();
461    ProductsA += (*h)->GetA_asInt();
462    ProductsZ += (*h)->GetZ_asInt();
463  }
464 
465  if (ProductsA != theInitialState.GetA_asInt()) {
466    G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl;
467    G4cout << "G4ExcitationHandler.cc: Barionic Number Conservation test for deexcitation fragments" 
468           << G4endl; 
469    G4cout << "Initial A = " << theInitialState.GetA_asInt() 
470           << "   Fragments A = " << ProductsA << "   Diference --> " 
471           << theInitialState.GetA() - ProductsA << G4endl;
472  }
473  if (ProductsZ != theInitialState.GetZ_asInt()) {
474    G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
475    G4cout << "G4ExcitationHandler.cc: Charge Conservation test for deexcitation fragments" 
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()) > 1.0*keV) {
482    G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
483    G4cout << "G4ExcitationHandler.cc: Energy Conservation test for deexcitation fragments" 
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()) > 1.0*keV || 
490      std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV ||
491      std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) {
492    G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
493    G4cout << "G4ExcitationHandler.cc: Momentum Conservation test for deexcitation fragments" 
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#endif
502
503
504
505
Note: See TracBrowser for help on using the repository browser.