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

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

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