Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4CascadeInterface.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: G4CascadeInterface.cc 71719 2013-06-21 00:01:54Z mkelsey $
27 //
28 // 20100114 M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
29 // 20100413 M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
30 // 20100414 M. Kelsey -- Check for K0L/K0S before using G4InuclElemPart::type
31 // 20100418 M. Kelsey -- Reference output particle lists via const-ref, use
32 // const_iterator for both.
33 // 20100428 M. Kelsey -- Use G4InuclParticleNames enum
34 // 20100429 M. Kelsey -- Change "case gamma:" to "case photon:"
35 // 20100517 M. Kelsey -- Follow new ctors for G4*Collider family.
36 // 20100520 M. Kelsey -- Simplify collision loop, move momentum rotations to
37 // G4CollisionOutput, copy G4DynamicParticle directly from
38 // G4InuclParticle, no switch-block required.
39 // 20100615 M. Kelsey -- Bug fix: For K0's need ekin in GEANT4 units
40 // 20100617 M. Kelsey -- Rename "debug_" preprocessor flag to G4CASCADE_DEBUG,
41 // and "BERTDEV" to "G4CASCADE_COULOMB_DEV"
42 // 20100617 M. Kelsey -- Make G4InuclCollider a local data member
43 // 20100618 M. Kelsey -- Deploy energy-conservation test on final state, with
44 // preprocessor flag G4CASCADE_SKIP_ECONS to skip test.
45 // 20100620 M. Kelsey -- Use new energy-conservation pseudo-collider
46 // 20100621 M. Kelsey -- Fix compiler warning from GCC 4.5
47 // 20100624 M. Kelsey -- Fix cascade loop to check nTries every time (had
48 // allowed for infinite loop on E-violation); dump event data
49 // to output if E-violation exceeds maxTries; use CheckBalance
50 // for baryon and charge conservation.
51 // 20100701 M. Kelsey -- Pass verbosity through to G4CollisionOutput
52 // 20100714 M. Kelsey -- Report number of iterations before success
53 // 20100720 M. Kelsey -- Use G4CASCADE_SKIP_ECONS flag for reporting
54 // 20100723 M. Kelsey -- Move G4CollisionOutput to .hh file for reuse
55 // 20100914 M. Kelsey -- Migrate to integer A and Z
56 // 20100916 M. Kelsey -- Simplify ApplyYourself() by encapsulating code blocks
57 // into numerous functions; make data-member colliders pointers;
58 // provide support for projectile nucleus
59 // 20100919 M. Kelsey -- Fix incorrect logic in retryInelasticNucleus()
60 // 20100922 M. Kelsey -- Add functions to select de-excitation method
61 // 20100924 M. Kelsey -- Migrate to "OutgoingNuclei" names in CollisionOutput
62 // 20110224 M. Kelsey -- Add createTarget() for use with Propagate(); split
63 // conservation law messages to separate function; begin to add
64 // infrastructure code to Propagate. Move verbose
65 // setting to .cc file, and apply to all member objects.
66 // 20110301 M. Kelsey -- Add copyPreviousCascade() for use with Propagate()
67 // along with new buffers and related particle-conversion
68 // functions. Encapulate buffer deletion in clear(). Add some
69 // diagnostic messages.
70 // 20110302 M. Kelsey -- Redo diagnostics inside G4CASCADE_DEBUG_INTERFACE
71 // 20110304 M. Kelsey -- Drop conversion of Propagate() arguments; pass
72 // directly to collider for processing. Rename makeReactionProduct
73 // to makeDynamicParticle.
74 // 20110316 M. Kelsey -- Move kaon-mixing for type-code into G4InuclParticle;
75 // add placeholders to capture and use bullet in Propagte
76 // 20110327 G. Folger -- Set up for E/p checking by G4HadronicProcess in ctor
77 // 20110328 M. Kelsey -- Modify balance() initialization to match Gunter's
78 // 20110404 M. Kelsey -- Get primary projectile from base class (ref-03)
79 // 20110502 M. Kelsey -- Add interface to capture random seeds for user
80 // 20110719 M. Kelsey -- Use trivialise() in case maximum retries are reached
81 // 20110720 M. Kelsey -- Discard elastic-cut array (no longer needed),
82 // discard local "theFinalState" (in base as "theParticleChange"),
83 // Modify createBullet() to set null pointer if bullet unusable,
84 // return empty final-state on failures.
85 // Fix charge violation report before throwing exception.
86 // 20110722 M. Kelsey -- In makeDynamicParticle(), allow invalid type codes
87 // in order to process, e.g., resonances from Propagate() input
88 // 20110728 M. Kelsey -- Per V.Ivantchenko, change NoInteraction to return
89 // zero particles, but set kinetic energy from projectile.
90 // 20110801 M. Kelsey -- Make bullet, target point to local buffers, no delete
91 // 20110802 M. Kelsey -- Use new decay handler for Propagate interface
92 // 20110922 M. Kelsey -- Follow migration of G4InuclParticle::print(), use
93 // G4ExceptionDescription for reporting before throwing exception
94 // 20120125 M. Kelsey -- In retryInelasticProton() check for empty output
95 // 20120525 M. Kelsey -- In NoInteraction, check for Ekin<0., set to zero;
96 // use SetEnergyChange(0.) explicitly for good final states.
97 // 20120822 M. Kelsey -- Move envvars to G4CascadeParameters.
98 // 20130508 D. Wright -- Add support for muon capture
99 // 20130804 M. Kelsey -- Fix bug #1513 -- "(Z=1)" in boolean expression
100 // 20140116 M. Kelsey -- Move statics to const data members to avoid weird
101 // interactions with MT.
102 // 20140929 M. Kelsey -- Explicitly call useCascadeDeexcitation() in ctor
103 // 20150506 M. Kelsey -- Call Initialize() in ctor for master thread only
104 // 20150608 M. Kelsey -- Label all while loops as terminating.
105 
106 #include <cmath>
107 #include <iostream>
108 
109 #include "G4CascadeInterface.hh"
110 #include "globals.hh"
111 #include "G4SystemOfUnits.hh"
112 #include "G4CascadeChannelTables.hh"
113 #include "G4CascadeCheckBalance.hh"
114 #include "G4CascadeParameters.hh"
115 #include "G4CollisionOutput.hh"
116 #include "G4DecayKineticTracks.hh"
117 #include "G4DynamicParticle.hh"
118 #include "G4HadronicException.hh"
119 #include "G4InuclCollider.hh"
121 #include "G4InuclNuclei.hh"
122 #include "G4InuclParticle.hh"
123 #include "G4InuclParticleNames.hh"
124 #include "G4KaonZeroLong.hh"
125 #include "G4KaonZeroShort.hh"
126 #include "G4KineticTrack.hh"
127 #include "G4KineticTrackVector.hh"
128 #include "G4Nucleus.hh"
129 #include "G4ParticleDefinition.hh"
131 #include "G4Threading.hh"
132 #include "G4Track.hh"
133 #include "G4V3DNucleus.hh"
134 #include "G4UnboundPN.hh"
135 #include "G4Dineutron.hh"
136 #include "G4Diproton.hh"
137 
138 using namespace G4InuclParticleNames;
139 
140 typedef std::vector<G4InuclElementaryParticle>::const_iterator particleIterator;
141 typedef std::vector<G4InuclNuclei>::const_iterator nucleiIterator;
142 
143 
144 // Constructor and destrutor
145 
148  randomFile(G4CascadeParameters::randomFile()),
149  maximumTries(20), numberOfTries(0),
150  collider(new G4InuclCollider), balance(new G4CascadeCheckBalance(name)),
151  bullet(0), target(0), output(new G4CollisionOutput) {
152  // Set up global objects for master thread or sequential build
154 
156  balance->setLimits(5*perCent, 10*MeV/GeV); // Bertini internal units
158 
161  else
163 }
164 
166  clear();
167  delete collider; collider=0;
168  delete balance; balance=0;
169  delete output; output=0;
170 }
171 
172 void G4CascadeInterface::ModelDescription(std::ostream& outFile) const
173 {
174  outFile << "The Bertini-style cascade implements the inelastic scattering\n"
175  << "of hadrons by nuclei. Nucleons, pions, kaons and hyperons\n"
176  << "from 0 to 15 GeV may be used as projectiles in this model.\n"
177  << "Final state hadrons are produced by a classical cascade\n"
178  << "consisting of individual hadron-nucleon scatterings which use\n"
179  << "free-space partial cross sections, corrected for various\n"
180  << "nuclear medium effects. The target nucleus is modeled as a\n"
181  << "set of 1, 3 or 6 spherical shells, in which scattered hadrons\n"
182  << "travel in straight lines until they are reflected from or\n"
183  << "transmitted through shell boundaries.\n";
184 }
185 
186 void G4CascadeInterface::DumpConfiguration(std::ostream& outFile) const {
188 }
189 
191  bullet=0;
192  target=0;
193 }
194 
195 
196 // Initialize shared objects for use across multiple threads
197 
203  if (!ch || !pn || !nn || !pp) return; // Avoid "unused variables"
204 }
205 
206 
207 // Select post-cascade processing (default will be CascadeDeexcitation)
208 // NOTE: Currently just calls through to Collider, in future will do something
209 
211  collider->useCascadeDeexcitation();
212 }
213 
215  collider->usePreCompoundDeexcitation();
216 }
217 
218 
219 // Apply verbosity to all member objects (overrides base class version)
220 
223  collider->setVerboseLevel(verbose);
224  balance->setVerboseLevel(verbose);
225  output->setVerboseLevel(verbose);
226 }
227 
228 
229 // Test whether inputs are valid for this model
230 
232  G4Nucleus& /* theNucleus */) {
233  return IsApplicable(aTrack.GetDefinition());
234 }
235 
237  if (aPD->GetAtomicMass() > 1) return true; // Nuclei are okay
238 
239  // Valid particle and have interactions available
241  return (G4CascadeChannelTables::GetTable(type));
242 }
243 
244 
245 // Main Actions
246 
249  G4Nucleus& theNucleus) {
250  if (verboseLevel)
251  G4cout << " >>> G4CascadeInterface::ApplyYourself" << G4endl;
252 
253  if (aTrack.GetKineticEnergy() < 0.) {
254  G4cerr << " >>> G4CascadeInterface got negative-energy track: "
255  << aTrack.GetDefinition()->GetParticleName() << " Ekin = "
256  << aTrack.GetKineticEnergy() << G4endl;
257  }
258 
259 #ifdef G4CASCADE_DEBUG_INTERFACE
260  static G4int counter(0);
261  counter++;
262  G4cerr << "Reaction number "<< counter << " "
263  << aTrack.GetDefinition()->GetParticleName() << " "
264  << aTrack.GetKineticEnergy() << " MeV" << G4endl;
265 #endif
266 
267  if (!randomFile.empty()) { // User requested random-seed capture
268  if (verboseLevel>1)
269  G4cout << " Saving random engine state to " << randomFile << G4endl;
271  }
272 
274  clear();
275 
276  // Abort processing if no interaction is possible
277  if (!IsApplicable(aTrack, theNucleus)) {
278  if (verboseLevel) G4cerr << " No interaction possible " << G4endl;
279  return NoInteraction(aTrack, theNucleus);
280  }
281 
282  // Make conversion between native Geant4 and Bertini cascade classes.
283  if (!createBullet(aTrack)) {
284  if (verboseLevel) G4cerr << " Unable to create usable bullet" << G4endl;
285  return NoInteraction(aTrack, theNucleus);
286  }
287 
288  if (!createTarget(theNucleus)) {
289  if (verboseLevel) G4cerr << " Unable to create usable target" << G4endl;
290  return NoInteraction(aTrack, theNucleus);
291  }
292 
293  // Different retry conditions for proton target vs. nucleus
294  const G4bool isHydrogen = (theNucleus.GetA_asInt() == 1);
295 
296  numberOfTries = 0;
297  do { // we try to create inelastic interaction
298  if (verboseLevel > 1)
299  G4cout << " Generating cascade attempt " << numberOfTries << G4endl;
300 
301  output->reset();
302  collider->collide(bullet, target, *output);
303  balance->collide(bullet, target, *output);
304 
305  numberOfTries++;
306  /* Loop checking 08.06.2015 MHK */
307  } while ( isHydrogen ? retryInelasticProton() : retryInelasticNucleus() );
308 
309  // Null event if unsuccessful
310  if (numberOfTries >= maximumTries) {
311  if (verboseLevel)
312  G4cout << " Cascade aborted after trials " << numberOfTries << G4endl;
313  return NoInteraction(aTrack, theNucleus);
314  }
315 
316  // Abort job if energy or momentum are not conserved
317  if (!balance->okay()) {
319  return NoInteraction(aTrack, theNucleus);
320  }
321 
322  // Successful cascade -- clean up and return
323  if (verboseLevel) {
324  G4cout << " Cascade output after trials " << numberOfTries << G4endl;
325  if (verboseLevel > 1) output->printCollisionOutput();
326  }
327 
328  // Rotate event to put Z axis along original projectile direction
329  output->rotateEvent(bulletInLabFrame);
330 
332 
333  // Report violations of conservation laws in original frame
335 
336  // Clean up and return final result;
337  clear();
338 /*
339  G4int nSec = theParticleChange.GetNumberOfSecondaries();
340  for (G4int i = 0; i < nSec; i++) {
341  G4HadSecondary* sec = theParticleChange.GetSecondary(i);
342  G4DynamicParticle* dp = sec->GetParticle();
343  if (dp->GetDefinition()->GetParticleName() == "neutron")
344  G4cout << dp->GetDefinition()->GetParticleName() << " has "
345  << dp->GetKineticEnergy()/MeV << " MeV " << G4endl;
346  }
347 */
348  return &theParticleChange;
349 }
350 
353  G4V3DNucleus* theNucleus) {
354  if (verboseLevel) G4cout << " >>> G4CascadeInterface::Propagate" << G4endl;
355 
356 #ifdef G4CASCADE_DEBUG_INTERFACE
357  if (verboseLevel>1) {
358  G4cout << " G4V3DNucleus A " << theNucleus->GetMassNumber()
359  << " Z " << theNucleus->GetCharge()
360  << "\n " << theSecondaries->size() << " secondaries:" << G4endl;
361  for (size_t i=0; i<theSecondaries->size(); i++) {
362  G4KineticTrack* kt = (*theSecondaries)[i];
363  G4cout << " " << i << ": " << kt->GetDefinition()->GetParticleName()
364  << " p " << kt->Get4Momentum() << " @ " << kt->GetPosition()
365  << " t " << kt->GetFormationTime() << G4endl;
366  }
367  }
368 #endif
369 
370  if (!randomFile.empty()) { // User requested random-seed capture
371  if (verboseLevel>1)
372  G4cout << " Saving random engine state to " << randomFile << G4endl;
374  }
375 
377  clear();
378 
379  // Process input secondaries list to eliminate resonances
380  G4DecayKineticTracks decay(theSecondaries);
381 
382  // NOTE: Requires 9.4-ref-03 mods to base class and G4TheoFSGenerator
383  const G4HadProjectile* projectile = GetPrimaryProjectile();
384  if (projectile) createBullet(*projectile);
385 
386  if (!createTarget(theNucleus)) {
387  if (verboseLevel) G4cerr << " Unable to create usable target" << G4endl;
388  return 0; // FIXME: This will cause a segfault later
389  }
390 
391  numberOfTries = 0;
392  do {
393  if (verboseLevel > 1)
394  G4cout << " Generating rescatter attempt " << numberOfTries << G4endl;
395 
396  output->reset();
397  collider->rescatter(bullet, theSecondaries, theNucleus, *output);
398  balance->collide(bullet, target, *output);
399 
400  numberOfTries++;
401  // FIXME: retry checks will SEGFAULT until we can define the bullet!
402  } while (retryInelasticNucleus()); /* Loop checking 08.06.2015 MHK */
403 
404  // Check whether repeated attempts have all failed; report and exit
405  if (numberOfTries >= maximumTries && !balance->okay()) {
406  throwNonConservationFailure(); // This terminates the job
407  }
408 
409  // Successful cascade -- clean up and return
410  if (verboseLevel) {
411  G4cout << " Cascade rescatter after trials " << numberOfTries << G4endl;
412  if (verboseLevel > 1) output->printCollisionOutput();
413  }
414 
415  // Does calling code take ownership? I hope so!
417 
418  // Clean up and and return final result
419  clear();
420  return propResult;
421 }
422 
423 
424 // Replicate input particles onto output
425 
428  G4Nucleus& /*theNucleus*/) {
429  if (verboseLevel)
430  G4cout << " >>> G4CascadeInterface::NoInteraction" << G4endl;
431 
434 
435  G4double ekin = aTrack.GetKineticEnergy()>0. ? aTrack.GetKineticEnergy() : 0.;
436  theParticleChange.SetEnergyChange(ekin); // Protect against rounding
437 
438  return &theParticleChange;
439 }
440 
441 
442 // Convert input projectile to Bertini internal object
443 
445  const G4ParticleDefinition* trkDef = aTrack.GetDefinition();
446  G4int bulletType = 0; // For elementary particles
447  G4int bulletA = 0, bulletZ = 0; // For nucleus projectile
448 
449  if (trkDef->GetAtomicMass() <= 1) {
450  bulletType = G4InuclElementaryParticle::type(trkDef);
451  } else {
452  bulletA = trkDef->GetAtomicMass();
453  bulletZ = trkDef->GetAtomicNumber();
454  }
455 
456  if (0 == bulletType && 0 == bulletA*bulletZ) {
457  if (verboseLevel) {
458  G4cerr << " G4CascadeInterface: " << trkDef->GetParticleName()
459  << " not usable as bullet." << G4endl;
460  }
461  bullet = 0;
462  return false;
463  }
464 
465  // Code momentum and energy -- Bertini wants z-axis and GeV units
466  G4LorentzVector projectileMomentum = aTrack.Get4Momentum()/GeV;
467 
468  // Rotation/boost to get from z-axis back to original frame
469  bulletInLabFrame = G4LorentzRotation::IDENTITY; // Initialize
470  bulletInLabFrame.rotateZ(-projectileMomentum.phi());
471  bulletInLabFrame.rotateY(-projectileMomentum.theta());
472  bulletInLabFrame.invert();
473 
474  G4LorentzVector momentumBullet(0., 0., projectileMomentum.rho(),
475  projectileMomentum.e());
476 
477  if (G4InuclElementaryParticle::valid(bulletType)) {
478  hadronBullet.fill(momentumBullet, bulletType);
479  bullet = &hadronBullet;
480  } else {
481  nucleusBullet.fill(momentumBullet, bulletA, bulletZ);
482  bullet = &nucleusBullet;
483  }
484 
485  if (verboseLevel > 2) G4cout << "Bullet: \n" << *bullet << G4endl;
486 
487  return true;
488 }
489 
490 
491 // Convert input nuclear target to Bertini internal object
492 
494  return createTarget(theNucleus.GetA_asInt(), theNucleus.GetZ_asInt());
495 }
496 
498  return createTarget(theNucleus->GetMassNumber(), theNucleus->GetCharge());
499 }
500 
502  if (A > 1) {
503  nucleusTarget.fill(A, Z);
504  target = &nucleusTarget;
505  } else {
506  hadronTarget.fill(0., (Z==1?proton:neutron));
507  target = &hadronTarget;
508  }
509 
510  if (verboseLevel > 2) G4cout << "Target: \n" << *target << G4endl;
511 
512  return true; // Right now, target never fails
513 }
514 
515 
516 // Convert Bertini particle to output (G4DynamicParticle)
517 
520  G4int outgoingType = iep.type();
521 
522  if (iep.quasi_deutron()) {
523  G4cerr << " ERROR: G4CascadeInterface incompatible particle type "
524  << outgoingType << G4endl;
525  return 0;
526  }
527 
528  // Copy local G4DynPart to public output (handle kaon mixing specially)
529  if (outgoingType == kaonZero || outgoingType == kaonZeroBar) {
530  G4ThreeVector momDir = iep.getMomentum().vect().unit();
531  G4double ekin = iep.getKineticEnergy()*GeV; // Bertini -> G4 units
532 
534  if (G4UniformRand() > 0.5) pd = G4KaonZeroLong::Definition();
535 
536  return new G4DynamicParticle(pd, momDir, ekin);
537  } else {
538  return new G4DynamicParticle(iep.getDynamicParticle());
539  }
540 
541  return 0; // Should never get here!
542 }
543 
546  if (verboseLevel > 2) {
547  G4cout << " Nuclei fragment: \n" << inuc << G4endl;
548  }
549 
550  // Copy local G4DynPart to public output
551  return new G4DynamicParticle(inuc.getDynamicParticle());
552 }
553 
554 
555 // Transfer Bertini internal final state to hadronics interface
556 
558  if (verboseLevel > 1)
559  G4cout << " >>> G4CascadeInterface::copyOutputToHadronicResult" << G4endl;
560 
561  const std::vector<G4InuclNuclei>& outgoingNuclei = output->getOutgoingNuclei();
562  const std::vector<G4InuclElementaryParticle>& particles = output->getOutgoingParticles();
563 
566 
567  // Get outcoming particles
568  if (!particles.empty()) {
569  particleIterator ipart = particles.begin();
570  for (; ipart != particles.end(); ipart++) {
572  }
573  }
574 
575  // get nuclei fragments
576  if (!outgoingNuclei.empty()) {
577  nucleiIterator ifrag = outgoingNuclei.begin();
578  for (; ifrag != outgoingNuclei.end(); ifrag++) {
580  }
581  }
582 }
583 
585  if (verboseLevel > 1)
586  G4cout << " >>> G4CascadeInterface::copyOutputToReactionProducts" << G4endl;
587 
588  const std::vector<G4InuclElementaryParticle>& particles = output->getOutgoingParticles();
589  const std::vector<G4InuclNuclei>& fragments = output->getOutgoingNuclei();
590 
592 
593  G4ReactionProduct* rp = 0; // Buffers to create outgoing tracks
594  G4DynamicParticle* dp = 0;
595 
596  // Get outcoming particles
597  if (!particles.empty()) {
598  particleIterator ipart = particles.begin();
599  for (; ipart != particles.end(); ipart++) {
600  rp = new G4ReactionProduct;
601  dp = makeDynamicParticle(*ipart);
602  (*rp) = (*dp); // This does all the necessary copying
603  propResult->push_back(rp);
604  delete dp;
605  }
606  }
607 
608  // get nuclei fragments
609  if (!fragments.empty()) {
610  nucleiIterator ifrag = fragments.begin();
611  for (; ifrag != fragments.end(); ifrag++) {
612  rp = new G4ReactionProduct;
613  dp = makeDynamicParticle(*ifrag);
614  (*rp) = (*dp); // This does all the necessary copying
615  propResult->push_back(rp);
616  delete dp;
617  }
618  }
619 
620  return propResult;
621 }
622 
623 
624 // Report violations of conservation laws in original frame
625 
627  balance->collide(bullet, target, *output);
628 
629  if (verboseLevel > 2) {
630  if (!balance->baryonOkay()) {
631  G4cerr << "ERROR: no baryon number conservation, sum of baryons = "
632  << balance->deltaB() << G4endl;
633  }
634 
635  if (!balance->chargeOkay()) {
636  G4cerr << "ERROR: no charge conservation, sum of charges = "
637  << balance->deltaQ() << G4endl;
638  }
639 
640  if (std::abs(balance->deltaKE()) > 0.01 ) { // GeV
641  G4cerr << "Kinetic energy conservation violated by "
642  << balance->deltaKE() << " GeV" << G4endl;
643  }
644 
645  G4double eInit = bullet->getEnergy() + target->getEnergy();
646  G4double eFinal = eInit + balance->deltaE();
647 
648  G4cout << "Initial energy " << eInit << " final energy " << eFinal
649  << "\nTotal energy conservation at level "
650  << balance->deltaE() * GeV << " MeV" << G4endl;
651 
652  if (balance->deltaKE() > 5.0e-5 ) { // 0.05 MeV
653  G4cerr << "FATAL ERROR: kinetic energy created "
654  << balance->deltaKE() * GeV << " MeV" << G4endl;
655  }
656  }
657 }
658 
659 
660 // Evaluate whether any outgoing particles penetrated Coulomb barrier
661 
663  G4bool violated = false; // by default coulomb analysis is OK
664 
665  const G4double coulumbBarrier = 8.7 * MeV/GeV; // Bertini uses GeV
666 
667  const std::vector<G4InuclElementaryParticle>& p =
668  output->getOutgoingParticles();
669 
670  for (particleIterator ipart=p.begin(); ipart != p.end(); ipart++) {
671  if (ipart->type() == proton) {
672  violated |= (ipart->getKineticEnergy() < coulumbBarrier);
673  }
674  }
675 
676  return violated;
677 }
678 
679 // Check whether inelastic collision on proton failed
680 
682  const std::vector<G4InuclElementaryParticle>& out =
683  output->getOutgoingParticles();
684 
685 #ifdef G4CASCADE_DEBUG_INTERFACE
686  // Report on all retry conditions, in order of return logic
687  G4cout << " retryInelasticProton: number of Tries "
688  << ((numberOfTries < maximumTries) ? "RETRY (t)" : "EXIT (f)")
689  << "\n retryInelasticProton: AND collision type ";
690  if (out.empty()) G4cout << "FAILED" << G4endl;
691  else {
692  G4cout << (out.size() == 2 ? "ELASTIC (t)" : "INELASTIC (f)")
693  << "\n retryInelasticProton: AND Leading particles bullet "
694  << (out.size() >= 2 &&
695  (out[0].getDefinition() == bullet->getDefinition() ||
696  out[1].getDefinition() == bullet->getDefinition())
697  ? "YES (t)" : "NO (f)")
698  << G4endl;
699  }
700 #endif
701 
702  return ( (numberOfTries < maximumTries) &&
703  (out.empty() ||
704  (out.size() == 2 &&
705  (out[0].getDefinition() == bullet->getDefinition() ||
706  out[1].getDefinition() == bullet->getDefinition())))
707  );
708 }
709 
710 // Check whether generic inelastic collision failed
711 // NOTE: some conditions are set by compiler flags
712 
714  // Quantities necessary for conditional block below
715  G4int npart = output->numberOfOutgoingParticles();
716  G4int nfrag = output->numberOfOutgoingNuclei();
717 
718  const G4ParticleDefinition* firstOut = (npart == 0) ? 0 :
719  output->getOutgoingParticles().begin()->getDefinition();
720 
721 #ifdef G4CASCADE_DEBUG_INTERFACE
722  // Report on all retry conditions, in order of return logic
723  G4cout << " retryInelasticNucleus: numberOfTries "
724  << ((numberOfTries < maximumTries) ? "RETRY (t)" : "EXIT (f)")
725  << "\n retryInelasticNucleus: AND outputParticles "
726  << ((npart != 0) ? "NON-ZERO (t)" : "EMPTY (f)")
727 #ifdef G4CASCADE_COULOMB_DEV
728  << "\n retryInelasticNucleus: AND coulombBarrier (COULOMB_DEV) "
729  << (coulombBarrierViolation() ? "VIOLATED (t)" : "PASSED (f)")
730  << "\n retryInelasticNucleus: AND collision type (COULOMB_DEV) "
731  << ((npart+nfrag > 2) ? "INELASTIC (t)" : "ELASTIC (f)")
732 #else
733  << "\n retryInelasticNucleus: AND collsion type "
734  << ((npart+nfrag < 3) ? "ELASTIC (t)" : "INELASTIC (f)")
735  << "\n retryInelasticNucleus: AND Leading particle bullet "
736  << ((firstOut == bullet->getDefinition()) ? "YES (t)" : "NO (f)")
737 #endif
738  << "\n retryInelasticNucleus: OR conservation "
739  << (!balance->okay() ? "FAILED (t)" : "PASSED (f)")
740  << G4endl;
741 #endif
742 
743  return ( (numberOfTries < maximumTries) &&
744  ( ((npart != 0) &&
745 #ifdef G4CASCADE_COULOMB_DEV
746  (coulombBarrierViolation() && npart+nfrag > 2)
747 #else
748  (npart+nfrag < 3 && firstOut == bullet->getDefinition())
749 #endif
750  )
751 #ifndef G4CASCADE_SKIP_ECONS
752  || (!balance->okay())
753 #endif
754  )
755  );
756 }
757 
758 
759 // Terminate job in case of persistent non-conservation
760 // FIXME: Need to migrate to G4ExceptionDescription
761 
763  // NOTE: Once G4HadronicException is changed, use the following line!
764  // G4ExceptionDescription errInfo;
765  std::ostream& errInfo = G4cerr;
766 
767  errInfo << " >>> G4CascadeInterface has non-conserving"
768  << " cascade after " << numberOfTries << " attempts." << G4endl;
769 
770  G4String throwMsg = "G4CascadeInterface - ";
771  if (!balance->energyOkay()) {
772  throwMsg += "Energy";
773  errInfo << " Energy conservation violated by " << balance->deltaE()
774  << " GeV (" << balance->relativeE() << ")" << G4endl;
775  }
776 
777  if (!balance->momentumOkay()) {
778  throwMsg += "Momentum";
779  errInfo << " Momentum conservation violated by " << balance->deltaP()
780  << " GeV/c (" << balance->relativeP() << ")" << G4endl;
781  }
782 
783  if (!balance->baryonOkay()) {
784  throwMsg += "Baryon number";
785  errInfo << " Baryon number violated by " << balance->deltaB() << G4endl;
786  }
787 
788  if (!balance->chargeOkay()) {
789  throwMsg += "Charge";
790  errInfo << " Charge conservation violated by " << balance->deltaQ()
791  << G4endl;
792  }
793 
794  errInfo << " Final event output, for debugging:\n"
795  << " Bullet: \n" << *bullet << G4endl
796  << " Target: \n" << *target << G4endl;
797  output->printCollisionOutput(errInfo);
798 
799  throwMsg += " non-conservation. More info in output.";
800  throw G4HadronicException(__FILE__, __LINE__, throwMsg); // Job ends here!
801 }
void fill(G4int a, G4int z, G4double exc=0., Model model=DefaultModel)
G4bool retryInelasticNucleus() const
void rescatter(G4InuclParticle *bullet, G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus, G4CollisionOutput &globalOutput)
void setVerboseLevel(G4int verbose)
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
const XML_Char * name
Definition: expat.h:151
G4CascadeInterface(const G4String &name="BertiniCascade")
static const G4CascadeChannel * GetTable(G4int initialState)
G4bool coulombBarrierViolation() const
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
const XML_Char * target
Definition: expat.h:268
virtual G4int GetCharge()=0
const G4DynamicParticle & getDynamicParticle() const
G4LorentzVector getMomentum() const
const G4HadProjectile * GetPrimaryProjectile() const
const G4ThreeVector & GetPosition() const
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
static constexpr double perCent
Definition: G4SIunits.hh:332
const char * p
Definition: xmltok.h:285
void printCollisionOutput(std::ostream &os=G4cout) const
virtual G4int GetMassNumber()=0
const G4ParticleDefinition * getDefinition() const
G4double getEnergy() const
void setLimits(G4double relative, G4double absolute)
static G4Diproton * Definition()
Definition: G4Diproton.cc:68
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
HepLorentzRotation & rotateY(double delta)
virtual void setVerboseLevel(G4int verbose=0)
G4int GetAtomicNumber() const
double phi() const
G4double getKineticEnergy() const
void SetStatusChange(G4HadFinalStateStatus aS)
std::vector< G4ReactionProduct * > G4ReactionProductVector
void useCascadeDeexcitation()
virtual void DumpConfiguration(std::ostream &outFile) const
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
const G4ParticleDefinition * GetDefinition() const
double theta() const
G4double GetFormationTime() const
G4ReactionProductVector * Propagate(G4KineticTrackVector *theSecondaries, G4V3DNucleus *theNucleus)
bool G4bool
Definition: G4Types.hh:79
void setVerboseLevel(G4int verbose=0)
G4double GetKineticEnergy() const
G4bool createTarget(G4Nucleus &theNucleus)
G4int numberOfOutgoingParticles() const
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &globalOutput)
double rho() const
G4ReactionProductVector * copyOutputToReactionProducts()
G4DynamicParticle * makeDynamicParticle(const G4InuclElementaryParticle &iep) const
void rotateEvent(const G4LorentzRotation &rotate)
static G4KaonZeroLong * Definition()
const G4LorentzVector & Get4Momentum() const
G4int GetAtomicMass() const
HepLorentzRotation & rotateZ(double delta)
HepLorentzRotation & invert()
static G4Dineutron * Definition()
Definition: G4Dineutron.cc:68
G4int numberOfOutgoingNuclei() const
void SetEnergyChange(G4double anEnergy)
static void saveEngineStatus(const char filename[]="Config.conf")
Definition: Random.cc:275
Hep3Vector unit() const
const std::vector< G4InuclNuclei > & getOutgoingNuclei() const
static G4KaonZeroShort * Definition()
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
static constexpr double GeV
Definition: G4SIunits.hh:217
void fill(G4int ityp, Model model=DefaultModel)
std::vector< G4InuclNuclei >::const_iterator nucleiIterator
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:65
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
virtual void ModelDescription(std::ostream &outFile) const
G4bool IsMasterThread()
Definition: G4Threading.cc:146
static DLL_API const HepLorentzRotation IDENTITY
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
double G4double
Definition: G4Types.hh:76
G4bool createBullet(const G4HadProjectile &aTrack)
void SetVerboseLevel(G4int verbose)
G4bool retryInelasticProton() const
const G4LorentzVector & Get4Momentum() const
void SetVerboseLevel(G4int value)
void SetEnergyMomentumCheckLevels(G4double relativeLevel, G4double absoluteLevel)
const G4ParticleDefinition * GetDefinition() const
void usePreCompoundDeexcitation()
static void DumpConfiguration(std::ostream &os)
static G4UnboundPN * Definition()
Definition: G4UnboundPN.cc:67
G4GLOB_DLL std::ostream G4cerr
G4bool IsApplicable(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
G4HadFinalState * NoInteraction(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
static G4bool usePreCompound()