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

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

geant4 tag 9.4
File size: 17.5 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.26 2010/11/23 18:10:10 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
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
54G4Evaporation::G4Evaporation() 
55{
56  //theChannelFactory = new G4EvaporationFactory();
57  theChannelFactory = new G4EvaporationDefaultGEMFactory();
58  InitialiseEvaporation();
59}
60
61G4Evaporation::G4Evaporation(std::vector<G4VEvaporationChannel*> * aChannelsVector) 
62  : theChannels(aChannelsVector), theChannelFactory(0), nChannels(0)
63{
64  InitialiseEvaporation();
65}
66
67G4Evaporation::~G4Evaporation()
68{
69  if (theChannels != 0) { theChannels = 0; }
70  if (theChannelFactory != 0) { delete theChannelFactory; }
71}
72
73void 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
83void 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
96void G4Evaporation::SetDefaultChannel()
97{
98  if (theChannelFactory != 0) delete theChannelFactory;
99  theChannelFactory = new G4EvaporationFactory();
100  InitialiseEvaporation();
101}
102
103void G4Evaporation::SetGEMChannel()
104{
105  if (theChannelFactory != 0) delete theChannelFactory;
106  theChannelFactory = new G4EvaporationGEMFactory();
107  InitialiseEvaporation();
108}
109
110void G4Evaporation::SetCombinedChannel()
111{
112  if (theChannelFactory != 0) delete theChannelFactory;
113  theChannelFactory = new G4EvaporationDefaultGEMFactory();
114  InitialiseEvaporation();
115}
116
117G4FragmentVector * 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 <= totprob*1.e-8 && oldprob <= totprob*1.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      //G4cout << "Start gamma evaporation" << G4endl;
187      theTempResult = (*theChannels)[0]->BreakUpFragment(theResidualNucleus);
188      if(theTempResult) {
189        size_t nsec = theTempResult->size();
190        for(size_t j=0; j<nsec; ++j) {
191          theResult->push_back((*theTempResult)[j]);
192        }
193        delete theTempResult;
194      }
195      totprob = 0.0;
196    }
197
198    // stable fragnent - evaporation is finished
199    if(0.0 == totprob) {
200
201      // if fragment is exotic, then try to decay it
202      if(0.0 == abun && Z < 20) {
203        //G4cout << "$$$ Decay exotic fragment" << G4endl;
204        theTempResult = unstableBreakUp.BreakUpFragment(theResidualNucleus);
205        if(theTempResult) {
206          size_t nsec = theTempResult->size();
207          for(size_t j=0; j<nsec; ++j) {
208            theResult->push_back((*theTempResult)[j]);
209          }
210          delete theTempResult;
211        }
212      }
213
214      // save residual fragment
215      theResult->push_back(theResidualNucleus);
216      return theResult;
217    }
218
219
220    // select channel
221    totprob *= G4UniformRand();
222    // loop over evaporation channels
223    for(i=0; i<maxchannel; ++i) { if(probabilities[i] >= totprob) { break; } }
224
225    // this should not happen
226    if(i >= nChannels) { i = nChannels - 1; }
227
228
229    // single photon evaporation, primary pointer is kept
230    if(0 == i) {
231      //G4cout << "Single gamma" << G4endl;
232      G4Fragment* gamma = (*theChannels)[0]->EmittedFragment(theResidualNucleus);
233      if(gamma) { theResult->push_back(gamma); }
234
235      // fission, return results to the main loop if fission is succesful
236    } else if(1 == i) {
237      //G4cout << "Fission" << G4endl;
238      theTempResult = (*theChannels)[1]->BreakUp(*theResidualNucleus);
239      if(theTempResult) {
240        size_t nsec = theTempResult->size();
241        G4bool deletePrimary = true;
242        for(size_t j=0; j<nsec; ++j) {
243          if(theResidualNucleus == (*theTempResult)[j]) { deletePrimary = false; }
244          theResult->push_back((*theTempResult)[j]);
245        }
246        if(deletePrimary) { delete theResidualNucleus; }
247        delete theTempResult;
248        return theResult;
249      }
250
251      // other channels
252    } else {
253      //G4cout << "Channel # " << i << G4endl;
254      theTempResult = (*theChannels)[i]->BreakUp(*theResidualNucleus);
255      if(theTempResult) {
256        size_t nsec = theTempResult->size();
257        if(nsec > 0) {
258          --nsec;
259          for(size_t j=0; j<nsec; ++j) {
260            theResult->push_back((*theTempResult)[j]);
261          }
262          // if the residual change its pointer
263          // then delete previous residual fragment and update to the new
264          if(theResidualNucleus != (*theTempResult)[nsec] ) { 
265            delete theResidualNucleus; 
266            theResidualNucleus = (*theTempResult)[nsec];
267          }
268        }
269        delete theTempResult;
270      }
271    }
272  }
273
274  // loop is stopped, save residual
275  theResult->push_back(theResidualNucleus);
276 
277#ifdef debug
278  G4cout << "======== Evaporation Conservation Test ===========\n"
279         << "==================================================\n";
280  CheckConservation(theNucleus,theResult);
281  G4cout << "==================================================\n";
282#endif
283  return theResult;
284}
285
286/*
287G4FragmentVector * G4Evaporation::BreakItUp(const G4Fragment &theNucleus)
288{
289  G4FragmentVector * theResult = new G4FragmentVector;
290
291    // CHECK that Excitation Energy != 0
292    if (theNucleus.GetExcitationEnergy() <= 0.0) {
293        theResult->push_back(new G4Fragment(theNucleus));
294        return theResult;
295    }
296
297    // The residual nucleus (after evaporation of each fragment)
298    G4Fragment theResidualNucleus = theNucleus;
299
300    // Number of channels
301    G4int TotNumberOfChannels = theChannels->size(); 
302
303    // Starts loop over evaporated particles
304    for (;;)
305 
306   {   
307        // loop over evaporation channels
308        std::vector<G4VEvaporationChannel*>::iterator i;
309        for (i=theChannels->begin(); i != theChannels->end(); i++)
310          {
311  // for inverse cross section choice
312            (*i)->SetOPTxs(OPTxs);
313  // for superimposed Coulomb Barrier for inverse cross sections
314            (*i)->UseSICB(useSICB);
315
316            (*i)->Initialize(theResidualNucleus);
317          }
318        // Can't use this form beacuse Initialize is a non const member function
319        //      for_each(theChannels->begin(),theChannels->end(),
320        //               bind2nd(mem_fun(&G4VEvaporationChannel::Initialize),theResidualNucleus));
321        // Work out total decay probability by summing over channels
322        G4double TotalProbability =  std::accumulate(theChannels->begin(),
323                                                       theChannels->end(),
324                                                       0.0,SumProbabilities());
325       
326        if (TotalProbability <= 0.0)
327          {
328            // Will be no evaporation more
329            // write information about residual nucleus
330            theResult->push_back(new G4Fragment(theResidualNucleus));
331            break;
332          }
333        else
334          {
335            // Selection of evaporation channel, fission or gamma
336            // G4double * EmissionProbChannel = new G4double(TotNumberOfChannels);
337            std::vector<G4double> EmissionProbChannel;
338            EmissionProbChannel.reserve(theChannels->size());
339           
340
341            // EmissionProbChannel[0] = theChannels->at(0)->GetEmissionProbability();
342
343
344            G4double first = theChannels->front()->GetEmissionProbability();
345
346            EmissionProbChannel.push_back(first >0 ? first : 0); // index 0
347
348
349            //      EmissionProbChannel.push_back(theChannels->front()->GetEmissionProbability()); // index 0
350           
351            for (i= (theChannels->begin()+1); i != theChannels->end(); i++)
352              {
353                // EmissionProbChannel[i] = EmissionProbChannel[i-1] +
354                // theChannels->at(i)->GetEmissionProbability();
355                //              EmissionProbChannel.push_back(EmissionProbChannel.back() + (*i)->GetEmissionProbability());
356                first = (*i)->GetEmissionProbability();
357                EmissionProbChannel.push_back(first> 0? EmissionProbChannel.back() + first : EmissionProbChannel.back());
358              }
359
360            G4double shoot = G4UniformRand() * TotalProbability;
361            G4int j;
362            for (j=0; j < TotNumberOfChannels; j++)
363              {
364                // if (shoot < EmissionProbChannel[i])
365                if (shoot < EmissionProbChannel[j])
366                  break;
367              }
368           
369            // delete [] EmissionProbChannel;
370            EmissionProbChannel.clear();
371                       
372            if( j >= TotNumberOfChannels )
373              {
374                G4cerr << " Residual A: " << theResidualNucleus.GetA() << " Residual Z: " << theResidualNucleus.GetZ() << " Excitation Energy: " << theResidualNucleus.GetExcitationEnergy() << G4endl;
375                G4cerr << " j has not chosen a channel, j = " << j << " TotNumberOfChannels " << TotNumberOfChannels << " Total Probability: " << TotalProbability << G4endl;
376                for (j=0; j < TotNumberOfChannels; j++)
377                  {
378                    G4cerr << " j: " << j << " EmissionProbChannel: " << EmissionProbChannel[j] << " and shoot: " << shoot << " (<ProbChannel?) " << G4endl;
379                  }             
380                throw G4HadronicException(__FILE__, __LINE__,  "G4Evaporation::BreakItUp: Can't define emission probability of the channels!" );
381              }
382            else
383              {
384                // Perform break-up
385                G4FragmentVector * theEvaporationResult = (*theChannels)[j]->BreakUp(theResidualNucleus);
386
387#ifdef debug
388                G4cout <<           "-----------------------------------------------------------\n";
389                G4cout << G4endl << " After the evaporation of a particle, testing conservation \n";
390                CheckConservation(theResidualNucleus,theEvaporationResult);
391                G4cout << G4endl
392                       <<           "------------------------------------------------------------\n";
393#endif 
394
395                // Check if chosen channel is fission (there are only two EXCITED fragments)
396                // or the channel could not evaporate anything
397                if ( theEvaporationResult->size() == 1 ||
398                     ((*(theEvaporationResult->begin()))->GetExcitationEnergy() > 0.0 &&
399                      (*(theEvaporationResult->end()-1))->GetExcitationEnergy() > 0.0) ) {
400                  // FISSION
401                  for (G4FragmentVector::iterator i = theEvaporationResult->begin();
402                       i != theEvaporationResult->end(); ++i)
403                    {
404                      theResult->push_back(*(i));
405                    }
406                  delete theEvaporationResult;
407                  break;
408                } else {
409                  // EVAPORATION
410                  for (G4FragmentVector::iterator i = theEvaporationResult->begin();
411                       i != theEvaporationResult->end()-1; ++i)
412                    {
413#ifdef PRECOMPOUND_TEST
414                      if ((*i)->GetA() == 0) (*i)->SetCreatorModel(G4String("G4PhotonEvaporation"));
415#endif
416                      theResult->push_back(*(i));
417                    }
418                  theResidualNucleus = *(theEvaporationResult->back());
419                  delete theEvaporationResult->back();
420                  delete theEvaporationResult;
421#ifdef PRECOMPOUND_TEST
422                  theResidualNucleus.SetCreatorModel(G4String("ResidualNucleus"));
423#endif
424
425                }
426              }
427          }
428      }
429   
430#ifdef debug
431    G4cout << "======== Evaporation Conservation Test ===========\n"
432           << "==================================================\n";
433    CheckConservation(theNucleus,theResult);
434    G4cout << "==================================================\n";
435#endif
436    return theResult;
437}
438*/
439
440
441#ifdef debug
442void G4Evaporation::CheckConservation(const G4Fragment & theInitialState,
443                                      G4FragmentVector * Result) const
444{
445    G4double ProductsEnergy =0;
446    G4ThreeVector ProductsMomentum;
447    G4int ProductsA = 0;
448    G4int ProductsZ = 0;
449    for (G4FragmentVector::iterator h = Result->begin(); h != Result->end(); h++) {
450        G4LorentzVector tmp = (*h)->GetMomentum();
451        ProductsEnergy += tmp.e();
452        ProductsMomentum += tmp.vect();
453        ProductsA += static_cast<G4int>((*h)->GetA());
454        ProductsZ += static_cast<G4int>((*h)->GetZ());
455    }
456
457    if (ProductsA != theInitialState.GetA()) {
458        G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl;
459        G4cout << "G4Evaporation.cc: Barionic Number Conservation test for evaporation fragments" 
460               << G4endl; 
461        G4cout << "Initial A = " << theInitialState.GetA() 
462               << "   Fragments A = " << ProductsA << "   Diference --> " 
463               << theInitialState.GetA() - ProductsA << G4endl;
464    }
465    if (ProductsZ != theInitialState.GetZ()) {
466        G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
467        G4cout << "G4Evaporation.cc: Charge Conservation test for evaporation fragments" 
468               << G4endl; 
469        G4cout << "Initial Z = " << theInitialState.GetZ() 
470               << "   Fragments Z = " << ProductsZ << "   Diference --> " 
471               << theInitialState.GetZ() - ProductsZ << G4endl;
472    }
473    if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) {
474        G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
475        G4cout << "G4Evaporation.cc: Energy Conservation test for evaporation fragments" 
476               << G4endl; 
477        G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
478               << "   Fragments E = " << ProductsEnergy/MeV  << " MeV   Diference --> " 
479               << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl;
480    } 
481    if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV || 
482        std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV ||
483        std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) {
484        G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
485        G4cout << "G4Evaporation.cc: Momentum Conservation test for evaporation fragments" 
486               << G4endl; 
487        G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
488               << "   Fragments P = " << ProductsMomentum  << " MeV   Diference --> " 
489               << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl;
490    }
491    return;
492}
493#endif
494
495
496
497
Note: See TracBrowser for help on using the repository browser.