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