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