Geant4  10.01.p02
G4IntraNucleiCascader.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: G4IntraNucleiCascader.cc 67746 2013-03-05 21:11:14Z mkelsey $
27 //
28 // 20100114 M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
29 // 20100307 M. Kelsey -- Bug fix: momentum_out[0] should be momentum_out.e()
30 // 20100309 M. Kelsey -- Eliminate some unnecessary std::pow()
31 // 20100407 M. Kelsey -- Pass "all_particles" as argument to initializeCascad,
32 // following recent change to G4NucleiModel.
33 // 20100413 M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
34 // 20100517 M. Kelsey -- Inherit from common base class, make other colliders
35 // simple data members
36 // 20100616 M. Kelsey -- Add reporting of final residual particle
37 // 20100617 M. Kelsey -- Remove "RUN" preprocessor flag and all "#else" code,
38 // pass verbosity to collider. Make G4NucleiModel a data member,
39 // instead of creating and deleting on every cycle.
40 // 20100620 M. Kelsey -- Improved diagnostic messages. Simplify kinematics
41 // of recoil nucleus.
42 // 20100622 M. Kelsey -- Use local "bindingEnergy()" to call through.
43 // 20100623 M. Kelsey -- Undo G4NucleiModel change from 0617. Does not work
44 // properly across multiple interactions.
45 // 20100627 M. Kelsey -- Protect recoil nucleus energy from floating roundoff
46 // by setting small +ve or -ve values to zero.
47 // 20100701 M. Kelsey -- Let excitation energy be handled by G4InuclNuclei,
48 // allow for ground-state recoil (goodCase == true for Eex==0.)
49 // 20100702 M. Kelsey -- Negative energy recoil should be rejected
50 // 20100706 D. Wright -- Copy "abandoned" cparticles to output list, copy
51 // mesonic "excitons" to output list; should be absorbed, fix up
52 // diagnostic messages.
53 // 20100713 M. Kelsey -- Add more diagnostics for Dennis' changes.
54 // 20100714 M. Kelsey -- Switch to new G4CascadeColliderBase class, remove
55 // sanity check on afin/zfin (not valid).
56 // 20100715 M. Kelsey -- Add diagnostic for ekin_in vs. actual ekin; reduce
57 // KE across Coulomb barrier. Rearrange end-of-loop if blocks,
58 // add conservation check at end.
59 // 20100716 M. Kelsey -- Eliminate inter_case; use base-class functionality.
60 // Add minimum-fragment requirement for recoil, in order to
61 // allow for momentum balancing
62 // 20100720 M. Kelsey -- Make EPCollider pointer member
63 // 20100721 M. Kelsey -- Turn on conservation checks unconditionally (override
64 // new G4CASCADE_CHECK_ECONS setting
65 // 20100722 M. Kelsey -- Move cascade output buffers to .hh file
66 // 20100728 M. Kelsey -- Make G4NucleiModel data member for persistence,
67 // delete colliders in destructor
68 // 20100906 M. Kelsey -- Hide "non-physical fragment" behind verbose flag
69 // 20100907 M. Kelsey -- Add makeResidualFragment function to create object
70 // 20100909 M. Kelsey -- Remove all local "fragment" stuff, use RecoilMaker.
71 // move goodCase() to RecoilMaker.
72 // 20100910 M. Kelsey -- Use RecoilMaker::makeRecoilFragment().
73 // 20100915 M. Kelsey -- Define functions to deal with trapped particles,
74 // move the exciton container to a data member
75 // 20100916 M. Kelsey -- Put decay photons directly onto output list
76 // 20100921 M. Kelsey -- Migrate to RecoilMaker::makeRecoilNuclei().
77 // 20100924 M. Kelsey -- Minor shuffling of post-cascade recoil building.
78 // Create G4Fragment for recoil and store in output.
79 // 20110131 M. Kelsey -- Move "momentum_in" calculation inside verbosity
80 // 20110214 M. Kelsey -- Follow G4InuclParticle::Model enumerator migration
81 // 20110224 M. Kelsey -- Add ::rescatter() function which takes a list of
82 // pre-existing secondaries as input. Split ::collide() into
83 // separate utility functions. Move cascade parameters to static
84 // data members. Add setVerboseLevel().
85 // 20110302 M. Kelsey -- Move G4NucleiModel::printModel() call to G4NucleiModel
86 // 20110303 M. Kelsey -- Add more cascade functions to support rescattering
87 // 20110304 M. Kelsey -- Get original Propagate() arguments here in rescatter()
88 // and convert to particles, nuclei and G4NucleiModel state.
89 // 20110308 M. Kelsey -- Don't put recoiling fragment onto output list any more
90 // 20110308 M. Kelsey -- Decay unstable hadrons from pre-cascade, use daughters
91 // 20110324 M. Kelsey -- Get locations of hit nuclei in ::rescatter(), pass
92 // to G4NucleiModel::reset().
93 // 20110404 M. Kelsey -- Reduce maximum number of retries to 100, reflection
94 // cut to 50.
95 // 20110721 M. Kelsey -- Put unusable pre-cascade particles directly on output,
96 // do not decay.
97 // 20110722 M. Kelsey -- Deprecate "output_particles" list in favor of using
98 // output directly (will help with pre-cascade issues).
99 // 20110801 M. Kelsey -- Use G4Inucl(Particle)::fill() functions to reduce
100 // creation of temporaries. Add local target buffer for
101 // rescattering, to avoid memory leak.
102 // 20110808 M. Kelsey -- Pass buffer to generateParticleFate() to avoid copy
103 // 20110919 M. Kelsey -- Add optional final-state clustering, controlled (for
104 // now) with compiler flag G4CASCADE_DO_COALESCENCE
105 // 20110922 M. Kelsey -- Follow migrations G4InuclParticle::print(ostream&)
106 // and G4CascadParticle::print(ostream&); drop Q,B printing
107 // 20110926 M. Kelsey -- Replace compiler flag with one-time envvar in ctor
108 // for final-state clustering.
109 // 20111003 M. Kelsey -- Prepare for gamma-N interactions by checking for
110 // final-state tables instead of particle "isPhoton()"
111 // 20120521 A. Ribon -- Specify mass when decay trapped particle.
112 // 20120822 M. Kelsey -- Move envvars to G4CascadeParameters.
113 // 20121205 M. Kelsey -- In processSecondary(), set generation to 1, as these
114 // particles are not true projectiles, but already embedded.
115 // 20130304 M. Kelsey -- Use new G4CascadeHistory to dump cascade structure
116 // 20140310 M. Kelsey -- (Bug #1584) Release memory allocated by DecayIt()
117 // 20140409 M. Kelsey -- Use const G4ParticleDefinition* everywhere
118 
119 #include <algorithm>
120 
121 #include "G4IntraNucleiCascader.hh"
122 #include "G4SystemOfUnits.hh"
123 #include "G4CascadeChannelTables.hh"
124 #include "G4CascadeCoalescence.hh"
125 #include "G4CascadeHistory.hh"
126 #include "G4CascadeParameters.hh"
127 #include "G4CascadeRecoilMaker.hh"
128 #include "G4CascadParticle.hh"
129 #include "G4CollisionOutput.hh"
130 #include "G4DecayProducts.hh"
131 #include "G4DecayTable.hh"
133 #include "G4ExitonConfiguration.hh"
134 #include "G4HadTmpUtil.hh"
136 #include "G4InuclNuclei.hh"
137 #include "G4InuclParticleNames.hh"
139 #include "G4KineticTrack.hh"
140 #include "G4KineticTrackVector.hh"
141 #include "G4LorentzConvertor.hh"
142 #include "G4Neutron.hh"
143 #include "G4NucleiModel.hh"
144 #include "G4ParticleLargerEkin.hh"
145 #include "G4Proton.hh"
146 #include "G4V3DNucleus.hh"
147 #include "Randomize.hh"
148 
149 using namespace G4InuclParticleNames;
150 using namespace G4InuclSpecialFunctions;
151 
152 
153 // Configuration parameters for cascade production
158 
159 
160 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
161 
163  : G4CascadeColliderBase("G4IntraNucleiCascader"), model(new G4NucleiModel),
164  theElementaryParticleCollider(new G4ElementaryParticleCollider),
165  theRecoilMaker(new G4CascadeRecoilMaker), theClusterMaker(0),
166  theCascadeHistory(0), tnuclei(0), bnuclei(0), bparticle(0),
167  minimum_recoil_A(0.), coulombBarrier(0.),
168  nucleusTarget(new G4InuclNuclei),
169  protonTarget(new G4InuclElementaryParticle) {
172 
175 }
176 
178  delete model;
180  delete theRecoilMaker;
181  delete theClusterMaker;
182  delete theCascadeHistory;
183  delete nucleusTarget;
184  delete protonTarget;
185 }
186 
189  model->setVerboseLevel(verbose);
192 
193  // Optional functionality
196 }
197 
198 
199 
201  G4InuclParticle* target,
202  G4CollisionOutput& globalOutput) {
203  if (verboseLevel) G4cout << " >>> G4IntraNucleiCascader::collide " << G4endl;
204 
205  if (!initialize(bullet, target)) return; // Load buffers and drivers
206 
207  G4int itry = 0;
208  do {
209  newCascade(++itry);
210  setupCascade();
211  generateCascade();
212  } while (!finishCascade() && itry<itry_max);
213 
214  // Report full structure of final cascade if requested
216 
217  finalize(itry, bullet, target, globalOutput);
218 }
219 
220 // For use with Propagate to preload a set of secondaries
221 // FIXME: So far, we don't have any bullet information from Propagate!
222 
224  G4KineticTrackVector* theSecondaries,
225  G4V3DNucleus* theNucleus,
226  G4CollisionOutput& globalOutput) {
227  if (verboseLevel)
228  G4cout << " >>> G4IntraNucleiCascader::rescatter " << G4endl;
229 
230  G4InuclParticle* target = createTarget(theNucleus);
231  if (!initialize(bullet, target)) return; // Load buffers and drivers
232 
233  G4int itry = 0;
234  do {
235  newCascade(++itry);
236  preloadCascade(theNucleus, theSecondaries);
237  generateCascade();
238  } while (!finishCascade() && itry<itry_max);
239 
240  // Report full structure of final cascade if requested
242 
243  finalize(itry, bullet, target, globalOutput);
244 }
245 
246 
248  G4InuclParticle* target) {
249  if (verboseLevel>1)
250  G4cout << " >>> G4IntraNucleiCascader::initialize " << G4endl;
251 
252  // Configure processing modules
254 
255  interCase.set(bullet,target); // Classify collision type
256 
257  if (verboseLevel > 3) {
258  G4cout << *interCase.getBullet() << G4endl
259  << *interCase.getTarget() << G4endl;
260  }
261 
262  // Bullet may be nucleus or simple particle
263  bnuclei = dynamic_cast<G4InuclNuclei*>(interCase.getBullet());
265 
266  if (!bnuclei && !bparticle) {
267  G4cerr << " G4IntraNucleiCascader: projectile is not a valid particle."
268  << G4endl;
269  return false;
270  }
271 
272  // Target _must_ be nucleus
273  tnuclei = dynamic_cast<G4InuclNuclei*>(interCase.getTarget());
274  if (!tnuclei) {
275  if (verboseLevel)
276  G4cerr << " Target is not a nucleus. Abandoning." << G4endl;
277  return false;
278  }
279 
281  coulombBarrier = 0.00126*tnuclei->getZ() / (1.+G4cbrt(tnuclei->getA()));
282 
283  // Energy/momentum conservation usually requires a recoiling nuclear fragment
284  // This cut will be increased on each "itry" if momentum could not balance.
285  minimum_recoil_A = 0.;
286 
287  if (verboseLevel > 3) {
288  G4LorentzVector momentum_in = bullet->getMomentum() + target->getMomentum();
289  G4cout << " intitial momentum E " << momentum_in.e() << " Px "
290  << momentum_in.x() << " Py " << momentum_in.y() << " Pz "
291  << momentum_in.z() << G4endl;
292  }
293 
294  return true;
295 }
296 
297 // Re-initialize buffers for new attempt at cascade
298 
300  if (verboseLevel > 1) {
301  G4cout << " IntraNucleiCascader itry " << itry << " inter_case "
302  << interCase.code() << G4endl;
303  }
304 
305  model->reset(); // Start new cascade process
306  output.reset();
307  new_cascad_particles.clear();
309  cascad_particles.clear(); // List of initial secondaries
310 
312 }
313 
314 
315 // Load initial cascade using nuclear-model calculations
316 
318  if (verboseLevel > 1)
319  G4cout << " >>> G4IntraNucleiCascader::setupCascade" << G4endl;
320 
321  if (interCase.hadNucleus()) { // particle with nuclei
322  if (verboseLevel > 3)
323  G4cout << " bparticle charge " << bparticle->getCharge()
324  << " baryon number " << bparticle->baryon() << G4endl;
325 
327  } else { // nuclei with nuclei
328  G4int ab = bnuclei->getA();
329  G4int zb = bnuclei->getZ();
330 
331  G4NucleiModel::modelLists all_particles; // Buffer to receive lists
332  model->initializeCascad(bnuclei, tnuclei, all_particles);
333 
334  cascad_particles = all_particles.first;
335  output.addOutgoingParticles(all_particles.second);
336 
337  if (cascad_particles.size() == 0) { // compound nuclei
338  G4int i;
339 
340  for (i = 0; i < ab; i++) {
341  G4int knd = i < zb ? 1 : 2;
343  };
344 
345  G4int ihn = G4int(2 * (ab-zb) * inuclRndm() + 0.5);
346  G4int ihz = G4int(2 * zb * inuclRndm() + 0.5);
347 
348  for (i = 0; i < ihn; i++) theExitonConfiguration.incrementHoles(2);
349  for (i = 0; i < ihz; i++) theExitonConfiguration.incrementHoles(1);
350  }
351  } // if (interCase ...
352 }
353 
354 
355 // Generate one possible cascade (all secondaries, etc.)
356 
358  if (verboseLevel>1) G4cout << " generateCascade " << G4endl;
359 
360  G4int iloop = 0;
361  while (!cascad_particles.empty() && !model->empty()) {
362  iloop++;
363 
364  if (verboseLevel > 2) {
365  G4cout << " Iteration " << iloop << ": Number of cparticles "
366  << cascad_particles.size() << " last one: \n"
367  << cascad_particles.back() << G4endl;
368  }
369 
370  // Record incident particle first, to get history ID
371  if (theCascadeHistory) {
373  if (verboseLevel > 2) {
374  G4cout << " active cparticle got history ID "
375  << cascad_particles.back().getHistoryId() << G4endl;
376  }
377  }
378 
382 
383  // Record interaction for later reporting (if desired)
384  if (theCascadeHistory && new_cascad_particles.size()>1)
386 
387  if (verboseLevel > 2) {
388  G4cout << " After generate fate: New particles "
389  << new_cascad_particles.size() << G4endl
390  << " Discarding last cparticle from list " << G4endl;
391  }
392 
393  cascad_particles.pop_back();
394 
395  // handle the result of a new step
396 
397  if (new_cascad_particles.size() == 1) { // last particle goes without interaction
398  const G4CascadParticle& currentCParticle = new_cascad_particles[0];
399 
400  if (model->stillInside(currentCParticle)) {
401  if (verboseLevel > 3)
402  G4cout << " particle still inside nucleus " << G4endl;
403 
404  if (currentCParticle.getNumberOfReflections() < reflection_cut &&
405  model->worthToPropagate(currentCParticle)) {
406  if (verboseLevel > 3) G4cout << " continue reflections " << G4endl;
407  cascad_particles.push_back(currentCParticle);
408  } else {
409  processTrappedParticle(currentCParticle);
410  } // reflection or exciton
411 
412  } else { // particle about to leave nucleus - check for Coulomb barrier
413  if (verboseLevel > 3) G4cout << " possible escape " << G4endl;
414 
415  const G4InuclElementaryParticle& currentParticle =
416  currentCParticle.getParticle();
417 
418  G4double KE = currentParticle.getKineticEnergy();
419  G4double mass = currentParticle.getMass();
420  G4double Q = currentParticle.getCharge();
421 
422  if (verboseLevel > 3)
423  G4cout << " KE " << KE << " barrier " << Q*coulombBarrier << G4endl;
424 
425  if (KE < Q*coulombBarrier) {
426  // Calculate barrier penetration
427  G4double CBP = 0.0;
428 
429  // if (KE > 0.0001) CBP = std::exp(-0.00126*tnuclei->getZ()*0.25*
430  // (1./KE - 1./coulombBarrier));
431  if (KE > 0.0001) CBP = std::exp(-0.0181*0.5*tnuclei->getZ()*
432  (1./KE - 1./coulombBarrier)*
433  std::sqrt(mass*(coulombBarrier-KE)) );
434 
435  if (G4UniformRand() < CBP) {
436  if (verboseLevel > 3)
437  G4cout << " tunneled\n" << currentParticle << G4endl;
438 
439  // Tunnelling through barrier leaves KE unchanged
440  output.addOutgoingParticle(currentParticle);
441  } else {
442  processTrappedParticle(currentCParticle);
443  }
444  } else {
445  output.addOutgoingParticle(currentParticle);
446 
447  if (verboseLevel > 3)
448  G4cout << " Goes out\n" << output.getOutgoingParticles().back()
449  << G4endl;
450  }
451  }
452  } else { // interaction
453  if (verboseLevel > 3)
454  G4cout << " interacted, adding new to list " << G4endl;
455 
456  cascad_particles.insert(cascad_particles.end(),
457  new_cascad_particles.begin(),
458  new_cascad_particles.end());
459 
460  std::pair<G4int, G4int> holes = model->getTypesOfNucleonsInvolved();
461  if (verboseLevel > 3)
462  G4cout << " adding new exciton holes " << holes.first << ","
463  << holes.second << G4endl;
464 
466 
467  if (holes.second > 0)
469  } // if (new_cascad_particles ...
470 
471  // Evaluate nuclear residue
474 
475  G4double aresid = theRecoilMaker->getRecoilA();
476  if (verboseLevel > 2) {
477  G4cout << " cparticles remaining " << cascad_particles.size()
478  << " nucleus (model) has "
479  << model->getNumberOfNeutrons() << " n, "
480  << model->getNumberOfProtons() << " p "
481  << " residual fragment A " << aresid << G4endl;
482  }
483 
484  if (aresid <= minimum_recoil_A) return; // Must have minimum fragment
485  } // while cascade-list and model
486 }
487 
488 
489 // Conslidate results of cascade and evaluate success
490 
492  if (verboseLevel > 1)
493  G4cout << " >>> G4IntraNucleiCascader::finishCascade ?" << G4endl;
494 
495  // Add left-over cascade particles to output
497  cascad_particles.clear();
498 
499  // Cascade is finished. Check if it's OK.
500  if (verboseLevel>3) {
501  G4cout << " G4IntraNucleiCascader finished" << G4endl;
503  }
504 
505  // Apply cluster coalesence model to produce light ions
506  if (theClusterMaker) {
509 
510  // Update recoil fragment after generating light ions
511  if (verboseLevel>3) G4cout << " Recomputing recoil fragment" << G4endl;
513  output);
514  if (verboseLevel>3) {
515  G4cout << " After cluster coalescence" << G4endl;
517  }
518  }
519 
520  // Use last created recoil fragment instead of re-constructing
521  G4int afin = theRecoilMaker->getRecoilA();
522  G4int zfin = theRecoilMaker->getRecoilZ();
523 
524  // FIXME: Should we deal with unbalanced (0,0) case before rejecting?
525 
526  // Sanity check before proceeding
528  if (verboseLevel > 1)
529  G4cerr << " Recoil nucleus is not physical: A=" << afin << " Z="
530  << zfin << G4endl;
531  return false; // Discard event and try again
532  }
533 
535 
536  if (verboseLevel > 1) {
537  G4cout << " afin " << afin << " zfin " << zfin << G4endl;
538  }
539 
540  if (afin == 0) return true; // Whole event fragmented, exit
541 
542  if (afin == 1) { // Add bare nucleon to particle list
543  G4int last_type = (zfin==1) ? 1 : 2; // proton=1, neutron=2
544 
546  G4double mres = presid.m();
547 
548  // Check for sensible kinematics
549  if (mres-mass < -small_ekin) { // Insufficient recoil energy
550  if (verboseLevel > 2) G4cerr << " unphysical recoil nucleon" << G4endl;
551  return false;
552  }
553 
554  if (mres-mass > small_ekin) { // Too much extra energy
555  if (verboseLevel > 2)
556  G4cerr << " extra energy with recoil nucleon" << G4endl;
557 
558  // FIXME: For now, we add the nucleon as unbalanced, and let
559  // "SetOnShell" fudge things. This should be abandoned.
560  }
561 
562  G4InuclElementaryParticle last_particle(presid, last_type,
564 
565  if (verboseLevel > 3) {
566  G4cout << " adding recoiling nucleon to output list\n"
567  << last_particle << G4endl;
568  }
569 
570  output.addOutgoingParticle(last_particle);
571 
572  // Update recoil to include residual nucleon
574  output);
575  }
576 
577  // Process recoil fragment for consistency, exit or reject
578  if (output.numberOfOutgoingParticles() == 1) {
580  if (std::abs(Eex) < quasielast_cut) {
581  if (verboseLevel > 3) {
582  G4cout << " quasi-elastic scatter with " << Eex << " MeV recoil"
583  << G4endl;
584  }
585 
587  if (verboseLevel > 3) {
588  G4cout << " Eex reset to " << theRecoilMaker->getRecoilExcitation()
589  << G4endl;
590  }
591  }
592  }
593 
594  if (theRecoilMaker->goodNucleus()) {
596 
598  if (!recoilFrag) {
599  G4cerr << "Got null pointer for recoil fragment!" << G4endl;
600  return false;
601  }
602 
603  if (verboseLevel > 2)
604  G4cout << " adding recoil fragment to output list" << G4endl;
605 
606  output.addRecoilFragment(*recoilFrag);
607  }
608 
609  // Put final-state particles in "leading order" for return
610  std::vector<G4InuclElementaryParticle>& opart = output.getOutgoingParticles();
611  std::sort(opart.begin(), opart.end(), G4ParticleLargerEkin());
612 
613  // Adjust final state to balance momentum and energy if necessary
618 
619  if (output.acceptable()) return true;
620  else if (verboseLevel>2) G4cerr << " Cascade setOnShell failed." << G4endl;
621  }
622 
623  // Cascade not physically reasonable
624  if (afin <= minimum_recoil_A && minimum_recoil_A < tnuclei->getA()) {
626  if (verboseLevel > 3) {
627  G4cout << " minimum recoil fragment increased to A " << minimum_recoil_A
628  << G4endl;
629  }
630  }
631 
632  if (verboseLevel>2) G4cerr << " Cascade failed. Retrying..." << G4endl;
633  return false;
634 }
635 
636 
637 // Transfer finished cascade to return buffer
638 
639 void
641  G4InuclParticle* target,
642  G4CollisionOutput& globalOutput) {
643  if (itry >= itry_max) {
644  if (verboseLevel) {
645  G4cout << " IntraNucleiCascader-> no inelastic interaction after "
646  << itry << " attempts " << G4endl;
647  }
648 
649  output.trivialise(bullet, target);
650  } else if (verboseLevel) {
651  G4cout << " IntraNucleiCascader output after trials " << itry << G4endl;
652  }
653 
654  // Copy final generated cascade to output buffer for return
655  globalOutput.add(output);
656 }
657 
658 
659 // Create simple nucleus from rescattering target
660 
663  G4int theNucleusA = theNucleus->GetMassNumber();
664  G4int theNucleusZ = theNucleus->GetCharge();
665 
666  if (theNucleusA > 1) {
667  if (!nucleusTarget) nucleusTarget = new G4InuclNuclei; // Just in case
668  nucleusTarget->fill(0., theNucleusA, theNucleusZ, 0.);
669  return nucleusTarget;
670  } else {
672  protonTarget->fill(0., (theNucleusZ==1)?proton:neutron);
673  return protonTarget;
674  }
675 
676  return 0; // Can never actually get here
677 }
678 
679 // Copy existing (rescattering) cascade for propagation
680 
681 void
683  G4KineticTrackVector* theSecondaries) {
684  if (verboseLevel > 1)
685  G4cout << " >>> G4IntraNucleiCascader::preloadCascade" << G4endl;
686 
687  copyWoundedNucleus(theNucleus); // Update interacted nucleon counts
688  copySecondaries(theSecondaries); // Copy original to internal list
689 }
690 
692  if (verboseLevel > 1)
693  G4cout << " >>> G4IntraNucleiCascader::copyWoundedNucleus" << G4endl;
694 
695  // Loop over nucleons and count hits as exciton holes
697  hitNucleons.clear();
698  if (theNucleus->StartLoop()) {
699  G4Nucleon* nucl = 0;
700  G4int nuclType = 0;
701  while ((nucl = theNucleus->GetNextNucleon())) {
702  if (nucl->AreYouHit()) { // Found previously interacted nucleon
705  hitNucleons.push_back(nucl->GetPosition());
706  }
707  }
708  }
709 
710  if (verboseLevel > 3)
711  G4cout << " nucleus has " << theExitonConfiguration.neutronHoles
712  << " neutrons hit, " << theExitonConfiguration.protonHoles
713  << " protons hit" << G4endl;
714 
715  // Preload nuclear model with confirmed hits, including locations
718 }
719 
720 void
722  if (verboseLevel > 1)
723  G4cout << " >>> G4IntraNucleiCascader::copySecondaries" << G4endl;
724 
725  for (size_t i=0; i<secondaries->size(); i++) {
726  if (verboseLevel > 3) G4cout << " processing secondary " << i << G4endl;
727 
728  processSecondary((*secondaries)[i]); // Copy to cascade or to output
729  }
730 
731  // Sort list of secondaries to put leading particle first
732  std::sort(cascad_particles.begin(), cascad_particles.end(),
734 
735  if (verboseLevel > 2) {
736  G4cout << " Original list of " << secondaries->size() << " secondaries"
737  << " produced " << cascad_particles.size() << " cascade, "
738  << output.numberOfOutgoingParticles() << " released particles, "
739  << output.numberOfOutgoingNuclei() << " fragments" << G4endl;
740  }
741 }
742 
743 
744 // Convert from pre-cascade secondary to local version
745 
747  if (!ktrack) return; // Sanity check
748 
749  // Get particle type to determine whether to keep or release
750  const G4ParticleDefinition* kpd = ktrack->GetDefinition();
751  if (!kpd) return;
752 
754  if (!ktype) {
755  releaseSecondary(ktrack);
756  return;
757  }
758 
759  if (verboseLevel > 1) {
760  G4cout << " >>> G4IntraNucleiCascader::processSecondary "
761  << kpd->GetParticleName() << G4endl;
762  }
763 
764  // Allocate next local particle in buffer and fill
765  cascad_particles.resize(cascad_particles.size()+1); // Like push_back();
766  G4CascadParticle& cpart = cascad_particles.back();
767 
768  // Convert momentum to Bertini internal units
769  cpart.getParticle().fill(ktrack->Get4Momentum()/GeV, ktype);
770  cpart.setGeneration(1);
771  cpart.setMovingInsideNuclei();
772  cpart.initializePath(0);
773 
774  // Convert position units to Bertini's internal scale
775  G4ThreeVector cpos = ktrack->GetPosition()/model->getRadiusUnits();
776 
777  cpart.updatePosition(cpos);
778  cpart.updateZone(model->getZone(cpos.mag()));
779 
780  if (verboseLevel > 2)
781  G4cout << " Created cascade particle \n" << cpart << G4endl;
782 }
783 
784 
785 // Transfer unusable pre-cascade secondaries directly to output
786 
788  const G4ParticleDefinition* kpd = ktrack->GetDefinition();
789 
790  if (verboseLevel > 1) {
791  G4cout << " >>> G4IntraNucleiCascader::releaseSecondary "
792  << kpd->GetParticleName() << G4endl;
793  }
794 
795  // Convert light ion into nucleus on fragment list
796  if (dynamic_cast<const G4Ions*>(kpd)) {
797  // Use resize() and fill() to avoid memory churn
799  G4InuclNuclei& inucl = output.getOutgoingNuclei().back();
800 
801  inucl.fill(ktrack->Get4Momentum()/GeV,
802  kpd->GetAtomicMass(), kpd->GetAtomicNumber());
803  if (verboseLevel > 2)
804  G4cout << " Created pre-cascade fragment\n" << inucl << G4endl;
805  } else {
806  // Use resize() and fill() to avoid memory churn
809 
810  // SPECIAL: Use G4PartDef directly, allowing unknown type code
811  ipart.fill(ktrack->Get4Momentum()/GeV, ktrack->GetDefinition());
812  if (verboseLevel > 2)
813  G4cout << " Created invalid pre-cascade particle\n" << ipart << G4endl;
814  }
815 }
816 
817 
818 // Convert particles which cannot escape into excitons (or eject/decay them)
819 
822  const G4InuclElementaryParticle& trappedP = trapped.getParticle();
823 
824  G4int xtype = trappedP.type();
825  if (verboseLevel > 3) G4cout << " exciton of type " << xtype << G4endl;
826 
827  if (trappedP.nucleon()) { // normal exciton (proton or neutron)
830  return;
831  }
832 
833  if (trappedP.hyperon()) { // Not nucleon, so must be hyperon
834  decayTrappedParticle(trapped);
836  return;
837  }
838 
839  // non-standard exciton; release it
840  // FIXME: this is a meson, so need to absorb it
841  if (verboseLevel > 3) {
842  G4cout << " non-standard should be absorbed, now released\n"
843  << trapped << G4endl;
844  }
845 
846  output.addOutgoingParticle(trappedP);
847 }
848 
849 
850 // Decay unstable trapped particles, and add secondaries to processing list
851 
854  if (verboseLevel > 3)
855  G4cout << " unstable must be decayed in flight" << G4endl;
856 
857  const G4InuclElementaryParticle& trappedP = trapped.getParticle();
858 
859  G4DecayTable* unstable = trappedP.getDefinition()->GetDecayTable();
860  if (!unstable) { // No decay table; cannot decay!
861  if (verboseLevel > 3)
862  G4cerr << " no decay table! Releasing trapped particle" << G4endl;
863 
864  output.addOutgoingParticle(trappedP);
865  return;
866  }
867 
868  // Get secondaries from decay in particle's rest frame
869  G4DecayProducts* daughters = unstable->SelectADecayChannel()->DecayIt( trappedP.getDefinition()->GetPDGMass() );
870  if (!daughters) { // No final state; cannot decay!
871  if (verboseLevel > 3)
872  G4cerr << " no daughters! Releasing trapped particle" << G4endl;
873 
874  output.addOutgoingParticle(trappedP);
875  return;
876  }
877 
878  if (verboseLevel > 3)
879  G4cout << " " << daughters->entries() << " decay daughters" << G4endl;
880 
881  // Convert secondaries to lab frame
882  G4double decayEnergy = trappedP.getEnergy();
883  G4ThreeVector decayDir = trappedP.getMomentum().vect().unit();
884  daughters->Boost(decayEnergy, decayDir);
885 
886  // Put all the secondaries onto the list for propagation
887  const G4ThreeVector& decayPos = trapped.getPosition();
888  G4int zone = trapped.getCurrentZone();
889  G4int gen = trapped.getGeneration()+1;
890 
891  for (G4int i=0; i<daughters->entries(); i++) {
892  G4DynamicParticle* idaug = (*daughters)[i];
893 
895 
896  // Propagate hadronic secondaries with known interactions (tables)
897  if (G4CascadeChannelTables::GetTable(idaugEP.type())) {
898  if (verboseLevel > 3) G4cout << " propagating " << idaugEP << G4endl;
899  cascad_particles.push_back(G4CascadParticle(idaugEP,decayPos,zone,0.,gen));
900  } else {
901  if (verboseLevel > 3) G4cout << " releasing " << idaugEP << G4endl;
902  output.addOutgoingParticle(idaugEP);
903  }
904  }
905 
906  delete daughters; // Clean up memory created by DecayIt()
907 }
G4bool hadNucleus() const
void reset(G4int nHitNeutrons=0, G4int nHitProtons=0, const std::vector< G4ThreeVector > *hitPoints=0)
void fill(G4int a, G4int z, G4double exc=0., Model model=DefaultModel)
void trivialise(G4InuclParticle *bullet, G4InuclParticle *target)
void setVerboseLevel(G4int verbose)
G4InuclElementaryParticle * protonTarget
void processTrappedParticle(const G4CascadParticle &trapped)
std::pair< std::vector< G4CascadParticle >, std::vector< G4InuclElementaryParticle > > modelLists
std::vector< G4CascadParticle > cascad_particles
static const double MeV
Definition: G4SIunits.hh:193
const G4ParticleDefinition * GetParticleType() const
Definition: G4Nucleon.hh:84
static const G4CascadeChannel * GetTable(G4int initialState)
G4int getZ() const
static const G4int reflection_cut
G4CascadeHistory * theCascadeHistory
std::pair< G4int, G4int > getTypesOfNucleonsInvolved() const
static G4bool showHistory()
virtual G4int GetCharge()=0
CLHEP::Hep3Vector G4ThreeVector
G4LorentzVector getMomentum() const
void setVerboseLevel(G4int verbose)
const G4ThreeVector & GetPosition() const
G4int getGeneration() const
void initializePath(G4double npath)
virtual G4bool StartLoop()=0
void printCollisionOutput(std::ostream &os=G4cout) const
static G4bool doCoalescence()
virtual G4int GetMassNumber()=0
void Print(std::ostream &os) const
const G4ParticleDefinition * getDefinition() const
G4ExitonConfiguration theExitonConfiguration
static const G4double quasielast_cut
G4bool stillInside(const G4CascadParticle &cparticle)
G4int code() const
void addOutgoingParticle(const G4InuclElementaryParticle &particle)
G4double getEnergy() const
void updatePosition(const G4ThreeVector &pos)
std::vector< G4ThreeVector > hitNucleons
const G4LorentzVector & getRecoilMomentum() const
virtual const G4ThreeVector & GetPosition() const
Definition: G4Nucleon.hh:68
static G4double getParticleMass(G4int type)
void Boost(G4double totalEnergy, const G4ThreeVector &momentumDirection)
G4CascadParticle initializeCascad(G4InuclElementaryParticle *particle)
void setGeneration(G4int gen)
int G4int
Definition: G4Types.hh:78
G4double getRecoilExcitation() const
G4InuclParticle * createTarget(G4V3DNucleus *theNucleus)
const G4InuclElementaryParticle & getParticle() const
const G4String & GetParticleName() const
G4bool acceptable() const
virtual void setVerboseLevel(G4int verbose=0)
G4int GetAtomicNumber() const
G4double getKineticEnergy() const
std::vector< G4InuclElementaryParticle >::iterator particleIterator
G4DecayTable * GetDecayTable() const
void updateZone(G4int izone)
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
void add(const G4CollisionOutput &right)
G4CascadeCoalescence * theClusterMaker
G4bool initialize(G4InuclParticle *bullet, G4InuclParticle *target)
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
G4int getA() const
void decayTrappedParticle(const G4CascadParticle &trapped)
G4CascadeRecoilMaker * theRecoilMaker
bool G4bool
Definition: G4Types.hh:79
G4bool AreYouHit() const
Definition: G4Nucleon.hh:97
void rescatter(G4InuclParticle *bullet, G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus, G4CollisionOutput &globalOutput)
G4int numberOfOutgoingParticles() const
static const double GeV
Definition: G4SIunits.hh:196
void generateModel(G4InuclNuclei *nuclei)
void DropEntry(const G4CascadParticle &cpart)
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &globalOutput)
G4int GetAtomicMass() const
G4VDecayChannel * SelectADecayChannel()
Definition: G4DecayTable.cc:81
void processSecondary(const G4KineticTrack *aSecondary)
void generateParticleFate(G4CascadParticle &cparticle, G4ElementaryParticleCollider *theEPCollider, std::vector< G4CascadParticle > &cascade)
G4int AddEntry(G4CascadParticle &cpart)
G4int numberOfOutgoingNuclei() const
void addOutgoingParticles(const std::vector< G4InuclElementaryParticle > &particles)
static const G4double small_ekin
G4int getNumberOfProtons() const
void set(G4InuclParticle *part1, G4InuclParticle *part2)
static const G4double ab
void copyWoundedNucleus(G4V3DNucleus *theNucleus)
G4bool empty() const
G4double GetPDGMass() const
void setMovingInsideNuclei(G4bool isMovingIn=true)
void setTolerance(G4double tolerance)
void copySecondaries(G4KineticTrackVector *theSecondaries)
G4bool worthToPropagate(const G4CascadParticle &cparticle) const
const std::vector< G4InuclNuclei > & getOutgoingNuclei() const
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
G4int getCurrentZone() const
G4int getNumberOfNeutrons() const
void fill(G4int ityp, Model model=DefaultModel)
G4double getCharge() const
void setVerboseLevel(G4int verbose)
G4InuclParticle * getBullet() const
#define G4endl
Definition: G4ios.hh:61
virtual void setVerboseLevel(G4int verbose=0)
void setVerboseLevel(G4int verbose=0)
void setVerboseLevel(G4int verbose=0)
virtual G4Nucleon * GetNextNucleon()=0
std::vector< G4CascadParticle > new_cascad_particles
G4int entries() const
const G4ThreeVector & getPosition() const
virtual G4DecayProducts * DecayIt(G4double parentMass=-1.0)=0
void addRecoilFragment(const G4Fragment *aFragment)
void FindClusters(G4CollisionOutput &finalState)
void setRecoilExcitation(G4double Eexc)
double G4double
Definition: G4Types.hh:76
void finalize(G4int itry, G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &globalOutput)
G4InuclParticle * getTarget() const
const G4LorentzVector & Get4Momentum() const
G4ElementaryParticleCollider * theElementaryParticleCollider
G4double getRadiusUnits() const
const G4ParticleDefinition * GetDefinition() const
G4int getNumberOfReflections() const
G4int AddVertex(G4CascadParticle &cpart, std::vector< G4CascadParticle > &daug)
G4int getZone(G4double r) const
void releaseSecondary(const G4KineticTrack *aSecondary)
void addExcitonConfiguration(const G4ExitonConfiguration exciton)
void preloadCascade(G4V3DNucleus *theNucleus, G4KineticTrackVector *theSecondaries)
G4Fragment * makeRecoilFragment()
void setOnShell(G4InuclParticle *bullet, G4InuclParticle *target)
G4GLOB_DLL std::ostream G4cerr
G4InuclElementaryParticle * bparticle
G4double getMass() const
CLHEP::HepLorentzVector G4LorentzVector