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