Geant4_10
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 
66 #include "G4InuclCollider.hh"
68 #include "G4CascadeCheckBalance.hh"
69 #include "G4CascadeDeexcitation.hh"
70 #include "G4CollisionOutput.hh"
72 #include "G4IntraNucleiCascader.hh"
74 #include "G4InuclNuclei.hh"
75 #include "G4LorentzConvertor.hh"
77 
78 
80  : G4CascadeColliderBase("G4InuclCollider"),
81  theElementaryParticleCollider(new G4ElementaryParticleCollider),
82  theIntraNucleiCascader(new G4IntraNucleiCascader),
83  theDeexcitation(new G4CascadeDeexcitation) {}
84 
86  delete theElementaryParticleCollider;
87  delete theIntraNucleiCascader;
88  delete theDeexcitation;
89 }
90 
91 
92 // Set verbosity and pass on to member objects
95 
96  theElementaryParticleCollider->setVerboseLevel(verboseLevel);
97  theIntraNucleiCascader->setVerboseLevel(verboseLevel);
98  theDeexcitation->setVerboseLevel(verboseLevel);
99 
101  DEXoutput.setVerboseLevel(verboseLevel);
102 }
103 
104 
105 // Select post-cascade processing (default will be CascadeDeexcitation)
106 
108  delete theDeexcitation;
109  theDeexcitation = new G4CascadeDeexcitation;
110  theDeexcitation->setVerboseLevel(verboseLevel);
111 }
112 
114  delete theDeexcitation;
115  theDeexcitation = new G4PreCompoundDeexcitation;
116  theDeexcitation->setVerboseLevel(verboseLevel);
117 }
118 
119 
120 // Main action
121 
123  G4CollisionOutput& globalOutput) {
124  if (verboseLevel) G4cout << " >>> G4InuclCollider::collide" << G4endl;
125 
126  const G4int itry_max = 100;
127 
128  // Particle-on-particle collision; no nucleus involved
129  if (useEPCollider(bullet,target)) {
130  if (verboseLevel > 2)
131  G4cout << " InuclCollider -> particle on particle collision" << G4endl;
132 
133  theElementaryParticleCollider->collide(bullet, target, globalOutput);
134  return;
135  }
136 
137  interCase.set(bullet,target); // Classify collision type
138  if (verboseLevel > 2) {
139  G4cout << " InuclCollider -> inter case " << interCase.code() << G4endl;
140  }
141 
142  if (!interCase.valid()) {
143  if (verboseLevel > 1)
144  G4cerr << " InuclCollider -> no collision possible " << G4endl;
145 
146  globalOutput.trivialise(bullet, target);
147  return;
148  }
149 
150  // Target must be a nucleus
151  G4InuclNuclei* ntarget = dynamic_cast<G4InuclNuclei*>(interCase.getTarget());
152  if (!ntarget) {
153  G4cerr << " InuclCollider -> ERROR target is not a nucleus " << G4endl;
154 
155  globalOutput.trivialise(bullet, target);
156  return;
157  }
158 
159  G4int btype = 0;
160  G4int ab = 0;
161  G4int zb = 0;
162 
163  if (interCase.hadNucleus()) { // particle with nuclei
164  G4InuclElementaryParticle* pbullet =
166 
167  if (!pbullet) {
168  G4cerr << " InuclCollider -> ERROR bullet is not a hadron " << G4endl;
169  globalOutput.trivialise(bullet, target);
170  return;
171  }
172 
173  if (!G4CascadeChannelTables::GetTable(pbullet->type())) {
174  G4cerr << " InuclCollider -> ERROR can not collide with "
175  << pbullet->getDefinition()->GetParticleName() << G4endl;
176  globalOutput.trivialise(bullet, target);
177  return;
178  }
179 
180  btype = pbullet->type();
181  } else { // nuclei with nuclei
182  G4InuclNuclei* nbullet =
183  dynamic_cast<G4InuclNuclei*>(interCase.getBullet());
184  if (!nbullet) {
185  G4cerr << " InuclCollider -> ERROR bullet is not a nucleus " << G4endl;
186  globalOutput.trivialise(bullet, target);
187  return;
188  }
189 
190  ab = nbullet->getA();
191  zb = nbullet->getZ();
192  }
193 
194  G4LorentzConvertor convertToTargetRestFrame(bullet, ntarget);
195  G4double ekin = convertToTargetRestFrame.getKinEnergyInTheTRS();
196 
197  if (verboseLevel > 3) G4cout << " ekin in trs " << ekin << G4endl;
198 
199  if (!inelasticInteractionPossible(bullet, target, ekin)) {
200  if (verboseLevel > 3)
201  G4cout << " InuclCollider -> inelastic interaction is impossible\n"
202  << " due to the coulomb barirer " << G4endl;
203 
204  globalOutput.trivialise(bullet, target);
205  return;
206  }
207 
208  // Generate interaction secondaries in rest frame of target nucleus
209  convertToTargetRestFrame.toTheTargetRestFrame();
210  if (verboseLevel > 3) {
211  G4cout << " degenerated? " << convertToTargetRestFrame.trivial()
212  << G4endl;
213  }
214 
215  G4LorentzVector bmom; // Bullet is along local Z
216  bmom.setZ(convertToTargetRestFrame.getTRSMomentum());
217 
218  // Need to make copy of bullet with momentum realigned
219  G4InuclParticle* zbullet = 0;
220  if (interCase.hadNucleus())
221  zbullet = new G4InuclElementaryParticle(bmom, btype);
222  else
223  zbullet = new G4InuclNuclei(bmom, ab, zb);
224 
225  G4int itry = 0;
226  while (itry < itry_max) {
227  itry++;
228  if (verboseLevel > 2)
229  G4cout << " InuclCollider itry " << itry << G4endl;
230 
231  globalOutput.reset(); // Clear buffers for this attempt
232  output.reset();
233 
234  theIntraNucleiCascader->collide(zbullet, target, output);
235 
236  if (verboseLevel > 1) G4cout << " After Cascade " << G4endl;
237 
238  deexcite(output.getRecoilFragment(), output);
239  output.removeRecoilFragment();
240 
241  if (verboseLevel > 2)
242  G4cout << " itry " << itry << " finished, moving to lab frame" << G4endl;
243 
244  // convert to the LAB frame and add to final result
245  output.boostToLabFrame(convertToTargetRestFrame);
246 
247  globalOutput.add(output);
248 
249  // Adjust final state particles to balance momentum and energy
250  // FIXME: This should no longer be necessary!
251  globalOutput.setOnShell(bullet, target);
252  if (globalOutput.acceptable()) {
253  if (verboseLevel)
254  G4cout << " InuclCollider output after trials " << itry << G4endl;
255  delete zbullet;
256  return;
257  } else {
258  if (verboseLevel>2)
259  G4cerr << " InuclCollider setOnShell failed." << G4endl;
260  }
261  } // while (itry < itry_max)
262 
263  if (verboseLevel) {
264  G4cout << " InuclCollider -> can not generate acceptable inter. after "
265  << itry_max << " attempts " << G4endl;
266  }
267 
268  globalOutput.trivialise(bullet, target);
269 
270  delete zbullet;
271  return;
272 }
273 
274 
275 // For use with Propagate to preload a set of secondaries
276 
278  G4KineticTrackVector* theSecondaries,
279  G4V3DNucleus* theNucleus,
280  G4CollisionOutput& globalOutput) {
281  if (verboseLevel) G4cout << " >>> G4InuclCollider::rescatter" << G4endl;
282 
283  G4int itry=1; // For diagnostic post-processing only
284  if (verboseLevel > 2) G4cout << " InuclCollider itry " << itry << G4endl;
285 
286  globalOutput.reset(); // Clear buffers for this attempt
287  output.reset();
288 
289  theIntraNucleiCascader->rescatter(bullet, theSecondaries, theNucleus,
290  output);
291 
292  if (verboseLevel > 1) G4cout << " After Rescatter" << G4endl;
293 
294  deexcite(output.getRecoilFragment(), output);
295  output.removeRecoilFragment();
296 
297  globalOutput.add(output); // Add local results to global output
298 
299  if (verboseLevel)
300  G4cout << " InuclCollider output after trials " << itry << G4endl;
301 }
302 
303 
304 // De-excite nuclear fragment to ground state
305 
307  G4CollisionOutput& globalOutput) {
308  if (fragment.GetA() <= 1) return; // Nothing real to be de-excited
309 
310  if (verboseLevel) G4cout << " >>> G4InuclCollider::deexcite" << G4endl;
311 
312  const G4int itry_max = 10; // Maximum number of attempts
313  G4int itry = 0;
314  do {
315  if (verboseLevel > 2) G4cout << " deexcite itry " << itry << G4endl;
316 
317  DEXoutput.reset();
318  theDeexcitation->deExcite(fragment, DEXoutput);
319  } while (!validateOutput(fragment, DEXoutput) && (++itry < itry_max));
320 
321  // Add de-excitation products to output buffer
322  globalOutput.add(DEXoutput);
323 }
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
virtual G4bool useEPCollider(G4InuclParticle *bullet, G4InuclParticle *target) const
G4double GetA() const
Definition: G4Fragment.hh:283
void removeRecoilFragment(G4int index=-1)
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
G4bool valid() const
G4int code() const
virtual G4bool inelasticInteractionPossible(G4InuclParticle *bullet, G4InuclParticle *target, G4double ekin) const
virtual ~G4InuclCollider()
void deexcite(const G4Fragment &fragment, G4CollisionOutput &globalOutput)
const XML_Char * target
Definition: expat.h:268
int G4int
Definition: G4Types.hh:78
G4bool trivial() const
const G4String & GetParticleName() const
G4bool acceptable() const
virtual void setVerboseLevel(G4int verbose=0)
virtual void deExcite(const G4Fragment &fragment, G4CollisionOutput &output)=0
void add(const G4CollisionOutput &right)
void useCascadeDeexcitation()
G4GLOB_DLL std::ostream G4cout
G4int getA() const
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
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &globalOutput)
void set(G4InuclParticle *part1, G4InuclParticle *part2)
virtual G4bool validateOutput(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
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
G4ParticleDefinition * getDefinition() const
void usePreCompoundDeexcitation()
void setOnShell(G4InuclParticle *bullet, G4InuclParticle *target)
G4GLOB_DLL std::ostream G4cerr