Context Navigation

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

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

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

Download in other formats: