Geant4  10.00.p03
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 
117 #include <algorithm>
118 
119 #include "G4IntraNucleiCascader.hh"
120 #include "G4SystemOfUnits.hh"
121 #include "G4CascadeChannelTables.hh"
122 #include "G4CascadeCoalescence.hh"
123 #include "G4CascadeHistory.hh"
124 #include "G4CascadeParameters.hh"
125 #include "G4CascadeRecoilMaker.hh"
126 #include "G4CascadParticle.hh"
127 #include "G4CollisionOutput.hh"
128 #include "G4DecayProducts.hh"
129 #include "G4DecayTable.hh"
131 #include "G4ExitonConfiguration.hh"
132 #include "G4HadTmpUtil.hh"
134 #include "G4InuclNuclei.hh"
135 #include "G4InuclParticleNames.hh"
137 #include "G4KineticTrack.hh"
138 #include "G4KineticTrackVector.hh"
139 #include "G4LorentzConvertor.hh"
140 #include "G4Neutron.hh"
141 #include "G4NucleiModel.hh"
142 #include "G4ParticleLargerEkin.hh"
143 #include "G4Proton.hh"
144 #include "G4V3DNucleus.hh"
145 #include "Randomize.hh"
146 
147 using namespace G4InuclParticleNames;
148 using namespace G4InuclSpecialFunctions;
149 
150 
151 // Configuration parameters for cascade production
156 
157 
158 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
159 
161  : G4CascadeColliderBase("G4IntraNucleiCascader"), model(new G4NucleiModel),
162  theElementaryParticleCollider(new G4ElementaryParticleCollider),
163  theRecoilMaker(new G4CascadeRecoilMaker), theClusterMaker(0),
164  theCascadeHistory(0), tnuclei(0), bnuclei(0), bparticle(0),
165  minimum_recoil_A(0.), coulombBarrier(0.),
166  nucleusTarget(new G4InuclNuclei),
167  protonTarget(new G4InuclElementaryParticle) {
170 
173 }
174 
176  delete model;
178  delete theRecoilMaker;
179  delete theClusterMaker;
180  delete theCascadeHistory;
181  delete nucleusTarget;
182  delete protonTarget;
183 }
184 
187  model->setVerboseLevel(verbose);
190 
191  // Optional functionality
194 }
195 
196 
197 
199  G4InuclParticle* target,
200  G4CollisionOutput& globalOutput) {
201  if (verboseLevel) G4cout << " >>> G4IntraNucleiCascader::collide " << G4endl;
202 
203  if (!initialize(bullet, target)) return; // Load buffers and drivers
204 
205  G4int itry = 0;
206  do {
207  newCascade(++itry);
208  setupCascade();
209  generateCascade();
210  } while (!finishCascade() && itry<itry_max);
211 
212  // Report full structure of final cascade if requested
214 
215  finalize(itry, bullet, target, globalOutput);
216 }
217 
218 // For use with Propagate to preload a set of secondaries
219 // FIXME: So far, we don't have any bullet information from Propagate!
220 
222  G4KineticTrackVector* theSecondaries,
223  G4V3DNucleus* theNucleus,
224  G4CollisionOutput& globalOutput) {
225  if (verboseLevel)
226  G4cout << " >>> G4IntraNucleiCascader::rescatter " << G4endl;
227 
228  G4InuclParticle* target = createTarget(theNucleus);
229  if (!initialize(bullet, target)) return; // Load buffers and drivers
230 
231  G4int itry = 0;
232  do {
233  newCascade(++itry);
234  preloadCascade(theNucleus, theSecondaries);
235  generateCascade();
236  } while (!finishCascade() && itry<itry_max);
237 
238  // Report full structure of final cascade if requested
240 
241  finalize(itry, bullet, target, globalOutput);
242 }
243 
244 
246  G4InuclParticle* target) {
247  if (verboseLevel>1)
248  G4cout << " >>> G4IntraNucleiCascader::initialize " << G4endl;
249 
250  // Configure processing modules
252 
253  interCase.set(bullet,target); // Classify collision type
254 
255  if (verboseLevel > 3) {
256  G4cout << *interCase.getBullet() << G4endl
257  << *interCase.getTarget() << G4endl;
258  }
259 
260  // Bullet may be nucleus or simple particle
261  bnuclei = dynamic_cast<G4InuclNuclei*>(interCase.getBullet());
263 
264  if (!bnuclei && !bparticle) {
265  G4cerr << " G4IntraNucleiCascader: projectile is not a valid particle."
266  << G4endl;
267  return false;
268  }
269 
270  // Target _must_ be nucleus
271  tnuclei = dynamic_cast<G4InuclNuclei*>(interCase.getTarget());
272  if (!tnuclei) {
273  if (verboseLevel)
274  G4cerr << " Target is not a nucleus. Abandoning." << G4endl;
275  return false;
276  }
277 
279  coulombBarrier = 0.00126*tnuclei->getZ() / (1.+G4cbrt(tnuclei->getA()));
280 
281  // Energy/momentum conservation usually requires a recoiling nuclear fragment
282  // This cut will be increased on each "itry" if momentum could not balance.
283  minimum_recoil_A = 0.;
284 
285  if (verboseLevel > 3) {
286  G4LorentzVector momentum_in = bullet->getMomentum() + target->getMomentum();
287  G4cout << " intitial momentum E " << momentum_in.e() << " Px "
288  << momentum_in.x() << " Py " << momentum_in.y() << " Pz "
289  << momentum_in.z() << G4endl;
290  }
291 
292  return true;
293 }
294 
295 // Re-initialize buffers for new attempt at cascade
296 
298  if (verboseLevel > 1) {
299  G4cout << " IntraNucleiCascader itry " << itry << " inter_case "
300  << interCase.code() << G4endl;
301  }
302 
303  model->reset(); // Start new cascade process
304  output.reset();
305  new_cascad_particles.clear();
307  cascad_particles.clear(); // List of initial secondaries
308 
310 }
311 
312 
313 // Load initial cascade using nuclear-model calculations
314 
316  if (verboseLevel > 1)
317  G4cout << " >>> G4IntraNucleiCascader::setupCascade" << G4endl;
318 
319  if (interCase.hadNucleus()) { // particle with nuclei
320  if (verboseLevel > 3)
321  G4cout << " bparticle charge " << bparticle->getCharge()
322  << " baryon number " << bparticle->baryon() << G4endl;
323 
325  } else { // nuclei with nuclei
326  G4int ab = bnuclei->getA();
327  G4int zb = bnuclei->getZ();
328 
329  G4NucleiModel::modelLists all_particles; // Buffer to receive lists
330  model->initializeCascad(bnuclei, tnuclei, all_particles);
331 
332  cascad_particles = all_particles.first;
333  output.addOutgoingParticles(all_particles.second);
334 
335  if (cascad_particles.size() == 0) { // compound nuclei
336  G4int i;
337 
338  for (i = 0; i < ab; i++) {
339  G4int knd = i < zb ? 1 : 2;
341  };
342 
343  G4int ihn = G4int(2 * (ab-zb) * inuclRndm() + 0.5);
344  G4int ihz = G4int(2 * zb * inuclRndm() + 0.5);
345 
346  for (i = 0; i < ihn; i++) theExitonConfiguration.incrementHoles(2);
347  for (i = 0; i < ihz; i++) theExitonConfiguration.incrementHoles(1);
348  }
349  } // if (interCase ...
350 }
351 
352 
353 // Generate one possible cascade (all secondaries, etc.)
354 
356  if (verboseLevel>1) G4cout << " generateCascade " << G4endl;
357 
358  G4int iloop = 0;
359  while (!cascad_particles.empty() && !model->empty()) {
360  iloop++;
361 
362  if (verboseLevel > 2) {
363  G4cout << " Iteration " << iloop << ": Number of cparticles "
364  << cascad_particles.size() << " last one: \n"
365  << cascad_particles.back() << G4endl;
366  }
367 
368  // Record incident particle first, to get history ID
369  if (theCascadeHistory) {
371  if (verboseLevel > 2) {
372  G4cout << " active cparticle got history ID "
373  << cascad_particles.back().getHistoryId() << G4endl;
374  }
375  }
376 
380 
381  // Record interaction for later reporting (if desired)
382  if (theCascadeHistory && new_cascad_particles.size()>1)
384 
385  if (verboseLevel > 2) {
386  G4cout << " After generate fate: New particles "
387  << new_cascad_particles.size() << G4endl
388  << " Discarding last cparticle from list " << G4endl;
389  }
390 
391  cascad_particles.pop_back();
392 
393  // handle the result of a new step
394 
395  if (new_cascad_particles.size() == 1) { // last particle goes without interaction
396  const G4CascadParticle& currentCParticle = new_cascad_particles[0];
397 
398  if (model->stillInside(currentCParticle)) {
399  if (verboseLevel > 3)
400  G4cout << " particle still inside nucleus " << G4endl;
401 
402  if (currentCParticle.getNumberOfReflections() < reflection_cut &&
403  model->worthToPropagate(currentCParticle)) {
404  if (verboseLevel > 3) G4cout << " continue reflections " << G4endl;
405  cascad_particles.push_back(currentCParticle);
406  } else {
407  processTrappedParticle(currentCParticle);
408  } // reflection or exciton
409 
410  } else { // particle about to leave nucleus - check for Coulomb barrier
411  if (verboseLevel > 3) G4cout << " possible escape " << G4endl;
412 
413  const G4InuclElementaryParticle& currentParticle =
414  currentCParticle.getParticle();
415 
416  G4double KE = currentParticle.getKineticEnergy();
417  G4double mass = currentParticle.getMass();
418  G4double Q = currentParticle.getCharge();
419 
420  if (verboseLevel > 3)
421  G4cout << " KE " << KE << " barrier " << Q*coulombBarrier << G4endl;
422 
423  if (KE < Q*coulombBarrier) {
424  // Calculate barrier penetration
425  G4double CBP = 0.0;
426 
427  // if (KE > 0.0001) CBP = std::exp(-0.00126*tnuclei->getZ()*0.25*
428  // (1./KE - 1./coulombBarrier));
429  if (KE > 0.0001) CBP = std::exp(-0.0181*0.5*tnuclei->getZ()*
430  (1./KE - 1./coulombBarrier)*
431  std::sqrt(mass*(coulombBarrier-KE)) );
432 
433  if (G4UniformRand() < CBP) {
434  if (verboseLevel > 3)
435  G4cout << " tunneled\n" << currentParticle << G4endl;
436 
437  // Tunnelling through barrier leaves KE unchanged
438  output.addOutgoingParticle(currentParticle);
439  } else {
440  processTrappedParticle(currentCParticle);
441  }
442  } else {
443  output.addOutgoingParticle(currentParticle);
444 
445  if (verboseLevel > 3)
446  G4cout << " Goes out\n" << output.getOutgoingParticles().back()
447  << G4endl;
448  }
449  }
450  } else { // interaction
451  if (verboseLevel > 3)
452  G4cout << " interacted, adding new to list " << G4endl;
453 
454  cascad_particles.insert(cascad_particles.end(),
455  new_cascad_particles.begin(),
456  new_cascad_particles.end());
457 
458  std::pair<G4int, G4int> holes = model->getTypesOfNucleonsInvolved();
459  if (verboseLevel > 3)
460  G4cout << " adding new exciton holes " << holes.first << ","
461  << holes.second << G4endl;
462 
464 
465  if (holes.second > 0)
467  } // if (new_cascad_particles ...
468 
469  // Evaluate nuclear residue
472 
473  G4double aresid = theRecoilMaker->getRecoilA();
474  if (verboseLevel > 2) {
475  G4cout << " cparticles remaining " << cascad_particles.size()
476  << " nucleus (model) has "
477  << model->getNumberOfNeutrons() << " n, "
478  << model->getNumberOfProtons() << " p "
479  << " residual fragment A " << aresid << G4endl;
480  }
481 
482  if (aresid <= minimum_recoil_A) return; // Must have minimum fragment
483  } // while cascade-list and model
484 }
485 
486 
487 // Conslidate results of cascade and evaluate success
488 
490  if (verboseLevel > 1)
491  G4cout << " >>> G4IntraNucleiCascader::finishCascade ?" << G4endl;
492 
493  // Add left-over cascade particles to output
495  cascad_particles.clear();
496 
497  // Cascade is finished. Check if it's OK.
498  if (verboseLevel>3) {
499  G4cout << " G4IntraNucleiCascader finished" << G4endl;
501  }
502 
503  // Apply cluster coalesence model to produce light ions
504  if (theClusterMaker) {
507 
508  // Update recoil fragment after generating light ions
509  if (verboseLevel>3) G4cout << " Recomputing recoil fragment" << G4endl;
511  output);
512  if (verboseLevel>3) {
513  G4cout << " After cluster coalescence" << G4endl;
515  }
516  }
517 
518  // Use last created recoil fragment instead of re-constructing
519  G4int afin = theRecoilMaker->getRecoilA();
520  G4int zfin = theRecoilMaker->getRecoilZ();
521 
522  // FIXME: Should we deal with unbalanced (0,0) case before rejecting?
523 
524  // Sanity check before proceeding
526  if (verboseLevel > 1)
527  G4cerr << " Recoil nucleus is not physical: A=" << afin << " Z="
528  << zfin << G4endl;
529  return false; // Discard event and try again
530  }
531 
533 
534  if (verboseLevel > 1) {
535  G4cout << " afin " << afin << " zfin " << zfin << G4endl;
536  }
537 
538  if (afin == 0) return true; // Whole event fragmented, exit
539 
540  if (afin == 1) { // Add bare nucleon to particle list
541  G4int last_type = (zfin==1) ? 1 : 2; // proton=1, neutron=2
542 
544  G4double mres = presid.m();
545 
546  // Check for sensible kinematics
547  if (mres-mass < -small_ekin) { // Insufficient recoil energy
548  if (verboseLevel > 2) G4cerr << " unphysical recoil nucleon" << G4endl;
549  return false;
550  }
551 
552  if (mres-mass > small_ekin) { // Too much extra energy
553  if (verboseLevel > 2)
554  G4cerr << " extra energy with recoil nucleon" << G4endl;
555 
556  // FIXME: For now, we add the nucleon as unbalanced, and let
557  // "SetOnShell" fudge things. This should be abandoned.
558  }
559 
560  G4InuclElementaryParticle last_particle(presid, last_type,
562 
563  if (verboseLevel > 3) {
564  G4cout << " adding recoiling nucleon to output list\n"
565  << last_particle << G4endl;
566  }
567 
568  output.addOutgoingParticle(last_particle);
569 
570  // Update recoil to include residual nucleon
572  output);
573  }
574 
575  // Process recoil fragment for consistency, exit or reject
576  if (output.numberOfOutgoingParticles() == 1) {
578  if (std::abs(Eex) < quasielast_cut) {
579  if (verboseLevel > 3) {
580  G4cout << " quasi-elastic scatter with " << Eex << " MeV recoil"
581  << G4endl;
582  }
583 
585  if (verboseLevel > 3) {
586  G4cout << " Eex reset to " << theRecoilMaker->getRecoilExcitation()
587  << G4endl;
588  }
589  }
590  }
591 
592  if (theRecoilMaker->goodNucleus()) {
594 
596  if (!recoilFrag) {
597  G4cerr << "Got null pointer for recoil fragment!" << G4endl;
598  return false;
599  }
600 
601  if (verboseLevel > 2)
602  G4cout << " adding recoil fragment to output list" << G4endl;
603 
604  output.addRecoilFragment(*recoilFrag);
605  }
606 
607  // Put final-state particles in "leading order" for return
608  std::vector<G4InuclElementaryParticle>& opart = output.getOutgoingParticles();
609  std::sort(opart.begin(), opart.end(), G4ParticleLargerEkin());
610 
611  // Adjust final state to balance momentum and energy if necessary
616 
617  if (output.acceptable()) return true;
618  else if (verboseLevel>2) G4cerr << " Cascade setOnShell failed." << G4endl;
619  }
620 
621  // Cascade not physically reasonable
622  if (afin <= minimum_recoil_A && minimum_recoil_A < tnuclei->getA()) {
624  if (verboseLevel > 3) {
625  G4cout << " minimum recoil fragment increased to A " << minimum_recoil_A
626  << G4endl;
627  }
628  }
629 
630  if (verboseLevel>2) G4cerr << " Cascade failed. Retrying..." << G4endl;
631  return false;
632 }
633 
634 
635 // Transfer finished cascade to return buffer
636 
637 void
639  G4InuclParticle* target,
640  G4CollisionOutput& globalOutput) {
641  if (itry >= itry_max) {
642  if (verboseLevel) {
643  G4cout << " IntraNucleiCascader-> no inelastic interaction after "
644  << itry << " attempts " << G4endl;
645  }
646 
647  output.trivialise(bullet, target);
648  } else if (verboseLevel) {
649  G4cout << " IntraNucleiCascader output after trials " << itry << G4endl;
650  }
651 
652  // Copy final generated cascade to output buffer for return
653  globalOutput.add(output);
654 }
655 
656 
657 // Create simple nucleus from rescattering target
658 
661  G4int theNucleusA = theNucleus->GetMassNumber();
662  G4int theNucleusZ = theNucleus->GetCharge();
663 
664  if (theNucleusA > 1) {
665  if (!nucleusTarget) nucleusTarget = new G4InuclNuclei; // Just in case
666  nucleusTarget->fill(0., theNucleusA, theNucleusZ, 0.);
667  return nucleusTarget;
668  } else {
670  protonTarget->fill(0., (theNucleusZ==1)?proton:neutron);
671  return protonTarget;
672  }
673 
674  return 0; // Can never actually get here
675 }
676 
677 // Copy existing (rescattering) cascade for propagation
678 
679 void
681  G4KineticTrackVector* theSecondaries) {
682  if (verboseLevel > 1)
683  G4cout << " >>> G4IntraNucleiCascader::preloadCascade" << G4endl;
684 
685  copyWoundedNucleus(theNucleus); // Update interacted nucleon counts
686  copySecondaries(theSecondaries); // Copy original to internal list
687 }
688 
690  if (verboseLevel > 1)
691  G4cout << " >>> G4IntraNucleiCascader::copyWoundedNucleus" << G4endl;
692 
693  // Loop over nucleons and count hits as exciton holes
695  hitNucleons.clear();
696  if (theNucleus->StartLoop()) {
697  G4Nucleon* nucl = 0;
698  G4int nuclType = 0;
699  while ((nucl = theNucleus->GetNextNucleon())) {
700  if (nucl->AreYouHit()) { // Found previously interacted nucleon
703  hitNucleons.push_back(nucl->GetPosition());
704  }
705  }
706  }
707 
708  if (verboseLevel > 3)
709  G4cout << " nucleus has " << theExitonConfiguration.neutronHoles
710  << " neutrons hit, " << theExitonConfiguration.protonHoles
711  << " protons hit" << G4endl;
712 
713  // Preload nuclear model with confirmed hits, including locations
716 }
717 
718 void
720  if (verboseLevel > 1)
721  G4cout << " >>> G4IntraNucleiCascader::copySecondaries" << G4endl;
722 
723  for (size_t i=0; i<secondaries->size(); i++) {
724  if (verboseLevel > 3) G4cout << " processing secondary " << i << G4endl;
725 
726  processSecondary((*secondaries)[i]); // Copy to cascade or to output
727  }
728 
729  // Sort list of secondaries to put leading particle first
730  std::sort(cascad_particles.begin(), cascad_particles.end(),
732 
733  if (verboseLevel > 2) {
734  G4cout << " Original list of " << secondaries->size() << " secondaries"
735  << " produced " << cascad_particles.size() << " cascade, "
736  << output.numberOfOutgoingParticles() << " released particles, "
737  << output.numberOfOutgoingNuclei() << " fragments" << G4endl;
738  }
739 }
740 
741 
742 // Convert from pre-cascade secondary to local version
743 
745  if (!ktrack) return; // Sanity check
746 
747  // Get particle type to determine whether to keep or release
748  G4ParticleDefinition* kpd = ktrack->GetDefinition();
749  if (!kpd) return;
750 
752  if (!ktype) {
753  releaseSecondary(ktrack);
754  return;
755  }
756 
757  if (verboseLevel > 1) {
758  G4cout << " >>> G4IntraNucleiCascader::processSecondary "
759  << kpd->GetParticleName() << G4endl;
760  }
761 
762  // Allocate next local particle in buffer and fill
763  cascad_particles.resize(cascad_particles.size()+1); // Like push_back();
764  G4CascadParticle& cpart = cascad_particles.back();
765 
766  // Convert momentum to Bertini internal units
767  cpart.getParticle().fill(ktrack->Get4Momentum()/GeV, ktype);
768  cpart.setGeneration(1);
769  cpart.setMovingInsideNuclei();
770  cpart.initializePath(0);
771 
772  // Convert position units to Bertini's internal scale
773  G4ThreeVector cpos = ktrack->GetPosition()/model->getRadiusUnits();
774 
775  cpart.updatePosition(cpos);
776  cpart.updateZone(model->getZone(cpos.mag()));
777 
778  if (verboseLevel > 2)
779  G4cout << " Created cascade particle \n" << cpart << G4endl;
780 }
781 
782 
783 // Transfer unusable pre-cascade secondaries directly to output
784 
786  G4ParticleDefinition* kpd = ktrack->GetDefinition();
787 
788  if (verboseLevel > 1) {
789  G4cout << " >>> G4IntraNucleiCascader::releaseSecondary "
790  << kpd->GetParticleName() << G4endl;
791  }
792 
793  // Convert light ion into nucleus on fragment list
794  if (dynamic_cast<G4Ions*>(kpd)) {
795  // Use resize() and fill() to avoid memory churn
797  G4InuclNuclei& inucl = output.getOutgoingNuclei().back();
798 
799  inucl.fill(ktrack->Get4Momentum()/GeV,
800  kpd->GetAtomicMass(), kpd->GetAtomicNumber());
801  if (verboseLevel > 2)
802  G4cout << " Created pre-cascade fragment\n" << inucl << G4endl;
803  } else {
804  // Use resize() and fill() to avoid memory churn
807 
808  // SPECIAL: Use G4PartDef directly, allowing unknown type code
809  ipart.fill(ktrack->Get4Momentum()/GeV, ktrack->GetDefinition());
810  if (verboseLevel > 2)
811  G4cout << " Created invalid pre-cascade particle\n" << ipart << G4endl;
812  }
813 }
814 
815 
816 // Convert particles which cannot escape into excitons (or eject/decay them)
817 
820  const G4InuclElementaryParticle& trappedP = trapped.getParticle();
821 
822  G4int xtype = trappedP.type();
823  if (verboseLevel > 3) G4cout << " exciton of type " << xtype << G4endl;
824 
825  if (trappedP.nucleon()) { // normal exciton (proton or neutron)
828  return;
829  }
830 
831  if (trappedP.hyperon()) { // Not nucleon, so must be hyperon
832  decayTrappedParticle(trapped);
834  return;
835  }
836 
837  // non-standard exciton; release it
838  // FIXME: this is a meson, so need to absorb it
839  if (verboseLevel > 3) {
840  G4cout << " non-standard should be absorbed, now released\n"
841  << trapped << G4endl;
842  }
843 
844  output.addOutgoingParticle(trappedP);
845 }
846 
847 
848 // Decay unstable trapped particles, and add secondaries to processing list
849 
852  if (verboseLevel > 3)
853  G4cout << " unstable must be decayed in flight" << G4endl;
854 
855  const G4InuclElementaryParticle& trappedP = trapped.getParticle();
856 
857  G4DecayTable* unstable = trappedP.getDefinition()->GetDecayTable();
858  if (!unstable) { // No decay table; cannot decay!
859  if (verboseLevel > 3)
860  G4cerr << " no decay table! Releasing trapped particle" << G4endl;
861 
862  output.addOutgoingParticle(trappedP);
863  return;
864  }
865 
866  // Get secondaries from decay in particle's rest frame
867  G4DecayProducts* daughters = unstable->SelectADecayChannel()->DecayIt( trappedP.getDefinition()->GetPDGMass() );
868  if (!daughters) { // No final state; cannot decay!
869  if (verboseLevel > 3)
870  G4cerr << " no daughters! Releasing trapped particle" << G4endl;
871 
872  output.addOutgoingParticle(trappedP);
873  return;
874  }
875 
876  if (verboseLevel > 3)
877  G4cout << " " << daughters->entries() << " decay daughters" << G4endl;
878 
879  // Convert secondaries to lab frame
880  G4double decayEnergy = trappedP.getEnergy();
881  G4ThreeVector decayDir = trappedP.getMomentum().vect().unit();
882  daughters->Boost(decayEnergy, decayDir);
883 
884  // Put all the secondaries onto the list for propagation
885  const G4ThreeVector& decayPos = trapped.getPosition();
886  G4int zone = trapped.getCurrentZone();
887  G4int gen = trapped.getGeneration()+1;
888 
889  for (G4int i=0; i<daughters->entries(); i++) {
890  G4DynamicParticle* idaug = (*daughters)[i];
891 
893 
894  // Propagate hadronic secondaries with known interactions (tables)
895  if (G4CascadeChannelTables::GetTable(idaugEP.type())) {
896  if (verboseLevel > 3) G4cout << " propagating " << idaugEP << G4endl;
897  cascad_particles.push_back(G4CascadParticle(idaugEP,decayPos,zone,0.,gen));
898  } else {
899  if (verboseLevel > 3) G4cout << " releasing " << idaugEP << G4endl;
900  output.addOutgoingParticle(idaugEP);
901  }
902  }
903 
904  delete daughters; // Clean up memory created by DecayIt()
905 }
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
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
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
G4ParticleDefinition * GetDefinition() const
G4bool initialize(G4InuclParticle *bullet, G4InuclParticle *target)
#define G4UniformRand()
Definition: Randomize.hh:87
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)
G4ParticleDefinition * GetParticleType() const
Definition: G4Nucleon.hh:84
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
G4int getNumberOfReflections() const
G4int AddVertex(G4CascadParticle &cpart, std::vector< G4CascadParticle > &daug)
G4ParticleDefinition * getDefinition() const
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