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