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

Last change on this file since 836 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 19.2 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// Hadronic Process: Nuclear De-excitations
28// by V. Lara (May 1998)
29// Modif (30 June 1998) by V. Lara:
30//      -Modified the Transform method for use G4ParticleTable and
31//       therefore G4IonTable. It makes possible to convert all kind
32//       of fragments (G4Fragment) produced in deexcitation to
33//       G4DynamicParticle
34//      -It uses default algorithms for:
35//              Evaporation: G4StatEvaporation
36//              MultiFragmentation: G4DummyMF (a dummy one)
37//              Fermi Breakup model: G4StatFermiBreakUp
38
39
40#include "G4ExcitationHandler.hh"
41#include <list>
42
43//#define debugphoton
44
45
46G4ExcitationHandler::G4ExcitationHandler():
47  maxZForFermiBreakUp(1),maxAForFermiBreakUp(1),minEForMultiFrag(4.0*GeV),
48  MyOwnEvaporationClass(true), MyOwnMultiFragmentationClass(true),MyOwnFermiBreakUpClass(true),
49  MyOwnPhotonEvaporationClass(true)
50{                                                                         
51  theTableOfParticles = G4ParticleTable::GetParticleTable();
52 
53  theEvaporation = new G4Evaporation;
54  theMultiFragmentation = new G4StatMF;
55  theFermiModel = new G4FermiBreakUp;
56  thePhotonEvaporation = new G4PhotonEvaporation;
57}
58
59G4ExcitationHandler::G4ExcitationHandler(const G4ExcitationHandler &)
60{
61  throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::copy_constructor: is meant to not be accessable! ");
62}
63
64
65G4ExcitationHandler::~G4ExcitationHandler()
66{
67  if (MyOwnEvaporationClass) delete theEvaporation;
68  if (MyOwnMultiFragmentationClass) delete theMultiFragmentation;
69  if (MyOwnFermiBreakUpClass) delete theFermiModel;
70  if (MyOwnPhotonEvaporationClass) delete thePhotonEvaporation;
71}
72
73
74const G4ExcitationHandler & G4ExcitationHandler::operator=(const G4ExcitationHandler &)
75{
76  throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator=: is meant to not be accessable! ");
77
78  return *this;
79}
80
81
82G4bool G4ExcitationHandler::operator==(const G4ExcitationHandler &) const
83{
84  throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator==: is meant to not be accessable! ");
85  return false;
86} 
87
88G4bool G4ExcitationHandler::operator!=(const G4ExcitationHandler &) const
89{
90  throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator!=: is meant to not be accessable! ");
91  return true;
92}
93
94
95G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment &theInitialState) const
96{
97
98  G4FragmentVector * theResult = 0; 
99  G4double exEnergy = theInitialState.GetExcitationEnergy();
100  //  G4cout << " first exEnergy in MeV: " << exEnergy/MeV << G4endl;
101  G4double A = theInitialState.GetA();
102  G4int Z = static_cast<G4int>(theInitialState.GetZ());
103  G4FragmentVector* theTempResult = 0; 
104  G4Fragment theExcitedNucleus;
105
106  // Test applicability
107  if (A > 4) 
108    {
109      // Initial State De-Excitation
110      if(A<GetMaxA()&&Z<GetMaxZ()) 
111        //     && exEnergy>G4NucleiPropertiesTable::GetBindingEnergy(Z,A)) {
112        {
113          theResult = theFermiModel->BreakItUp(theInitialState);
114        }
115      else   if (exEnergy>GetMinE()*A) 
116        {
117          theResult = theMultiFragmentation->BreakItUp(theInitialState);
118        }
119      else 
120        {
121          theResult = theEvaporation->BreakItUp(theInitialState);
122        }
123   
124
125
126
127      // De-Excitation loop
128      // ------------------
129      // Check if there are excited fragments
130      std::list<G4Fragment*> theResultList;
131      G4FragmentVector::iterator j;
132      std::list<G4Fragment*>::iterator i;
133      for (j = theResult->begin(); j != theResult->end();j++) 
134        {
135          theResultList.push_back(*j);
136        }
137      theResult->clear();
138      for (i = theResultList.begin(); i != theResultList.end(); i++) 
139        {
140          exEnergy = (*i)->GetExcitationEnergy();
141          //      G4cout << " exEnergy in MeV: " << exEnergy/MeV << G4endl;
142          if (exEnergy > 0.0) 
143            {
144              A = (*i)->GetA();
145              Z = static_cast<G4int>((*i)->GetZ());
146              theExcitedNucleus = *(*i);
147              // try to de-excite this fragment
148              if( A < GetMaxA() && Z < GetMaxZ() )
149                // && exEnergy>G4NucleiPropertiesTable::GetBindingEnergy(Z,A))
150                {
151                  // Fermi Breakup not now called for for exotic fragments for good reasons...
152                  // theTempResult = theFermiModel->BreakItUp(theExcitedNucleus);
153                  //if (theTempResult->size() == 1)
154                  //  {
155                  //    std::for_each(theTempResult->begin(),theTempResult->end(), G4Delete());
156                  //    delete theTempResult;
157                  //  }
158                  theTempResult = theEvaporation->BreakItUp(theExcitedNucleus);
159                } 
160              else 
161                {
162                  // Evaporation
163                  theTempResult = theEvaporation->BreakItUp(theExcitedNucleus);
164                }
165              // The Nucleus has been fragmented?
166              if (theTempResult->size() > 1) 
167                // If so :
168                {
169                  // Remove excited fragment from the result
170                  //    delete theResult->removeAt(i--);
171                  delete (*i);
172                  i = theResultList.erase(i);
173                  // and add theTempResult elements to theResult
174                  for (G4FragmentVector::reverse_iterator ri = theTempResult->rbegin();
175                       ri != theTempResult->rend(); ++ri)
176                    {
177                      theResultList.push_back(*ri);
178                    }
179                  delete theTempResult;
180                } 
181              else 
182                // If not :
183                {
184                  // it doesn't matter, we Follow with the next fragment but
185                  // I have to clean up
186                  std::for_each(theTempResult->begin(),theTempResult->end(), G4Delete());
187                  delete theTempResult;
188                }
189            }
190        }
191      for (i = theResultList.begin(); i != theResultList.end(); i++)
192        {
193          theResult->push_back(*i);
194        }
195      theResultList.clear();
196    }
197  else   // if A > 4
198    {
199      theResult = new G4FragmentVector();
200      theResult->push_back(new G4Fragment(theInitialState));
201    }
202
203  // Now we try to deexcite by means of PhotonEvaporation those fragments
204  // which are excited.
205 
206  theTempResult = 0;
207  std::list<G4Fragment*> theFinalResultList;
208  //AHtest  std::list<G4Fragment*> theFinalPhotonResultList;
209  std::list<G4Fragment*> theResultList;
210  std::list<G4Fragment*>::iterator j;
211  G4FragmentVector::iterator i;
212  for (i = theResult->begin(); i != theResult->end();i++) 
213    {
214      theResultList.push_back(*i);
215      //      G4cout << " Before loop list energy in MeV: " << ((*i)->GetExcitationEnergy())/MeV << G4endl;
216    }
217  theResult->clear();
218
219  for (j = theResultList.begin(); j != theResultList.end(); j++) {
220    //    G4cout << " Test loop list: " << (*j)->GetExcitationEnergy() << " size: " << theResultList.size() << G4endl;
221  }
222
223  //  for (j = theResultList.begin(); j != theResultList.end(); j++)
224  j = theResultList.begin();  //AH
225  while (j != theResultList.end()) //AH
226    {
227      if ((*j)->GetA() > 1 && (*j)->GetExcitationEnergy() > 0.1*eV) 
228        {
229          theExcitedNucleus = *(*j);
230          theTempResult = thePhotonEvaporation->BreakItUp(theExcitedNucleus);
231          // If Gamma Evaporation has succeed then
232          if (theTempResult->size() > 1) 
233            {
234              // Remove excited fragment from the result
235              //              delete (*j);
236              //              theResultList.erase(j--);
237              //              theResultList.erase(j); don't delete as there's no push back...
238              // and add theTempResult elements to theResult
239              for (G4FragmentVector::reverse_iterator ri = theTempResult->rbegin();
240                   ri != theTempResult->rend(); ++ri)
241                {
242#ifdef PRECOMPOUND_TEST
243                  if ((*ri)->GetA() == 0)
244                    (*ri)->SetCreatorModel(G4String("G4PhotonEvaporation"));
245                  else
246                    (*ri)->SetCreatorModel(G4String("ResidualNucleus"));
247#endif
248                  theResultList.push_back(*ri);
249                  //AHtest                theFinalPhotonResultList.push_back(*ri);
250                  //              theFinalResultList.push_back(*ri); don't add to final result as they'll go through the loop
251                }
252              delete theTempResult;
253            }
254          // In other case, just clean theTempResult and continue
255          else 
256            {
257              std::for_each(theTempResult->begin(), theTempResult->end(), DeleteFragment());
258              delete theTempResult;
259#ifdef debugphoton
260              G4cout << "G4ExcitationHandler: Gamma Evaporation could not deexcite the nucleus: \n"
261                     << "-----------------------------------------------------------------------\n"
262                     << theExcitedNucleus << '\n'
263                     << "-----------------------------------------------------------------------\n";
264#endif
265              G4double GammaEnergy = (*j)->GetExcitationEnergy();
266              G4double cosTheta = 1. - 2. * G4UniformRand();
267              G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
268              G4double phi = twopi * G4UniformRand();
269              G4ThreeVector GammaP(GammaEnergy * sinTheta * std::cos(phi),
270                                   GammaEnergy * sinTheta * std::sin(phi),
271                                   GammaEnergy * cosTheta );
272              G4LorentzVector Gamma4P(GammaP,GammaEnergy);
273              G4Fragment * theHandlerPhoton = new G4Fragment(Gamma4P,G4Gamma::GammaDefinition());
274             
275             
276             
277              G4double Mass = (*j)->GetGroundStateMass();
278              G4ThreeVector ResidualP((*j)->GetMomentum().vect() - GammaP);
279              G4double ResidualE = std::sqrt(ResidualP*ResidualP + Mass*Mass);
280              G4LorentzVector Residual4P(ResidualP,ResidualE);
281              (*j)->SetMomentum(Residual4P);
282             
283             
284       
285#ifdef PRECOMPOUND_TEST
286              theHandlerPhoton->SetCreatorModel("G4ExcitationHandler");
287#endif
288//            theFinalPhotonResultList.push_back( theHandlerPhoton );
289//            G4cout << " adding photon fragment " << G4endl;
290              theResultList.push_back( theHandlerPhoton );
291              //              theFinalResultList.push_back( theHandlerPhoton );
292              theFinalResultList.push_back(*j);
293#ifdef debugphoton
294              G4cout << "Emmited photon:\n"
295                     << theResultList.back() << '\n'
296                     << "Residual nucleus after photon emission:\n"
297                     << *(*j) << '\n'
298                     << "-----------------------------------------------------------------------\n";
299#endif
300              //test          j++; // AH only increment if not erased:
301            }   
302        } else {
303          //test          j++; // AH increment iterator if a proton or excitation energy small
304          theFinalResultList.push_back(*j);
305        }
306//      G4cout << " Inside loop list: " << (*j)->GetExcitationEnergy() << " size: " << theFinalResultList.size() << G4endl;
307      j++;
308    }
309  //  for (j = theResultList.begin(); j != theResultList.end(); j++)
310  for (j = theFinalResultList.begin(); j != theFinalResultList.end(); j++)
311    {
312      theResult->push_back(*j);
313    }
314
315//AHtest   for (j = theFinalPhotonResultList.begin(); j != theFinalPhotonResultList.end(); j++)
316//AHtest     {
317//AHtest       theResult->push_back(*j);
318//AHtest       number_results++;
319//AHtest     }
320
321
322  theResultList.clear();
323  theFinalResultList.clear();
324  //AHtest  theFinalPhotonResultList.clear();
325 
326 
327#ifdef debug
328  CheckConservation(theInitialState,theResult);
329#endif
330  // Change G4FragmentVector by G4DynamicParticle
331  return Transform(theResult);
332}
333
334G4ReactionProductVector * 
335G4ExcitationHandler::Transform(G4FragmentVector * theFragmentVector) const
336{
337  if (theFragmentVector == 0) return 0;
338 
339  // Conversion from G4FragmentVector to G4ReactionProductVector
340  G4ParticleDefinition *theGamma = G4Gamma::GammaDefinition();
341  G4ParticleDefinition *theNeutron = G4Neutron::NeutronDefinition();
342  G4ParticleDefinition *theProton = G4Proton::ProtonDefinition();   
343  G4ParticleDefinition *theDeuteron = G4Deuteron::DeuteronDefinition();
344  G4ParticleDefinition *theTriton = G4Triton::TritonDefinition();
345  G4ParticleDefinition *theHelium3 = G4He3::He3Definition();
346  G4ParticleDefinition *theAlpha = G4Alpha::AlphaDefinition();
347  G4ParticleDefinition *theKindOfFragment = 0;
348  theNeutron->SetVerboseLevel(2);
349  G4ReactionProductVector * theReactionProductVector = new G4ReactionProductVector;
350  G4int theFragmentA, theFragmentZ;
351  G4LorentzVector theFragmentMomentum;
352
353  G4FragmentVector::iterator i;
354  for (i = theFragmentVector->begin(); i != theFragmentVector->end(); i++) {
355    //    std::cout << (*i) <<'\n';
356    theFragmentA = static_cast<G4int>((*i)->GetA());
357    theFragmentZ = static_cast<G4int>((*i)->GetZ());
358    theFragmentMomentum = (*i)->GetMomentum();
359    theKindOfFragment = 0;
360    if (theFragmentA == 0 && theFragmentZ == 0) {       // photon
361      theKindOfFragment = theGamma;     
362    } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
363      theKindOfFragment = theNeutron;
364    } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
365      theKindOfFragment = theProton;
366    } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
367      theKindOfFragment = theDeuteron;
368    } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
369      theKindOfFragment = theTriton;
370    } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
371      theKindOfFragment = theHelium3;
372    } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha
373      theKindOfFragment = theAlpha;
374    } else {
375      theKindOfFragment = theTableOfParticles->FindIon(theFragmentZ,theFragmentA,0,theFragmentZ);
376    }
377    if (theKindOfFragment != 0) 
378      {
379        G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
380        theNew->SetMomentum(theFragmentMomentum.vect());
381        theNew->SetTotalEnergy(theFragmentMomentum.e());
382        theNew->SetFormationTime((*i)->GetCreationTime());
383#ifdef PRECOMPOUND_TEST
384        theNew->SetCreatorModel((*i)->GetCreatorModel());
385#endif
386        theReactionProductVector->push_back(theNew);
387      }
388  }
389  if (theFragmentVector != 0)
390    { 
391      std::for_each(theFragmentVector->begin(), theFragmentVector->end(), DeleteFragment());
392      delete theFragmentVector;
393    }
394  G4ReactionProductVector::iterator debugit;
395  for(debugit=theReactionProductVector->begin(); 
396      debugit!=theReactionProductVector->end(); debugit++)
397    {
398    if((*debugit)->GetTotalEnergy()<1.*eV)
399      {
400        if(getenv("G4DebugPhotonevaporationData"))
401          {
402            G4cerr << "G4ExcitationHandler: Warning: Photonevaporation data not exact."<<G4endl;
403            G4cerr << "G4ExcitationHandler: Warning: Found gamma with energy = "
404                   << (*debugit)->GetTotalEnergy()/MeV << "MeV"
405                   << G4endl;
406          }
407        delete (*debugit);
408        *debugit = 0;
409      }
410  }
411  G4ReactionProduct* tmpPtr=0;
412  theReactionProductVector->erase(std::remove_if(theReactionProductVector->begin(),
413                                                 theReactionProductVector->end(),
414                                                 std::bind2nd(std::equal_to<G4ReactionProduct*>(),
415                                                              tmpPtr)),
416                                  theReactionProductVector->end());
417  return theReactionProductVector;
418}
419
420
421#ifdef debug
422void G4ExcitationHandler::CheckConservation(const G4Fragment & theInitialState,
423                                            G4FragmentVector * Result) const
424{
425  G4double ProductsEnergy =0;
426  G4ThreeVector ProductsMomentum;
427  G4int ProductsA = 0;
428  G4int ProductsZ = 0;
429  G4FragmentVector::iterator h;
430  for (h = Result->begin(); h != Result->end(); h++) {
431    G4LorentzVector tmp = (*h)->GetMomentum();
432    ProductsEnergy += tmp.e();
433    ProductsMomentum += tmp.vect();
434    ProductsA += static_cast<G4int>((*h)->GetA());
435    ProductsZ += static_cast<G4int>((*h)->GetZ());
436  }
437 
438  if (ProductsA != theInitialState.GetA()) {
439    G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl;
440    G4cout << "G4ExcitationHandler.cc: Barionic Number Conservation test for deexcitation fragments" 
441           << G4endl; 
442    G4cout << "Initial A = " << theInitialState.GetA() 
443           << "   Fragments A = " << ProductsA << "   Diference --> " 
444           << theInitialState.GetA() - ProductsA << G4endl;
445  }
446  if (ProductsZ != theInitialState.GetZ()) {
447    G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
448    G4cout << "G4ExcitationHandler.cc: Charge Conservation test for deexcitation fragments" 
449           << G4endl; 
450    G4cout << "Initial Z = " << theInitialState.GetZ() 
451           << "   Fragments Z = " << ProductsZ << "   Diference --> " 
452           << theInitialState.GetZ() - ProductsZ << G4endl;
453  }
454  if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) {
455    G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
456    G4cout << "G4ExcitationHandler.cc: Energy Conservation test for deexcitation fragments" 
457           << G4endl; 
458    G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
459           << "   Fragments E = " << ProductsEnergy/MeV  << " MeV   Diference --> " 
460           << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl;
461  } 
462  if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV || 
463      std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV ||
464      std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) {
465    G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
466    G4cout << "G4ExcitationHandler.cc: Momentum Conservation test for deexcitation fragments" 
467           << G4endl; 
468    G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
469           << "   Fragments P = " << ProductsMomentum  << " MeV   Diference --> " 
470           << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl;
471  }
472  return;
473}
474#endif
475
476
477
478
Note: See TracBrowser for help on using the repository browser.