00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099 #include <cmath>
00100 #include <iostream>
00101
00102 #include "G4CascadeInterface.hh"
00103 #include "globals.hh"
00104 #include "G4SystemOfUnits.hh"
00105 #include "G4CascadeChannelTables.hh"
00106 #include "G4CascadeCheckBalance.hh"
00107 #include "G4CascadeParameters.hh"
00108 #include "G4CollisionOutput.hh"
00109 #include "G4DecayKineticTracks.hh"
00110 #include "G4DynamicParticle.hh"
00111 #include "G4HadronicException.hh"
00112 #include "G4InuclCollider.hh"
00113 #include "G4InuclElementaryParticle.hh"
00114 #include "G4InuclNuclei.hh"
00115 #include "G4InuclParticle.hh"
00116 #include "G4InuclParticleNames.hh"
00117 #include "G4KaonZeroLong.hh"
00118 #include "G4KaonZeroShort.hh"
00119 #include "G4KineticTrack.hh"
00120 #include "G4KineticTrackVector.hh"
00121 #include "G4Nucleus.hh"
00122 #include "G4ParticleDefinition.hh"
00123 #include "G4ReactionProductVector.hh"
00124 #include "G4Track.hh"
00125 #include "G4V3DNucleus.hh"
00126
00127 using namespace G4InuclParticleNames;
00128
00129 typedef std::vector<G4InuclElementaryParticle>::const_iterator particleIterator;
00130 typedef std::vector<G4InuclNuclei>::const_iterator nucleiIterator;
00131
00132
00133
00134
00135 const G4String G4CascadeInterface::randomFile = G4CascadeParameters::randomFile();
00136
00137
00138
00139 const G4int G4CascadeInterface::maximumTries = 20;
00140
00141
00142
00143
00144 G4CascadeInterface::G4CascadeInterface(const G4String& name)
00145 : G4VIntraNuclearTransportModel(name), numberOfTries(0),
00146 collider(new G4InuclCollider), balance(new G4CascadeCheckBalance(name)),
00147 bullet(0), target(0), output(new G4CollisionOutput) {
00148 SetEnergyMomentumCheckLevels(5*perCent, 10*MeV);
00149 balance->setLimits(5*perCent, 10*MeV/GeV);
00150
00151
00152 this->SetVerboseLevel(G4CascadeParameters::verbose());
00153 if (G4CascadeParameters::usePreCompound()) usePreCompoundDeexcitation();
00154 }
00155
00156
00157 G4CascadeInterface::~G4CascadeInterface() {
00158 clear();
00159 delete collider; collider=0;
00160 delete balance; balance=0;
00161 delete output; output=0;
00162 }
00163
00164 void G4CascadeInterface::ModelDescription(std::ostream& outFile) const
00165 {
00166 outFile << "The Bertini-style cascade implements the inelastic scattering\n"
00167 << "of hadrons by nuclei. Nucleons, pions, kaons and hyperons\n"
00168 << "from 0 to 15 GeV may be used as projectiles in this model.\n"
00169 << "Final state hadrons are produced by a classical cascade\n"
00170 << "consisting of individual hadron-nucleon scatterings which use\n"
00171 << "free-space partial cross sections, corrected for various\n"
00172 << "nuclear medium effects. The target nucleus is modeled as a\n"
00173 << "set of 1, 3 or 6 spherical shells, in which scattered hadrons\n"
00174 << "travel in straight lines until they are reflected from or\n"
00175 << "transmitted through shell boundaries.\n";
00176 }
00177
00178 void G4CascadeInterface::DumpConfiguration(std::ostream& outFile) const {
00179 G4CascadeParameters::DumpConfiguration(outFile);
00180 }
00181
00182
00183 void G4CascadeInterface::clear() {
00184 bullet=0;
00185 target=0;
00186 }
00187
00188
00189
00190
00191
00192 void G4CascadeInterface::useCascadeDeexcitation() {
00193 collider->useCascadeDeexcitation();
00194 }
00195
00196 void G4CascadeInterface::usePreCompoundDeexcitation() {
00197 collider->usePreCompoundDeexcitation();
00198 }
00199
00200
00201
00202
00203 void G4CascadeInterface::SetVerboseLevel(G4int verbose) {
00204 G4HadronicInteraction::SetVerboseLevel(verbose);
00205 collider->setVerboseLevel(verbose);
00206 balance->setVerboseLevel(verbose);
00207 output->setVerboseLevel(verbose);
00208 }
00209
00210
00211
00212
00213 G4bool G4CascadeInterface::IsApplicable(const G4HadProjectile& aTrack,
00214 G4Nucleus& ) {
00215 return IsApplicable(aTrack.GetDefinition());
00216 }
00217
00218 G4bool G4CascadeInterface::IsApplicable(const G4ParticleDefinition* aPD) const {
00219 if (aPD->GetAtomicMass() > 1) return true;
00220
00221
00222 G4int type = G4InuclElementaryParticle::type(aPD);
00223 return (type>0 && G4CascadeChannelTables::GetTable(type));
00224 }
00225
00226
00227
00228
00229 G4HadFinalState*
00230 G4CascadeInterface::ApplyYourself(const G4HadProjectile& aTrack,
00231 G4Nucleus& theNucleus) {
00232 if (verboseLevel)
00233 G4cout << " >>> G4CascadeInterface::ApplyYourself" << G4endl;
00234
00235 if (aTrack.GetKineticEnergy() < 0.) {
00236 G4cerr << " >>> G4CascadeInterface got negative-energy track: "
00237 << aTrack.GetDefinition()->GetParticleName() << " Ekin = "
00238 << aTrack.GetKineticEnergy() << G4endl;
00239 }
00240
00241 #ifdef G4CASCADE_DEBUG_INTERFACE
00242 static G4int counter(0);
00243 counter++;
00244 G4cerr << "Reaction number "<< counter << " "
00245 << aTrack.GetDefinition()->GetParticleName() << " "
00246 << aTrack.GetKineticEnergy() << " MeV" << G4endl;
00247 #endif
00248
00249 if (!randomFile.empty()) {
00250 if (verboseLevel>1)
00251 G4cout << " Saving random engine state to " << randomFile << G4endl;
00252 CLHEP::HepRandom::saveEngineStatus(randomFile);
00253 }
00254
00255 theParticleChange.Clear();
00256 clear();
00257
00258
00259 if (!IsApplicable(aTrack, theNucleus)) {
00260 if (verboseLevel) G4cerr << " No interaction possible " << G4endl;
00261 return NoInteraction(aTrack, theNucleus);
00262 }
00263
00264
00265 if (!createBullet(aTrack)) {
00266 if (verboseLevel) G4cerr << " Unable to create usable bullet" << G4endl;
00267 return NoInteraction(aTrack, theNucleus);
00268 }
00269
00270 if (!createTarget(theNucleus)) {
00271 if (verboseLevel) G4cerr << " Unable to create usable target" << G4endl;
00272 return NoInteraction(aTrack, theNucleus);
00273 }
00274
00275
00276 const G4bool isHydrogen = (theNucleus.GetA_asInt() == 1);
00277
00278 numberOfTries = 0;
00279 do {
00280 if (verboseLevel > 1)
00281 G4cout << " Generating cascade attempt " << numberOfTries << G4endl;
00282
00283 output->reset();
00284 collider->collide(bullet, target, *output);
00285 balance->collide(bullet, target, *output);
00286
00287 numberOfTries++;
00288 } while ( isHydrogen ? retryInelasticProton() : retryInelasticNucleus() );
00289
00290
00291 if (numberOfTries >= maximumTries) {
00292 if (verboseLevel)
00293 G4cout << " Cascade aborted after trials " << numberOfTries << G4endl;
00294 return NoInteraction(aTrack, theNucleus);
00295 }
00296
00297
00298 if (!balance->okay()) {
00299 throwNonConservationFailure();
00300 return NoInteraction(aTrack, theNucleus);
00301 }
00302
00303
00304 if (verboseLevel) {
00305 G4cout << " Cascade output after trials " << numberOfTries << G4endl;
00306 if (verboseLevel > 1) output->printCollisionOutput();
00307 }
00308
00309
00310 output->rotateEvent(bulletInLabFrame);
00311
00312 copyOutputToHadronicResult();
00313
00314
00315 checkFinalResult();
00316
00317
00318 clear();
00319 return &theParticleChange;
00320 }
00321
00322 G4ReactionProductVector*
00323 G4CascadeInterface::Propagate(G4KineticTrackVector* theSecondaries,
00324 G4V3DNucleus* theNucleus) {
00325 if (verboseLevel) G4cout << " >>> G4CascadeInterface::Propagate" << G4endl;
00326
00327 #ifdef G4CASCADE_DEBUG_INTERFACE
00328 if (verboseLevel>1) {
00329 G4cout << " G4V3DNucleus A " << theNucleus->GetMassNumber()
00330 << " Z " << theNucleus->GetCharge()
00331 << "\n " << theSecondaries->size() << " secondaries:" << G4endl;
00332 for (size_t i=0; i<theSecondaries->size(); i++) {
00333 G4KineticTrack* kt = (*theSecondaries)[i];
00334 G4cout << " " << i << ": " << kt->GetDefinition()->GetParticleName()
00335 << " p " << kt->Get4Momentum() << " @ " << kt->GetPosition()
00336 << " t " << kt->GetFormationTime() << G4endl;
00337 }
00338 }
00339 #endif
00340
00341 if (!randomFile.empty()) {
00342 if (verboseLevel>1)
00343 G4cout << " Saving random engine state to " << randomFile << G4endl;
00344 CLHEP::HepRandom::saveEngineStatus(randomFile);
00345 }
00346
00347 theParticleChange.Clear();
00348 clear();
00349
00350
00351 G4DecayKineticTracks decay(theSecondaries);
00352
00353
00354 const G4HadProjectile* projectile = GetPrimaryProjectile();
00355 if (projectile) createBullet(*projectile);
00356
00357 if (!createTarget(theNucleus)) {
00358 if (verboseLevel) G4cerr << " Unable to create usable target" << G4endl;
00359 return 0;
00360 }
00361
00362 numberOfTries = 0;
00363 do {
00364 if (verboseLevel > 1)
00365 G4cout << " Generating rescatter attempt " << numberOfTries << G4endl;
00366
00367 output->reset();
00368 collider->rescatter(bullet, theSecondaries, theNucleus, *output);
00369 balance->collide(bullet, target, *output);
00370
00371 numberOfTries++;
00372
00373 } while (retryInelasticNucleus());
00374
00375
00376 if (numberOfTries >= maximumTries && !balance->okay()) {
00377 throwNonConservationFailure();
00378 }
00379
00380
00381 if (verboseLevel) {
00382 G4cout << " Cascade rescatter after trials " << numberOfTries << G4endl;
00383 if (verboseLevel > 1) output->printCollisionOutput();
00384 }
00385
00386
00387 G4ReactionProductVector* propResult = copyOutputToReactionProducts();
00388
00389
00390 clear();
00391 return propResult;
00392 }
00393
00394
00395
00396
00397 G4HadFinalState*
00398 G4CascadeInterface::NoInteraction(const G4HadProjectile& aTrack,
00399 G4Nucleus& ) {
00400 if (verboseLevel)
00401 G4cout << " >>> G4CascadeInterface::NoInteraction" << G4endl;
00402
00403 theParticleChange.Clear();
00404 theParticleChange.SetStatusChange(isAlive);
00405
00406 G4double ekin = aTrack.GetKineticEnergy()>0. ? aTrack.GetKineticEnergy() : 0.;
00407 theParticleChange.SetEnergyChange(ekin);
00408
00409 return &theParticleChange;
00410 }
00411
00412
00413
00414
00415 G4bool G4CascadeInterface::createBullet(const G4HadProjectile& aTrack) {
00416 const G4ParticleDefinition* trkDef = aTrack.GetDefinition();
00417
00418 G4int bulletType = 0;
00419 G4int bulletA = 0, bulletZ = 0;
00420
00421 if (trkDef->GetAtomicMass() <= 1) {
00422 bulletType = G4InuclElementaryParticle::type(trkDef);
00423 } else {
00424 bulletA = trkDef->GetAtomicMass();
00425 bulletZ = trkDef->GetAtomicNumber();
00426 }
00427
00428 if (0 == bulletType && 0 == bulletA*bulletZ) {
00429 if (verboseLevel) {
00430 G4cerr << " G4CascadeInterface: " << trkDef->GetParticleName()
00431 << " not usable as bullet." << G4endl;
00432 }
00433 bullet = 0;
00434 return false;
00435 }
00436
00437
00438 G4LorentzVector projectileMomentum = aTrack.Get4Momentum()/GeV;
00439
00440
00441 bulletInLabFrame = G4LorentzRotation::IDENTITY;
00442 bulletInLabFrame.rotateZ(-projectileMomentum.phi());
00443 bulletInLabFrame.rotateY(-projectileMomentum.theta());
00444 bulletInLabFrame.invert();
00445
00446 G4LorentzVector momentumBullet(0., 0., projectileMomentum.rho(),
00447 projectileMomentum.e());
00448
00449 if (bulletType > 0) {
00450 hadronBullet.fill(momentumBullet, bulletType);
00451 bullet = &hadronBullet;
00452 } else {
00453 nucleusBullet.fill(momentumBullet, bulletA, bulletZ);
00454 bullet = &nucleusBullet;
00455 }
00456
00457 if (verboseLevel > 2) G4cout << "Bullet: \n" << *bullet << G4endl;
00458
00459 return true;
00460 }
00461
00462
00463
00464
00465 G4bool G4CascadeInterface::createTarget(G4Nucleus& theNucleus) {
00466 return createTarget(theNucleus.GetA_asInt(), theNucleus.GetZ_asInt());
00467 }
00468
00469 G4bool G4CascadeInterface::createTarget(G4V3DNucleus* theNucleus) {
00470 return createTarget(theNucleus->GetMassNumber(), theNucleus->GetCharge());
00471 }
00472
00473 G4bool G4CascadeInterface::createTarget(G4int A, G4int Z) {
00474 if (A > 1) {
00475 nucleusTarget.fill(A, Z);
00476 target = &nucleusTarget;
00477 } else {
00478 hadronTarget.fill(0., (Z=1?proton:neutron));
00479 target = &hadronTarget;
00480 }
00481
00482 if (verboseLevel > 2) G4cout << "Target: \n" << *target << G4endl;
00483
00484 return true;
00485 }
00486
00487
00488
00489
00490 G4DynamicParticle* G4CascadeInterface::
00491 makeDynamicParticle(const G4InuclElementaryParticle& iep) const {
00492 G4int outgoingType = iep.type();
00493
00494 if (iep.quasi_deutron()) {
00495 G4cerr << " ERROR: G4CascadeInterface incompatible particle type "
00496 << outgoingType << G4endl;
00497 return 0;
00498 }
00499
00500
00501 if (outgoingType == kaonZero || outgoingType == kaonZeroBar) {
00502 G4ThreeVector momDir = iep.getMomentum().vect().unit();
00503 G4double ekin = iep.getKineticEnergy()*GeV;
00504
00505 G4ParticleDefinition* pd = G4KaonZeroShort::Definition();
00506 if (G4UniformRand() > 0.5) pd = G4KaonZeroLong::Definition();
00507
00508 return new G4DynamicParticle(pd, momDir, ekin);
00509 } else {
00510 return new G4DynamicParticle(iep.getDynamicParticle());
00511 }
00512
00513 return 0;
00514 }
00515
00516 G4DynamicParticle*
00517 G4CascadeInterface::makeDynamicParticle(const G4InuclNuclei& inuc) const {
00518 if (verboseLevel > 2) {
00519 G4cout << " Nuclei fragment: \n" << inuc << G4endl;
00520 }
00521
00522
00523 return new G4DynamicParticle(inuc.getDynamicParticle());
00524 }
00525
00526
00527
00528
00529 void G4CascadeInterface::copyOutputToHadronicResult() {
00530 if (verboseLevel > 1)
00531 G4cout << " >>> G4CascadeInterface::copyOutputToHadronicResult" << G4endl;
00532
00533 const std::vector<G4InuclNuclei>& outgoingNuclei = output->getOutgoingNuclei();
00534 const std::vector<G4InuclElementaryParticle>& particles = output->getOutgoingParticles();
00535
00536 theParticleChange.SetStatusChange(stopAndKill);
00537 theParticleChange.SetEnergyChange(0.);
00538
00539
00540 if (!particles.empty()) {
00541 particleIterator ipart = particles.begin();
00542 for (; ipart != particles.end(); ipart++) {
00543 theParticleChange.AddSecondary(makeDynamicParticle(*ipart));
00544 }
00545 }
00546
00547
00548 if (!outgoingNuclei.empty()) {
00549 nucleiIterator ifrag = outgoingNuclei.begin();
00550 for (; ifrag != outgoingNuclei.end(); ifrag++) {
00551 theParticleChange.AddSecondary(makeDynamicParticle(*ifrag));
00552 }
00553 }
00554 }
00555
00556 G4ReactionProductVector* G4CascadeInterface::copyOutputToReactionProducts() {
00557 if (verboseLevel > 1)
00558 G4cout << " >>> G4CascadeInterface::copyOutputToReactionProducts" << G4endl;
00559
00560 const std::vector<G4InuclElementaryParticle>& particles = output->getOutgoingParticles();
00561 const std::vector<G4InuclNuclei>& fragments = output->getOutgoingNuclei();
00562
00563 G4ReactionProductVector* propResult = new G4ReactionProductVector;
00564
00565 G4ReactionProduct* rp = 0;
00566 G4DynamicParticle* dp = 0;
00567
00568
00569 if (!particles.empty()) {
00570 particleIterator ipart = particles.begin();
00571 for (; ipart != particles.end(); ipart++) {
00572 rp = new G4ReactionProduct;
00573 dp = makeDynamicParticle(*ipart);
00574 (*rp) = (*dp);
00575 propResult->push_back(rp);
00576 delete dp;
00577 }
00578 }
00579
00580
00581 if (!fragments.empty()) {
00582 nucleiIterator ifrag = fragments.begin();
00583 for (; ifrag != fragments.end(); ifrag++) {
00584 rp = new G4ReactionProduct;
00585 dp = makeDynamicParticle(*ifrag);
00586 (*rp) = (*dp);
00587 propResult->push_back(rp);
00588 delete dp;
00589 }
00590 }
00591
00592 return propResult;
00593 }
00594
00595
00596
00597
00598 void G4CascadeInterface::checkFinalResult() {
00599 balance->collide(bullet, target, *output);
00600
00601 if (verboseLevel > 2) {
00602 if (!balance->baryonOkay()) {
00603 G4cerr << "ERROR: no baryon number conservation, sum of baryons = "
00604 << balance->deltaB() << G4endl;
00605 }
00606
00607 if (!balance->chargeOkay()) {
00608 G4cerr << "ERROR: no charge conservation, sum of charges = "
00609 << balance->deltaQ() << G4endl;
00610 }
00611
00612 if (std::abs(balance->deltaKE()) > 0.01 ) {
00613 G4cerr << "Kinetic energy conservation violated by "
00614 << balance->deltaKE() << " GeV" << G4endl;
00615 }
00616
00617 G4double eInit = bullet->getEnergy() + target->getEnergy();
00618 G4double eFinal = eInit + balance->deltaE();
00619
00620 G4cout << "Initial energy " << eInit << " final energy " << eFinal
00621 << "\nTotal energy conservation at level "
00622 << balance->deltaE() * GeV << " MeV" << G4endl;
00623
00624 if (balance->deltaKE() > 5.0e-5 ) {
00625 G4cerr << "FATAL ERROR: kinetic energy created "
00626 << balance->deltaKE() * GeV << " MeV" << G4endl;
00627 }
00628 }
00629 }
00630
00631
00632
00633
00634 G4bool G4CascadeInterface::coulombBarrierViolation() const {
00635 G4bool violated = false;
00636
00637 const G4double coulumbBarrier = 8.7 * MeV/GeV;
00638
00639 const std::vector<G4InuclElementaryParticle>& p =
00640 output->getOutgoingParticles();
00641
00642 for (particleIterator ipart=p.begin(); ipart != p.end(); ipart++) {
00643 if (ipart->type() == proton) {
00644 violated |= (ipart->getKineticEnergy() < coulumbBarrier);
00645 }
00646 }
00647
00648 return violated;
00649 }
00650
00651
00652
00653 G4bool G4CascadeInterface::retryInelasticProton() const {
00654 const std::vector<G4InuclElementaryParticle>& out =
00655 output->getOutgoingParticles();
00656
00657 #ifdef G4CASCADE_DEBUG_INTERFACE
00658
00659 G4cout << " retryInelasticProton: number of Tries "
00660 << ((numberOfTries < maximumTries) ? "RETRY (t)" : "EXIT (f)")
00661 << "\n retryInelasticProton: AND collision type ";
00662 if (out.empty()) G4cout << "FAILED" << G4endl;
00663 else {
00664 G4cout << (out.size() == 2 ? "ELASTIC (t)" : "INELASTIC (f)")
00665 << "\n retryInelasticProton: AND Leading particles bullet "
00666 << (out.size() >= 2 &&
00667 (out[0].getDefinition() == bullet->getDefinition() ||
00668 out[1].getDefinition() == bullet->getDefinition())
00669 ? "YES (t)" : "NO (f)")
00670 << G4endl;
00671 }
00672 #endif
00673
00674 return ( (numberOfTries < maximumTries) &&
00675 (out.empty() ||
00676 (out.size() == 2 &&
00677 (out[0].getDefinition() == bullet->getDefinition() ||
00678 out[1].getDefinition() == bullet->getDefinition())))
00679 );
00680 }
00681
00682
00683
00684
00685 G4bool G4CascadeInterface::retryInelasticNucleus() const {
00686
00687 G4int npart = output->numberOfOutgoingParticles();
00688 G4int nfrag = output->numberOfOutgoingNuclei();
00689
00690 const G4ParticleDefinition* firstOut = (npart == 0) ? 0 :
00691 output->getOutgoingParticles().begin()->getDefinition();
00692
00693 #ifdef G4CASCADE_DEBUG_INTERFACE
00694
00695 G4cout << " retryInelasticNucleus: numberOfTries "
00696 << ((numberOfTries < maximumTries) ? "RETRY (t)" : "EXIT (f)")
00697 << "\n retryInelasticNucleus: AND outputParticles "
00698 << ((npart != 0) ? "NON-ZERO (t)" : "EMPTY (f)")
00699 #ifdef G4CASCADE_COULOMB_DEV
00700 << "\n retryInelasticNucleus: AND coulombBarrier (COULOMB_DEV) "
00701 << (coulombBarrierViolation() ? "VIOLATED (t)" : "PASSED (f)")
00702 << "\n retryInelasticNucleus: AND collision type (COULOMB_DEV) "
00703 << ((npart+nfrag > 2) ? "INELASTIC (t)" : "ELASTIC (f)")
00704 #else
00705 << "\n retryInelasticNucleus: AND collsion type "
00706 << ((npart+nfrag < 3) ? "ELASTIC (t)" : "INELASTIC (f)")
00707 << "\n retryInelasticNucleus: AND Leading particle bullet "
00708 << ((firstOut == bullet->getDefinition()) ? "YES (t)" : "NO (f)")
00709 #endif
00710 << "\n retryInelasticNucleus: OR conservation "
00711 << (!balance->okay() ? "FAILED (t)" : "PASSED (f)")
00712 << G4endl;
00713 #endif
00714
00715 return ( ((numberOfTries < maximumTries) &&
00716 (npart != 0) &&
00717 #ifdef G4CASCADE_COULOMB_DEV
00718 (coulombBarrierViolation() && npart+nfrag > 2)
00719 #else
00720 (npart+nfrag < 3 && firstOut == bullet->getDefinition())
00721 #endif
00722 )
00723 #ifndef G4CASCADE_SKIP_ECONS
00724 || (!balance->okay())
00725 #endif
00726 );
00727 }
00728
00729
00730
00731
00732
00733 void G4CascadeInterface::throwNonConservationFailure() {
00734
00735
00736 std::ostream& errInfo = G4cerr;
00737
00738 errInfo << " >>> G4CascadeInterface has non-conserving"
00739 << " cascade after " << numberOfTries << " attempts." << G4endl;
00740
00741 G4String throwMsg = "G4CascadeInterface - ";
00742 if (!balance->energyOkay()) {
00743 throwMsg += "Energy";
00744 errInfo << " Energy conservation violated by " << balance->deltaE()
00745 << " GeV (" << balance->relativeE() << ")" << G4endl;
00746 }
00747
00748 if (!balance->momentumOkay()) {
00749 throwMsg += "Momentum";
00750 errInfo << " Momentum conservation violated by " << balance->deltaP()
00751 << " GeV/c (" << balance->relativeP() << ")" << G4endl;
00752 }
00753
00754 if (!balance->baryonOkay()) {
00755 throwMsg += "Baryon number";
00756 errInfo << " Baryon number violated by " << balance->deltaB() << G4endl;
00757 }
00758
00759 if (!balance->chargeOkay()) {
00760 throwMsg += "Charge";
00761 errInfo << " Charge conservation violated by " << balance->deltaQ()
00762 << G4endl;
00763 }
00764
00765 errInfo << " Final event output, for debugging:\n"
00766 << " Bullet: \n" << *bullet << G4endl
00767 << " Target: \n" << *target << G4endl;
00768 output->printCollisionOutput(errInfo);
00769
00770 throwMsg += " non-conservation. More info in output.";
00771 throw G4HadronicException(__FILE__, __LINE__, throwMsg);
00772 }