Context Navigation

source: trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4Evaporation.cc @ 1337

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago
tag geant4.9.4 beta 1 + modifs locales
File size: 17.3 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	// $Id: G4Evaporation.cc,v 1.25 2010/06/09 11:56:47 vnivanch Exp $
28	// GEANT4 tag $Name: geant4-09-04-beta-01 $
29	//
30	// Hadronic Process: Nuclear De-excitations
31	// by V. Lara (Oct 1998)
32	//
33	// Alex Howard - added protection for negative probabilities in the sum, 14/2/07
34	//
35	// Modif (03 September 2008) by J. M. Quesada for external choice of inverse
36	// cross section option
37	// JMQ (06 September 2008) Also external choices have been added for
38	// superimposed Coulomb barrier (if useSICBis set true, by default is false)
39	//
40	// V.Ivanchenko (27 July 2009) added G4EvaporationDefaultGEMFactory option
41	// V.Ivanchenko (10 May 2010) rewrited BreakItUp method: do not make new/delete
42	// photon channel first, fission second,
43	// added G4UnstableFragmentBreakUp to decay
44	// unphysical fragments (like 2n or 2p),
45	// use Z and A integer
46
47	#include "G4Evaporation.hh"
48	#include "G4EvaporationFactory.hh"
49	#include "G4EvaporationGEMFactory.hh"
50	#include "G4EvaporationDefaultGEMFactory.hh"
51	#include "G4HadronicException.hh"
52	#include "G4NistManager.hh"
53
54	G4Evaporation::G4Evaporation()
55	{
56	//theChannelFactory = new G4EvaporationFactory();
57	theChannelFactory = new G4EvaporationDefaultGEMFactory();
58	InitialiseEvaporation();
59	}
60
61	G4Evaporation::G4Evaporation(std::vector<G4VEvaporationChannel> aChannelsVector)
62	: theChannels(aChannelsVector), theChannelFactory(0), nChannels(0)
63	{
64	InitialiseEvaporation();
65	}
66
67	G4Evaporation::~G4Evaporation()
68	{
69	if (theChannels != 0) { theChannels = 0; }
70	if (theChannelFactory != 0) { delete theChannelFactory; }
71	}
72
73	void G4Evaporation::InitialiseEvaporation()
74	{
75	nist = G4NistManager::Instance();
76	minExcitation = CLHEP::keV;
77	if(theChannelFactory) { theChannels = theChannelFactory->GetChannel(); }
78	nChannels = theChannels->size();
79	probabilities.resize(nChannels, 0.0);
80	Initialise();
81	}
82
83	void G4Evaporation::Initialise()
84	{
85	// loop over evaporation channels
86	std::vector<G4VEvaporationChannel*>::iterator i;
87	for (i=theChannels->begin(); i != theChannels->end(); ++i)
88	{
89	// for inverse cross section choice
90	(*i)->SetOPTxs(OPTxs);
91	// for superimposed Coulomb Barrier for inverse cross sections
92	(*i)->UseSICB(useSICB);
93	}
94	}
95
96	void G4Evaporation::SetDefaultChannel()
97	{
98	if (theChannelFactory != 0) delete theChannelFactory;
99	theChannelFactory = new G4EvaporationFactory();
100	InitialiseEvaporation();
101	}
102
103	void G4Evaporation::SetGEMChannel()
104	{
105	if (theChannelFactory != 0) delete theChannelFactory;
106	theChannelFactory = new G4EvaporationGEMFactory();
107	InitialiseEvaporation();
108	}
109
110	void G4Evaporation::SetCombinedChannel()
111	{
112	if (theChannelFactory != 0) delete theChannelFactory;
113	theChannelFactory = new G4EvaporationDefaultGEMFactory();
114	InitialiseEvaporation();
115	}
116
117	G4FragmentVector * G4Evaporation::BreakItUp(const G4Fragment &theNucleus)
118	{
119	G4FragmentVector * theResult = new G4FragmentVector;
120	G4FragmentVector * theTempResult;
121
122	// The residual nucleus (after evaporation of each fragment)
123	G4Fragment* theResidualNucleus = new G4Fragment(theNucleus);
124
125	G4double totprob, prob, oldprob = 0.0;
126	G4int maxchannel, i;
127
128	G4int Amax = theResidualNucleus->GetA_asInt();
129
130	// Starts loop over evaporated particles, loop is limited by number
131	// of nucleons
132	for(G4int ia=0; ia<Amax; ++ia) {
133
134	// g,n,p - evaporation is finished
135	G4int A = theResidualNucleus->GetA_asInt();
136	if(1 >= A) {
137	theResult->push_back(theResidualNucleus);
138	return theResult;
139	}
140
141	// check if it is stable, then finish evaporation
142	G4int Z = theResidualNucleus->GetZ_asInt();
143	G4double abun = nist->GetIsotopeAbundance(Z, A);
144	/*
145	G4cout << "### G4Evaporation::BreakItUp step " << ia << " Z= " << Z
146	<< " A= " << A << " Eex(MeV)= "
147	<< theResidualNucleus->GetExcitationEnergy()
148	<< " aban= " << abun << G4endl;
149	*/
150	if(theResidualNucleus->GetExcitationEnergy() <= minExcitation &&
151	(abun > 0.0))
152	{
153	theResult->push_back(theResidualNucleus);
154	return theResult;
155	}
156
157	totprob = 0.0;
158	maxchannel = nChannels;
159
160	//G4cout << "### Evaporation loop #" << ia
161	// << " Fragment: " << theResidualNucleus << G4endl;
162
163	// loop over evaporation channels
164	for(i=0; i<nChannels; ++i) {
165	(theChannels)[i]->Initialize(theResidualNucleus);
166	prob = (*theChannels)[i]->GetEmissionProbability();
167	//G4cout << " Channel# " << i << " prob= " << prob << G4endl;
168
169	//if(0 == i && 0.0 == abun) { prob = 0.0; }
170
171	totprob += prob;
172	probabilities[i] = totprob;
173	// if two recent probabilities are near zero stop computations
174	if(i>=8) {
175	if(prob <= totprob1.e-8 && oldprob <= totprob1.e-8) {
176	maxchannel = i+1;
177	break;
178	}
179	}
180	oldprob = prob;
181	}
182
183	// photon evaporation in the case of no other channels available
184	// do evaporation chain and reset total probability
185	if(0.0 < totprob && probabilities[0] == totprob) {
186	theTempResult = (*theChannels)[0]->BreakUpFragment(theResidualNucleus);
187	if(theTempResult) {
188	size_t nsec = theTempResult->size();
189	for(size_t j=0; j<nsec; ++j) {
190	theResult->push_back((*theTempResult)[j]);
191	}
192	delete theTempResult;
193	}
194	totprob = 0.0;
195	}
196
197	// stable fragnent - evaporation is finished
198	if(0.0 == totprob) {
199
200	// if fragment is exotic, then try to decay it
201	if(0.0 == abun && Z < 20) {
202	//G4cout << "$$$ Decay exotic fragment" << G4endl;
203	theTempResult = unstableBreakUp.BreakUpFragment(theResidualNucleus);
204	if(theTempResult) {
205	size_t nsec = theTempResult->size();
206	for(size_t j=0; j<nsec; ++j) {
207	theResult->push_back((*theTempResult)[j]);
208	}
209	delete theTempResult;
210	}
211	}
212
213	// save residual fragment
214	theResult->push_back(theResidualNucleus);
215	return theResult;
216	}
217
218
219	// select channel
220	totprob *= G4UniformRand();
221	// loop over evaporation channels
222	for(i=0; i<maxchannel; ++i) { if(probabilities[i] >= totprob) { break; } }
223
224	// this should not happen
225	if(i >= nChannels) { i = nChannels - 1; }
226
227
228	// single photon evaporation, primary pointer is kept
229	if(0 == i) {
230	G4Fragment* gamma = (*theChannels)[0]->EmittedFragment(theResidualNucleus);
231	if(gamma) { theResult->push_back(gamma); }
232
233	// fission, return results to the main loop if fission is succesful
234	} else if(1 == i) {
235	theTempResult = (theChannels)[1]->BreakUp(theResidualNucleus);
236	if(theTempResult) {
237	size_t nsec = theTempResult->size();
238	G4bool deletePrimary = true;
239	for(size_t j=0; j<nsec; ++j) {
240	if(theResidualNucleus == (*theTempResult)[j]) { deletePrimary = false; }
241	theResult->push_back((*theTempResult)[j]);
242	}
243	if(deletePrimary) { delete theResidualNucleus; }
244	delete theTempResult;
245	return theResult;
246	}
247
248	// other channels
249	} else {
250	theTempResult = (theChannels)[i]->BreakUp(theResidualNucleus);
251	if(theTempResult) {
252	size_t nsec = theTempResult->size();
253	if(nsec > 0) {
254	--nsec;
255	for(size_t j=0; j<nsec; ++j) {
256	theResult->push_back((*theTempResult)[j]);
257	}
258	// if the residual change its pointer
259	// then delete previous residual fragment and update to the new
260	if(theResidualNucleus != (*theTempResult)[nsec] ) {
261	delete theResidualNucleus;
262	theResidualNucleus = (*theTempResult)[nsec];
263	}
264	}
265	delete theTempResult;
266	}
267	}
268	}
269
270	// loop is stopped, save residual
271	theResult->push_back(theResidualNucleus);
272
273	#ifdef debug
274	G4cout << "======== Evaporation Conservation Test ===========\n"
275	<< "==================================================\n";
276	CheckConservation(theNucleus,theResult);
277	G4cout << "==================================================\n";
278	#endif
279	return theResult;
280	}
281
282	/*
283	G4FragmentVector * G4Evaporation::BreakItUp(const G4Fragment &theNucleus)
284	{
285	G4FragmentVector * theResult = new G4FragmentVector;
286
287	// CHECK that Excitation Energy != 0
288	if (theNucleus.GetExcitationEnergy() <= 0.0) {
289	theResult->push_back(new G4Fragment(theNucleus));
290	return theResult;
291	}
292
293	// The residual nucleus (after evaporation of each fragment)
294	G4Fragment theResidualNucleus = theNucleus;
295
296	// Number of channels
297	G4int TotNumberOfChannels = theChannels->size();
298
299	// Starts loop over evaporated particles
300	for (;;)
301
302	{
303	// loop over evaporation channels
304	std::vector<G4VEvaporationChannel*>::iterator i;
305	for (i=theChannels->begin(); i != theChannels->end(); i++)
306	{
307	// for inverse cross section choice
308	(*i)->SetOPTxs(OPTxs);
309	// for superimposed Coulomb Barrier for inverse cross sections
310	(*i)->UseSICB(useSICB);
311
312	(*i)->Initialize(theResidualNucleus);
313	}
314	// Can't use this form beacuse Initialize is a non const member function
315	// for_each(theChannels->begin(),theChannels->end(),
316	// bind2nd(mem_fun(&G4VEvaporationChannel::Initialize),theResidualNucleus));
317	// Work out total decay probability by summing over channels
318	G4double TotalProbability = std::accumulate(theChannels->begin(),
319	theChannels->end(),
320	0.0,SumProbabilities());
321
322	if (TotalProbability <= 0.0)
323	{
324	// Will be no evaporation more
325	// write information about residual nucleus
326	theResult->push_back(new G4Fragment(theResidualNucleus));
327	break;
328	}
329	else
330	{
331	// Selection of evaporation channel, fission or gamma
332	// G4double * EmissionProbChannel = new G4double(TotNumberOfChannels);
333	std::vector<G4double> EmissionProbChannel;
334	EmissionProbChannel.reserve(theChannels->size());
335
336
337	// EmissionProbChannel[0] = theChannels->at(0)->GetEmissionProbability();
338
339
340	G4double first = theChannels->front()->GetEmissionProbability();
341
342	EmissionProbChannel.push_back(first >0 ? first : 0); // index 0
343
344
345	// EmissionProbChannel.push_back(theChannels->front()->GetEmissionProbability()); // index 0
346
347	for (i= (theChannels->begin()+1); i != theChannels->end(); i++)
348	{
349	// EmissionProbChannel[i] = EmissionProbChannel[i-1] +
350	// theChannels->at(i)->GetEmissionProbability();
351	// EmissionProbChannel.push_back(EmissionProbChannel.back() + (*i)->GetEmissionProbability());
352	first = (*i)->GetEmissionProbability();
353	EmissionProbChannel.push_back(first> 0? EmissionProbChannel.back() + first : EmissionProbChannel.back());
354	}
355
356	G4double shoot = G4UniformRand() * TotalProbability;
357	G4int j;
358	for (j=0; j < TotNumberOfChannels; j++)
359	{
360	// if (shoot < EmissionProbChannel[i])
361	if (shoot < EmissionProbChannel[j])
362	break;
363	}
364
365	// delete [] EmissionProbChannel;
366	EmissionProbChannel.clear();
367
368	if( j >= TotNumberOfChannels )
369	{
370	G4cerr << " Residual A: " << theResidualNucleus.GetA() << " Residual Z: " << theResidualNucleus.GetZ() << " Excitation Energy: " << theResidualNucleus.GetExcitationEnergy() << G4endl;
371	G4cerr << " j has not chosen a channel, j = " << j << " TotNumberOfChannels " << TotNumberOfChannels << " Total Probability: " << TotalProbability << G4endl;
372	for (j=0; j < TotNumberOfChannels; j++)
373	{
374	G4cerr << " j: " << j << " EmissionProbChannel: " << EmissionProbChannel[j] << " and shoot: " << shoot << " (<ProbChannel?) " << G4endl;
375	}
376	throw G4HadronicException(__FILE__, __LINE__, "G4Evaporation::BreakItUp: Can't define emission probability of the channels!" );
377	}
378	else
379	{
380	// Perform break-up
381	G4FragmentVector * theEvaporationResult = (*theChannels)[j]->BreakUp(theResidualNucleus);
382
383	#ifdef debug
384	G4cout << "-----------------------------------------------------------\n";
385	G4cout << G4endl << " After the evaporation of a particle, testing conservation \n";
386	CheckConservation(theResidualNucleus,theEvaporationResult);
387	G4cout << G4endl
388	<< "------------------------------------------------------------\n";
389	#endif
390
391	// Check if chosen channel is fission (there are only two EXCITED fragments)
392	// or the channel could not evaporate anything
393	if ( theEvaporationResult->size() == 1 \|\|
394	((*(theEvaporationResult->begin()))->GetExcitationEnergy() > 0.0 &&
395	(*(theEvaporationResult->end()-1))->GetExcitationEnergy() > 0.0) ) {
396	// FISSION
397	for (G4FragmentVector::iterator i = theEvaporationResult->begin();
398	i != theEvaporationResult->end(); ++i)
399	{
400	theResult->push_back(*(i));
401	}
402	delete theEvaporationResult;
403	break;
404	} else {
405	// EVAPORATION
406	for (G4FragmentVector::iterator i = theEvaporationResult->begin();
407	i != theEvaporationResult->end()-1; ++i)
408	{
409	#ifdef PRECOMPOUND_TEST
410	if ((i)->GetA() == 0) (i)->SetCreatorModel(G4String("G4PhotonEvaporation"));
411	#endif
412	theResult->push_back(*(i));
413	}
414	theResidualNucleus = *(theEvaporationResult->back());
415	delete theEvaporationResult->back();
416	delete theEvaporationResult;
417	#ifdef PRECOMPOUND_TEST
418	theResidualNucleus.SetCreatorModel(G4String("ResidualNucleus"));
419	#endif
420
421	}
422	}
423	}
424	}
425
426	#ifdef debug
427	G4cout << "======== Evaporation Conservation Test ===========\n"
428	<< "==================================================\n";
429	CheckConservation(theNucleus,theResult);
430	G4cout << "==================================================\n";
431	#endif
432	return theResult;
433	}
434	*/
435
436
437	#ifdef debug
438	void G4Evaporation::CheckConservation(const G4Fragment & theInitialState,
439	G4FragmentVector * Result) const
440	{
441	G4double ProductsEnergy =0;
442	G4ThreeVector ProductsMomentum;
443	G4int ProductsA = 0;
444	G4int ProductsZ = 0;
445	for (G4FragmentVector::iterator h = Result->begin(); h != Result->end(); h++) {
446	G4LorentzVector tmp = (*h)->GetMomentum();
447	ProductsEnergy += tmp.e();
448	ProductsMomentum += tmp.vect();
449	ProductsA += static_cast<G4int>((*h)->GetA());
450	ProductsZ += static_cast<G4int>((*h)->GetZ());
451	}
452
453	if (ProductsA != theInitialState.GetA()) {
454	G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl;
455	G4cout << "G4Evaporation.cc: Barionic Number Conservation test for evaporation fragments"
456	<< G4endl;
457	G4cout << "Initial A = " << theInitialState.GetA()
458	<< " Fragments A = " << ProductsA << " Diference --> "
459	<< theInitialState.GetA() - ProductsA << G4endl;
460	}
461	if (ProductsZ != theInitialState.GetZ()) {
462	G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
463	G4cout << "G4Evaporation.cc: Charge Conservation test for evaporation fragments"
464	<< G4endl;
465	G4cout << "Initial Z = " << theInitialState.GetZ()
466	<< " Fragments Z = " << ProductsZ << " Diference --> "
467	<< theInitialState.GetZ() - ProductsZ << G4endl;
468	}
469	if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) {
470	G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
471	G4cout << "G4Evaporation.cc: Energy Conservation test for evaporation fragments"
472	<< G4endl;
473	G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
474	<< " Fragments E = " << ProductsEnergy/MeV << " MeV Diference --> "
475	<< (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl;
476	}
477	if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV \|\|
478	std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV \|\|
479	std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) {
480	G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
481	G4cout << "G4Evaporation.cc: Momentum Conservation test for evaporation fragments"
482	<< G4endl;
483	G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
484	<< " Fragments P = " << ProductsMomentum << " MeV Diference --> "
485	<< theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl;
486	}
487	return;
488	}
489	#endif
490
491
492
493

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

Download in other formats: