G4CascadeCheckBalance.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // Verify and report four-momentum conservation for collision output; uses
00029 // same interface as collision generators.
00030 //
00031 // 20100624  M. Kelsey -- Add baryon number, charge, and kinetic energy
00032 // 20100624  M. Kelsey -- Bug fix:  All checks should be |delta| !
00033 // 20100627  M. Kelsey -- Always report violations on cerr, regardless of
00034 //              verbosity.
00035 // 20100628  M. Kelsey -- Add interface to take list of particles directly,
00036 //              bug fix reporting of charge conservation error.
00037 // 20100630  M. Kelsey -- for nuclei, include excitation energies in total.
00038 // 20100701  M. Kelsey -- Undo previous change, handled by G4InuclNuclei.
00039 // 20100711  M. Kelsey -- Use name of parent collider for reporting messages
00040 // 20100713  M. kelsey -- Hide conservation errors behind verbosity
00041 // 20100715  M. Kelsey -- Use new G4CollisionOutput totals instead of loops,
00042 //              move temporary buffer to be data member
00043 // 20100719  M. Kelsey -- Change zero tolerance to 10 keV instead of 1 keV.
00044 // 20100909  M. Kelsey -- Add interface for both kinds of particle lists
00045 // 20101019  M. Kelsey -- CoVerity report: unitialized constructor
00046 // 20110328  M. Kelsey -- Add default ctor and explicit limit setting
00047 // 20110722  M. Kelsey -- For IntraNucleiCascader, take G4CollOut as argument
00048 // 20120525  M. Kelsey -- Follow process-level checking: allow _either_ rel.
00049 //              or abs. to pass (instead of requiring both)
00050 // 20121002  M. Kelsey -- Add strangeness check (useful for Omega- beam)
00051 
00052 #include "G4CascadeCheckBalance.hh"
00053 #include "globals.hh"
00054 #include "G4CascadParticle.hh"
00055 #include "G4InuclElementaryParticle.hh"
00056 #include "G4InuclNuclei.hh"
00057 #include "G4InuclParticle.hh"
00058 #include "G4CollisionOutput.hh"
00059 #include "G4LorentzVector.hh"
00060 #include <vector>
00061 
00062 
00063 // Constructor sets acceptance limits
00064 
00065 const G4double G4CascadeCheckBalance::tolerance = 1e-6; // How small is zero?
00066 
00067 G4CascadeCheckBalance::G4CascadeCheckBalance(const char* owner)
00068   : G4VCascadeCollider(owner), relativeLimit(G4CascadeCheckBalance::tolerance),
00069     absoluteLimit(G4CascadeCheckBalance::tolerance), initialBaryon(0),
00070     finalBaryon(0), initialCharge(0), finalCharge(0), initialStrange(0),
00071     finalStrange(0) {}
00072 
00073 G4CascadeCheckBalance::G4CascadeCheckBalance(G4double relative,
00074                                              G4double absolute,
00075                                              const char* owner)
00076   : G4VCascadeCollider(owner), relativeLimit(relative),
00077     absoluteLimit(absolute), initialBaryon(0), finalBaryon(0),
00078     initialCharge(0), finalCharge(0), initialStrange(0),
00079     finalStrange(0) {}
00080 
00081 
00082 // Pseudo-collision just computes input and output four-vectors
00083 
00084 void G4CascadeCheckBalance::collide(G4InuclParticle* bullet,
00085                                     G4InuclParticle* target,
00086                                     G4CollisionOutput& output) {
00087   if (verboseLevel)
00088     G4cout << " >>> G4CascadeCheckBalance(" << theName << ")::collide"
00089            << G4endl;
00090 
00091   initial *= 0.;        // Fast reset; some colliders only have one pointer
00092   if (bullet) initial += bullet->getMomentum();
00093   if (target) initial += target->getMomentum();
00094 
00095   // Baryon number, charge and strangeness must be computed "by hand"
00096   initialCharge = 0;
00097   if (bullet) initialCharge += G4int(bullet->getCharge());
00098   if (target) initialCharge += G4int(target->getCharge());
00099 
00100   G4InuclElementaryParticle* pbullet =
00101     dynamic_cast<G4InuclElementaryParticle*>(bullet);
00102   G4InuclElementaryParticle* ptarget =
00103     dynamic_cast<G4InuclElementaryParticle*>(target);
00104 
00105   G4InuclNuclei* nbullet = dynamic_cast<G4InuclNuclei*>(bullet);
00106   G4InuclNuclei* ntarget = dynamic_cast<G4InuclNuclei*>(target);
00107 
00108   initialBaryon =
00109     ((pbullet ? pbullet->baryon() : nbullet ? nbullet->getA() : 0) +
00110      (ptarget ? ptarget->baryon() : ntarget ? ntarget->getA() : 0) );
00111 
00112   // NOTE:  Currently we ignore possibility of hypernucleus target
00113   initialStrange = 0;
00114   if (pbullet) initialStrange += pbullet->getStrangeness();
00115   if (ptarget) initialStrange += ptarget->getStrangeness();
00116 
00117   // Final state totals are computed for us
00118   final = output.getTotalOutputMomentum();
00119   finalBaryon = output.getTotalBaryonNumber();
00120   finalCharge = output.getTotalCharge();
00121   finalStrange = output.getTotalStrangeness();
00122 
00123   // Report results
00124   if (verboseLevel) {
00125     G4cout << " initial px " << initial.px() << " py " << initial.py()
00126            << " pz " << initial.pz() << " E " << initial.e()
00127            << " baryon " << initialBaryon << " charge " << initialCharge
00128            << " strange " << initialStrange << G4endl
00129            << "   final px " << final.px() << " py " << final.py()
00130            << " pz " << final.pz() << " E " << final.e()
00131            << " baryon " << finalBaryon << " charge " << finalCharge
00132            << " strange " << finalStrange << G4endl;
00133   }
00134 }
00135 
00136 // Take list of output particles directly (e.g., from G4EPCollider internals)
00137 
00138 void G4CascadeCheckBalance::collide(G4InuclParticle* bullet, 
00139                                     G4InuclParticle* target,
00140     const std::vector<G4InuclElementaryParticle>& particles) {
00141   if (verboseLevel)
00142     G4cout << " >>> G4CascadeCheckBalance(" << theName << ")::collide(<vector>)"
00143            << G4endl;
00144 
00145   tempOutput.reset();                   // Buffer for processing
00146   tempOutput.addOutgoingParticles(particles);
00147   collide(bullet, target, tempOutput);
00148 }
00149 
00150 
00151 // Take list of nuclear fragments directly (e.g., from G4Fissioner internals)
00152 
00153 void G4CascadeCheckBalance::collide(G4InuclParticle* bullet, 
00154                                     G4InuclParticle* target,
00155     const std::vector<G4InuclNuclei>& fragments) {
00156   if (verboseLevel)
00157     G4cout << " >>> G4CascadeCheckBalance(" << theName << ")::collide(<vector>)"
00158            << G4endl;
00159 
00160   tempOutput.reset();                   // Buffer for processing
00161   tempOutput.addOutgoingNuclei(fragments);
00162   collide(bullet, target, tempOutput);
00163 }
00164 
00165 
00166 // Take list of "cparticles" (e.g., from G4NucleiModel internals)
00167 
00168 void G4CascadeCheckBalance::collide(G4InuclParticle* bullet,
00169                                     G4InuclParticle* target,
00170                     const std::vector<G4CascadParticle>& particles) {
00171   if (verboseLevel)
00172     G4cout << " >>> G4CascadeCheckBalance(" << theName
00173            << ")::collide(<cparticles>)" << G4endl;
00174 
00175   tempOutput.reset();                   // Buffer for processing
00176   tempOutput.addOutgoingParticles(particles);
00177   collide(bullet, target, tempOutput);
00178 }
00179 
00180 
00181 // Take lists of both G4InuclEP & G4CP (e.g., from G4IntraNucleiCascader)
00182 
00183 void G4CascadeCheckBalance::collide(G4InuclParticle* bullet,
00184                                    G4InuclParticle* target,
00185                                     G4CollisionOutput& output,
00186              const std::vector<G4CascadParticle>& cparticles) {
00187   if (verboseLevel)
00188     G4cout << " >>> G4CascadeCheckBalance(" << theName
00189            << ")::collide(<EP>,<CP>)" << G4endl;
00190 
00191   tempOutput.reset();                   // Buffer for processing
00192   tempOutput.add(output);
00193   tempOutput.addOutgoingParticles(cparticles);
00194   collide(bullet, target, tempOutput);
00195 }
00196 
00197 
00198 // Compare relative and absolute violations to limits, and report
00199 
00200 G4bool G4CascadeCheckBalance::energyOkay() const {
00201   G4bool relokay = (std::abs(relativeE()) < relativeLimit);
00202   G4bool absokay = (std::abs(deltaE()) < absoluteLimit);
00203 
00204   if (verboseLevel && !(relokay || absokay)) {
00205     G4cerr << theName << ": Energy conservation: relative " << relativeE()
00206            << (relokay ? " conserved" : " VIOLATED")
00207            << " absolute " << deltaE()
00208            << (absokay ? " conserved" : " VIOLATED") << G4endl;
00209   } else if (verboseLevel > 1) {
00210     G4cout << theName << ": Energy conservation: relative " << relativeE()
00211            << " conserved absolute " << deltaE() << " conserved" << G4endl;
00212   }
00213 
00214   return (relokay && absokay);
00215 }
00216 
00217 G4bool G4CascadeCheckBalance::ekinOkay() const {
00218   G4bool relokay = (std::abs(relativeKE()) < relativeLimit);
00219   G4bool absokay = (std::abs(deltaKE()) < absoluteLimit);
00220 
00221   if (verboseLevel && !(relokay || absokay)) {
00222     G4cerr << theName << ": Kinetic energy balance: relative "
00223            << relativeKE() << (relokay ? " conserved" : " VIOLATED")
00224            << " absolute " << deltaKE()
00225            << (absokay ? " conserved" : " VIOLATED") << G4endl;
00226   } else if (verboseLevel > 1) {
00227     G4cout << theName << ": Kinetic energy balance: relative "
00228            << relativeKE() << " conserved absolute " << deltaKE()
00229            << " conserved" << G4endl;
00230   }
00231 
00232   return (relokay && absokay);
00233 }
00234 
00235 G4bool G4CascadeCheckBalance::momentumOkay() const {
00236   G4bool relokay = (std::abs(relativeP()) < relativeLimit);
00237   G4bool absokay = (std::abs(deltaP()) < absoluteLimit);
00238 
00239   if (verboseLevel && !(relokay || absokay)) {
00240     G4cerr << theName << ": Momentum conservation: relative " << relativeP()
00241            << (relokay ? " conserved" : " VIOLATED")
00242            << " absolute " << deltaP()
00243            << (absokay ? " conserved" : " VIOLATED") << G4endl;
00244   } else if (verboseLevel > 1) {
00245     G4cout << theName << ": Momentum conservation: relative " << relativeP()
00246            << " conserved absolute " << deltaP() << " conserved" << G4endl;
00247   }
00248 
00249   return (relokay && absokay);
00250 }
00251 
00252 G4bool G4CascadeCheckBalance::baryonOkay() const {
00253   G4bool bokay = (deltaB() == 0);       // Must be perfect!
00254 
00255   if (verboseLevel && !bokay)
00256     G4cerr << theName << ": Baryon number VIOLATED " << deltaB() << G4endl;
00257 
00258   return bokay;
00259 }
00260 
00261 G4bool G4CascadeCheckBalance::chargeOkay() const {
00262   G4bool qokay = (deltaQ() == 0);       // Must be perfect!
00263 
00264   if (verboseLevel && !qokay)
00265     G4cerr << theName << ": Charge conservation VIOLATED " << deltaQ()
00266            << G4endl;
00267 
00268   return qokay;
00269 }
00270 
00271 G4bool G4CascadeCheckBalance::strangeOkay() const {
00272   G4bool sokay = (deltaS() == 0);       // Must be perfect!
00273 
00274   if (verboseLevel && !sokay)
00275     G4cerr << theName << ": Strangeness conservation VIOLATED " << deltaS()
00276            << G4endl;
00277 
00278   return sokay;
00279 }

Generated on Mon May 27 17:47:49 2013 for Geant4 by  doxygen 1.4.7