Bug Summary

File:scratch/gluex/scan-build-work/hdgeant4/src/G4fixes/G4ExcitationHandler.cc
Location:line 141, column 26
Description:Value stored to 'param' during its initialization is never read

Annotated Source Code

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// $Id: G4ExcitationHandler.cc 104984 2017-07-03 15:13:37Z gcosmo $
27//
28// Hadronic Process: Nuclear De-excitations
29// by V. Lara (May 1998)
30//
31//
32// Modified:
33// 30 June 1998 by V. Lara:
34// -Modified the Transform method for use G4ParticleTable and
35// therefore G4IonTable. It makes possible to convert all kind
36// of fragments (G4Fragment) produced in deexcitation to
37// G4DynamicParticle
38// -It uses default algorithms for:
39// Evaporation: G4Evaporation
40// MultiFragmentation: G4StatMF
41// Fermi Breakup model: G4FermiBreakUp
42// 24 Jul 2008 by M. A. Cortes Giraldo:
43// -Max Z,A for Fermi Break-Up turns to 9,17 by default
44// -BreakItUp() reorganised and bug in Evaporation loop fixed
45// -Transform() optimised
46// (September 2008) by J. M. Quesada. External choices have been added for :
47// -inverse cross section option (default OPTxs=3)
48// -superimposed Coulomb barrier (if useSICB is set true, by default it is false)
49// September 2009 by J. M. Quesada:
50// -according to Igor Pshenichnov, SMM will be applied (just in case) only once.
51// 27 Nov 2009 by V.Ivanchenko:
52// -cleanup the logic, reduce number internal vectors, fixed memory leak.
53// 11 May 2010 by V.Ivanchenko:
54// -FermiBreakUp activated, used integer Z and A, used BreakUpFragment method for
55// final photon deexcitation; used check on adundance of a fragment, decay
56// unstable fragments with A <5
57// 22 March 2011 by V.Ivanchenko: general cleanup and addition of a condition:
58// products of Fermi Break Up cannot be further deexcited by this model
59// 30 March 2011 by V.Ivanchenko removed private inline methods, moved Set methods
60// to the source
61// 23 January 2012 by V.Ivanchenko general cleanup including destruction of
62// objects, propagate G4PhotonEvaporation pointer to G4Evaporation class and
63// not delete it here
64
65#include "G4ExcitationHandler.hh"
66#include "G4SystemOfUnits.hh"
67#include "G4LorentzVector.hh"
68#include "G4NistManager.hh"
69#include "G4ParticleTable.hh"
70#include "G4ParticleTypes.hh"
71#include "G4Ions.hh"
72#include "G4Electron.hh"
73
74#include "G4VMultiFragmentation.hh"
75#include "G4VFermiBreakUp.hh"
76
77#include "G4VEvaporation.hh"
78#include "G4VEvaporationChannel.hh"
79#include "G4Evaporation.hh"
80#include "G4PhotonEvaporation.hh"
81#include "G4StatMF.hh"
82#include "G4FermiBreakUpVI.hh"
83#include "G4NuclearLevelData.hh"
84#include "G4Pow.hh"
85
86G4ExcitationHandler::G4ExcitationHandler()
87 : maxZForFermiBreakUp(9),maxAForFermiBreakUp(17),
88 fVerbose(0),isInitialised(false),isEvapLocal(true)
89{
90 theTableOfIons = G4ParticleTable::GetParticleTable()->GetIonTable();
91 nist = G4NistManager::Instance();
92
93 theMultiFragmentation = nullptr;
94 theFermiModel = nullptr;
95 G4Pow::GetInstance();
96 theEvaporation = new G4Evaporation();
97 thePhotonEvaporation = theEvaporation->GetPhotonEvaporation();
98 theResults.reserve(60);
99 results.reserve(30);
100 theEvapList.reserve(30);
101 thePhotoEvapList.reserve(10);
102
103 SetParameters();
104 electron = G4Electron::Electron();
105
106 if(fVerbose > 0) { G4cout(*G4cout_p) << "### New handler " << this << G4endlstd::endl; }
107}
108
109G4ExcitationHandler::~G4ExcitationHandler()
110{
111 //G4cout << "### Delete handler " << this << G4endl;
112 delete theMultiFragmentation;
113 delete theFermiModel;
114 if(isEvapLocal) { delete theEvaporation; }
115}
116
117void G4ExcitationHandler::SetParameters()
118{
119 G4DeexPrecoParameters* param =
120 G4NuclearLevelData::GetInstance()->GetParameters();
121 isActive = true;
122 if(fDummy == param->GetDeexChannelsType()) { isActive = false; }
123 minEForMultiFrag = param->GetMinExPerNucleounForMF();
124 minExcitation = param->GetMinExcitation();
125 icID = param->GetInternalConversionID();
126 if(isActive) {
127 if(!thePhotonEvaporation) {
128 SetPhotonEvaporation(new G4PhotonEvaporation());
129 }
130 if(!theFermiModel) { SetFermiModel(new G4FermiBreakUpVI()); }
131 if(!theMultiFragmentation) { SetMultiFragmentation(new G4StatMF()); }
132 }
133}
134
135void G4ExcitationHandler::Initialise()
136{
137 if(isInitialised) { return; }
138 if(fVerbose > 0) {
139 G4cout(*G4cout_p) << "G4ExcitationHandler::Initialise() started " << this << G4endlstd::endl;
140 }
141 G4DeexPrecoParameters* param =
Value stored to 'param' during its initialization is never read
142 G4NuclearLevelData::GetInstance()->GetParameters();
143 isInitialised = true;
144 SetParameters();
145 if(isActive) {
146 theFermiModel->Initialise();
147 theEvaporation->InitialiseChannels();
148 }
149 if (1==0)
150 param->Dump();
151}
152
153void G4ExcitationHandler::SetEvaporation(G4VEvaporation* ptr, G4bool flag)
154{
155 if(ptr && ptr != theEvaporation) {
156 delete theEvaporation;
157 theEvaporation = ptr;
158 thePhotonEvaporation = ptr->GetPhotonEvaporation();
159 theEvaporation->SetFermiBreakUp(theFermiModel);
160 isEvapLocal = flag;
161 }
162}
163
164void
165G4ExcitationHandler::SetMultiFragmentation(G4VMultiFragmentation* ptr)
166{
167 if(ptr && ptr != theMultiFragmentation) {
168 delete theMultiFragmentation;
169 theMultiFragmentation = ptr;
170 }
171}
172
173void G4ExcitationHandler::SetFermiModel(G4VFermiBreakUp* ptr)
174{
175 if(ptr && ptr != theFermiModel) {
176 delete theFermiModel;
177 theFermiModel = ptr;
178 theEvaporation->SetFermiBreakUp(theFermiModel);
179 }
180}
181
182void
183G4ExcitationHandler::SetPhotonEvaporation(G4VEvaporationChannel* ptr)
184{
185 if(ptr && ptr != thePhotonEvaporation) {
186 thePhotonEvaporation = ptr;
187 theEvaporation->SetPhotonEvaporation(ptr);
188 }
189}
190
191void G4ExcitationHandler::SetDeexChannelsType(G4DeexChannelType val)
192{
193 G4Evaporation* evap = static_cast<G4Evaporation*>(theEvaporation);
194 if(val == fDummy) {
195 isActive = false;
196 return;
197 }
198 if(!evap) { return; }
199 if(val == fEvaporation) {
200 evap->SetDefaultChannel();
201 } else if(val == fCombined) {
202 evap->SetCombinedChannel();
203 } else if(val == fGEM) {
204 evap->SetGEMChannel();
205 }
206 evap->InitialiseChannels();
207 if(G4Threading::IsMasterThread()) {
208 G4cout(*G4cout_p) << "Number of de-excitation channels is changed to "
209 << theEvaporation->GetNumberOfChannels();
210 if(fVerbose > 0) { G4cout(*G4cout_p) << " " << this; }
211 G4cout(*G4cout_p) << G4endlstd::endl;
212 }
213}
214
215G4ReactionProductVector *
216G4ExcitationHandler::BreakItUp(const G4Fragment & theInitialState)
217{
218 // Variables existing until end of method
219 G4Fragment * theInitialStatePtr = new G4Fragment(theInitialState);
220 if(fVerbose > 1) {
221 G4cout(*G4cout_p) << "@@@@@@@@@@ Start G4Excitation Handler @@@@@@@@@@@@@ " << G4endlstd::endl;
222 G4cout(*G4cout_p) << theInitialState << G4endlstd::endl;
223 }
224 if(!isInitialised) { Initialise(); }
225
226 // pointer to fragment vector which receives temporal results
227 G4FragmentVector * theTempResult = nullptr;
228
229 theResults.clear();
230 thePhotoEvapList.clear();
231 theEvapList.clear();
232
233 // Variables to describe the excited configuration
234 G4double exEnergy = theInitialState.GetExcitationEnergy();
235 G4int A = theInitialState.GetA_asInt();
236 G4int Z = theInitialState.GetZ_asInt();
237
238 // In case A <= 1 the fragment will not perform any nucleon emission
239 if (A <= 1 || !isActive) {
240 theResults.push_back( theInitialStatePtr );
241
242 // check if a fragment is stable
243 } else if(exEnergy < minExcitation && nist->GetIsotopeAbundance(Z, A) > 0.0) {
244 theResults.push_back( theInitialStatePtr );
245
246 // JMQ 150909: first step in de-excitation is treated separately
247 // Fragments after the first step are stored in theEvapList
248 } else {
249 if((A<maxAForFermiBreakUp && Z<maxZForFermiBreakUp)
250 || exEnergy <= minEForMultiFrag*A) {
251 theEvapList.push_back(theInitialStatePtr);
252
253 // Statistical Multifragmentation will take place only once
254 } else {
255 theTempResult = theMultiFragmentation->BreakItUp(theInitialState);
256 if(!theTempResult) {
257 theEvapList.push_back(theInitialStatePtr);
258 } else {
259 size_t nsec = theTempResult->size();
260
261 // no fragmentation
262 if(0 == nsec) {
263 theEvapList.push_back(theInitialStatePtr);
264
265 // secondary are produced - sort out secondary fragments
266 } else {
267 G4bool deletePrimary = true;
268 G4FragmentVector::iterator j;
269 for (j = theTempResult->begin(); j != theTempResult->end(); ++j) {
270
271 if((*j) == theInitialStatePtr) { deletePrimary = false; }
272 A = (*j)->GetA_asInt();
273
274 // gamma, p, n
275 if(A <= 1) {
276 theResults.push_back(*j);
277
278 // Analyse fragment A > 1
279 } else {
280 G4double exEnergy1 = (*j)->GetExcitationEnergy();
281
282 // cold fragments
283 if(exEnergy1 < minExcitation) {
284 Z = (*j)->GetZ_asInt();
285 if(nist->GetIsotopeAbundance(Z, A) > 0.0) {
286 theResults.push_back(*j); // stable fragment
287 } else {
288 theEvapList.push_back(*j);
289 }
290 // hot fragments are unstable
291 } else {
292 theEvapList.push_back(*j);
293 }
294 }
295 }
296 if( deletePrimary ) { delete theInitialStatePtr; }
297 }
298 delete theTempResult; // end multifragmentation
299 }
300 }
301 }
302 if(fVerbose > 2) {
303 G4cout(*G4cout_p) << "## After first step " << theEvapList.size() << " for evap; "
304 << thePhotoEvapList.size() << " for photo-evap; "
305 << theResults.size() << " results. " << G4endlstd::endl;
306 }
307 // -----------------------------------
308 // FermiBreakUp and De-excitation loop
309 // -----------------------------------
310
311 static const G4int countmax = 1000;
312 G4Fragment* frag;
313 size_t kk;
314 for (kk=0; kk<theEvapList.size(); ++kk) {
315 frag = theEvapList[kk];
316 if(fVerbose > 2) {
317 G4cout(*G4cout_p) << "Next evaporate: " << G4endlstd::endl;
318 G4cout(*G4cout_p) << *frag << G4endlstd::endl;
319 }
320 if(kk >= countmax) {
321 G4ExceptionDescription ed;
322 ed << "Infinite loop in the de-excitation module: " << kk
323 << " iterations \n"
324 << " Initial fragment: \n" << theInitialState
325 << "\n Current fragment: \n" << *frag;
326 G4Exception("G4ExcitationHandler::BreakItUp","had0333",FatalException,
327 ed,"Stop execution");
328
329 }
330 A = frag->GetA_asInt();
331 Z = frag->GetZ_asInt();
332 results.clear();
333
334 // Fermi Break-Up
335 if(theFermiModel->IsApplicable(Z, A, frag->GetExcitationEnergy())) {
336 theFermiModel->BreakFragment(&results, frag);
337 size_t nsec = results.size();
338 if(fVerbose > 2) { G4cout(*G4cout_p) << "FermiBreakUp Nsec= " << nsec << G4endlstd::endl; }
339
340 // FBU takes care to delete input fragment or add it to the results
341 // The secondary may be excited - photo-evaporation should be applied
342 for(size_t j=0; j<nsec; ++j) {
343 exEnergy = results[j]->GetExcitationEnergy();
344 if(exEnergy < minExcitation) { theResults.push_back(results[j]); }
345 else { thePhotoEvapList.push_back(results[j]); }
346 }
347 continue;
348 }
349 // apply Evaporation, residual nucleus is always added to the results
350 theEvaporation->BreakFragment(&results, frag);
351 size_t nsec = results.size();
352 if(fVerbose > 2) { G4cout(*G4cout_p) << "Evaporation Nsec= " << nsec << G4endlstd::endl; }
353
354 // no evaporation
355 if(1 >= nsec) {
356 theResults.push_back(frag);
357 continue;
358 }
359
360 // Sort out secondary fragments
361 for (size_t j = 0; j<nsec; ++j) {
362 if(fVerbose > 3) {
363 G4cout(*G4cout_p) << "Evaporated product #" << j << G4endlstd::endl;
364 G4cout(*G4cout_p) << results[j] << G4endlstd::endl;
365 }
366 A = results[j]->GetA_asInt();
367 if(A <= 1) {
368 theResults.push_back(results[j]); // gamma, p, n
369 continue;
370 }
371 exEnergy = results[j]->GetExcitationEnergy();
372
373 // hot fragment
374 if(exEnergy >= minExcitation) {
375 theEvapList.push_back(results[j]);
376
377 // cold fragment
378 } else {
379 Z = results[j]->GetZ_asInt();
380
381 // natural isotope
382 if(nist->GetIsotopeAbundance(Z, A) > 0.0) {
383 theResults.push_back(results[j]); // stable fragment
384
385 } else {
386 theEvapList.push_back(results[j]);
387 }
388 }
389 } // end of loop on secondary
390 } // end of the loop over theEvapList
391 if(fVerbose > 2) {
392 G4cout(*G4cout_p) << "## After 2nd step " << theEvapList.size() << " was evap; "
393 << thePhotoEvapList.size() << " for photo-evap; "
394 << theResults.size() << " results. " << G4endlstd::endl;
395 }
396 // -----------------------
397 // Photon-Evaporation loop
398 // -----------------------
399
400 // at this point only photon evaporation is possible
401 size_t kkmax = thePhotoEvapList.size();
402 for (kk=0; kk<kkmax; ++kk) {
403 frag = thePhotoEvapList[kk];
404 if(fVerbose > 2) {
405 G4cout(*G4cout_p) << "Next photon evaporate: " << thePhotonEvaporation << G4endlstd::endl;
406 G4cout(*G4cout_p) << *frag << G4endlstd::endl;
407 }
408 exEnergy = frag->GetExcitationEnergy();
409
410 // photon de-excitation only for hot fragments
411 if(exEnergy > minExcitation) {
412 thePhotonEvaporation->BreakUpChain(&theResults, frag);
413 }
414
415 // primary fragment is kept
416 theResults.push_back(frag);
417
418 } // end of photon-evaporation loop
419 if(fVerbose > 2) {
420 G4cout(*G4cout_p) << "## After 3d step " << theEvapList.size() << " was evap; "
421 << thePhotoEvapList.size() << " was photo-evap; "
422 << theResults.size() << " results. " << G4endlstd::endl;
423 }
424 G4ReactionProductVector * theReactionProductVector =
425 new G4ReactionProductVector();
426
427 // MAC (24/07/08)
428 // To optimise the storing speed, we reserve space in memory for the vector
429 theReactionProductVector->reserve( theResults.size() );
430
431 G4int theFragmentA, theFragmentZ;
432
433 if(fVerbose > 1) {
434 G4cout(*G4cout_p) << "### ExcitationHandler provides " << theResults.size()
435 << " evaporated products:" << G4endlstd::endl;
436 }
437 kkmax = theResults.size();
438 for (kk=0; kk<kkmax; ++kk) {
439 frag = theResults[kk];
440
441 // in the case of dummy de-excitation, excitation energy is transfered
442 // into kinetic energy
443 if(!isActive && 0 == kk) {
444 G4double mass = frag->GetGroundStateMass();
445 G4double ptot = (frag->GetMomentum()).vect().mag();
446 G4double etot = (frag->GetMomentum()).e();
447 G4double fac = (etot <= mass || 0.0 == ptot) ? 0.0
448 : std::sqrt((etot - mass)*(etot + mass))/ptot;
449 G4LorentzVector lv((frag->GetMomentum()).px()*fac,
450 (frag->GetMomentum()).py()*fac,
451 (frag->GetMomentum()).pz()*fac, etot);
452 frag->SetMomentum(lv);
453 }
454 if(fVerbose > 1) {
455 G4cout(*G4cout_p) << kk << "-th fragment " << frag;
456 if(frag->NuclearPolarization()) {
457 G4cout(*G4cout_p) << " " << frag->NuclearPolarization();
458 }
459 G4cout(*G4cout_p) << G4endlstd::endl;
460 G4cout(*G4cout_p) << *frag << G4endlstd::endl;
461 }
462
463 theFragmentA = frag->GetA_asInt();
464 theFragmentZ = frag->GetZ_asInt();
465 G4double etot= frag->GetMomentum().e();
466 G4double eexc = 0.0;
467 const G4ParticleDefinition* theKindOfFragment = nullptr;
468 if (theFragmentA == 0) { // photon or e-
469 theKindOfFragment = frag->GetParticleDefinition();
470 } else if (theFragmentA == 1 && theFragmentZ == 0) { // neutron
471 theKindOfFragment = G4Neutron::NeutronDefinition();
472 } else if (theFragmentA == 1 && theFragmentZ == 1) { // proton
473 theKindOfFragment = G4Proton::ProtonDefinition();
474 } else if (theFragmentA == 2 && theFragmentZ == 1) { // deuteron
475 theKindOfFragment = G4Deuteron::DeuteronDefinition();
476 } else if (theFragmentA == 3 && theFragmentZ == 1) { // triton
477 theKindOfFragment = G4Triton::TritonDefinition();
478 } else if (theFragmentA == 3 && theFragmentZ == 2) { // helium3
479 theKindOfFragment = G4He3::He3Definition();
480 } else if (theFragmentA == 4 && theFragmentZ == 2) { // alpha
481 theKindOfFragment = G4Alpha::AlphaDefinition();;
482 } else {
483
484 // fragment
485 eexc = frag->GetExcitationEnergy();
486 G4int idxf = frag->GetFloatingLevelNumber();
487 if(eexc < minExcitation) {
488 eexc = 0.0;
489 idxf = 0;
490 }
491
492 theKindOfFragment = theTableOfIons->GetIon(theFragmentZ,theFragmentA,eexc,
493 G4Ions::FloatLevelBase(idxf));
494 if(fVerbose > 1) {
495 G4cout(*G4cout_p) << "### EXCH: Find ion Z= " << theFragmentZ << " A= " << theFragmentA
496 << " Eexc(MeV)= " << eexc/MeV << " idx= " << idxf
497 << " " << theKindOfFragment << G4endlstd::endl;
498 }
499 }
500 // fragment identified
501 if(theKindOfFragment) {
502 G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
503 theNew->SetMomentum(frag->GetMomentum().vect());
504 theNew->SetTotalEnergy(etot);
505 theNew->SetFormationTime(frag->GetCreationTime());
506 if(theKindOfFragment == electron) { theNew->SetCreatorModel(icID); }
507 theReactionProductVector->push_back(theNew);
508
509 // fragment not found out ground state is created
510 } else {
511 theKindOfFragment =
512 theTableOfIons->GetIon(theFragmentZ,theFragmentA,0.0,noFloatG4Ions::G4FloatLevelBase::no_Float,0);
513 if(theKindOfFragment) {
514 G4ThreeVector mom(0.0,0.0,0.0);
515 G4double ionmass = theKindOfFragment->GetPDGMass();
516 if(etot <= ionmass) {
517 etot = ionmass;
518 } else {
519 G4double ptot = std::sqrt((etot - ionmass)*(etot + ionmass));
520 mom = (frag->GetMomentum().vect().unit())*ptot;
521 }
522 G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
523 theNew->SetMomentum(mom);
524 theNew->SetTotalEnergy(etot);
525 theNew->SetFormationTime(frag->GetCreationTime());
526 theReactionProductVector->push_back(theNew);
527 if(fVerbose > 2) {
528 G4cout(*G4cout_p) << "### Find ion Z= " << theFragmentZ << " A= " << theFragmentA
529 << " ground state, energy corrected E(MeV)= " << etot << G4endlstd::endl;
530 }
531 }
532 }
533 delete frag;
534 if(fVerbose > 1) { G4cout(*G4cout_p) << "G4Fragment #" << kk << " is deleted" << G4endlstd::endl; }
535 }
536 if(fVerbose > 2) {
537 G4cout(*G4cout_p) << "@@@@@@@@@@ End G4Excitation Handler "<< G4endlstd::endl;
538 }
539 return theReactionProductVector;
540}
541
542void G4ExcitationHandler::ModelDescription(std::ostream& outFile) const
543{
544 outFile << "G4ExcitationHandler description\n"
545 << "This class samples de-excitation of excited nucleus using\n"
546 << "Fermi Break-up model for light fragments (Z < 9, A < 17), "
547 << "evaporation, fission, and photo-evaporation models. Evaporated\n"
548 << "particle may be proton, neutron, and other light fragment \n"
549 << "(Z < 13, A < 29). During photon evaporation produced gamma \n"
550 << "or electrons due to internal conversion \n";
551}
552
553
554
555
556
557