Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4InuclCollider.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: G4InuclCollider.cc 71769 2013-06-21 21:23:50Z mkelsey $
27 //
28 // 20100114 M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
29 // 20100309 M. Kelsey -- Eliminate some unnecessary std::pow()
30 // 20100413 M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
31 // 20100418 M. Kelsey -- Move lab-frame transformation code to G4CollisonOutput
32 // 20100429 M. Kelsey -- Change "photon()" to "isPhoton()"
33 // 20100517 M. Kelsey -- Inherit from common base class, make other colliders
34 // simple data members, consolidate code
35 // 20100620 M. Kelsey -- Reorganize top level if-blocks to reduce nesting,
36 // use new four-vector conservation check.
37 // 20100701 M. Kelsey -- Bug fix energy-conservation after equilibrium evap,
38 // pass verbosity through to G4CollisionOutput
39 // 20100714 M. Kelsey -- Move conservation checking to base class, report
40 // number of iterations at end
41 // 20100715 M. Kelsey -- Remove all the "Before xxx" and "After xxx"
42 // conservation checks, as the colliders now all do so. Move
43 // local buffers outside while() loop, use new "combined add()"
44 // interface for copying local buffers to global.
45 // 20100716 M. Kelsey -- Drop G4IntraNucleiCascader::setInteractionCase()
46 // 20100720 M. Kelsey -- Make all the collders pointer members (to reducde
47 // external compile dependences).
48 // 20100915 M. Kelsey -- Move post-cascade colliders to G4CascadeDeexcitation,
49 // simplify operational code somewhat
50 // 20100922 M. Kelsey -- Add functions to select de-excitation method;
51 // default is G4CascadeDeexcitation (i.e., built-in modules)
52 // 20100924 M. Kelsey -- Migrate to integer A and Z
53 // 20101019 M. Kelsey -- CoVerity report: check dynamic_cast<> for null
54 // 20110224 M. Kelsey -- Add ::rescatter() function which takes a list of
55 // pre-existing secondaries as input. Add setVerboseLevel().
56 // 20110301 M. Kelsey -- Pass verbosity to new or changed de-excitation
57 // 20110304 M. Kelsey -- Modify rescatter to use original Propagate() input
58 // 20110308 M. Kelsey -- Separate de-excitation block from collide(); check
59 // for single-nucleon "fragment", rather than for null fragment
60 // 20110413 M. Kelsey -- Modify diagnostic messages in ::rescatter() to be
61 // equivalent to those from ::collide().
62 // 20111003 M. Kelsey -- Prepare for gamma-N interactions by checking for
63 // final-state tables instead of particle "isPhoton()"
64 // 20130621 M. Kelsey -- Pass G4Fragment to de-excitation modules directly
65 // 20140929 M. Kelsey -- Make PreCompound the default de-excitation
66 // 20141111 M. Kelsey -- Revert default use of PreCompound; replace
67 // G4Fragment::GetA() call with GetA_asInt().
68 // 20150205 M. Kelsey -- New photonuclearOkay() filter to reject events
69 // around giant dipole resonance with no hadronic secondaries.
70 // Addresses bug #1680.
71 // 20150220 M. Kelsey -- Improve photonuclearOkay() filter by just checking
72 // final-state nucleus vs. target, rather than all secondaries.
73 // 20150608 M. Kelsey -- Label all while loops as terminating.
74 
75 #include "G4InuclCollider.hh"
77 #include "G4CascadeCheckBalance.hh"
78 #include "G4CascadeDeexcitation.hh"
79 #include "G4CollisionOutput.hh"
81 #include "G4IntraNucleiCascader.hh"
83 #include "G4InuclNuclei.hh"
84 #include "G4LorentzConvertor.hh"
86 
87 
89  : G4CascadeColliderBase("G4InuclCollider"),
90  theElementaryParticleCollider(new G4ElementaryParticleCollider),
91  theIntraNucleiCascader(new G4IntraNucleiCascader),
92  theDeexcitation(new G4PreCompoundDeexcitation) {}
93 
95  delete theElementaryParticleCollider;
96  delete theIntraNucleiCascader;
97  delete theDeexcitation;
98 }
99 
100 
101 // Set verbosity and pass on to member objects
104 
105  theElementaryParticleCollider->setVerboseLevel(verboseLevel);
106  theIntraNucleiCascader->setVerboseLevel(verboseLevel);
107  theDeexcitation->setVerboseLevel(verboseLevel);
108 
110  DEXoutput.setVerboseLevel(verboseLevel);
111 }
112 
113 
114 // Select post-cascade processing (default will be CascadeDeexcitation)
115 
117  delete theDeexcitation;
118  theDeexcitation = new G4CascadeDeexcitation;
119  theDeexcitation->setVerboseLevel(verboseLevel);
120 }
121 
123  delete theDeexcitation;
124  theDeexcitation = new G4PreCompoundDeexcitation;
125  theDeexcitation->setVerboseLevel(verboseLevel);
126 }
127 
128 
129 // Main action
130 
132  G4CollisionOutput& globalOutput) {
133  if (verboseLevel) G4cout << " >>> G4InuclCollider::collide" << G4endl;
134 
135  const G4int itry_max = 100;
136 
137  // Particle-on-particle collision; no nucleus involved
138  if (useEPCollider(bullet,target)) {
139  if (verboseLevel > 2)
140  G4cout << " InuclCollider -> particle on particle collision" << G4endl;
141 
142  theElementaryParticleCollider->collide(bullet, target, globalOutput);
143  return;
144  }
145 
146  interCase.set(bullet,target); // Classify collision type
147  if (verboseLevel > 2) {
148  G4cout << " InuclCollider -> inter case " << interCase.code() << G4endl;
149  }
150 
151  if (!interCase.valid()) {
152  if (verboseLevel > 1)
153  G4cerr << " InuclCollider -> no collision possible " << G4endl;
154 
155  globalOutput.trivialise(bullet, target);
156  return;
157  }
158 
159  // Target must be a nucleus
160  G4InuclNuclei* ntarget = dynamic_cast<G4InuclNuclei*>(interCase.getTarget());
161  if (!ntarget) {
162  G4cerr << " InuclCollider -> ERROR target is not a nucleus " << G4endl;
163 
164  globalOutput.trivialise(bullet, target);
165  return;
166  }
167 
168  G4int btype = 0;
169  G4int ab = 0;
170  G4int zb = 0;
171 
172  if (interCase.hadNucleus()) { // particle with nuclei
173  G4InuclElementaryParticle* pbullet =
175 
176  if (!pbullet) {
177  G4cerr << " InuclCollider -> ERROR bullet is not a hadron " << G4endl;
178  globalOutput.trivialise(bullet, target);
179  return;
180  }
181 
182  if (!G4CascadeChannelTables::GetTable(pbullet->type())) {
183  G4cerr << " InuclCollider -> ERROR can not collide with "
184  << pbullet->getDefinition()->GetParticleName() << G4endl;
185  globalOutput.trivialise(bullet, target);
186  return;
187  }
188 
189  btype = pbullet->type();
190  } else { // nuclei with nuclei
191  G4InuclNuclei* nbullet =
192  dynamic_cast<G4InuclNuclei*>(interCase.getBullet());
193  if (!nbullet) {
194  G4cerr << " InuclCollider -> ERROR bullet is not a nucleus " << G4endl;
195  globalOutput.trivialise(bullet, target);
196  return;
197  }
198 
199  ab = nbullet->getA();
200  zb = nbullet->getZ();
201  }
202 
203  G4LorentzConvertor convertToTargetRestFrame(bullet, ntarget);
204  G4double ekin = convertToTargetRestFrame.getKinEnergyInTheTRS();
205 
206  if (verboseLevel > 3) G4cout << " ekin in trs " << ekin << G4endl;
207 
208  if (!inelasticInteractionPossible(bullet, target, ekin)) {
209  if (verboseLevel > 3)
210  G4cout << " InuclCollider -> inelastic interaction is impossible\n"
211  << " due to the coulomb barirer " << G4endl;
212 
213  globalOutput.trivialise(bullet, target);
214  return;
215  }
216 
217  // Generate interaction secondaries in rest frame of target nucleus
218  convertToTargetRestFrame.toTheTargetRestFrame();
219  if (verboseLevel > 3) {
220  G4cout << " degenerated? " << convertToTargetRestFrame.trivial()
221  << G4endl;
222  }
223 
224  G4LorentzVector bmom; // Bullet is along local Z
225  bmom.setZ(convertToTargetRestFrame.getTRSMomentum());
226 
227  // Need to make copy of bullet with momentum realigned
228  G4InuclParticle* zbullet = 0;
229  if (interCase.hadNucleus())
230  zbullet = new G4InuclElementaryParticle(bmom, btype);
231  else
232  zbullet = new G4InuclNuclei(bmom, ab, zb);
233 
234  G4int itry = 0;
235  while (itry < itry_max) { /* Loop checking 08.06.2015 MHK */
236  itry++;
237  if (verboseLevel > 2)
238  G4cout << " InuclCollider itry " << itry << G4endl;
239 
240  globalOutput.reset(); // Clear buffers for this attempt
241  output.reset();
242 
243  theIntraNucleiCascader->collide(zbullet, target, output);
244 
245  if (verboseLevel > 1) G4cout << " After Cascade " << G4endl;
246 
247  deexcite(output.getRecoilFragment(), output);
248  output.removeRecoilFragment();
249 
250  //*** TEMPORARY, USE ENVVAR TO ENABLE/DISABLE THIS TEST ***
251  if (getenv("G4CASCADE_CHECK_PHOTONUCLEAR"))
252  if (!photonuclearOkay(output)) continue;
253 
254  if (verboseLevel > 2)
255  G4cout << " itry " << itry << " finished, moving to lab frame" << G4endl;
256 
257  // convert to the LAB frame and add to final result
258  output.boostToLabFrame(convertToTargetRestFrame);
259 
260  globalOutput.add(output);
261 
262  // Adjust final state particles to balance momentum and energy
263  // FIXME: This should no longer be necessary!
264  globalOutput.setOnShell(bullet, target);
265  if (globalOutput.acceptable()) {
266  if (verboseLevel)
267  G4cout << " InuclCollider output after trials " << itry << G4endl;
268  delete zbullet;
269  return;
270  } else {
271  if (verboseLevel>2)
272  G4cerr << " InuclCollider setOnShell failed." << G4endl;
273  }
274  } // while (itry < itry_max)
275 
276  if (verboseLevel) {
277  G4cout << " InuclCollider -> can not generate acceptable inter. after "
278  << itry_max << " attempts " << G4endl;
279  }
280 
281  globalOutput.trivialise(bullet, target);
282 
283  delete zbullet;
284  return;
285 }
286 
287 
288 // For use with Propagate to preload a set of secondaries
289 
291  G4KineticTrackVector* theSecondaries,
292  G4V3DNucleus* theNucleus,
293  G4CollisionOutput& globalOutput) {
294  if (verboseLevel) G4cout << " >>> G4InuclCollider::rescatter" << G4endl;
295 
296  G4int itry=1; // For diagnostic post-processing only
297  if (verboseLevel > 2) G4cout << " InuclCollider itry " << itry << G4endl;
298 
299  globalOutput.reset(); // Clear buffers for this attempt
300  output.reset();
301 
302  theIntraNucleiCascader->rescatter(bullet, theSecondaries, theNucleus,
303  output);
304 
305  if (verboseLevel > 1) G4cout << " After Rescatter" << G4endl;
306 
307  deexcite(output.getRecoilFragment(), output);
308  output.removeRecoilFragment();
309 
310  globalOutput.add(output); // Add local results to global output
311 
312  if (verboseLevel)
313  G4cout << " InuclCollider output after trials " << itry << G4endl;
314 }
315 
316 
317 // De-excite nuclear fragment to ground state
318 
320  G4CollisionOutput& globalOutput) {
321  if (fragment.GetA_asInt() <= 1) return; // Nothing real to be de-excited
322 
323  if (verboseLevel) G4cout << " >>> G4InuclCollider::deexcite" << G4endl;
324 
325  const G4int itry_max = 10; // Maximum number of attempts
326  G4int itry = 0;
327  do { /* Loop checking 08.06.2015 MHK */
328  if (verboseLevel > 2) G4cout << " deexcite itry " << itry << G4endl;
329 
330  DEXoutput.reset();
331  theDeexcitation->deExcite(fragment, DEXoutput);
332  } while (!validateOutput(fragment, DEXoutput) && (++itry < itry_max));
333 
334  // Add de-excitation products to output buffer
335  globalOutput.add(DEXoutput);
336 }
337 
338 
339 // Looks for non-gamma final state in photonuclear or leptonuclear
340 
342  if (interCase.twoNuclei()) return true; // A-A is not photonuclear
343 
344  G4InuclElementaryParticle* bullet =
346  if (!bullet || !(bullet->isPhoton() || bullet->isElectron())) return true;
347 
348  if (verboseLevel>1)
349  G4cout << " >>> G4InuclCollider::photonuclearOkay" << G4endl;
350 
351  if (bullet->getKineticEnergy() > 0.050) return true;
352 
353  if (verboseLevel>2) {
354  if (checkOutput.numberOfOutgoingNuclei() > 0) {
355  G4cout << " comparing final nucleus with initial target:\n"
356  << checkOutput.getOutgoingNuclei()[0] << G4endl
357  << *(interCase.getTarget()) << G4endl;
358  } else {
359  G4cout << " no final nucleus remains when target was "
360  << *(interCase.getTarget()) << G4endl;
361  }
362  }
363 
364  // Hadron production changes target nucleus
365  G4double mfinalNuc = 0.0;
366  if (checkOutput.numberOfOutgoingNuclei() > 0)
367  mfinalNuc = checkOutput.getOutgoingNuclei()[0].getMass();
368  G4double mtargetNuc = interCase.getTarget()->getMass();
369  if (mfinalNuc != mtargetNuc) return true; // Mass from G4Ions is fixed
370 
371  if (verboseLevel>2)
372  G4cout << " photonuclear produced only gammas. Try again." << G4endl;
373 
374  return false; // Final state is entirely de-excitation photons
375 }
G4bool hadNucleus() const
void trivialise(G4InuclParticle *bullet, G4InuclParticle *target)
void rescatter(G4InuclParticle *bullet, G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus, G4CollisionOutput &globalOutput)
void setVerboseLevel(G4int verbose)
static const G4CascadeChannel * GetTable(G4int initialState)
G4int getZ() const
G4double getTRSMomentum() const
const XML_Char * target
Definition: expat.h:268
virtual G4bool useEPCollider(G4InuclParticle *bullet, G4InuclParticle *target) const
void removeRecoilFragment(G4int index=-1)
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
G4bool valid() const
const G4ParticleDefinition * getDefinition() const
G4int code() const
virtual G4bool inelasticInteractionPossible(G4InuclParticle *bullet, G4InuclParticle *target, G4double ekin) const
virtual ~G4InuclCollider()
G4bool photonuclearOkay(G4CollisionOutput &checkOutput) const
void deexcite(const G4Fragment &fragment, G4CollisionOutput &globalOutput)
int G4int
Definition: G4Types.hh:78
G4bool trivial() const
const G4String & GetParticleName() const
G4bool acceptable() const
virtual void setVerboseLevel(G4int verbose=0)
G4double getKineticEnergy() const
virtual void deExcite(const G4Fragment &fragment, G4CollisionOutput &output)=0
void add(const G4CollisionOutput &right)
void useCascadeDeexcitation()
G4GLOB_DLL std::ostream G4cout
G4int GetA_asInt() const
Definition: G4Fragment.hh:266
G4int getA() const
bool G4bool
Definition: G4Types.hh:79
void setVerboseLevel(G4int verbose=0)
void rescatter(G4InuclParticle *bullet, G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus, G4CollisionOutput &globalOutput)
void boostToLabFrame(const G4LorentzConvertor &convertor)
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &globalOutput)
G4double getKinEnergyInTheTRS() const
G4bool twoNuclei() const
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &globalOutput)
G4int numberOfOutgoingNuclei() const
void set(G4InuclParticle *part1, G4InuclParticle *part2)
static const G4double ab
virtual G4bool validateOutput(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
const std::vector< G4InuclNuclei > & getOutgoingNuclei() const
G4InuclParticle * getBullet() const
#define G4endl
Definition: G4ios.hh:61
virtual void setVerboseLevel(G4int verbose=0)
void setVerboseLevel(G4int verbose=0)
const G4Fragment & getRecoilFragment(G4int index=0) const
double G4double
Definition: G4Types.hh:76
G4InuclParticle * getTarget() const
void usePreCompoundDeexcitation()
void setOnShell(G4InuclParticle *bullet, G4InuclParticle *target)
G4GLOB_DLL std::ostream G4cerr
G4double getMass() const