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