Context Navigation

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

Last change on this file since 1007 was 962, checked in by garnier, 15 years ago
update processes
File size: 20.0 KB

Line
1	//
2	// ********************************************************************
3	// * License and Disclaimer *
4	// * *
5	// * The Geant4 software is copyright of the Copyright Holders of *
6	// * the Geant4 Collaboration. It is provided under the terms and *
7	// * conditions of the Geant4 Software License, included in the file *
8	// * LICENSE and available at http://cern.ch/geant4/license . These *
9	// * include a list of copyright holders. *
10	// * *
11	// * Neither the authors of this software system, nor their employing *
12	// * institutes,nor the agencies providing financial support for this *
13	// * work make any representation or warranty, express or implied, *
14	// * regarding this software system or assume any liability for its *
15	// * use. Please see the license in the file LICENSE and URL above *
16	// * for the full disclaimer and the limitation of liability. *
17	// * *
18	// * This code implementation is the result of the scientific and *
19	// * technical work of the GEANT4 collaboration. *
20	// * By using, copying, modifying or distributing the software (or *
21	// * any work based on the software) you agree to acknowledge its *
22	// * use in resulting scientific publications, and indicate your *
23	// * acceptance of all terms of the Geant4 Software license. *
24	// ********************************************************************
25	//
26	//
27	// Hadronic Process: Nuclear De-excitations
28	// by V. Lara (May 1998)
29	//
30	// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
31	// cross section option
32	//
33	// Modif (24 Jul 2008) by M. A. Cortes Giraldo:
34	// -Max Z,A for Fermi Break-Up turns to 9,17 by default
35	// -BreakItUp() reorganised and bug in Evaporation loop fixed
36	// -Transform() optimised
37	// Modif (30 June 1998) by V. Lara:
38	// -Modified the Transform method for use G4ParticleTable and
39	// therefore G4IonTable. It makes possible to convert all kind
40	// of fragments (G4Fragment) produced in deexcitation to
41	// G4DynamicParticle
42	// -It uses default algorithms for:
43	// Evaporation: G4Evaporation
44	// MultiFragmentation: G4StatMF
45	// Fermi Breakup model: G4FermiBreakUp
46	//
47	// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
48	// cross section option (default OPTxs=3)
49	// JMQ (06 September 2008) Also external choices have been added for
50	// superimposed Coulomb barrier (if useSICBis set true, by default is false)
51
52	#include "G4ExcitationHandler.hh"
53	#include <list>
54
55	//#define debugphoton
56
57
58	G4ExcitationHandler::G4ExcitationHandler():
59	// Fermi BreakUp is on and MultiFrag is off by default
60	// maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),minEForMultiFrag(4.0*GeV),
61	maxZForFermiBreakUp(1),maxAForFermiBreakUp(1),minEForMultiFrag(4.0*GeV),
62	MyOwnEvaporationClass(true), MyOwnMultiFragmentationClass(true),MyOwnFermiBreakUpClass(true),
63	MyOwnPhotonEvaporationClass(true),OPTxs(3),useSICB(false)
64	{
65	theTableOfParticles = G4ParticleTable::GetParticleTable();
66
67	theEvaporation = new G4Evaporation;
68	theMultiFragmentation = new G4StatMF;
69	theFermiModel = new G4FermiBreakUp;
70	thePhotonEvaporation = new G4PhotonEvaporation;
71	}
72
73	G4ExcitationHandler::G4ExcitationHandler(const G4ExcitationHandler &)
74	{
75	throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::copy_constructor: is meant to not be accessable! ");
76	}
77
78
79	G4ExcitationHandler::~G4ExcitationHandler()
80	{
81	if (MyOwnEvaporationClass) delete theEvaporation;
82	if (MyOwnMultiFragmentationClass) delete theMultiFragmentation;
83	if (MyOwnFermiBreakUpClass) delete theFermiModel;
84	if (MyOwnPhotonEvaporationClass) delete thePhotonEvaporation;
85	}
86
87
88	const G4ExcitationHandler & G4ExcitationHandler::operator=(const G4ExcitationHandler &)
89	{
90	throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator=: is meant to not be accessable! ");
91
92	return *this;
93	}
94
95
96	G4bool G4ExcitationHandler::operator==(const G4ExcitationHandler &) const
97	{
98	throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator==: is meant to not be accessable! ");
99	return false;
100	}
101
102	G4bool G4ExcitationHandler::operator!=(const G4ExcitationHandler &) const
103	{
104	throw G4HadronicException(__FILE__, __LINE__, "G4ExcitationHandler::operator!=: is meant to not be accessable! ");
105	return true;
106	}
107
108	////////////////////////////////////////////////////////////////////////////////////////////////
109	/// 25/07/08 16:45 Proposed by MAC ////////////////////////////////////////////////////////////
110	////////////////////////////////////////////////////////////////////////////////////////////////
111
112	G4ReactionProductVector * G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState) const
113	{
114	//for inverse cross section choice
115	theEvaporation->SetOPTxs(OPTxs);
116	//for the choice of superimposed Coulomb Barrier for inverse cross sections
117	theEvaporation->UseSICB(useSICB);
118
119	// Pointer which will be used to return the final production vector
120	G4FragmentVector * theResult = 0;
121
122	// Variables existing until end of method
123	G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
124	G4Fragment theExcitedNucleus; // object to be passed in BreakItUp methods
125	G4FragmentVector * theTempResult = 0; // pointer which receives temporal results
126	std::list<G4Fragment*> theEvapList; // list to apply Evaporation, SMF or Fermi Break-Up
127	std::list<G4Fragment*> theEvapStableList; // list to apply PhotonEvaporation
128	std::list<G4Fragment*> theFinalStableList; // list to store final result
129	std::list<G4Fragment*>::iterator iList;
130
131	// Variables to describe the excited configuration
132	G4double exEnergy = theInitialState.GetExcitationEnergy();
133	G4int A = static_cast<G4int>( theInitialState.GetA() +0.5 );
134	G4int Z = static_cast<G4int>( theInitialState.GetZ() +0.5 );
135
136	// In case A <= 4 the fragment will not perform any nucleon emission
137	if (A <= 4)
138	{
139	// I store G4Fragment* in theEvapStableList to apply thePhotonEvaporation later
140
141	theEvapStableList.push_back( theInitialStatePtr );
142	}
143
144	else // If A > 4 we try to apply theFermiModel, theMultiFragmentation or theEvaporation
145	{
146
147	// Store original state in theEvapList
148	theEvapList.push_back( theInitialStatePtr );
149
150	// ------------------------------
151	// De-excitation loop
152	// ------------------------------
153
154	for (iList = theEvapList.begin(); iList != theEvapList.end(); iList++)
155	{
156	exEnergy = (*iList)->GetExcitationEnergy();
157
158	if (exEnergy > 0.0)
159	{
160	// Check conditions for each model
161	theExcitedNucleus = (iList);
162	A = static_cast<G4int>((*iList)->GetA()+0.5); // +0.5 to avoid bad truncation
163	Z = static_cast<G4int>((*iList)->GetZ()+0.5);
164
165	if ( A < GetMaxA() && Z < GetMaxZ() ) // if satisfied apply Fermi Break-Up
166	{
167	theTempResult = theFermiModel->BreakItUp(theExcitedNucleus);
168	}
169	else if (exEnergy > GetMinE()*A) // if satisfied apply SMF
170	{
171	theTempResult = theMultiFragmentation->BreakItUp(theExcitedNucleus);
172	}
173	else // apply Evaporation in another case
174	{
175	theTempResult = theEvaporation->BreakItUp(theExcitedNucleus);
176	}
177
178	// New configuration is stored in theTempResult, so we can free
179	// the memory where the previous configuration is
180
181	delete (*iList);
182
183	// And now the theTempResult->size() tells us if the configuration has changed
184
185	if ( theTempResult->size() > 1 )
186	{
187	// push_back the result to the end of theEvapList (this same list)
188
189	for (G4FragmentVector::iterator j = theTempResult->begin();
190	j != theTempResult->end(); j++)
191	{
192	theEvapList.push_back(*j);
193	}
194	}
195	else
196	{
197	// push_back the result to theEvapStableList, because
198	// is still excited, but cannot emmit more nucleons
199
200	for (G4FragmentVector::iterator j = theTempResult->begin();
201	j != theTempResult->end(); j++)
202	{
203	theEvapStableList.push_back(*j);
204	}
205	}
206
207	// after working with theTempResult, clear and delete it
208	theTempResult->clear();
209	delete theTempResult;
210
211	}
212	else // exEnergy = 0.0
213	{
214	// if this fragment is at ground state,
215	// store it in theFinalStableList
216
217	theFinalStableList.push_back(*iList);
218
219	} // endif (exEnergy > 0.0)
220
221	} // end of the loop over theEvapList
222
223	theEvapList.clear(); // clear all the list and do not free memory pointed by
224	// each element because this have been done before!
225	}
226
227	// Now we try to deexcite by means of PhotonEvaporation those fragments
228	// which are still excited.
229
230
231	// -----------------------
232	// Photon-Evaporation loop
233	// -----------------------
234
235	for (iList = theEvapStableList.begin(); iList != theEvapStableList.end(); iList++)
236	{
237	A = static_cast<G4int>((*iList)->GetA()+0.5);
238	exEnergy = (*iList)->GetExcitationEnergy();
239
240	if ( A > 1 && exEnergy > 0.1*eV ) // if so, photon-evaporation is applied
241	{
242	theExcitedNucleus = (iList);
243	theTempResult = thePhotonEvaporation->BreakItUp(theExcitedNucleus);
244
245	// if there is a gamma emission then
246	if (theTempResult->size() > 1)
247	{
248	// first free the memory occupied by the previous state
249	delete (*iList);
250
251	// and now add the final state from gamma emission to the end of
252	// theEvapStableList
253	for (G4FragmentVector::reverse_iterator ri = theTempResult->rbegin();
254	ri != theTempResult->rend(); ++ri)
255	// reversed is applied in order to have residual nucleus in first position
256	{
257
258	#ifdef PRECOMPOUND_TEST
259	if ((*ri)->GetA() == 0)
260	(*ri)->SetCreatorModel(G4String("G4PhotonEvaporation"));
261	else
262	(*ri)->SetCreatorModel(G4String("ResidualNucleus"));
263	#endif
264
265	theEvapStableList.push_back(*ri);
266	}
267
268	// now we clean and remove the temporal vector
269	theTempResult->clear();
270	delete theTempResult;
271
272	}
273	else // if theTempResult->size() = 1
274	{
275	// if there is not any gamma emission from this excited fragment
276	// we have to emmit a gamma which forces the deexcitation
277
278	// First I clean completely theTempResult
279
280	for (G4FragmentVector::iterator j = theTempResult->begin();
281	j != theTempResult->end(); j++)
282	{
283	delete (*j);
284	}
285
286	theTempResult->clear();
287	delete theTempResult;
288
289	#ifdef debugphoton
290	G4cout << "G4ExcitationHandler: Gamma Evaporation could not deexcite the nucleus: \n"
291	<< "-----------------------------------------------------------------------\n"
292	<< theExcitedNucleus << '\n'
293	<< "-----------------------------------------------------------------------\n";
294	#endif
295
296	// Let's create a G4Fragment pointer representing the gamma emmited
297	G4double GammaEnergy = (*iList)->GetExcitationEnergy();
298	G4double cosTheta = 1. - 2. * G4UniformRand();
299	G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta);
300	G4double phi = twopi * G4UniformRand();
301	G4ThreeVector GammaP(GammaEnergy * sinTheta * std::cos(phi),
302	GammaEnergy * sinTheta * std::sin(phi),
303	GammaEnergy * cosTheta );
304	G4LorentzVector Gamma4P(GammaP,GammaEnergy);
305	G4Fragment * theHandlerPhoton = new G4Fragment(Gamma4P,G4Gamma::GammaDefinition());
306
307	// And now we update momentum and energy for the nucleus
308	G4double Mass = (*iList)->GetGroundStateMass();
309	G4ThreeVector ResidualP((*iList)->GetMomentum().vect() - GammaP);
310	G4double ResidualE = std::sqrt(ResidualPResidualP + MassMass);
311	G4LorentzVector Residual4P(ResidualP,ResidualE);
312	(*iList)->SetMomentum(Residual4P); // Now this fragment has been deexcited!
313
314	// we store the deexcited fragment in theFinalStableList
315	theFinalStableList.push_back(*iList);
316
317	#ifdef PRECOMPOUND_TEST
318	theHandlerPhoton->SetCreatorModel("G4ExcitationHandler");
319	#endif
320
321	// Finally, we add theHandlerPhoton to theFinalStableList
322	theFinalStableList.push_back(theHandlerPhoton);
323
324	#ifdef debugphoton
325	G4cout << "Emmited photon:\n"
326	<< theFinalStableList.back() << '\n'
327	<< "Residual nucleus after photon emission:\n"
328	<< (iList) << '\n'
329	<< "-----------------------------------------------------------------------\n";
330	#endif
331
332	}
333	}
334	else // case of a nucleon, gamma or very small excitation energy
335	{
336	// we don't have to do anything, just store the fragment in theFinalStableList
337
338	theFinalStableList.push_back(*iList);
339
340	} // A > 1 && exEnergy > 0.1*eV
341
342	} // end of photon-evaporation loop
343
344
345	// The deexcitation from fragments inside theEvapStableList has been finished, so...
346	theEvapStableList.clear();
347
348	// Now the final state is in theFinalStableList, and we have to send it to theResult vector
349
350	theResult = new G4FragmentVector;
351	theResult->reserve( theFinalStableList.size() * sizeof(G4Fragment*) );
352	// We reserve enough memory to optimise the storing speed
353
354	for (iList = theFinalStableList.begin(); iList != theFinalStableList.end(); iList++)
355	{
356	theResult->push_back(*iList);
357	}
358
359	// After storing the final state , we can clear theFinalStableList
360	theFinalStableList.clear();
361
362	#ifdef debug
363	CheckConservation(theInitialState,theResult);
364	#endif
365
366	// Change G4FragmentVector* to G4ReactionProductVector*
367	return Transform(theResult);
368	}
369
370
371
372	G4ReactionProductVector *
373	G4ExcitationHandler::Transform(G4FragmentVector * theFragmentVector) const
374	{
375	if (theFragmentVector == 0) return 0;
376
377	// Conversion from G4FragmentVector to G4ReactionProductVector
378	G4ParticleDefinition *theGamma = G4Gamma::GammaDefinition();
379	G4ParticleDefinition *theNeutron = G4Neutron::NeutronDefinition();
380	G4ParticleDefinition *theProton = G4Proton::ProtonDefinition();
381	G4ParticleDefinition *theDeuteron = G4Deuteron::DeuteronDefinition();
382	G4ParticleDefinition *theTriton = G4Triton::TritonDefinition();
383	G4ParticleDefinition *theHelium3 = G4He3::He3Definition();
384	G4ParticleDefinition *theAlpha = G4Alpha::AlphaDefinition();
385	G4ParticleDefinition *theKindOfFragment = 0;
386	theNeutron->SetVerboseLevel(2);
387	G4ReactionProductVector * theReactionProductVector = new G4ReactionProductVector;
388
389	// MAC (24/07/08)
390	// To optimise the storing speed, we reserve space in memory for the vector
391	theReactionProductVector->reserve( theFragmentVector->size() * sizeof(G4ReactionProduct*) );
392
393	G4int theFragmentA, theFragmentZ;
394	G4LorentzVector theFragmentMomentum;
395
396	G4FragmentVector::iterator i;
397	for (i = theFragmentVector->begin(); i != theFragmentVector->end(); i++) {
398	// std::cout << (*i) <<'\n';
399	theFragmentA = static_cast<G4int>((*i)->GetA());
400	theFragmentZ = static_cast<G4int>((*i)->GetZ());
401	theFragmentMomentum = (*i)->GetMomentum();
402	theKindOfFragment = 0;
403	if (theFragmentA == 0 && theFragmentZ == 0) { // photon
404	theKindOfFragment = theGamma;
405	} else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
406	theKindOfFragment = theNeutron;
407	} else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
408	theKindOfFragment = theProton;
409	} else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
410	theKindOfFragment = theDeuteron;
411	} else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
412	theKindOfFragment = theTriton;
413	} else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
414	theKindOfFragment = theHelium3;
415	} else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha
416	theKindOfFragment = theAlpha;
417	} else {
418	theKindOfFragment = theTableOfParticles->FindIon(theFragmentZ,theFragmentA,0,theFragmentZ);
419	}
420	if (theKindOfFragment != 0)
421	{
422	G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
423	theNew->SetMomentum(theFragmentMomentum.vect());
424	theNew->SetTotalEnergy(theFragmentMomentum.e());
425	theNew->SetFormationTime((*i)->GetCreationTime());
426	#ifdef PRECOMPOUND_TEST
427	theNew->SetCreatorModel((*i)->GetCreatorModel());
428	#endif
429	theReactionProductVector->push_back(theNew);
430	}
431	}
432	if (theFragmentVector != 0)
433	{
434	std::for_each(theFragmentVector->begin(), theFragmentVector->end(), DeleteFragment());
435	delete theFragmentVector;
436	}
437	G4ReactionProductVector::iterator debugit;
438	for(debugit=theReactionProductVector->begin();
439	debugit!=theReactionProductVector->end(); debugit++)
440	{
441	if((debugit)->GetTotalEnergy()<1.eV)
442	{
443	if(getenv("G4DebugPhotonevaporationData"))
444	{
445	G4cerr << "G4ExcitationHandler: Warning: Photonevaporation data not exact."<<G4endl;
446	G4cerr << "G4ExcitationHandler: Warning: Found gamma with energy = "
447	<< (*debugit)->GetTotalEnergy()/MeV << "MeV"
448	<< G4endl;
449	}
450	delete (*debugit);
451	*debugit = 0;
452	}
453	}
454	G4ReactionProduct* tmpPtr=0;
455	theReactionProductVector->erase(std::remove_if(theReactionProductVector->begin(),
456	theReactionProductVector->end(),
457	std::bind2nd(std::equal_to<G4ReactionProduct*>(),
458	tmpPtr)),
459	theReactionProductVector->end());
460	return theReactionProductVector;
461	}
462
463
464	#ifdef debug
465	void G4ExcitationHandler::CheckConservation(const G4Fragment & theInitialState,
466	G4FragmentVector * Result) const
467	{
468	G4double ProductsEnergy =0;
469	G4ThreeVector ProductsMomentum;
470	G4int ProductsA = 0;
471	G4int ProductsZ = 0;
472	G4FragmentVector::iterator h;
473	for (h = Result->begin(); h != Result->end(); h++) {
474	G4LorentzVector tmp = (*h)->GetMomentum();
475	ProductsEnergy += tmp.e();
476	ProductsMomentum += tmp.vect();
477	ProductsA += static_cast<G4int>((*h)->GetA());
478	ProductsZ += static_cast<G4int>((*h)->GetZ());
479	}
480
481	if (ProductsA != theInitialState.GetA()) {
482	G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl;
483	G4cout << "G4ExcitationHandler.cc: Barionic Number Conservation test for deexcitation fragments"
484	<< G4endl;
485	G4cout << "Initial A = " << theInitialState.GetA()
486	<< " Fragments A = " << ProductsA << " Diference --> "
487	<< theInitialState.GetA() - ProductsA << G4endl;
488	}
489	if (ProductsZ != theInitialState.GetZ()) {
490	G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
491	G4cout << "G4ExcitationHandler.cc: Charge Conservation test for deexcitation fragments"
492	<< G4endl;
493	G4cout << "Initial Z = " << theInitialState.GetZ()
494	<< " Fragments Z = " << ProductsZ << " Diference --> "
495	<< theInitialState.GetZ() - ProductsZ << G4endl;
496	}
497	if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) {
498	G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
499	G4cout << "G4ExcitationHandler.cc: Energy Conservation test for deexcitation fragments"
500	<< G4endl;
501	G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
502	<< " Fragments E = " << ProductsEnergy/MeV << " MeV Diference --> "
503	<< (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl;
504	}
505	if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV \|\|
506	std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV \|\|
507	std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) {
508	G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
509	G4cout << "G4ExcitationHandler.cc: Momentum Conservation test for deexcitation fragments"
510	<< G4endl;
511	G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
512	<< " Fragments P = " << ProductsMomentum << " MeV Diference --> "
513	<< theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl;
514	}
515	return;
516	}
517	#endif
518
519
520
521

Note: See TracBrowser for help on using the repository browser.

Download in other formats: