Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/de_excitation/handler/src/G4ExcitationHandler.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/de_excitation/handler/src/G4ExcitationHandler.cc (Version 11.3.0) and /processes/hadronic/models/de_excitation/handler/src/G4ExcitationHandler.cc (Version 9.0.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 //
 26 // Hadronic Process: Nuclear De-excitations        27 // Hadronic Process: Nuclear De-excitations
 27 // by V. Lara (May 1998)                           28 // by V. Lara (May 1998)
 28 //                                             <<  29 // Modif (30 June 1998) by V. Lara:
 29 //                                             << 
 30 // Modified:                                   << 
 31 // 30 June 1998 by V. Lara:                    << 
 32 //      -Modified the Transform method for use     30 //      -Modified the Transform method for use G4ParticleTable and 
 33 //       therefore G4IonTable. It makes possib     31 //       therefore G4IonTable. It makes possible to convert all kind 
 34 //       of fragments (G4Fragment) produced in     32 //       of fragments (G4Fragment) produced in deexcitation to 
 35 //       G4DynamicParticle                         33 //       G4DynamicParticle
 36 //      -It uses default algorithms for:           34 //      -It uses default algorithms for:
 37 //              Evaporation: G4Evaporation     <<  35 //              Evaporation: G4StatEvaporation
 38 //              MultiFragmentation: G4StatMF   <<  36 //              MultiFragmentation: G4DummyMF (a dummy one)
 39 //              Fermi Breakup model: G4FermiBr <<  37 //              Fermi Breakup model: G4StatFermiBreakUp
 40 // 24 Jul 2008 by M. A. Cortes Giraldo:        <<  38 
 41 //      -Max Z,A for Fermi Break-Up turns to 9 << 
 42 //      -BreakItUp() reorganised and bug in Ev << 
 43 //      -Transform() optimised                 << 
 44 // (September 2008) by J. M. Quesada. External << 
 45 //      -inverse cross section option (default << 
 46 //      -superimposed Coulomb barrier (if useS << 
 47 // September 2009 by J. M. Quesada:            << 
 48 //      -according to Igor Pshenichnov, SMM wi << 
 49 // 27 Nov 2009 by V.Ivanchenko:                << 
 50 //      -cleanup the logic, reduce number inte << 
 51 // 11 May 2010 by V.Ivanchenko:                << 
 52 //      -FermiBreakUp activated, used integer  << 
 53 //       final photon deexcitation; used check << 
 54 //       unstable fragments with A <5          << 
 55 // 22 March 2011 by V.Ivanchenko: general clea << 
 56 //       products of Fermi Break Up cannot be  << 
 57 // 30 March 2011 by V.Ivanchenko removed priva << 
 58 //       to the source                         << 
 59 // 23 January 2012 by V.Ivanchenko general cle << 
 60 //    objects, propagate G4PhotonEvaporation p << 
 61 //    not delete it here                       << 
 62                                                    39 
 63 #include "G4ExcitationHandler.hh"                  40 #include "G4ExcitationHandler.hh"
 64 #include "G4SystemOfUnits.hh"                  <<  41 #include <list>
 65 #include "G4LorentzVector.hh"                  << 
 66 #include "G4ThreeVector.hh"                    << 
 67 #include "G4ParticleTable.hh"                  << 
 68 #include "G4ParticleTypes.hh"                  << 
 69 #include "G4Ions.hh"                           << 
 70 #include "G4Electron.hh"                       << 
 71 #include "G4Lambda.hh"                         << 
 72                                                << 
 73 #include "G4VMultiFragmentation.hh"            << 
 74 #include "G4VFermiBreakUp.hh"                  << 
 75 #include "G4Element.hh"                        << 
 76 #include "G4ElementTable.hh"                   << 
 77                                                << 
 78 #include "G4VEvaporation.hh"                   << 
 79 #include "G4VEvaporationChannel.hh"            << 
 80 #include "G4Evaporation.hh"                    << 
 81 #include "G4PhotonEvaporation.hh"              << 
 82 #include "G4StatMF.hh"                         << 
 83 #include "G4FermiBreakUpVI.hh"                 << 
 84 #include "G4NuclearLevelData.hh"               << 
 85 #include "G4PhysicsModelCatalog.hh"            << 
 86                                                << 
 87 G4ExcitationHandler::G4ExcitationHandler()     << 
 88   : minEForMultiFrag(1.*CLHEP::TeV), minExcita << 
 89     maxExcitation(100.*CLHEP::MeV)             << 
 90 {                                              << 
 91   thePartTable = G4ParticleTable::GetParticleT << 
 92   theTableOfIons = thePartTable->GetIonTable() << 
 93   nist = G4NistManager::Instance();            << 
 94                                                << 
 95   theMultiFragmentation = new G4StatMF();      << 
 96   theFermiModel = new G4FermiBreakUpVI();      << 
 97   thePhotonEvaporation = new G4PhotonEvaporati << 
 98   SetEvaporation(new G4Evaporation(thePhotonEv << 
 99   theResults.reserve(60);                      << 
100   results.reserve(30);                         << 
101   theEvapList.reserve(30);                     << 
102                                                << 
103   theElectron = G4Electron::Electron();        << 
104   theNeutron = G4Neutron::NeutronDefinition(); << 
105   theProton = G4Proton::ProtonDefinition();    << 
106   theDeuteron = G4Deuteron::DeuteronDefinition << 
107   theTriton = G4Triton::TritonDefinition();    << 
108   theHe3 = G4He3::He3Definition();             << 
109   theAlpha = G4Alpha::AlphaDefinition();       << 
110   theLambda = G4Lambda::Lambda();              << 
111                                                    42 
112   fLambdaMass = theLambda->GetPDGMass();       <<  43 //#define debugphoton
113                                                    44 
114   if(fVerbose > 1) { G4cout << "### New handle << 
115 }                                              << 
116                                                    45 
117 G4ExcitationHandler::~G4ExcitationHandler()    <<  46 G4ExcitationHandler::G4ExcitationHandler():
118 {                                              <<  47   maxZForFermiBreakUp(1),maxAForFermiBreakUp(1),minEForMultiFrag(4.0*GeV),
119   delete theMultiFragmentation;                <<  48   MyOwnEvaporationClass(true), MyOwnMultiFragmentationClass(true),MyOwnFermiBreakUpClass(true),
120   delete theFermiModel;                        <<  49   MyOwnPhotonEvaporationClass(true)
121   if(isEvapLocal) { delete theEvaporation; }   <<  50 {                                                                          
                                                   >>  51   theTableOfParticles = G4ParticleTable::GetParticleTable();
                                                   >>  52   
                                                   >>  53   theEvaporation = new G4Evaporation;
                                                   >>  54   theMultiFragmentation = new G4StatMF;
                                                   >>  55   theFermiModel = new G4FermiBreakUp;
                                                   >>  56   thePhotonEvaporation = new G4PhotonEvaporation;
122 }                                                  57 }
123                                                    58 
124 void G4ExcitationHandler::SetParameters()      <<  59 G4ExcitationHandler::G4ExcitationHandler(const G4ExcitationHandler &)
125 {                                                  60 {
126   G4NuclearLevelData* ndata = G4NuclearLevelDa <<  61   throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::copy_constructor: is meant to not be accessable! ");
127   auto param = ndata->GetParameters();         << 
128   isActive = true;                             << 
129   // check if de-excitation is needed          << 
130   if (fDummy == param->GetDeexChannelsType())  << 
131     isActive = false;                          << 
132   } else {                                     << 
133     // upload data for elements used in geomet << 
134     G4int Zmax = 20;                           << 
135     const G4ElementTable* table = G4Element::G << 
136     for (auto const & elm : *table) { Zmax = s << 
137     ndata->UploadNuclearLevelData(Zmax+1);     << 
138   }                                            << 
139   minEForMultiFrag = param->GetMinExPerNucleou << 
140   minExcitation = param->GetMinExcitation();   << 
141   maxExcitation = param->GetPrecoHighEnergy(); << 
142                                                << 
143   // allowing local debug printout             << 
144   fVerbose = std::max(fVerbose, param->GetVerb << 
145   if (isActive) {                              << 
146     if (nullptr == thePhotonEvaporation) {     << 
147       SetPhotonEvaporation(new G4PhotonEvapora << 
148     }                                          << 
149     if (nullptr == theFermiModel) {            << 
150       SetFermiModel(new G4FermiBreakUpVI());   << 
151     }                                          << 
152     if (nullptr == theMultiFragmentation) {    << 
153       SetMultiFragmentation(new G4StatMF());   << 
154     }                                          << 
155     if (nullptr == theEvaporation) {           << 
156       SetEvaporation(new G4Evaporation(thePhot << 
157     }                                          << 
158   }                                            << 
159   theFermiModel->SetVerbose(fVerbose);         << 
160   if(fVerbose > 1) {                           << 
161     G4cout << "G4ExcitationHandler::SetParamet << 
162   }                                            << 
163 }                                                  62 }
164                                                    63 
165 void G4ExcitationHandler::Initialise()         << 
166 {                                              << 
167   if(isInitialised) { return; }                << 
168   if(fVerbose > 1) {                           << 
169     G4cout << "G4ExcitationHandler::Initialise << 
170   }                                            << 
171   G4DeexPrecoParameters* param =               << 
172     G4NuclearLevelData::GetInstance()->GetPara << 
173   isInitialised = true;                        << 
174   SetParameters();                             << 
175   if(isActive) {                               << 
176     theFermiModel->Initialise();               << 
177     theEvaporation->InitialiseChannels();      << 
178   }                                            << 
179   // dump level is controlled by parameter cla << 
180   param->Dump();                               << 
181 }                                              << 
182                                                    64 
183 void G4ExcitationHandler::SetEvaporation(G4VEv <<  65 G4ExcitationHandler::~G4ExcitationHandler()
184 {                                                  66 {
185   if(nullptr != ptr && ptr != theEvaporation)  <<  67   if (MyOwnEvaporationClass) delete theEvaporation;
186     theEvaporation = ptr;                      <<  68   if (MyOwnMultiFragmentationClass) delete theMultiFragmentation;
187     theEvaporation->SetPhotonEvaporation(thePh <<  69   if (MyOwnFermiBreakUpClass) delete theFermiModel;
188     theEvaporation->SetFermiBreakUp(theFermiMo <<  70   if (MyOwnPhotonEvaporationClass) delete thePhotonEvaporation;
189     isEvapLocal = flag;                        << 
190     if(fVerbose > 1) {                         << 
191       G4cout << "G4ExcitationHandler::SetEvapo << 
192     }                                          << 
193   }                                            << 
194 }                                                  71 }
195                                                    72 
196 void                                           << 
197 G4ExcitationHandler::SetMultiFragmentation(G4V << 
198 {                                              << 
199   if(nullptr != ptr && ptr != theMultiFragment << 
200     delete theMultiFragmentation;              << 
201     theMultiFragmentation = ptr;               << 
202   }                                            << 
203 }                                              << 
204                                                    73 
205 void G4ExcitationHandler::SetFermiModel(G4VFer <<  74 const G4ExcitationHandler & G4ExcitationHandler::operator=(const G4ExcitationHandler &)
206 {                                                  75 {
207   if(nullptr != ptr && ptr != theFermiModel) { <<  76   throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator=: is meant to not be accessable! ");
208     delete theFermiModel;                      << 
209     theFermiModel = ptr;                       << 
210     if(nullptr != theEvaporation) {            << 
211       theEvaporation->SetFermiBreakUp(theFermi << 
212     }                                          << 
213   }                                            << 
214 }                                              << 
215                                                    77 
216 void                                           <<  78   return *this;
217 G4ExcitationHandler::SetPhotonEvaporation(G4VE << 
218 {                                              << 
219   if(nullptr != ptr && ptr != thePhotonEvapora << 
220     delete thePhotonEvaporation;               << 
221     thePhotonEvaporation = ptr;                << 
222     if(nullptr != theEvaporation) {            << 
223       theEvaporation->SetPhotonEvaporation(ptr << 
224     }                                          << 
225     if(fVerbose > 1) {                         << 
226       G4cout << "G4ExcitationHandler::SetPhoto << 
227              << " for handler " << this << G4e << 
228     }                                          << 
229   }                                            << 
230 }                                                  79 }
231                                                    80 
232 void G4ExcitationHandler::SetDeexChannelsType( << 
233 {                                              << 
234   G4Evaporation* evap = static_cast<G4Evaporat << 
235   if(fVerbose > 1) {                           << 
236     G4cout << "G4ExcitationHandler::SetDeexCha << 
237      << " for " << this << G4endl;             << 
238   }                                            << 
239   if(val == fDummy) {                          << 
240     isActive = false;                          << 
241     return;                                    << 
242   }                                            << 
243   if (nullptr == evap) { return; }             << 
244   if (val == fEvaporation) {                   << 
245     evap->SetDefaultChannel();                 << 
246   } else if (val == fCombined) {               << 
247     evap->SetCombinedChannel();                << 
248   } else if (val == fGEM) {                    << 
249     evap->SetGEMChannel();                     << 
250   } else if (val == fGEMVI) {                  << 
251     evap->SetGEMVIChannel();                   << 
252   }                                            << 
253   evap->InitialiseChannels();                  << 
254   if (fVerbose > 1) {                          << 
255     if (G4Threading::IsMasterThread()) {       << 
256       G4cout << "Number of de-excitation chann << 
257        << theEvaporation->GetNumberOfChannels( << 
258       G4cout << " " << this;                   << 
259     }                                          << 
260     G4cout << G4endl;                          << 
261   }                                            << 
262 }                                              << 
263                                                    81 
264 G4VEvaporation* G4ExcitationHandler::GetEvapor <<  82 G4bool G4ExcitationHandler::operator==(const G4ExcitationHandler &) const
265 {                                                  83 {
266   if (nullptr != theEvaporation) { SetParamete <<  84   throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator==: is meant to not be accessable! ");
267   return theEvaporation;                       <<  85   return false;
268 }                                              <<  86 } 
269                                                    87 
270 G4VMultiFragmentation* G4ExcitationHandler::Ge <<  88 G4bool G4ExcitationHandler::operator!=(const G4ExcitationHandler &) const
271 {                                                  89 {
272   if (nullptr != theMultiFragmentation) { SetP <<  90   throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator!=: is meant to not be accessable! ");
273   return theMultiFragmentation;                <<  91   return true;
274 }                                                  92 }
275                                                    93 
276 G4VFermiBreakUp* G4ExcitationHandler::GetFermi << 
277 {                                              << 
278   if (nullptr != theFermiModel) { SetParameter << 
279   return theFermiModel;                        << 
280 }                                              << 
281                                                    94 
282 G4VEvaporationChannel* G4ExcitationHandler::Ge <<  95 G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment &theInitialState) const
283 {                                                  96 {
284   if(nullptr != thePhotonEvaporation) { SetPar << 
285   return thePhotonEvaporation;                 << 
286 }                                              << 
287                                                << 
288 G4ReactionProductVector *                      << 
289 G4ExcitationHandler::BreakItUp(const G4Fragmen << 
290 {                                              << 
291   // Variables existing until end of method    << 
292   G4Fragment * theInitialStatePtr = new G4Frag << 
293   if (fVerbose > 1) {                          << 
294     G4cout << "@@@@@@@@@@ Start G4Excitation H << 
295     G4cout << theInitialState << G4endl;       << 
296   }                                            << 
297   if (!isInitialised) { Initialise(); }        << 
298                                                    97 
299   // pointer to fragment vector which receives <<  98   G4FragmentVector * theResult = 0; 
300   G4FragmentVector * theTempResult = nullptr;  << 
301                                                << 
302   theResults.clear();                          << 
303   theEvapList.clear();                         << 
304                                                << 
305   // Variables to describe the excited configu << 
306   G4double exEnergy = theInitialState.GetExcit     99   G4double exEnergy = theInitialState.GetExcitationEnergy();
307   G4int A = theInitialState.GetA_asInt();      << 100   //  G4cout << " first exEnergy in MeV: " << exEnergy/MeV << G4endl;
308   G4int Z = theInitialState.GetZ_asInt();      << 101   G4double A = theInitialState.GetA();
309   G4int nL = theInitialState.GetNumberOfLambda << 102   G4int Z = static_cast<G4int>(theInitialState.GetZ());
310                                                << 103   G4FragmentVector* theTempResult = 0; 
311   // too much excitation                       << 104   G4Fragment theExcitedNucleus;
312   if (exEnergy > A*maxExcitation && A > 0) {   << 105 
313     ++fWarnings;                               << 106   // Test applicability
314     if(fWarnings < 0) {                        << 107   if (A > 4) 
315       G4ExceptionDescription ed;               << 108     {
316       ed << "High excitation Fragment Z= " <<  << 109       // Initial State De-Excitation 
317    << " Eex/A(MeV)= " << exEnergy/A;           << 110       if(A<GetMaxA()&&Z<GetMaxZ()) 
318       G4Exception("G4ExcitationHandler::BreakI << 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));
319     }                                             201     }
320   }                                            << 
321                                                   202 
322   // for hyper-nuclei subtract lambdas from th << 203   // Now we try to deexcite by means of PhotonEvaporation those fragments
323   G4double lambdaF = 0.0;                      << 204   // which are excited.
324   G4LorentzVector lambdaLV = theInitialStatePt << 205   
325   if (0 < nL) {                                << 206   theTempResult = 0;
326                                                << 207   std::list<G4Fragment*> theFinalResultList;
327     // is it a stable hyper-nuclei?            << 208   //AHtest  std::list<G4Fragment*> theFinalPhotonResultList;
328     if(A >= 3 && A <= 5 && nL <= 2) {          << 209   std::list<G4Fragment*> theResultList;
329       G4int pdg = 0;                           << 210   std::list<G4Fragment*>::iterator j;
330       if(3 == A && 1 == nL) {                  << 211   G4FragmentVector::iterator i;
331         pdg = 1010010030;                      << 212   for (i = theResult->begin(); i != theResult->end();i++) 
332       } else if(5 == A && 2 == Z && 1 == nL) { << 213     {
333         pdg = 1010020050;                      << 214       theResultList.push_back(*i);
334       } else if(4 == A) {                      << 215       //      G4cout << " Before loop list energy in MeV: " << ((*i)->GetExcitationEnergy())/MeV << G4endl;
335   if(1 == Z && 1 == nL) {                      << 216     }
336     pdg = 1010010040;                          << 217   theResult->clear();
337   } else if(2 == Z && 1 == nL) {               << 218 
338     pdg = 1010020040;                          << 219   for (j = theResultList.begin(); j != theResultList.end(); j++) {
339   } else if(0 == Z && 2 == nL) {               << 220     //    G4cout << " Test loop list: " << (*j)->GetExcitationEnergy() << " size: " << theResultList.size() << G4endl;
340     pdg = 1020000040;                          << 221   }
341   } else if(1 == Z && 2 == nL) {               << 222 
342     pdg = 1020010040;                          << 223   //  for (j = theResultList.begin(); j != theResultList.end(); j++) 
343   }                                            << 224   j = theResultList.begin();  //AH
344       }                                        << 225   while (j != theResultList.end()) //AH
345       // initial state is one of hyper-nuclei  << 226     {
346       if (0 < pdg) {                           << 227       if ((*j)->GetA() > 1 && (*j)->GetExcitationEnergy() > 0.1*eV) 
347   const G4ParticleDefinition* part = thePartTa << 228         {
348         if(nullptr != part) {                  << 229           theExcitedNucleus = *(*j);
349     G4ReactionProduct* theNew = new G4Reaction << 230           theTempResult = thePhotonEvaporation->BreakItUp(theExcitedNucleus);
350     G4ThreeVector dir = G4ThreeVector( 0.0, 0. << 231           // If Gamma Evaporation has succeed then
351           if ( lambdaLV.vect().mag() > CLHEP:: << 232           if (theTempResult->size() > 1) 
352       dir = lambdaLV.vect().unit();            << 233             {
353           }                                    << 234               // Remove excited fragment from the result 
354     G4double mass = part->GetPDGMass();        << 235         //        delete (*j);
355     G4double etot = std::max(lambdaLV.e(), mas << 236         //              theResultList.erase(j--);
356           dir *= std::sqrt((etot - mass)*(etot << 237         //        theResultList.erase(j); don't delete as there's no push back...
357     theNew->SetMomentum(dir);                  << 238               // and add theTempResult elements to theResult
358     theNew->SetTotalEnergy(etot);              << 239         for (G4FragmentVector::reverse_iterator ri = theTempResult->rbegin();
359     theNew->SetFormationTime(theInitialState.G << 240        ri != theTempResult->rend(); ++ri)
360     theNew->SetCreatorModelID(theInitialState. << 241                 {
361     G4ReactionProductVector* v = new G4Reactio << 242 #ifdef PRECOMPOUND_TEST
362           v->push_back(theNew);                << 243                   if ((*ri)->GetA() == 0)
363     return v;                                  << 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);
364   }                                               305   }
365       }                                        << 306 //      G4cout << " Inside loop list: " << (*j)->GetExcitationEnergy() << " size: " << theFinalResultList.size() << G4endl;
                                                   >> 307       j++;
366     }                                             308     }
367     G4double mass = theInitialStatePtr->GetGro << 309   //  for (j = theResultList.begin(); j != theResultList.end(); j++)
368     lambdaF = nL*(fLambdaMass - CLHEP::neutron << 310   for (j = theFinalResultList.begin(); j != theFinalResultList.end(); j++)
369                                                << 311     {
370     // de-excitation with neutrons instead of  << 312       theResult->push_back(*j);
371     theInitialStatePtr->SetZAandMomentum(lambd << 
372                                                << 
373     // 4-momentum not used in de-excitation    << 
374     lambdaLV *= lambdaF;                       << 
375   } else if (0 > nL) {                         << 
376     ++fWarnings;                               << 
377     if(fWarnings < 0) {                        << 
378       G4ExceptionDescription ed;               << 
379       ed << "Fragment with negative L: Z=" <<  << 
380    << " Eex/A(MeV)= " << exEnergy/A;           << 
381       G4Exception("G4ExcitationHandler::BreakI << 
382     }                                             313     }
383   }                                            << 
384                                                   314 
385   // In case A <= 1 the fragment will not perf << 315 //AHtest   for (j = theFinalPhotonResultList.begin(); j != theFinalPhotonResultList.end(); j++)
386   if (A <= 1 || !isActive) {                   << 316 //AHtest     {
387     theResults.push_back( theInitialStatePtr ) << 317 //AHtest       theResult->push_back(*j);
388                                                << 318 //AHtest       number_results++;
389     // check if a fragment is stable           << 319 //AHtest     }
390   } else if (exEnergy < minExcitation && nist- << 
391     theResults.push_back( theInitialStatePtr ) << 
392                                                << 
393     // JMQ 150909: first step in de-excitation << 
394     // Fragments after the first step are stor << 
395   } else {                                     << 
396     if ((A<maxAForFermiBreakUp && Z<maxZForFer << 
397   || exEnergy <= minEForMultiFrag*A) {         << 
398       theEvapList.push_back(theInitialStatePtr << 
399                                                   320 
400     // Statistical Multifragmentation will tak << 
401     } else {                                   << 
402       theTempResult = theMultiFragmentation->B << 
403       if (nullptr == theTempResult) {          << 
404   theEvapList.push_back(theInitialStatePtr);   << 
405       } else {                                 << 
406   std::size_t nsec = theTempResult->size();    << 
407                                                << 
408   // no fragmentation                          << 
409   if (0 == nsec) {                             << 
410     theEvapList.push_back(theInitialStatePtr); << 
411                                                << 
412     // secondary are produced - sort out secon << 
413   } else {                                     << 
414     G4bool deletePrimary = true;               << 
415     for (auto const & ptr : *theTempResult) {  << 
416       if (ptr == theInitialStatePtr) { deleteP << 
417       SortSecondaryFragment(ptr);              << 
418     }                                          << 
419     if (deletePrimary) { delete theInitialStat << 
420   }                                            << 
421   delete theTempResult; // end multifragmentat << 
422       }                                        << 
423     }                                          << 
424   }                                            << 
425   if (fVerbose > 2) {                          << 
426     G4cout << "## After first step of handler  << 
427            << " for evap;  "                   << 
428      << theResults.size() << " results. " << G << 
429   }                                            << 
430   // -----------------------------------       << 
431   // FermiBreakUp and De-excitation loop       << 
432   // -----------------------------------       << 
433                                                << 
434   static const G4int countmax = 1000;          << 
435   std::size_t kk;                              << 
436   for (kk=0; kk<theEvapList.size(); ++kk) {    << 
437     G4Fragment* frag = theEvapList[kk];        << 
438     if (fVerbose > 3) {                        << 
439       G4cout << "Next evaporate: " << G4endl;  << 
440       G4cout << *frag << G4endl;               << 
441     }                                          << 
442     if (kk >= countmax) {                      << 
443       G4ExceptionDescription ed;               << 
444       ed << "Infinite loop in the de-excitatio << 
445    << " iterations \n"                         << 
446    << "      Initial fragment: \n" << theIniti << 
447    << "\n      Current fragment: \n" << *frag; << 
448       G4Exception("G4ExcitationHandler::BreakI << 
449       ed,"Stop execution");                    << 
450                                                << 
451     }                                          << 
452     A = frag->GetA_asInt();                    << 
453     Z = frag->GetZ_asInt();                    << 
454     results.clear();                           << 
455     if (fVerbose > 2) {                        << 
456       G4cout << "G4ExcitationHandler# " << kk  << 
457              << " Eex(MeV)= " << frag->GetExci << 
458     }                                          << 
459     // Fermi Break-Up                          << 
460     if (theFermiModel->IsApplicable(Z, A, frag << 
461       theFermiModel->BreakFragment(&results, f << 
462       std::size_t nsec = results.size();       << 
463       if (fVerbose > 2) { G4cout << "FermiBrea << 
464                                                << 
465       // FBU takes care to delete input fragme << 
466       // The secondary may be excited - photo- << 
467       if (1 < nsec) {                          << 
468   for (auto const & res : results) {           << 
469     SortSecondaryFragment(res);                << 
470   }                                            << 
471   continue;                                    << 
472       }                                        << 
473       // evaporation will be applied           << 
474     }                                          << 
475     // apply Evaporation, residual nucleus is  << 
476     // photon evaporation is possible          << 
477     theEvaporation->BreakFragment(&results, fr << 
478     if (fVerbose > 3) {                        << 
479       G4cout << "Evaporation Nsec= " << result << 
480     }                                          << 
481     if (0 == results.size()) {                 << 
482       theResults.push_back(frag);              << 
483     } else {                                   << 
484       SortSecondaryFragment(frag);             << 
485     }                                          << 
486                                                << 
487     // Sort out secondary fragments            << 
488     for (auto const & res : results) {         << 
489       if(fVerbose > 4) {                       << 
490   G4cout << "Evaporated product #" << *res <<  << 
491       }                                        << 
492       SortSecondaryFragment(res);              << 
493     } // end of loop on secondary              << 
494   } // end of the loop over theEvapList        << 
495   if (fVerbose > 2) {                          << 
496     G4cout << "## After 2nd step of handler "  << 
497            << " was evap;  "                   << 
498      << theResults.size() << " results. " << G << 
499   }                                            << 
500   G4ReactionProductVector * theReactionProduct << 
501     new G4ReactionProductVector();             << 
502                                                   321 
503   // MAC (24/07/08)                            << 322   theResultList.clear();
504   // To optimise the storing speed, we reserve << 323   theFinalResultList.clear();
505   // in memory for the vector                  << 324   //AHtest  theFinalPhotonResultList.clear();
506   theReactionProductVector->reserve( theResult << 325   
507                                                << 326   
508   if (fVerbose > 1) {                          << 327 #ifdef debug
509     G4cout << "### ExcitationHandler provides  << 328   CheckConservation(theInitialState,theResult);
510      << " evaporated products:" << G4endl;     << 329 #endif
511   }                                            << 330   // Change G4FragmentVector by G4DynamicParticle
512   G4LorentzVector partOfLambdaLV;              << 331   return Transform(theResult);
513   if ( nL > 0 ) partOfLambdaLV = lambdaLV/(G4d << 332 }
514   for (auto const & frag : theResults) {       << 
515     G4LorentzVector lv0 = frag->GetMomentum(); << 
516     G4double etot = lv0.e();                   << 
517                                                << 
518     // in the case of dummy de-excitation, exc << 
519     // into kinetic energy of output ion       << 
520     if (!isActive) {                           << 
521       G4double mass = frag->GetGroundStateMass << 
522       G4double ptot = lv0.vect().mag();        << 
523       G4double fac  = (etot <= mass || 0.0 ==  << 
524   : std::sqrt((etot - mass)*(etot + mass))/pto << 
525       G4LorentzVector lv((frag->GetMomentum()) << 
526        (frag->GetMomentum()).py()*fac,         << 
527        (frag->GetMomentum()).pz()*fac, etot);  << 
528       frag->SetMomentum(lv);                   << 
529     }                                          << 
530     if (fVerbose > 3) {                        << 
531       G4cout << *frag;                         << 
532       if (frag->NuclearPolarization()) {       << 
533   G4cout << "  " << frag->NuclearPolarization( << 
534       }                                        << 
535       G4cout << G4endl;                        << 
536     }                                          << 
537                                                   333 
538     G4int fragmentA = frag->GetA_asInt();      << 334 G4ReactionProductVector * 
539     G4int fragmentZ = frag->GetZ_asInt();      << 335 G4ExcitationHandler::Transform(G4FragmentVector * theFragmentVector) const
540     G4double eexc = 0.0;                       << 336 {
541     const G4ParticleDefinition* theKindOfFragm << 337   if (theFragmentVector == 0) return 0;
542     G4bool isHyperN = false;                   << 338   
543     if (fragmentA == 0) {       // photon or e << 339   // Conversion from G4FragmentVector to G4ReactionProductVector
544       theKindOfFragment = frag->GetParticleDef << 340   G4ParticleDefinition *theGamma = G4Gamma::GammaDefinition();
545     } else if (fragmentA == 1 && fragmentZ ==  << 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
546       theKindOfFragment = theNeutron;             363       theKindOfFragment = theNeutron;
547     } else if (fragmentA == 1 && fragmentZ ==  << 364     } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
548       theKindOfFragment = theProton;              365       theKindOfFragment = theProton;
549     } else if (fragmentA == 2 && fragmentZ ==  << 366     } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
550       theKindOfFragment = theDeuteron;            367       theKindOfFragment = theDeuteron;
551     } else if (fragmentA == 3 && fragmentZ ==  << 368     } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
552       theKindOfFragment = theTriton;              369       theKindOfFragment = theTriton;
553       if(0 < nL) {                             << 370     } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
554         const G4ParticleDefinition* p = thePar << 371       theKindOfFragment = theHelium3;
555         if(nullptr != p) {                     << 372     } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha
556     theKindOfFragment = p;                     << 
557     isHyperN = true;                           << 
558     --nL;                                      << 
559   }                                            << 
560       }                                        << 
561     } else if (fragmentA == 3 && fragmentZ ==  << 
562       theKindOfFragment = theHe3;              << 
563     } else if (fragmentA == 4 && fragmentZ ==  << 
564       theKindOfFragment = theAlpha;               373       theKindOfFragment = theAlpha;
565       if (0 < nL) {                            << 
566         const G4ParticleDefinition* p = thePar << 
567         if(nullptr != p) {                     << 
568     theKindOfFragment = p;                     << 
569     isHyperN = true;                           << 
570     --nL;                                      << 
571   }                                            << 
572       }                                        << 
573     } else {                                      374     } else {
574                                                << 375       theKindOfFragment = theTableOfParticles->FindIon(theFragmentZ,theFragmentA,0,theFragmentZ);
575       // fragment                              << 
576       eexc = frag->GetExcitationEnergy();      << 
577       G4int idxf = frag->GetFloatingLevelNumbe << 
578       if (eexc < minExcitation) {              << 
579   eexc = 0.0;                                  << 
580         idxf = 0;                              << 
581       }                                        << 
582                                                << 
583       theKindOfFragment = theTableOfIons->GetI << 
584                                                << 
585       if (fVerbose > 3) {                      << 
586   G4cout << "### EXCH: Find ion Z= " << fragme << 
587                << " A= " << fragmentA          << 
588          << " Eexc(MeV)= " << eexc/MeV << " id << 
589          << " " << theKindOfFragment->GetParti << 
590          << G4endl;                            << 
591       }                                        << 
592     }                                             376     }
593     // fragment identified                     << 377     if (theKindOfFragment != 0) 
594     if (nullptr != theKindOfFragment) {        << 378       {
595       G4ReactionProduct * theNew = new G4React << 
596       if (isHyperN) {                          << 
597         G4LorentzVector lv = lv0 + partOfLambd << 
598   G4ThreeVector dir = lv.vect().unit();        << 
599         G4double mass = theKindOfFragment->Get << 
600         etot = std::max(lv.e(), mass);         << 
601         G4double ptot = std::sqrt((etot - mass << 
602         dir *= ptot;                           << 
603   theNew->SetMomentum(dir);                    << 
604   // remaining not compensated 4-momentum      << 
605         lambdaLV += (lv0 - G4LorentzVector(dir << 
606       } else {                                 << 
607   theNew->SetMomentum(lv0.vect());             << 
608       }                                        << 
609       theNew->SetTotalEnergy(etot);            << 
610       theNew->SetFormationTime(frag->GetCreati << 
611       theNew->SetCreatorModelID(frag->GetCreat << 
612       theReactionProductVector->push_back(theN << 
613                                                << 
614       // fragment not found out ground state i << 
615     } else {                                   << 
616       theKindOfFragment =                      << 
617   theTableOfIons->GetIon(fragmentZ,fragmentA,0 << 
618       if (theKindOfFragment) {                 << 
619   G4ThreeVector mom(0.0,0.0,0.0);              << 
620   G4double ionmass = theKindOfFragment->GetPDG << 
621   if (etot <= ionmass) {                       << 
622     etot = ionmass;                            << 
623   } else {                                     << 
624     G4double ptot = std::sqrt((etot - ionmass) << 
625     mom = (frag->GetMomentum().vect().unit())* << 
626   }                                            << 
627   G4ReactionProduct * theNew = new G4ReactionP    379   G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
628   theNew->SetMomentum(mom);                    << 380   theNew->SetMomentum(theFragmentMomentum.vect());
629   theNew->SetTotalEnergy(etot);                << 381   theNew->SetTotalEnergy(theFragmentMomentum.e());
630   theNew->SetFormationTime(frag->GetCreationTi << 382   theNew->SetFormationTime((*i)->GetCreationTime());
631   theNew->SetCreatorModelID(frag->GetCreatorMo << 383 #ifdef PRECOMPOUND_TEST
                                                   >> 384   theNew->SetCreatorModel((*i)->GetCreatorModel());
                                                   >> 385 #endif
632   theReactionProductVector->push_back(theNew);    386   theReactionProductVector->push_back(theNew);
633   if (fVerbose > 3) {                          << 
634     G4cout << "          ground state, energy  << 
635                  << etot << G4endl;            << 
636   }                                            << 
637       }                                           387       }
638     }                                          << 
639     delete frag;                               << 
640   }                                               388   }
641   // remaining lambdas are free; conserve quan << 389   if (theFragmentVector != 0)
642   // not 4-momentum                            << 390     { 
643   if (0 < nL) {                                << 391       std::for_each(theFragmentVector->begin(), theFragmentVector->end(), DeleteFragment());
644     G4ThreeVector dir = G4ThreeVector(0.0, 0.0 << 392       delete theFragmentVector;
645     if (lambdaLV.vect().mag() > CLHEP::eV) {   << 393     }
646       dir = lambdaLV.vect().unit();            << 394   G4ReactionProductVector::iterator debugit;
647     }                                          << 395   for(debugit=theReactionProductVector->begin(); 
648     G4double etot = std::max(lambdaLV.e()/(G4d << 396       debugit!=theReactionProductVector->end(); debugit++)
649     dir *= std::sqrt((etot - fLambdaMass)*(eto << 397     {
650     for (G4int i=0; i<nL; ++i) {               << 398     if((*debugit)->GetTotalEnergy()<1.*eV)
651       G4ReactionProduct* theNew = new G4Reacti << 399       {
652       theNew->SetMomentum(dir);                << 400   if(getenv("G4DebugPhotonevaporationData"))
653       theNew->SetTotalEnergy(etot);            << 401     {
654       theNew->SetFormationTime(theInitialState << 402       G4cerr << "G4ExcitationHandler: Warning: Photonevaporation data not exact."<<G4endl;
655       theNew->SetCreatorModelID(theInitialStat << 403       G4cerr << "G4ExcitationHandler: Warning: Found gamma with energy = "
656       theReactionProductVector->push_back(theN << 404        << (*debugit)->GetTotalEnergy()/MeV << "MeV"
657     }                                          << 405        << G4endl;
658   }                                            << 406     }
659   if (fVerbose > 3) {                          << 407   delete (*debugit);
660     G4cout << "@@@@@@@@@@ End G4Excitation Han << 408   *debugit = 0;
                                                   >> 409       }
661   }                                               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());
662   return theReactionProductVector;                417   return theReactionProductVector;
663 }                                                 418 }
664                                                   419 
665 void G4ExcitationHandler::ModelDescription(std << 420 
666 {                                              << 421 #ifdef debug
667   outFile << "G4ExcitationHandler description\ << 422 void G4ExcitationHandler::CheckConservation(const G4Fragment & theInitialState,
668     << "This class samples de-excitation of ex << 423               G4FragmentVector * Result) const
669     << "Fermi Break-up model for light fragmen << 424 {
670     << "evaporation, fission, and photo-evapor << 425   G4double ProductsEnergy =0;
671     << "particle may be proton, neutron, and o << 426   G4ThreeVector ProductsMomentum;
672     << "(Z < 13, A < 29). During photon evapor << 427   G4int ProductsA = 0;
673     << "or electrons due to internal conversio << 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;
674 }                                                 473 }
                                                   >> 474 #endif
                                                   >> 475 
675                                                   476 
676                                                   477 
677                                                   478 
678                                                   479