| Geant4 Cross Reference |
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 // $Id: G4ExcitationHandler.cc,v 1.26.2.1 2010/04/01 09:41:42 gcosmo Exp $
>> 27 // GEANT4 tag $Name: geant4-09-03-patch-02 $
>> 28 //
26 // Hadronic Process: Nuclear De-excitations 29 // Hadronic Process: Nuclear De-excitations
27 // by V. Lara (May 1998) 30 // by V. Lara (May 1998)
28 // 31 //
29 // 32 //
30 // Modified: 33 // Modified:
31 // 30 June 1998 by V. Lara: << 34 // (30 June 1998) by V. Lara:
32 // -Modified the Transform method for use 35 // -Modified the Transform method for use G4ParticleTable and
33 // therefore G4IonTable. It makes possib 36 // therefore G4IonTable. It makes possible to convert all kind
34 // of fragments (G4Fragment) produced in 37 // of fragments (G4Fragment) produced in deexcitation to
35 // G4DynamicParticle 38 // G4DynamicParticle
36 // -It uses default algorithms for: 39 // -It uses default algorithms for:
37 // Evaporation: G4Evaporation 40 // Evaporation: G4Evaporation
38 // MultiFragmentation: G4StatMF 41 // MultiFragmentation: G4StatMF
39 // Fermi Breakup model: G4FermiBr 42 // Fermi Breakup model: G4FermiBreakUp
40 // 24 Jul 2008 by M. A. Cortes Giraldo: << 43 // (24 Jul 2008) by M. A. Cortes Giraldo:
41 // -Max Z,A for Fermi Break-Up turns to 9 44 // -Max Z,A for Fermi Break-Up turns to 9,17 by default
42 // -BreakItUp() reorganised and bug in Ev 45 // -BreakItUp() reorganised and bug in Evaporation loop fixed
43 // -Transform() optimised 46 // -Transform() optimised
44 // (September 2008) by J. M. Quesada. External 47 // (September 2008) by J. M. Quesada. External choices have been added for :
45 // -inverse cross section option (default 48 // -inverse cross section option (default OPTxs=3)
46 // -superimposed Coulomb barrier (if useS 49 // -superimposed Coulomb barrier (if useSICB is set true, by default it is false)
47 // September 2009 by J. M. Quesada: << 50 // (September 2009) by J. M. Quesada:
48 // -according to Igor Pshenichnov, SMM wi 51 // -according to Igor Pshenichnov, SMM will be applied (just in case) only once.
49 // 27 Nov 2009 by V.Ivanchenko: << 52 // (27 Nov 2009) by V.Ivanchenko:
50 // -cleanup the logic, reduce number inte 53 // -cleanup the logic, reduce number internal vectors, fixed memory leak.
51 // 11 May 2010 by V.Ivanchenko: << 54 //
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 55
63 #include "G4ExcitationHandler.hh" 56 #include "G4ExcitationHandler.hh"
64 #include "G4SystemOfUnits.hh" << 57 #include "globals.hh"
65 #include "G4LorentzVector.hh" 58 #include "G4LorentzVector.hh"
66 #include "G4ThreeVector.hh" << 59 #include <list>
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 60
112 fLambdaMass = theLambda->GetPDGMass(); << 61 //#define debugphoton
113 62
114 if(fVerbose > 1) { G4cout << "### New handle <<
115 } <<
116 63
117 G4ExcitationHandler::~G4ExcitationHandler() << 64 G4ExcitationHandler::G4ExcitationHandler():
118 { << 65 // JMQ 160909 Fermi BreakUp & MultiFrag are on by default
119 delete theMultiFragmentation; << 66 // This is needed for activation of such models when G4BinaryLightIonReaction is used
120 delete theFermiModel; << 67 // since no interface (for external activation via macro input file) is still available.
121 if(isEvapLocal) { delete theEvaporation; } << 68 maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(4.0*GeV),
>> 69 // maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(3.0*MeV),
>> 70 //maxZForFermiBreakUp(1),maxAForFermiBreakUp(1),minEForMultiFrag(4.0*GeV),
>> 71 MyOwnEvaporationClass(true), MyOwnMultiFragmentationClass(true),MyOwnFermiBreakUpClass(true),
>> 72 MyOwnPhotonEvaporationClass(true),OPTxs(3),useSICB(false)
>> 73 {
>> 74 theTableOfParticles = G4ParticleTable::GetParticleTable();
>> 75
>> 76 theEvaporation = new G4Evaporation;
>> 77 theMultiFragmentation = new G4StatMF;
>> 78 theFermiModel = new G4FermiBreakUp;
>> 79 thePhotonEvaporation = new G4PhotonEvaporation;
122 } 80 }
123 81
124 void G4ExcitationHandler::SetParameters() << 82 G4ExcitationHandler::G4ExcitationHandler(const G4ExcitationHandler &)
125 { 83 {
126 G4NuclearLevelData* ndata = G4NuclearLevelDa << 84 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 } 85 }
164 86
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 <<
183 void G4ExcitationHandler::SetEvaporation(G4VEv <<
184 { <<
185 if(nullptr != ptr && ptr != theEvaporation) <<
186 theEvaporation = ptr; <<
187 theEvaporation->SetPhotonEvaporation(thePh <<
188 theEvaporation->SetFermiBreakUp(theFermiMo <<
189 isEvapLocal = flag; <<
190 if(fVerbose > 1) { <<
191 G4cout << "G4ExcitationHandler::SetEvapo <<
192 } <<
193 } <<
194 } <<
195 87
196 void << 88 G4ExcitationHandler::~G4ExcitationHandler()
197 G4ExcitationHandler::SetMultiFragmentation(G4V <<
198 { 89 {
199 if(nullptr != ptr && ptr != theMultiFragment << 90 if (MyOwnEvaporationClass) delete theEvaporation;
200 delete theMultiFragmentation; << 91 if (MyOwnMultiFragmentationClass) delete theMultiFragmentation;
201 theMultiFragmentation = ptr; << 92 if (MyOwnFermiBreakUpClass) delete theFermiModel;
202 } << 93 if (MyOwnPhotonEvaporationClass) delete thePhotonEvaporation;
203 } 94 }
204 95
205 void G4ExcitationHandler::SetFermiModel(G4VFer <<
206 { <<
207 if(nullptr != ptr && ptr != theFermiModel) { <<
208 delete theFermiModel; <<
209 theFermiModel = ptr; <<
210 if(nullptr != theEvaporation) { <<
211 theEvaporation->SetFermiBreakUp(theFermi <<
212 } <<
213 } <<
214 } <<
215 96
216 void << 97 const G4ExcitationHandler & G4ExcitationHandler::operator=(const G4ExcitationHandler &)
217 G4ExcitationHandler::SetPhotonEvaporation(G4VE <<
218 { 98 {
219 if(nullptr != ptr && ptr != thePhotonEvapora << 99 throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator=: is meant to not be accessable! ");
220 delete thePhotonEvaporation; << 100
221 thePhotonEvaporation = ptr; << 101 return *this;
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 } <<
231 <<
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 } 102 }
263 103
264 G4VEvaporation* G4ExcitationHandler::GetEvapor <<
265 { <<
266 if (nullptr != theEvaporation) { SetParamete <<
267 return theEvaporation; <<
268 } <<
269 104
270 G4VMultiFragmentation* G4ExcitationHandler::Ge << 105 G4bool G4ExcitationHandler::operator==(const G4ExcitationHandler &) const
271 { 106 {
272 if (nullptr != theMultiFragmentation) { SetP << 107 throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator==: is meant to not be accessable! ");
273 return theMultiFragmentation; << 108 return false;
274 } << 109 }
275 110
276 G4VFermiBreakUp* G4ExcitationHandler::GetFermi << 111 G4bool G4ExcitationHandler::operator!=(const G4ExcitationHandler &) const
277 { 112 {
278 if (nullptr != theFermiModel) { SetParameter << 113 throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator!=: is meant to not be accessable! ");
279 return theFermiModel; << 114 return true;
280 } 115 }
281 116
282 G4VEvaporationChannel* G4ExcitationHandler::Ge << 117 ////////////////////////////////////////////////////////////////////////////////////////////////
283 { << 118 /// 25/07/08 16:45 Proposed by MAC ////////////////////////////////////////////////////////////
284 if(nullptr != thePhotonEvaporation) { SetPar << 119 ////////////////////////////////////////////////////////////////////////////////////////////////
285 return thePhotonEvaporation; <<
286 } <<
287 120
288 G4ReactionProductVector * << 121 G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) const
289 G4ExcitationHandler::BreakItUp(const G4Fragmen << 122 {
290 { << 123 //for inverse cross section choice
>> 124 theEvaporation->SetOPTxs(OPTxs);
>> 125 //for the choice of superimposed Coulomb Barrier for inverse cross sections
>> 126 theEvaporation->UseSICB(useSICB);
>> 127
>> 128 // Pointer which will be used to return the final production vector
>> 129 //G4FragmentVector * theResult = new G4FragmentVector;
>> 130
291 // Variables existing until end of method 131 // Variables existing until end of method
>> 132 //G4Fragment * theInitialStatePtr = const_cast<G4Fragment*>(&theInitialState);
292 G4Fragment * theInitialStatePtr = new G4Frag 133 G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
293 if (fVerbose > 1) { << 134 G4FragmentVector * theTempResult = 0; // pointer which receives temporal results
294 G4cout << "@@@@@@@@@@ Start G4Excitation H << 135 std::list<G4Fragment*> theEvapList; // list to apply Evaporation, SMF or Fermi Break-Up
295 G4cout << theInitialState << G4endl; << 136 std::list<G4Fragment*> theEvapStableList; // list to apply PhotonEvaporation
296 } << 137 std::list<G4Fragment*> theResults; // list to store final result
297 if (!isInitialised) { Initialise(); } << 138 std::list<G4Fragment*>::iterator iList;
298 << 139 //
299 // pointer to fragment vector which receives << 140 //G4cout << "@@@@@@@@@@ Start G4Exitation Handler @@@@@@@@@@@@@" << G4endl;
300 G4FragmentVector * theTempResult = nullptr; << 141 //G4cout << theInitialState << G4endl;
301 << 142
302 theResults.clear(); <<
303 theEvapList.clear(); <<
304 <<
305 // Variables to describe the excited configu 143 // Variables to describe the excited configuration
306 G4double exEnergy = theInitialState.GetExcit 144 G4double exEnergy = theInitialState.GetExcitationEnergy();
307 G4int A = theInitialState.GetA_asInt(); << 145 G4int A = static_cast<G4int>( theInitialState.GetA() +0.5 );
308 G4int Z = theInitialState.GetZ_asInt(); << 146 G4int Z = static_cast<G4int>( theInitialState.GetZ() +0.5 );
309 G4int nL = theInitialState.GetNumberOfLambda << 147
310 << 148 // JMQ 150909: first step in de-excitation chain (SMM will be used only here)
311 // too much excitation << 149
312 if (exEnergy > A*maxExcitation && A > 0) { << 150 // In case A <= 4 the fragment will not perform any nucleon emission
313 ++fWarnings; << 151 if (A <= 4)
314 if(fWarnings < 0) { << 152 {
315 G4ExceptionDescription ed; << 153 // I store G4Fragment* in theEvapStableList to apply thePhotonEvaporation later
316 ed << "High excitation Fragment Z= " << << 154 theEvapStableList.push_back( theInitialStatePtr );
317 << " Eex/A(MeV)= " << exEnergy/A; <<
318 G4Exception("G4ExcitationHandler::BreakI <<
319 } 155 }
320 } << 156
321 << 157 else // If A > 4 we try to apply theFermiModel, theMultiFragmentation or theEvaporation
322 // for hyper-nuclei subtract lambdas from th << 158 {
323 G4double lambdaF = 0.0; << 159
324 G4LorentzVector lambdaLV = theInitialStatePt << 160 // JMQ 150909: first step in de-excitation is treated separately
325 if (0 < nL) { << 161 // Fragments after the first step are stored in theEvapList
326 << 162 // Statistical Multifragmentation will take place (just in case) only here
327 // is it a stable hyper-nuclei? << 163 //
328 if(A >= 3 && A <= 5 && nL <= 2) { << 164 // Test applicability
329 G4int pdg = 0; << 165 // Initial State De-Excitation
330 if(3 == A && 1 == nL) { << 166 if(A<GetMaxA()&&Z<GetMaxZ())
331 pdg = 1010010030; << 167 {
332 } else if(5 == A && 2 == Z && 1 == nL) { << 168 theTempResult = theFermiModel->BreakItUp(theInitialState);
333 pdg = 1010020050; << 169 }
334 } else if(4 == A) { << 170 else if (exEnergy>GetMinE()*A)
335 if(1 == Z && 1 == nL) { << 171 {
336 pdg = 1010010040; << 172 theTempResult = theMultiFragmentation->BreakItUp(theInitialState);
337 } else if(2 == Z && 1 == nL) { << 173 }
338 pdg = 1010020040; << 174 else
339 } else if(0 == Z && 2 == nL) { << 175 {
340 pdg = 1020000040; << 176 theTempResult = theEvaporation->BreakItUp(theInitialState);
341 } else if(1 == Z && 2 == nL) { << 177 }
342 pdg = 1020010040; << 178
343 } << 179 G4bool deletePrimary = true;
344 } << 180 if(theTempResult->size() > 0)
345 // initial state is one of hyper-nuclei << 181 {
346 if (0 < pdg) { << 182 // Store original state in theEvapList
347 const G4ParticleDefinition* part = thePartTa << 183 G4FragmentVector::iterator j;
348 if(nullptr != part) { << 184 for (j = theTempResult->begin(); j != theTempResult->end(); ++j)
349 G4ReactionProduct* theNew = new G4Reaction << 185 {
350 G4ThreeVector dir = G4ThreeVector( 0.0, 0. << 186 if((*j) == theInitialStatePtr) { deletePrimary = false; }
351 if ( lambdaLV.vect().mag() > CLHEP:: << 187 A = static_cast<G4int>((*j)->GetA()+0.5); // +0.5 to avoid bad truncation
352 dir = lambdaLV.vect().unit(); << 188
353 } << 189 if(A <= 1) { theResults.push_back(*j); } // gamma, p, n
354 G4double mass = part->GetPDGMass(); << 190 else if(A <= 4) { theEvapStableList.push_back(*j); } // evaporation is not possible
355 G4double etot = std::max(lambdaLV.e(), mas << 191 else { theEvapList.push_back(*j); } // evaporation is possible
356 dir *= std::sqrt((etot - mass)*(etot << 192 }
357 theNew->SetMomentum(dir); <<
358 theNew->SetTotalEnergy(etot); <<
359 theNew->SetFormationTime(theInitialState.G <<
360 theNew->SetCreatorModelID(theInitialState. <<
361 G4ReactionProductVector* v = new G4Reactio <<
362 v->push_back(theNew); <<
363 return v; <<
364 } 193 }
365 } << 194 if( deletePrimary ) { delete theInitialStatePtr; }
>> 195 delete theTempResult;
366 } 196 }
367 G4double mass = theInitialStatePtr->GetGro << 197 //
368 lambdaF = nL*(fLambdaMass - CLHEP::neutron << 198 // JMQ 150909: Further steps in de-excitation chain follow ..
>> 199
>> 200 //G4cout << "## After first step " << theEvapList.size() << " for evap; "
>> 201 // << theEvapStableList.size() << " for photo-evap; "
>> 202 // << theResults.size() << " results. " << G4endl;
>> 203
>> 204 // ------------------------------
>> 205 // De-excitation loop
>> 206 // ------------------------------
>> 207
>> 208 for (iList = theEvapList.begin(); iList != theEvapList.end(); ++iList)
>> 209 {
>> 210 A = static_cast<G4int>((*iList)->GetA()+0.5); // +0.5 to avoid bad truncation
>> 211 Z = static_cast<G4int>((*iList)->GetZ()+0.5);
>> 212
>> 213 // In case A <= 4 the fragment will not perform any nucleon emission
>> 214 if (A <= 4)
>> 215 {
>> 216 // storing G4Fragment* in theEvapStableList to apply thePhotonEvaporation later
>> 217 theEvapStableList.push_back(*iList );
>> 218 }
>> 219
>> 220 else // If A > 4 we try to apply theFermiModel or theEvaporation
>> 221 {
>> 222 // stable fragment
>> 223 if ((*iList)->GetExcitationEnergy() <= 0.1*eV)
>> 224 {
>> 225 theResults.push_back(*iList);
>> 226 }
>> 227 else
>> 228 {
>> 229 if ( A < GetMaxA() && Z < GetMaxZ() ) // if satisfied apply Fermi Break-Up
>> 230 {
>> 231 theTempResult = theFermiModel->BreakItUp(*(*iList));
>> 232 }
>> 233 else // apply Evaporation in another case
>> 234 {
>> 235 theTempResult = theEvaporation->BreakItUp(*(*iList));
>> 236 }
>> 237
>> 238 // New configuration is stored in theTempResult, so we can free
>> 239 // the memory where the previous configuration is
>> 240
>> 241 G4bool deletePrimary = true;
>> 242 G4int nsec = theTempResult->size();
>> 243
>> 244 // The number of secondaries tells us if the configuration has changed
>> 245 if ( nsec > 0 )
>> 246 {
>> 247 G4FragmentVector::iterator j;
>> 248 for (j = theTempResult->begin(); j != theTempResult->end(); ++j)
>> 249 {
>> 250 if((*j) == (*iList)) { deletePrimary = false; }
>> 251 A = static_cast<G4int>((*j)->GetA()+0.5); // +0.5 to avoid bad truncation
>> 252
>> 253 if(A <= 1) { theResults.push_back(*j); } // gamma, p, n
>> 254 else if(A <= 4 || 1 == nsec) { theEvapStableList.push_back(*j); } // no evaporation
>> 255 else { theEvapList.push_back(*j); }
>> 256 }
>> 257 }
>> 258 if( deletePrimary ) { delete (*iList); }
>> 259 delete theTempResult;
>> 260
>> 261 }
>> 262 } // endif (A <=4)
>> 263 } // end of the loop over theEvapList
>> 264
>> 265 //G4cout << "## After 2nd step " << theEvapList.size() << " was evap; "
>> 266 // << theEvapStableList.size() << " for photo-evap; "
>> 267 // << theResults.size() << " results. " << G4endl;
>> 268
>> 269 // -----------------------
>> 270 // Photon-Evaporation loop
>> 271 // -----------------------
>> 272
>> 273 for (iList = theEvapStableList.begin(); iList != theEvapStableList.end(); ++iList)
>> 274 {
>> 275 // take out stable particles and fragments
>> 276 A = static_cast<G4int>((*iList)->GetA()+0.5);
>> 277 if ( A <= 1 ) { theResults.push_back(*iList); }
>> 278 else if ((*iList)->GetExcitationEnergy() <= 0.1*eV) { theResults.push_back(*iList); }
>> 279
>> 280 else
>> 281 {
>> 282 // photon-evaporation is applied
>> 283 theTempResult = thePhotonEvaporation->BreakItUp(*(*iList));
>> 284
>> 285 G4bool deletePrimary = true;
>> 286 G4int nsec = theTempResult->size();
>> 287
>> 288 // if there is a gamma emission then
>> 289 if (nsec > 1)
>> 290 {
>> 291 G4FragmentVector::iterator j;
>> 292 for (j = theTempResult->begin(); j != theTempResult->end(); ++j)
>> 293 {
>> 294 if((*j) == (*iList)) { deletePrimary = false; }
>> 295 A = static_cast<G4int>((*j)->GetA()+0.5); // +0.5 to avoid bad truncation
>> 296
>> 297 if(A <= 1) { theResults.push_back(*j); } // gamma, p, n
>> 298 else if((*j)->GetExcitationEnergy() <= 0.1*eV) { theResults.push_back(*j); } // stable fragment
>> 299 else { theEvapStableList.push_back(*j); }
>> 300 }
>> 301 }
>> 302
>> 303 else if(1 == nsec)
>> 304 {
>> 305 G4FragmentVector::iterator j = theTempResult->begin();
>> 306 if((*j) == (*iList)) { deletePrimary = false; }
>> 307 // Let's create a G4Fragment pointer representing the gamma emmited
>> 308 G4LorentzVector lv = (*j)->GetMomentum();
>> 309 G4double Mass = (*j)->GetGroundStateMass();
>> 310 G4double Ecm = lv.m();
>> 311 if(Ecm - Mass > 0.1*eV)
>> 312 {
>> 313 G4ThreeVector bst = lv.boostVector();
>> 314 G4double GammaEnergy = 0.5*(Ecm - Mass)*(Ecm + Mass)/Ecm;
>> 315 G4double cosTheta = 1. - 2. * G4UniformRand();
>> 316 G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
>> 317 G4double phi = twopi * G4UniformRand();
>> 318 G4LorentzVector Gamma4P(GammaEnergy * sinTheta * std::cos(phi),
>> 319 GammaEnergy * sinTheta * std::sin(phi),
>> 320 GammaEnergy * cosTheta,
>> 321 GammaEnergy);
>> 322 Gamma4P.boost(bst);
>> 323 G4Fragment * theHandlerPhoton = new G4Fragment(Gamma4P,G4Gamma::GammaDefinition());
>> 324 theResults.push_back(theHandlerPhoton);
>> 325
>> 326 // And now we update momentum and energy for the nucleus
>> 327 lv -= Gamma4P;
>> 328 (*j)->SetMomentum(lv); // Now this fragment has been deexcited!
>> 329 }
>> 330 // we store the deexcited fragment
>> 331 theResults.push_back(*j);
>> 332 }
>> 333 if( deletePrimary ) { delete (*iList); }
>> 334 delete theTempResult;
>> 335 }
>> 336 } // end of photon-evaporation loop
>> 337
>> 338 //G4cout << "## After 3d step " << theEvapList.size() << " was evap; "
>> 339 // << theEvapStableList.size() << " was photo-evap; "
>> 340 // << theResults.size() << " results. " << G4endl;
>> 341
>> 342 #ifdef debug
>> 343 CheckConservation(theInitialState,*theResults);
>> 344 #endif
369 345
370 // de-excitation with neutrons instead of << 346 G4ReactionProductVector * theReactionProductVector = new G4ReactionProductVector;
371 theInitialStatePtr->SetZAandMomentum(lambd <<
372 347
373 // 4-momentum not used in de-excitation << 348 // MAC (24/07/08)
374 lambdaLV *= lambdaF; << 349 // To optimise the storing speed, we reserve space in memory for the vector
375 } else if (0 > nL) { << 350 theReactionProductVector->reserve( theResults.size() );
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 } <<
383 } <<
384 351
385 // In case A <= 1 the fragment will not perf << 352 G4int theFragmentA, theFragmentZ;
386 if (A <= 1 || !isActive) { << 353 G4LorentzVector theFragmentMomentum;
387 theResults.push_back( theInitialStatePtr ) <<
388 <<
389 // check if a fragment is stable <<
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 354
400 // Statistical Multifragmentation will tak << 355 std::list<G4Fragment*>::iterator i;
401 } else { << 356 for (i = theResults.begin(); i != theResults.end(); ++i)
402 theTempResult = theMultiFragmentation->B << 357 {
403 if (nullptr == theTempResult) { << 358 theFragmentA = static_cast<G4int>((*i)->GetA());
404 theEvapList.push_back(theInitialStatePtr); << 359 theFragmentZ = static_cast<G4int>((*i)->GetZ());
>> 360 theFragmentMomentum = (*i)->GetMomentum();
>> 361 G4ParticleDefinition* theKindOfFragment = 0;
>> 362 if (theFragmentA == 0 && theFragmentZ == 0) { // photon
>> 363 theKindOfFragment = G4Gamma::GammaDefinition();
>> 364 } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
>> 365 theKindOfFragment = G4Neutron::NeutronDefinition();
>> 366 } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
>> 367 theKindOfFragment = G4Proton::ProtonDefinition();
>> 368 } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
>> 369 theKindOfFragment = G4Deuteron::DeuteronDefinition();
>> 370 } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
>> 371 theKindOfFragment = G4Triton::TritonDefinition();
>> 372 } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
>> 373 theKindOfFragment = G4He3::He3Definition();
>> 374 } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha
>> 375 theKindOfFragment = G4Alpha::AlphaDefinition();;
405 } else { 376 } else {
406 std::size_t nsec = theTempResult->size(); << 377 theKindOfFragment = theTableOfParticles->FindIon(theFragmentZ,theFragmentA,0,theFragmentZ);
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 } 378 }
423 } << 379 if (theKindOfFragment != 0)
424 } << 380 {
425 if (fVerbose > 2) { << 381 G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
426 G4cout << "## After first step of handler << 382 theNew->SetMomentum(theFragmentMomentum.vect());
427 << " for evap; " << 383 theNew->SetTotalEnergy(theFragmentMomentum.e());
428 << theResults.size() << " results. " << G << 384 theNew->SetFormationTime((*i)->GetCreationTime());
429 } << 385 theReactionProductVector->push_back(theNew);
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 } 386 }
471 continue; << 387 delete (*i);
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 } 388 }
486 389
487 // Sort out secondary fragments << 390 return theReactionProductVector;
488 for (auto const & res : results) { << 391 }
489 if(fVerbose > 4) { << 392
490 G4cout << "Evaporated product #" << *res << << 393
491 } << 394
492 SortSecondaryFragment(res); << 395 G4ReactionProductVector *
493 } // end of loop on secondary << 396 G4ExcitationHandler::Transform(G4FragmentVector * theFragmentVector) const
494 } // end of the loop over theEvapList << 397 {
495 if (fVerbose > 2) { << 398 if (theFragmentVector == 0) return 0;
496 G4cout << "## After 2nd step of handler " << 399
497 << " was evap; " << 400 // Conversion from G4FragmentVector to G4ReactionProductVector
498 << theResults.size() << " results. " << G << 401 G4ParticleDefinition *theGamma = G4Gamma::GammaDefinition();
499 } << 402 G4ParticleDefinition *theNeutron = G4Neutron::NeutronDefinition();
500 G4ReactionProductVector * theReactionProduct << 403 G4ParticleDefinition *theProton = G4Proton::ProtonDefinition();
501 new G4ReactionProductVector(); << 404 G4ParticleDefinition *theDeuteron = G4Deuteron::DeuteronDefinition();
>> 405 G4ParticleDefinition *theTriton = G4Triton::TritonDefinition();
>> 406 G4ParticleDefinition *theHelium3 = G4He3::He3Definition();
>> 407 G4ParticleDefinition *theAlpha = G4Alpha::AlphaDefinition();
>> 408 G4ParticleDefinition *theKindOfFragment = 0;
>> 409 theNeutron->SetVerboseLevel(2);
>> 410 G4ReactionProductVector * theReactionProductVector = new G4ReactionProductVector;
502 411
503 // MAC (24/07/08) 412 // MAC (24/07/08)
504 // To optimise the storing speed, we reserve << 413 // To optimise the storing speed, we reserve space in memory for the vector
505 // in memory for the vector << 414 theReactionProductVector->reserve( theFragmentVector->size() * sizeof(G4ReactionProduct*) );
506 theReactionProductVector->reserve( theResult <<
507 415
508 if (fVerbose > 1) { << 416 G4int theFragmentA, theFragmentZ;
509 G4cout << "### ExcitationHandler provides << 417 G4LorentzVector theFragmentMomentum;
510 << " evaporated products:" << G4endl; <<
511 } <<
512 G4LorentzVector partOfLambdaLV; <<
513 if ( nL > 0 ) partOfLambdaLV = lambdaLV/(G4d <<
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 418
538 G4int fragmentA = frag->GetA_asInt(); << 419 G4FragmentVector::iterator i;
539 G4int fragmentZ = frag->GetZ_asInt(); << 420 for (i = theFragmentVector->begin(); i != theFragmentVector->end(); i++) {
540 G4double eexc = 0.0; << 421 // std::cout << (*i) <<'\n';
541 const G4ParticleDefinition* theKindOfFragm << 422 theFragmentA = static_cast<G4int>((*i)->GetA());
542 G4bool isHyperN = false; << 423 theFragmentZ = static_cast<G4int>((*i)->GetZ());
543 if (fragmentA == 0) { // photon or e << 424 theFragmentMomentum = (*i)->GetMomentum();
544 theKindOfFragment = frag->GetParticleDef << 425 theKindOfFragment = 0;
545 } else if (fragmentA == 1 && fragmentZ == << 426 if (theFragmentA == 0 && theFragmentZ == 0) { // photon
>> 427 theKindOfFragment = theGamma;
>> 428 } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
546 theKindOfFragment = theNeutron; 429 theKindOfFragment = theNeutron;
547 } else if (fragmentA == 1 && fragmentZ == << 430 } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
548 theKindOfFragment = theProton; 431 theKindOfFragment = theProton;
549 } else if (fragmentA == 2 && fragmentZ == << 432 } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
550 theKindOfFragment = theDeuteron; 433 theKindOfFragment = theDeuteron;
551 } else if (fragmentA == 3 && fragmentZ == << 434 } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
552 theKindOfFragment = theTriton; 435 theKindOfFragment = theTriton;
553 if(0 < nL) { << 436 } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
554 const G4ParticleDefinition* p = thePar << 437 theKindOfFragment = theHelium3;
555 if(nullptr != p) { << 438 } 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; 439 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 { 440 } else {
574 << 441 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 } 442 }
593 // fragment identified << 443 if (theKindOfFragment != 0)
594 if (nullptr != theKindOfFragment) { << 444 {
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 445 G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
628 theNew->SetMomentum(mom); << 446 theNew->SetMomentum(theFragmentMomentum.vect());
629 theNew->SetTotalEnergy(etot); << 447 theNew->SetTotalEnergy(theFragmentMomentum.e());
630 theNew->SetFormationTime(frag->GetCreationTi << 448 theNew->SetFormationTime((*i)->GetCreationTime());
631 theNew->SetCreatorModelID(frag->GetCreatorMo << 449 #ifdef PRECOMPOUND_TEST
>> 450 theNew->SetCreatorModel((*i)->GetCreatorModel());
>> 451 #endif
632 theReactionProductVector->push_back(theNew); 452 theReactionProductVector->push_back(theNew);
633 if (fVerbose > 3) { <<
634 G4cout << " ground state, energy <<
635 << etot << G4endl; <<
636 } <<
637 } 453 }
638 } <<
639 delete frag; <<
640 } <<
641 // remaining lambdas are free; conserve quan <<
642 // not 4-momentum <<
643 if (0 < nL) { <<
644 G4ThreeVector dir = G4ThreeVector(0.0, 0.0 <<
645 if (lambdaLV.vect().mag() > CLHEP::eV) { <<
646 dir = lambdaLV.vect().unit(); <<
647 } <<
648 G4double etot = std::max(lambdaLV.e()/(G4d <<
649 dir *= std::sqrt((etot - fLambdaMass)*(eto <<
650 for (G4int i=0; i<nL; ++i) { <<
651 G4ReactionProduct* theNew = new G4Reacti <<
652 theNew->SetMomentum(dir); <<
653 theNew->SetTotalEnergy(etot); <<
654 theNew->SetFormationTime(theInitialState <<
655 theNew->SetCreatorModelID(theInitialStat <<
656 theReactionProductVector->push_back(theN <<
657 } <<
658 } 454 }
659 if (fVerbose > 3) { << 455 if (theFragmentVector != 0)
660 G4cout << "@@@@@@@@@@ End G4Excitation Han << 456 {
>> 457 std::for_each(theFragmentVector->begin(), theFragmentVector->end(), DeleteFragment());
>> 458 delete theFragmentVector;
>> 459 }
>> 460 G4ReactionProductVector::iterator debugit;
>> 461 for(debugit=theReactionProductVector->begin();
>> 462 debugit!=theReactionProductVector->end(); debugit++)
>> 463 {
>> 464 if((*debugit)->GetTotalEnergy()<1.*eV)
>> 465 {
>> 466 if(getenv("G4DebugPhotonevaporationData"))
>> 467 {
>> 468 G4cerr << "G4ExcitationHandler: Warning: Photonevaporation data not exact."<<G4endl;
>> 469 G4cerr << "G4ExcitationHandler: Warning: Found gamma with energy = "
>> 470 << (*debugit)->GetTotalEnergy()/MeV << "MeV"
>> 471 << G4endl;
>> 472 }
>> 473 delete (*debugit);
>> 474 *debugit = 0;
>> 475 }
661 } 476 }
>> 477 G4ReactionProduct* tmpPtr=0;
>> 478 theReactionProductVector->erase(std::remove_if(theReactionProductVector->begin(),
>> 479 theReactionProductVector->end(),
>> 480 std::bind2nd(std::equal_to<G4ReactionProduct*>(),
>> 481 tmpPtr)),
>> 482 theReactionProductVector->end());
662 return theReactionProductVector; 483 return theReactionProductVector;
663 } 484 }
664 485
665 void G4ExcitationHandler::ModelDescription(std << 486
666 { << 487 #ifdef debug
667 outFile << "G4ExcitationHandler description\ << 488 void G4ExcitationHandler::CheckConservation(const G4Fragment & theInitialState,
668 << "This class samples de-excitation of ex << 489 G4FragmentVector * Result) const
669 << "Fermi Break-up model for light fragmen << 490 {
670 << "evaporation, fission, and photo-evapor << 491 G4double ProductsEnergy =0;
671 << "particle may be proton, neutron, and o << 492 G4ThreeVector ProductsMomentum;
672 << "(Z < 13, A < 29). During photon evapor << 493 G4int ProductsA = 0;
673 << "or electrons due to internal conversio << 494 G4int ProductsZ = 0;
>> 495 G4FragmentVector::iterator h;
>> 496 for (h = Result->begin(); h != Result->end(); h++) {
>> 497 G4LorentzVector tmp = (*h)->GetMomentum();
>> 498 ProductsEnergy += tmp.e();
>> 499 ProductsMomentum += tmp.vect();
>> 500 ProductsA += static_cast<G4int>((*h)->GetA());
>> 501 ProductsZ += static_cast<G4int>((*h)->GetZ());
>> 502 }
>> 503
>> 504 if (ProductsA != theInitialState.GetA()) {
>> 505 G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl;
>> 506 G4cout << "G4ExcitationHandler.cc: Barionic Number Conservation test for deexcitation fragments"
>> 507 << G4endl;
>> 508 G4cout << "Initial A = " << theInitialState.GetA()
>> 509 << " Fragments A = " << ProductsA << " Diference --> "
>> 510 << theInitialState.GetA() - ProductsA << G4endl;
>> 511 }
>> 512 if (ProductsZ != theInitialState.GetZ()) {
>> 513 G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
>> 514 G4cout << "G4ExcitationHandler.cc: Charge Conservation test for deexcitation fragments"
>> 515 << G4endl;
>> 516 G4cout << "Initial Z = " << theInitialState.GetZ()
>> 517 << " Fragments Z = " << ProductsZ << " Diference --> "
>> 518 << theInitialState.GetZ() - ProductsZ << G4endl;
>> 519 }
>> 520 if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) {
>> 521 G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
>> 522 G4cout << "G4ExcitationHandler.cc: Energy Conservation test for deexcitation fragments"
>> 523 << G4endl;
>> 524 G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
>> 525 << " Fragments E = " << ProductsEnergy/MeV << " MeV Diference --> "
>> 526 << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl;
>> 527 }
>> 528 if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV ||
>> 529 std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV ||
>> 530 std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) {
>> 531 G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
>> 532 G4cout << "G4ExcitationHandler.cc: Momentum Conservation test for deexcitation fragments"
>> 533 << G4endl;
>> 534 G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
>> 535 << " Fragments P = " << ProductsMomentum << " MeV Diference --> "
>> 536 << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl;
>> 537 }
>> 538 return;
674 } 539 }
>> 540 #endif
>> 541
675 542
676 543
677 544
678 545