ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4InuclNuclei.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4InuclNuclei.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 //
27 // 20100301 M. Kelsey -- Add function to create unphysical nuclei for use
28 // as temporary final-state fragments.
29 // 20100319 M. Kelsey -- Add information message to makeNuclearFragment().
30 // Use new GetBindingEnergy() function instead of bindingEnergy().
31 // 20100622 M. Kelsey -- Use local "bindingEnergy()" function to call through.
32 // 20100627 M. Kelsey -- Test for non-physical fragments and abort job.
33 // 20100630 M. Kelsey -- Use excitation energy in G4Ions
34 // 20100714 M. Kelsey -- Use G4DynamicParticle::theDynamicalMass to deal with
35 // excitation energy without instantianting "infinite" G4PartDefns.
36 // 20100719 M. Kelsey -- Change excitation energy without altering momentum
37 // 20100906 M. Kelsey -- Add fill() functions to rewrite contents
38 // 20100910 M. Kelsey -- Add clearExitonConfiguration() to fill() functions
39 // 20100914 M. Kelsey -- Make printout symmetric with G4InuclElemPart,
40 // migrate to integer A and Z
41 // 20100924 M. Kelsey -- Add constructor to copy G4Fragment input, and output
42 // functions to create G4Fragment
43 // 20110214 M. Kelsey -- Replace integer "model" with enum
44 // 20110308 M. Kelsey -- Follow new G4Fragment interface for hole types
45 // 20110427 M. Kelsey -- Remove PDG-code warning
46 // 20110721 M. Kelsey -- Follow base-class ctor change to pass model directly
47 // 20110829 M. Kelsey -- Add constructor to copy G4V3DNucleus input
48 // 20110919 M. Kelsey -- Special case: Allow fill(A=0,Z=0) to make dummy
49 // 20110922 M. Kelsey -- Add stream argument to printParticle() => print()
50 // 20121009 M. Kelsey -- Add report of excitons if non-empty
51 // 20130314 M. Kelsey -- Use G4IonList typedef for fragment map, encapsulate
52 // it in a static function with mutexes.
53 // 20130620 M. Kelsey -- Address Coverity #37503, check self in op=()
54 // 20140523 M. Kelsey -- Avoid FPE in setExcitationEnergy() for zero Ekin
55 // 20150608 M. Kelsey -- Label all while loops as terminating.
56 
57 #include "G4InuclNuclei.hh"
58 #include "G4AutoLock.hh"
59 #include "G4Fragment.hh"
60 #include "G4HadronicException.hh"
62 #include "G4IonTable.hh"
63 #include "G4Ions.hh"
64 #include "G4NucleiProperties.hh"
65 #include "G4Nucleon.hh"
66 #include "G4ParticleDefinition.hh"
67 #include "G4ParticleTable.hh"
68 #include "G4SystemOfUnits.hh"
69 #include "G4Threading.hh"
70 #include "G4V3DNucleus.hh"
71 
72 #include <assert.h>
73 #include <sstream>
74 #include <map>
75 
76 using namespace G4InuclSpecialFunctions;
77 
78 
79 // Convert contents from (via constructor) and to G4Fragment
80 
83  : G4InuclParticle() {
84  copy(aFragment, model);
85 }
86 
87 void G4InuclNuclei::copy(const G4Fragment& aFragment, Model model) {
88  fill(aFragment.GetMomentum()/GeV, aFragment.GetA_asInt(),
89  aFragment.GetZ_asInt(), aFragment.GetExcitationEnergy(), model);
90 
91  // Exciton configuration must be set by hand
93 
95  aFragment.GetNumberOfParticles() - aFragment.GetNumberOfCharged();
96 
98 
101 }
102 
103 
104 // FIXME: Should we have a local buffer and return by const-reference instead?
106  G4Fragment frag(getA(), getZ(), getMomentum()*GeV); // From Bertini units
107 
108  // Note: exciton configuration has to be set piece by piece
112 
116 
117  return frag;
118 }
119 
120 G4InuclNuclei::operator G4Fragment() const {
121  return makeG4Fragment();
122 }
123 
124 
125 // Convert contents from (via constructor) G4V3DNucleus
126 
129  : G4InuclParticle() {
130  copy(a3DNucleus, model);
131 }
132 
134  if (!a3DNucleus) return; // Null pointer means no action
135 
136  fill(0., a3DNucleus->GetMassNumber(), a3DNucleus->GetCharge(), 0., model);
137 
138  // Convert every hit nucleon into an exciton hole
139  if (a3DNucleus->StartLoop()) {
140  G4Nucleon* nucl = 0;
141 
142  /* Loop checking 08.06.2015 MHK */
143  while ((nucl = a3DNucleus->GetNextNucleon())) {
144  if (nucl->AreYouHit()) { // Found previously interacted nucleon
145  if (nucl->GetParticleType() == G4Proton::Definition())
147 
150  }
151  }
152  }
153 }
154 
155 
156 // Overwrite data structure (avoids creating/copying temporaries)
157 
161  setMomentum(mom);
162  setExitationEnergy(exc);
164  setModel(model);
165 }
166 
170  setKineticEnergy(ekin);
171  setExitationEnergy(exc);
173  setModel(model);
174 }
175 
177  setDefinition(0);
180 }
181 
182 
183 // Change excitation energy while keeping momentum vector constant
184 
186  G4double ekin = getKineticEnergy(); // Current kinetic energy
187 
188  G4double emass = getNucleiMass() + e*MeV/GeV; // From Bertini to G4 units
189 
190  // Safety check -- if zero energy, don't do computation
191  G4double ekin_new = (ekin == 0.) ? 0.
192  : std::sqrt(emass*emass + ekin*(2.*getMass()+ekin)) - emass;
193 
194  setMass(emass); // Momentum is computed from mass and Ekin
195  setKineticEnergy(ekin_new);
196 }
197 
198 
199 // Convert nuclear configuration to standard GEANT4 pointer
200 
201 // WARNING: Opposite conventions! G4InuclNuclei uses (A,Z) everywhere, while
202 // G4ParticleTable::GetIon() uses (Z,A)!
203 
205  // SPECIAL CASE: (0,0) means create dummy without definition
206  if (0 == a && 0 == z) return 0;
207 
209  G4ParticleDefinition *pd = pTable->GetIonTable()->GetIon(z, a, 0);
210 
211  // SPECIAL CASE: Non-physical nuclear fragment, for final-state return
212  if (!pd) pd = makeNuclearFragment(a,z);
213 
214  return pd; // This could return a null pointer if above fails
215 }
216 
217 
218 // Shared buffer of nuclear fragments created below, to avoid memory leaks
219 
220 namespace {
221  static std::map<G4int,G4ParticleDefinition*> fragmentList;
222  G4Mutex fragListMutex = G4MUTEX_INITIALIZER;
223 }
224 
225 // Creates a non-physical pseudo-nucleus, for return as final-state fragment
226 // from G4IntraNuclearCascader
227 
230  if (a<=0 || z<0 || a<z) {
231  G4cerr << " >>> G4InuclNuclei::makeNuclearFragment() called with"
232  << " impossible arguments A=" << a << " Z=" << z << G4endl;
233  throw G4HadronicException(__FILE__, __LINE__,
234  "G4InuclNuclei impossible A/Z arguments");
235  }
236 
238 
239  // Use local lookup table (see above) to maintain singletons
240  // NOTE: G4ParticleDefinitions don't need to be explicitly deleted
241  // (see comments in G4IonTable.cc::~G4IonTable)
242 
243  G4AutoLock fragListLock(&fragListMutex);
244  if (fragmentList.find(code) != fragmentList.end()) return fragmentList[code];
245  fragListLock.unlock();
246 
247  // Name string follows format in G4IonTable.cc::GetIonName(Z,A,E)
248  std::stringstream zstr, astr;
249  zstr << z;
250  astr << a;
251 
252  G4String name = "Z" + zstr.str() + "A" + astr.str();
253 
254  G4double mass = getNucleiMass(a,z) *GeV/MeV; // From Bertini to GEANT4 units
255 
256  // Arguments for constructor are as follows
257  // name mass width charge
258  // 2*spin parity C-conjugation
259  // 2*Isospin 2*Isospin3 G-parity
260  // type lepton number baryon number PDG encoding
261  // stable lifetime decay table
262  // shortlived subType anti_encoding Excitation-energy
263 
264  G4Ions* fragPD = new G4Ions(name, mass, 0., z*eplus,
265  0, +1, 0,
266  0, 0, 0,
267  "nucleus", 0, a, code,
268  true, 0., 0,
269  true, "generic", 0, 0.);
270  fragPD->SetAntiPDGEncoding(0);
271 
272  fragListLock.lock(); // Protect before saving new fragment
273  return (fragmentList[code] = fragPD); // Store in table for next lookup
274 }
275 
277  // Simple minded mass calculation use constants in CLHEP (all in MeV)
279 
280  return mass*MeV/GeV; // Convert from GEANT4 to Bertini units
281 }
282 
283 // Assignment operator for use with std::sort()
285  if (this != &right) {
288  }
289  return *this;
290 }
291 
292 // Dump particle properties for diagnostics
293 
294 void G4InuclNuclei::print(std::ostream& os) const {
296  os << G4endl << " Nucleus: " << getDefinition()->GetParticleName()
297  << " A " << getA() << " Z " << getZ() << " mass " << getMass()
298  << " Eex (MeV) " << getExitationEnergy();
299 
301  os << G4endl << " " << theExitonConfiguration;
302 }