00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058 #include "G4PhotonEvaporation.hh"
00059
00060 #include "globals.hh"
00061 #include "G4SystemOfUnits.hh"
00062 #include "Randomize.hh"
00063 #include "G4Gamma.hh"
00064 #include "G4LorentzVector.hh"
00065 #include "G4VGammaTransition.hh"
00066 #include "G4Fragment.hh"
00067 #include "G4FragmentVector.hh"
00068 #include "G4ContinuumGammaDeexcitation.hh"
00069 #include "G4DiscreteGammaDeexcitation.hh"
00070 #include "G4E1Probability.hh"
00071
00072 G4PhotonEvaporation::G4PhotonEvaporation()
00073 :_verbose(0),_myOwnProbAlgorithm (true),
00074 _eOccupancy(0), _vShellNumber(-1),_gammaE(0.)
00075 {
00076 _probAlgorithm = new G4E1Probability;
00077
00078 G4double timeLimit = DBL_MAX;
00079 char* env = getenv("G4AddTimeLimitToPhotonEvaporation");
00080 if(env) { timeLimit = 1.e-16*second; }
00081
00082 G4DiscreteGammaDeexcitation* p = new G4DiscreteGammaDeexcitation();
00083 p->SetICM(false);
00084 p->SetTimeLimit(timeLimit);
00085
00086 _discrDeexcitation = p;
00087 _contDeexcitation = new G4ContinuumGammaDeexcitation;
00088 _nucleus = 0;
00089 }
00090
00091 G4PhotonEvaporation::~G4PhotonEvaporation()
00092 {
00093 if(_myOwnProbAlgorithm) delete _probAlgorithm;
00094 delete _discrDeexcitation;
00095 delete _contDeexcitation;
00096 }
00097
00098 G4Fragment* G4PhotonEvaporation::EmittedFragment(G4Fragment* nucleus)
00099 {
00100
00101 _nucleus = nucleus;
00102
00103
00104 _contDeexcitation->SetNucleus(_nucleus);
00105 _contDeexcitation->Initialize();
00106
00107 if(_contDeexcitation->CanDoTransition()) {
00108 G4Fragment* gamma = _contDeexcitation->GenerateGamma();
00109 if(gamma) {
00110 if (_verbose > 0) {
00111 G4cout << "G4PhotonEvaporation::EmittedFragment continium deex: "
00112 << gamma << G4endl;
00113 G4cout << " Residual: " << nucleus << G4endl;
00114 }
00115 return gamma;
00116 }
00117 }
00118
00119
00120 _discrDeexcitation->SetNucleus(_nucleus);
00121 _discrDeexcitation->Initialize();
00122
00123 if(_discrDeexcitation->CanDoTransition()) {
00124 G4Fragment* gamma = _discrDeexcitation->GenerateGamma();
00125 if(gamma) {
00126 if (_verbose > 0) {
00127 G4cout << "G4PhotonEvaporation::EmittedFragment discrete deex: "
00128 << gamma << G4endl;
00129 G4cout << " Residual: " << nucleus << G4endl;
00130 }
00131 return gamma;
00132 }
00133 }
00134
00135 if (_verbose > 0) {
00136 G4cout << "G4PhotonEvaporation unable emit gamma: "
00137 << nucleus << G4endl;
00138 }
00139 return 0;
00140 }
00141
00142 G4FragmentVector* G4PhotonEvaporation::BreakUpFragment(G4Fragment* nucleus)
00143 {
00144
00145
00146 _nucleus = nucleus;
00147 _contDeexcitation->SetNucleus(_nucleus);
00148 _discrDeexcitation->SetNucleus(_nucleus);
00149
00150
00151 G4FragmentVector* products = _contDeexcitation->DoChain();
00152 if( !products ) { products = new G4FragmentVector(); }
00153
00154 if (_verbose > 0) {
00155 G4cout << "G4PhotonEvaporation::BreakUpFragment " << products->size()
00156 << " gammas from ContinuumDeexcitation " << G4endl;
00157 G4cout << " Residual: " << nucleus << G4endl;
00158 }
00159
00160 G4FragmentVector* discrProducts = _discrDeexcitation->DoChain();
00161 if(discrProducts) {
00162 _eOccupancy = _discrDeexcitation->GetEO();
00163 _vShellNumber = _discrDeexcitation->GetVacantSN();
00164
00165
00166 _discrDeexcitation->SetVaccantSN(-1);
00167
00168 if (_verbose > 0) {
00169 G4cout << "G4PhotonEvaporation::BreakUpFragment " << discrProducts->size()
00170 << " gammas from DiscreteDeexcitation " << G4endl;
00171 G4cout << " Residual: " << nucleus << G4endl;
00172 }
00173 G4FragmentVector::iterator i;
00174 for (i = discrProducts->begin(); i != discrProducts->end(); ++i)
00175 {
00176 products->push_back(*i);
00177 }
00178 delete discrProducts;
00179 }
00180
00181 if (_verbose > 0) {
00182 G4cout << "*-*-* Photon evaporation: " << products->size() << G4endl;
00183 }
00184 return products;
00185 }
00186
00187 G4FragmentVector* G4PhotonEvaporation::BreakUp(const G4Fragment& nucleus)
00188 {
00189
00190 _nucleus = new G4Fragment(nucleus);
00191
00192 _contDeexcitation->SetNucleus(_nucleus);
00193 _discrDeexcitation->SetNucleus(_nucleus);
00194
00195
00196
00197
00198
00199 G4FragmentVector* products = _contDeexcitation->DoTransition();
00200 if( !products ) { products = new G4FragmentVector(); }
00201 else if(_verbose > 0) {
00202 G4cout << "G4PhotonEvaporation::BreakUp " << products->size()
00203 << " gammas from ContinuesDeexcitation " << G4endl;
00204 G4cout << " Residual: " << nucleus << G4endl;
00205 }
00206
00207 if (0 == products->size())
00208 {
00209
00210 G4FragmentVector* discrProducts = _discrDeexcitation->DoTransition();
00211
00212 if (discrProducts) {
00213 _eOccupancy = _discrDeexcitation->GetEO();
00214 _vShellNumber = _discrDeexcitation->GetVacantSN();
00215
00216
00217 _discrDeexcitation->SetVaccantSN(-1);
00218
00219 if (_verbose > 0) {
00220 G4cout << " = BreakUp = " << discrProducts->size()
00221 << " gammas from DiscreteDeexcitation "
00222 << G4endl;
00223 G4cout << " Residual: " << nucleus << G4endl;
00224 }
00225 G4FragmentVector::iterator i;
00226 for (i = discrProducts->begin(); i != discrProducts->end(); ++i)
00227 {
00228 products->push_back(*i);
00229 }
00230 delete discrProducts;
00231 }
00232 }
00233
00234
00235 products->push_back(_nucleus);
00236
00237 if (_verbose > 0) {
00238 G4cout << "*-*-*-* Photon evaporation: " << products->size() << G4endl;
00239 }
00240
00241 return products;
00242 }
00243
00244 G4FragmentVector* G4PhotonEvaporation::BreakItUp(const G4Fragment& nucleus)
00245 {
00246
00247 _nucleus = new G4Fragment(nucleus);
00248 _contDeexcitation->SetNucleus(_nucleus);
00249 _discrDeexcitation->SetNucleus(_nucleus);
00250
00251
00252
00253
00254 G4FragmentVector* products = _contDeexcitation->DoChain();
00255 if( !products ) { products = new G4FragmentVector; }
00256
00257
00258 if (_verbose > 0) {
00259 G4cout << " = BreakItUp = " << products->size()
00260 << " gammas from ContinuumDeexcitation " << G4endl;
00261 }
00262
00263
00264 G4FragmentVector* discrProducts = _discrDeexcitation->DoChain();
00265 if(discrProducts) {
00266 _eOccupancy = _discrDeexcitation->GetEO();
00267 _vShellNumber = _discrDeexcitation->GetVacantSN();
00268
00269
00270 _discrDeexcitation->SetVaccantSN(-1);
00271
00272 if (_verbose > 0) {
00273 G4cout << " = BreakItUp = " << discrProducts->size()
00274 << " gammas from DiscreteDeexcitation " << G4endl;
00275 }
00276 G4FragmentVector::iterator i;
00277 for (i = discrProducts->begin(); i != discrProducts->end(); ++i)
00278 {
00279 products->push_back(*i);
00280 }
00281 delete discrProducts;
00282 }
00283
00284 products->push_back(_nucleus);
00285
00286 if (_verbose > 0) {
00287 G4cout << "*-*-* Photon evaporation: " << products->size() << G4endl;
00288 }
00289 return products;
00290 }
00291
00292 G4double
00293 G4PhotonEvaporation::GetEmissionProbability(G4Fragment* theNucleus)
00294 {
00295 _nucleus = theNucleus;
00296 G4double prob =
00297 _probAlgorithm->EmissionProbability(*_nucleus,_nucleus->GetExcitationEnergy());
00298 return prob;
00299 }
00300
00301
00302 void G4PhotonEvaporation::SetEmissionStrategy(G4VEmissionProbability * probAlgorithm)
00303 {
00304
00305
00306
00307 if(_myOwnProbAlgorithm) delete _probAlgorithm;
00308
00309 _probAlgorithm = probAlgorithm;
00310
00311 _myOwnProbAlgorithm = false;
00312 }
00313
00314
00315 void G4PhotonEvaporation::SetVerboseLevel(G4int verbose)
00316 {
00317 _verbose = verbose;
00318 _contDeexcitation->SetVerboseLevel(verbose);
00319 _discrDeexcitation->SetVerboseLevel(verbose);
00320 }
00321
00322 void G4PhotonEvaporation::SetICM(G4bool ic)
00323 {
00324 (static_cast <G4DiscreteGammaDeexcitation*> (_discrDeexcitation))->SetICM(ic);
00325 }
00326
00327 void G4PhotonEvaporation::SetMaxHalfLife(G4double hl)
00328 {
00329 (static_cast <G4DiscreteGammaDeexcitation*> (_discrDeexcitation))->SetHL(hl);
00330 }
00331
00332 void G4PhotonEvaporation::SetTimeLimit(G4double val)
00333 {
00334 _discrDeexcitation->SetTimeLimit(val);
00335 }
00336
00337 void G4PhotonEvaporation::RDMForced(G4bool fromRDM)
00338 {
00339 (static_cast <G4DiscreteGammaDeexcitation*> (_discrDeexcitation))->SetRDM(fromRDM);
00340 }
00341
00342 void G4PhotonEvaporation::SetEOccupancy(G4ElectronOccupancy eo)
00343 {
00344 _discrDeexcitation->SetEO(eo);
00345 }
00346
00347 #ifdef debug
00348 void G4PhotonEvaporation::CheckConservation(const G4Fragment & theInitialState,
00349 G4FragmentVector * Result) const
00350 {
00351 G4double ProductsEnergy =0;
00352 G4ThreeVector ProductsMomentum;
00353 G4int ProductsA = 0;
00354 G4int ProductsZ = 0;
00355 G4FragmentVector::iterator h;
00356 for (h = Result->begin(); h != Result->end(); h++) {
00357 G4LorentzVector tmp = (*h)->GetMomentum();
00358 ProductsEnergy += tmp.e();
00359 ProductsMomentum += tmp.vect();
00360 ProductsA += (*h)->GetA_asInt();
00361 ProductsZ += (*h)->GetZ_asInt();
00362 }
00363
00364 if (ProductsA != theInitialState.GetA_asInt()) {
00365 G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl;
00366 G4cout << "G4PhotonEvaporation.cc: Barionic Number Conservation test for evaporation fragments"
00367 << G4endl;
00368 G4cout << "Initial A = " << theInitialState.GetA_asInt()
00369 << " Fragments A = " << ProductsA << " Diference --> "
00370 << theInitialState.GetA_asInt() - ProductsA << G4endl;
00371 }
00372 if (ProductsZ != theInitialState.GetZ_asInt()) {
00373 G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
00374 G4cout << "G4PhotonEvaporation.cc: Charge Conservation test for evaporation fragments"
00375 << G4endl;
00376 G4cout << "Initial Z = " << theInitialState.GetZ_asInt()
00377 << " Fragments Z = " << ProductsZ << " Diference --> "
00378 << theInitialState.GetZ_asInt() - ProductsZ << G4endl;
00379 }
00380 if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 1.0*keV) {
00381 G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
00382 G4cout << "G4PhotonEvaporation.cc: Energy Conservation test for evaporation fragments"
00383 << G4endl;
00384 G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
00385 << " Fragments E = " << ProductsEnergy/MeV << " MeV Diference --> "
00386 << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl;
00387 }
00388 if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 1.0*keV ||
00389 std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 1.0*keV ||
00390 std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 1.0*keV) {
00391 G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
00392 G4cout << "G4PhotonEvaporation.cc: Momentum Conservation test for evaporation fragments"
00393 << G4endl;
00394 G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
00395 << " Fragments P = " << ProductsMomentum << " MeV Diference --> "
00396 << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl;
00397 }
00398 return;
00399 }
00400 #endif
00401
00402
00403