Geant4_10
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 
101 #include <cmath>
102 #include <iostream>
103 
104 #include "G4CascadeInterface.hh"
105 #include "globals.hh"
106 #include "G4SystemOfUnits.hh"
107 #include "G4CascadeChannelTables.hh"
108 #include "G4CascadeCheckBalance.hh"
109 #include "G4CascadeParameters.hh"
110 #include "G4CollisionOutput.hh"
111 #include "G4DecayKineticTracks.hh"
112 #include "G4DynamicParticle.hh"
113 #include "G4HadronicException.hh"
114 #include "G4InuclCollider.hh"
116 #include "G4InuclNuclei.hh"
117 #include "G4InuclParticle.hh"
118 #include "G4InuclParticleNames.hh"
119 #include "G4KaonZeroLong.hh"
120 #include "G4KaonZeroShort.hh"
121 #include "G4KineticTrack.hh"
122 #include "G4KineticTrackVector.hh"
123 #include "G4Nucleus.hh"
124 #include "G4ParticleDefinition.hh"
126 #include "G4Track.hh"
127 #include "G4V3DNucleus.hh"
128 #include "G4UnboundPN.hh"
129 #include "G4Dineutron.hh"
130 #include "G4Diproton.hh"
131 
132 using namespace G4InuclParticleNames;
133 
134 typedef std::vector<G4InuclElementaryParticle>::const_iterator particleIterator;
135 typedef std::vector<G4InuclNuclei>::const_iterator nucleiIterator;
136 
137 
138 // Filename to capture random seed, if specified by user at runtime
139 
140 const G4String G4CascadeInterface::randomFile = G4CascadeParameters::randomFile();
141 
142 // Maximum number of iterations allowed for inelastic collision attempts
143 
144 const G4int G4CascadeInterface::maximumTries = 20;
145 
146 
147 // Constructor and destrutor
148 
150  : G4VIntraNuclearTransportModel(name), numberOfTries(0),
151  collider(new G4InuclCollider), balance(new G4CascadeCheckBalance(name)),
152  bullet(0), target(0), output(new G4CollisionOutput) {
154  balance->setLimits(5*perCent, 10*MeV/GeV); // Bertini internal units
155  // Description(); // Model description
156 
159 }
160 
162  clear();
163  delete collider; collider=0;
164  delete balance; balance=0;
165  delete output; output=0;
166 }
167 
169 {
170  outFile << "The Bertini-style cascade implements the inelastic scattering\n"
171  << "of hadrons by nuclei. Nucleons, pions, kaons and hyperons\n"
172  << "from 0 to 15 GeV may be used as projectiles in this model.\n"
173  << "Final state hadrons are produced by a classical cascade\n"
174  << "consisting of individual hadron-nucleon scatterings which use\n"
175  << "free-space partial cross sections, corrected for various\n"
176  << "nuclear medium effects. The target nucleus is modeled as a\n"
177  << "set of 1, 3 or 6 spherical shells, in which scattered hadrons\n"
178  << "travel in straight lines until they are reflected from or\n"
179  << "transmitted through shell boundaries.\n";
180 }
181 
184 }
185 
187  bullet=0;
188  target=0;
189 }
190 
191 
192 // Initialize shared objects for use across multiple threads
193 
199  if (!ch || !pn || !nn || !pp) return; // Avoid "unused variables"
200 }
201 
202 
203 // Select post-cascade processing (default will be CascadeDeexcitation)
204 // NOTE: Currently just calls through to Collider, in future will do something
205 
207  collider->useCascadeDeexcitation();
208 }
209 
211  collider->usePreCompoundDeexcitation();
212 }
213 
214 
215 // Apply verbosity to all member objects (overrides base class version)
216 
219  collider->setVerboseLevel(verbose);
220  balance->setVerboseLevel(verbose);
221  output->setVerboseLevel(verbose);
222 }
223 
224 
225 // Test whether inputs are valid for this model
226 
228  G4Nucleus& /* theNucleus */) {
229  return IsApplicable(aTrack.GetDefinition());
230 }
231 
233  if (aPD->GetAtomicMass() > 1) return true; // Nuclei are okay
234 
235  // Valid particle and have interactions available
237  return (G4CascadeChannelTables::GetTable(type));
238 }
239 
240 
241 // Main Actions
242 
245  G4Nucleus& theNucleus) {
246  if (verboseLevel)
247  G4cout << " >>> G4CascadeInterface::ApplyYourself" << G4endl;
248 
249  if (aTrack.GetKineticEnergy() < 0.) {
250  G4cerr << " >>> G4CascadeInterface got negative-energy track: "
251  << aTrack.GetDefinition()->GetParticleName() << " Ekin = "
252  << aTrack.GetKineticEnergy() << G4endl;
253  }
254 
255 #ifdef G4CASCADE_DEBUG_INTERFACE
256  static G4int counter(0);
257  counter++;
258  G4cerr << "Reaction number "<< counter << " "
259  << aTrack.GetDefinition()->GetParticleName() << " "
260  << aTrack.GetKineticEnergy() << " MeV" << G4endl;
261 #endif
262 
263  if (!randomFile.empty()) { // User requested random-seed capture
264  if (verboseLevel>1)
265  G4cout << " Saving random engine state to " << randomFile << G4endl;
267  }
268 
270  clear();
271 
272  // Abort processing if no interaction is possible
273  if (!IsApplicable(aTrack, theNucleus)) {
274  if (verboseLevel) G4cerr << " No interaction possible " << G4endl;
275  return NoInteraction(aTrack, theNucleus);
276  }
277 
278  // Make conversion between native Geant4 and Bertini cascade classes.
279  if (!createBullet(aTrack)) {
280  if (verboseLevel) G4cerr << " Unable to create usable bullet" << G4endl;
281  return NoInteraction(aTrack, theNucleus);
282  }
283 
284  if (!createTarget(theNucleus)) {
285  if (verboseLevel) G4cerr << " Unable to create usable target" << G4endl;
286  return NoInteraction(aTrack, theNucleus);
287  }
288 
289  // Different retry conditions for proton target vs. nucleus
290  const G4bool isHydrogen = (theNucleus.GetA_asInt() == 1);
291 
292  numberOfTries = 0;
293  do { // we try to create inelastic interaction
294  if (verboseLevel > 1)
295  G4cout << " Generating cascade attempt " << numberOfTries << G4endl;
296 
297  output->reset();
298  collider->collide(bullet, target, *output);
299  balance->collide(bullet, target, *output);
300 
301  numberOfTries++;
302  } while ( isHydrogen ? retryInelasticProton() : retryInelasticNucleus() );
303 
304  // Null event if unsuccessful
305  if (numberOfTries >= maximumTries) {
306  if (verboseLevel)
307  G4cout << " Cascade aborted after trials " << numberOfTries << G4endl;
308  return NoInteraction(aTrack, theNucleus);
309  }
310 
311  // Abort job if energy or momentum are not conserved
312  if (!balance->okay()) {
314  return NoInteraction(aTrack, theNucleus);
315  }
316 
317  // Successful cascade -- clean up and return
318  if (verboseLevel) {
319  G4cout << " Cascade output after trials " << numberOfTries << G4endl;
320  if (verboseLevel > 1) output->printCollisionOutput();
321  }
322 
323  // Rotate event to put Z axis along original projectile direction
324  output->rotateEvent(bulletInLabFrame);
325 
327 
328  // Report violations of conservation laws in original frame
330 
331  // Clean up and return final result;
332  clear();
333 /*
334  G4int nSec = theParticleChange.GetNumberOfSecondaries();
335  for (G4int i = 0; i < nSec; i++) {
336  G4HadSecondary* sec = theParticleChange.GetSecondary(i);
337  G4DynamicParticle* dp = sec->GetParticle();
338  if (dp->GetDefinition()->GetParticleName() == "neutron")
339  G4cout << dp->GetDefinition()->GetParticleName() << " has "
340  << dp->GetKineticEnergy()/MeV << " MeV " << G4endl;
341  }
342 */
343  return &theParticleChange;
344 }
345 
348  G4V3DNucleus* theNucleus) {
349  if (verboseLevel) G4cout << " >>> G4CascadeInterface::Propagate" << G4endl;
350 
351 #ifdef G4CASCADE_DEBUG_INTERFACE
352  if (verboseLevel>1) {
353  G4cout << " G4V3DNucleus A " << theNucleus->GetMassNumber()
354  << " Z " << theNucleus->GetCharge()
355  << "\n " << theSecondaries->size() << " secondaries:" << G4endl;
356  for (size_t i=0; i<theSecondaries->size(); i++) {
357  G4KineticTrack* kt = (*theSecondaries)[i];
358  G4cout << " " << i << ": " << kt->GetDefinition()->GetParticleName()
359  << " p " << kt->Get4Momentum() << " @ " << kt->GetPosition()
360  << " t " << kt->GetFormationTime() << G4endl;
361  }
362  }
363 #endif
364 
365  if (!randomFile.empty()) { // User requested random-seed capture
366  if (verboseLevel>1)
367  G4cout << " Saving random engine state to " << randomFile << G4endl;
369  }
370 
372  clear();
373 
374  // Process input secondaries list to eliminate resonances
375  G4DecayKineticTracks decay(theSecondaries);
376 
377  // NOTE: Requires 9.4-ref-03 mods to base class and G4TheoFSGenerator
378  const G4HadProjectile* projectile = GetPrimaryProjectile();
379  if (projectile) createBullet(*projectile);
380 
381  if (!createTarget(theNucleus)) {
382  if (verboseLevel) G4cerr << " Unable to create usable target" << G4endl;
383  return 0; // FIXME: This will cause a segfault later
384  }
385 
386  numberOfTries = 0;
387  do {
388  if (verboseLevel > 1)
389  G4cout << " Generating rescatter attempt " << numberOfTries << G4endl;
390 
391  output->reset();
392  collider->rescatter(bullet, theSecondaries, theNucleus, *output);
393  balance->collide(bullet, target, *output);
394 
395  numberOfTries++;
396  // FIXME: retry checks will SEGFAULT until we can define the bullet!
397  } while (retryInelasticNucleus());
398 
399  // Check whether repeated attempts have all failed; report and exit
400  if (numberOfTries >= maximumTries && !balance->okay()) {
401  throwNonConservationFailure(); // This terminates the job
402  }
403 
404  // Successful cascade -- clean up and return
405  if (verboseLevel) {
406  G4cout << " Cascade rescatter after trials " << numberOfTries << G4endl;
407  if (verboseLevel > 1) output->printCollisionOutput();
408  }
409 
410  // Does calling code take ownership? I hope so!
412 
413  // Clean up and and return final result
414  clear();
415  return propResult;
416 }
417 
418 
419 // Replicate input particles onto output
420 
423  G4Nucleus& /*theNucleus*/) {
424  if (verboseLevel)
425  G4cout << " >>> G4CascadeInterface::NoInteraction" << G4endl;
426 
429 
430  G4double ekin = aTrack.GetKineticEnergy()>0. ? aTrack.GetKineticEnergy() : 0.;
431  theParticleChange.SetEnergyChange(ekin); // Protect against rounding
432 
433  return &theParticleChange;
434 }
435 
436 
437 // Convert input projectile to Bertini internal object
438 
440  const G4ParticleDefinition* trkDef = aTrack.GetDefinition();
441  G4int bulletType = 0; // For elementary particles
442  G4int bulletA = 0, bulletZ = 0; // For nucleus projectile
443 
444  if (trkDef->GetAtomicMass() <= 1) {
445  bulletType = G4InuclElementaryParticle::type(trkDef);
446  } else {
447  bulletA = trkDef->GetAtomicMass();
448  bulletZ = trkDef->GetAtomicNumber();
449  }
450 
451  if (0 == bulletType && 0 == bulletA*bulletZ) {
452  if (verboseLevel) {
453  G4cerr << " G4CascadeInterface: " << trkDef->GetParticleName()
454  << " not usable as bullet." << G4endl;
455  }
456  bullet = 0;
457  return false;
458  }
459 
460  // Code momentum and energy -- Bertini wants z-axis and GeV units
461  G4LorentzVector projectileMomentum = aTrack.Get4Momentum()/GeV;
462 
463  // Rotation/boost to get from z-axis back to original frame
464  bulletInLabFrame = G4LorentzRotation::IDENTITY; // Initialize
465  bulletInLabFrame.rotateZ(-projectileMomentum.phi());
466  bulletInLabFrame.rotateY(-projectileMomentum.theta());
467  bulletInLabFrame.invert();
468 
469  G4LorentzVector momentumBullet(0., 0., projectileMomentum.rho(),
470  projectileMomentum.e());
471 
472  if (G4InuclElementaryParticle::valid(bulletType)) {
473  hadronBullet.fill(momentumBullet, bulletType);
474  bullet = &hadronBullet;
475  } else {
476  nucleusBullet.fill(momentumBullet, bulletA, bulletZ);
477  bullet = &nucleusBullet;
478  }
479 
480  if (verboseLevel > 2) G4cout << "Bullet: \n" << *bullet << G4endl;
481 
482  return true;
483 }
484 
485 
486 // Convert input nuclear target to Bertini internal object
487 
489  return createTarget(theNucleus.GetA_asInt(), theNucleus.GetZ_asInt());
490 }
491 
493  return createTarget(theNucleus->GetMassNumber(), theNucleus->GetCharge());
494 }
495 
497  if (A > 1) {
498  nucleusTarget.fill(A, Z);
499  target = &nucleusTarget;
500  } else {
501  hadronTarget.fill(0., (Z==1?proton:neutron));
502  target = &hadronTarget;
503  }
504 
505  if (verboseLevel > 2) G4cout << "Target: \n" << *target << G4endl;
506 
507  return true; // Right now, target never fails
508 }
509 
510 
511 // Convert Bertini particle to output (G4DynamicParticle)
512 
515  G4int outgoingType = iep.type();
516 
517  if (iep.quasi_deutron()) {
518  G4cerr << " ERROR: G4CascadeInterface incompatible particle type "
519  << outgoingType << G4endl;
520  return 0;
521  }
522 
523  // Copy local G4DynPart to public output (handle kaon mixing specially)
524  if (outgoingType == kaonZero || outgoingType == kaonZeroBar) {
525  G4ThreeVector momDir = iep.getMomentum().vect().unit();
526  G4double ekin = iep.getKineticEnergy()*GeV; // Bertini -> G4 units
527 
529  if (G4UniformRand() > 0.5) pd = G4KaonZeroLong::Definition();
530 
531  return new G4DynamicParticle(pd, momDir, ekin);
532  } else {
533  return new G4DynamicParticle(iep.getDynamicParticle());
534  }
535 
536  return 0; // Should never get here!
537 }
538 
541  if (verboseLevel > 2) {
542  G4cout << " Nuclei fragment: \n" << inuc << G4endl;
543  }
544 
545  // Copy local G4DynPart to public output
546  return new G4DynamicParticle(inuc.getDynamicParticle());
547 }
548 
549 
550 // Transfer Bertini internal final state to hadronics interface
551 
553  if (verboseLevel > 1)
554  G4cout << " >>> G4CascadeInterface::copyOutputToHadronicResult" << G4endl;
555 
556  const std::vector<G4InuclNuclei>& outgoingNuclei = output->getOutgoingNuclei();
557  const std::vector<G4InuclElementaryParticle>& particles = output->getOutgoingParticles();
558 
561 
562  // Get outcoming particles
563  if (!particles.empty()) {
564  particleIterator ipart = particles.begin();
565  for (; ipart != particles.end(); ipart++) {
567  }
568  }
569 
570  // get nuclei fragments
571  if (!outgoingNuclei.empty()) {
572  nucleiIterator ifrag = outgoingNuclei.begin();
573  for (; ifrag != outgoingNuclei.end(); ifrag++) {
575  }
576  }
577 }
578 
580  if (verboseLevel > 1)
581  G4cout << " >>> G4CascadeInterface::copyOutputToReactionProducts" << G4endl;
582 
583  const std::vector<G4InuclElementaryParticle>& particles = output->getOutgoingParticles();
584  const std::vector<G4InuclNuclei>& fragments = output->getOutgoingNuclei();
585 
587 
588  G4ReactionProduct* rp = 0; // Buffers to create outgoing tracks
589  G4DynamicParticle* dp = 0;
590 
591  // Get outcoming particles
592  if (!particles.empty()) {
593  particleIterator ipart = particles.begin();
594  for (; ipart != particles.end(); ipart++) {
595  rp = new G4ReactionProduct;
596  dp = makeDynamicParticle(*ipart);
597  (*rp) = (*dp); // This does all the necessary copying
598  propResult->push_back(rp);
599  delete dp;
600  }
601  }
602 
603  // get nuclei fragments
604  if (!fragments.empty()) {
605  nucleiIterator ifrag = fragments.begin();
606  for (; ifrag != fragments.end(); ifrag++) {
607  rp = new G4ReactionProduct;
608  dp = makeDynamicParticle(*ifrag);
609  (*rp) = (*dp); // This does all the necessary copying
610  propResult->push_back(rp);
611  delete dp;
612  }
613  }
614 
615  return propResult;
616 }
617 
618 
619 // Report violations of conservation laws in original frame
620 
622  balance->collide(bullet, target, *output);
623 
624  if (verboseLevel > 2) {
625  if (!balance->baryonOkay()) {
626  G4cerr << "ERROR: no baryon number conservation, sum of baryons = "
627  << balance->deltaB() << G4endl;
628  }
629 
630  if (!balance->chargeOkay()) {
631  G4cerr << "ERROR: no charge conservation, sum of charges = "
632  << balance->deltaQ() << G4endl;
633  }
634 
635  if (std::abs(balance->deltaKE()) > 0.01 ) { // GeV
636  G4cerr << "Kinetic energy conservation violated by "
637  << balance->deltaKE() << " GeV" << G4endl;
638  }
639 
640  G4double eInit = bullet->getEnergy() + target->getEnergy();
641  G4double eFinal = eInit + balance->deltaE();
642 
643  G4cout << "Initial energy " << eInit << " final energy " << eFinal
644  << "\nTotal energy conservation at level "
645  << balance->deltaE() * GeV << " MeV" << G4endl;
646 
647  if (balance->deltaKE() > 5.0e-5 ) { // 0.05 MeV
648  G4cerr << "FATAL ERROR: kinetic energy created "
649  << balance->deltaKE() * GeV << " MeV" << G4endl;
650  }
651  }
652 }
653 
654 
655 // Evaluate whether any outgoing particles penetrated Coulomb barrier
656 
658  G4bool violated = false; // by default coulomb analysis is OK
659 
660  const G4double coulumbBarrier = 8.7 * MeV/GeV; // Bertini uses GeV
661 
662  const std::vector<G4InuclElementaryParticle>& p =
663  output->getOutgoingParticles();
664 
665  for (particleIterator ipart=p.begin(); ipart != p.end(); ipart++) {
666  if (ipart->type() == proton) {
667  violated |= (ipart->getKineticEnergy() < coulumbBarrier);
668  }
669  }
670 
671  return violated;
672 }
673 
674 // Check whether inelastic collision on proton failed
675 
677  const std::vector<G4InuclElementaryParticle>& out =
678  output->getOutgoingParticles();
679 
680 #ifdef G4CASCADE_DEBUG_INTERFACE
681  // Report on all retry conditions, in order of return logic
682  G4cout << " retryInelasticProton: number of Tries "
683  << ((numberOfTries < maximumTries) ? "RETRY (t)" : "EXIT (f)")
684  << "\n retryInelasticProton: AND collision type ";
685  if (out.empty()) G4cout << "FAILED" << G4endl;
686  else {
687  G4cout << (out.size() == 2 ? "ELASTIC (t)" : "INELASTIC (f)")
688  << "\n retryInelasticProton: AND Leading particles bullet "
689  << (out.size() >= 2 &&
690  (out[0].getDefinition() == bullet->getDefinition() ||
691  out[1].getDefinition() == bullet->getDefinition())
692  ? "YES (t)" : "NO (f)")
693  << G4endl;
694  }
695 #endif
696 
697  return ( (numberOfTries < maximumTries) &&
698  (out.empty() ||
699  (out.size() == 2 &&
700  (out[0].getDefinition() == bullet->getDefinition() ||
701  out[1].getDefinition() == bullet->getDefinition())))
702  );
703 }
704 
705 // Check whether generic inelastic collision failed
706 // NOTE: some conditions are set by compiler flags
707 
709  // Quantities necessary for conditional block below
711  G4int nfrag = output->numberOfOutgoingNuclei();
712 
713  const G4ParticleDefinition* firstOut = (npart == 0) ? 0 :
714  output->getOutgoingParticles().begin()->getDefinition();
715 
716 #ifdef G4CASCADE_DEBUG_INTERFACE
717  // Report on all retry conditions, in order of return logic
718  G4cout << " retryInelasticNucleus: numberOfTries "
719  << ((numberOfTries < maximumTries) ? "RETRY (t)" : "EXIT (f)")
720  << "\n retryInelasticNucleus: AND outputParticles "
721  << ((npart != 0) ? "NON-ZERO (t)" : "EMPTY (f)")
722 #ifdef G4CASCADE_COULOMB_DEV
723  << "\n retryInelasticNucleus: AND coulombBarrier (COULOMB_DEV) "
724  << (coulombBarrierViolation() ? "VIOLATED (t)" : "PASSED (f)")
725  << "\n retryInelasticNucleus: AND collision type (COULOMB_DEV) "
726  << ((npart+nfrag > 2) ? "INELASTIC (t)" : "ELASTIC (f)")
727 #else
728  << "\n retryInelasticNucleus: AND collsion type "
729  << ((npart+nfrag < 3) ? "ELASTIC (t)" : "INELASTIC (f)")
730  << "\n retryInelasticNucleus: AND Leading particle bullet "
731  << ((firstOut == bullet->getDefinition()) ? "YES (t)" : "NO (f)")
732 #endif
733  << "\n retryInelasticNucleus: OR conservation "
734  << (!balance->okay() ? "FAILED (t)" : "PASSED (f)")
735  << G4endl;
736 #endif
737 
738  return ( ((numberOfTries < maximumTries) &&
739  (npart != 0) &&
740 #ifdef G4CASCADE_COULOMB_DEV
741  (coulombBarrierViolation() && npart+nfrag > 2)
742 #else
743  (npart+nfrag < 3 && firstOut == bullet->getDefinition())
744 #endif
745  )
746 #ifndef G4CASCADE_SKIP_ECONS
747  || (!balance->okay())
748 #endif
749  );
750 }
751 
752 
753 // Terminate job in case of persistent non-conservation
754 // FIXME: Need to migrate to G4ExceptionDescription
755 
757  // NOTE: Once G4HadronicException is changed, use the following line!
758  // G4ExceptionDescription errInfo;
759  std::ostream& errInfo = G4cerr;
760 
761  errInfo << " >>> G4CascadeInterface has non-conserving"
762  << " cascade after " << numberOfTries << " attempts." << G4endl;
763 
764  G4String throwMsg = "G4CascadeInterface - ";
765  if (!balance->energyOkay()) {
766  throwMsg += "Energy";
767  errInfo << " Energy conservation violated by " << balance->deltaE()
768  << " GeV (" << balance->relativeE() << ")" << G4endl;
769  }
770 
771  if (!balance->momentumOkay()) {
772  throwMsg += "Momentum";
773  errInfo << " Momentum conservation violated by " << balance->deltaP()
774  << " GeV/c (" << balance->relativeP() << ")" << G4endl;
775  }
776 
777  if (!balance->baryonOkay()) {
778  throwMsg += "Baryon number";
779  errInfo << " Baryon number violated by " << balance->deltaB() << G4endl;
780  }
781 
782  if (!balance->chargeOkay()) {
783  throwMsg += "Charge";
784  errInfo << " Charge conservation violated by " << balance->deltaQ()
785  << G4endl;
786  }
787 
788  errInfo << " Final event output, for debugging:\n"
789  << " Bullet: \n" << *bullet << G4endl
790  << " Target: \n" << *target << G4endl;
791  output->printCollisionOutput(errInfo);
792 
793  throwMsg += " non-conservation. More info in output.";
794  throw G4HadronicException(__FILE__, __LINE__, throwMsg); // Job ends here!
795 }
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
G4CascadeInterface(const G4String &name="BertiniCascade")
static const G4CascadeChannel * GetTable(G4int initialState)
G4bool coulombBarrierViolation() const
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
virtual G4int GetCharge()=0
std::ofstream outFile
Definition: GammaRayTel.cc:68
const G4DynamicParticle & getDynamicParticle() const
G4LorentzVector getMomentum() const
const G4HadProjectile * GetPrimaryProjectile() const
const G4ThreeVector & GetPosition() const
void collide(G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &output)
Int_t ipart
Definition: Style.C:10
const char * p
Definition: xmltok.h:285
void printCollisionOutput(std::ostream &os=G4cout) const
virtual G4int GetMassNumber()=0
const XML_Char * name
Definition: expat.h:151
G4double getEnergy() const
const XML_Char * target
Definition: expat.h:268
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
static const G4String & randomFile()
void useCascadeDeexcitation()
G4ParticleDefinition * GetDefinition() const
virtual void DumpConfiguration(std::ostream &outFile) const
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
const G4ParticleDefinition * GetDefinition() const
double theta() const
Float_t Z
Definition: plot.C:39
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:175
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
virtual void SetVerboseLevel(G4int value)
void fill(G4int ityp, Model model=DefaultModel)
std::vector< G4InuclNuclei >::const_iterator nucleiIterator
float perCent
Definition: hepunit.py:239
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:64
Int_t npart
Definition: Style.C:7
#define G4endl
Definition: G4ios.hh:61
virtual void ModelDescription(std::ostream &outFile) const
static DLL_API const HepLorentzRotation IDENTITY
double G4double
Definition: G4Types.hh:76
G4bool createBullet(const G4HadProjectile &aTrack)
void SetVerboseLevel(G4int verbose)
G4bool retryInelasticProton() const
const G4LorentzVector & Get4Momentum() const
void SetEnergyMomentumCheckLevels(G4double relativeLevel, G4double absoluteLevel)
G4ParticleDefinition * getDefinition() const
void usePreCompoundDeexcitation()
void AddSecondary(G4DynamicParticle *aP)
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()