G4CascadeRecoilMaker.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 // Collects generated cascade data (using Collider::collide() interface)
00029 // and computes the nuclear recoil kinematics needed to balance the event.
00030 //
00031 // 20100909  M. Kelsey -- Inspired by G4CascadeCheckBalance
00032 // 20100909  M. Kelsey -- Move G4IntraNucleiCascader::goodCase() here, add
00033 //              tolerance for "almost zero" excitation energy
00034 // 20100910  M. Kelsey -- Drop getRecoilFragment() in favor of user calling
00035 //              makeRecoilFragment() with returned non-const pointer.  Drop
00036 //              handling of excitons.
00037 // 20100921  M. Kelsey -- Return G4InuclNuclei using "makeRecoilNuclei()".
00038 //              Repurpose "makeRecoilFragment()" to return G4Fragment.
00039 // 20100924  M. Kelsey -- Remove unusable G4Fragment::SetExcitationEnergy().
00040 //              Add deltaM to compute mass difference, use excitationEnergy
00041 //              to force G4Fragment four-vector to match.
00042 // 20110214  M. Kelsey -- Follow G4InuclParticle::Model enumerator migration
00043 // 20110303  M. Kelsey -- Add diagnostic messages to goodNucleus().
00044 // 20110308  M. Kelsey -- Follow new G4Fragment interface for hole types
00045 // 20110722  M. Kelsey -- For IntraNucleiCascader, take G4CollOut as argument
00046 
00047 #include <vector>
00048 
00049 #include "G4CascadeRecoilMaker.hh"
00050 #include "globals.hh"
00051 #include "G4SystemOfUnits.hh"
00052 #include "G4CascadParticle.hh"
00053 #include "G4CascadeCheckBalance.hh"
00054 #include "G4CollisionOutput.hh"
00055 #include "G4Fragment.hh"
00056 #include "G4InuclElementaryParticle.hh"
00057 #include "G4InuclNuclei.hh"
00058 #include "G4InuclParticle.hh"
00059 #include "G4InuclSpecialFunctions.hh"
00060 #include "G4LorentzVector.hh"
00061 
00062 using namespace G4InuclSpecialFunctions;
00063 
00064 
00065 // Constructor and destructor
00066 
00067 G4CascadeRecoilMaker::G4CascadeRecoilMaker(G4double tolerance)
00068   : G4VCascadeCollider("G4CascadeRecoilMaker"),
00069     excTolerance(tolerance), inputEkin(0.),
00070     recoilA(0), recoilZ(0), excitationEnergy(0.) {
00071   balance = new G4CascadeCheckBalance(tolerance, tolerance, theName);
00072 }
00073   
00074 G4CascadeRecoilMaker::~G4CascadeRecoilMaker() {
00075   delete balance;
00076 }
00077 
00078 
00079 // Standard Collider interface (non-const output "buffer")
00080 
00081 void G4CascadeRecoilMaker::collide(G4InuclParticle* bullet,
00082                                     G4InuclParticle* target,
00083                                     G4CollisionOutput& output) {
00084   if (verboseLevel > 1)
00085     G4cout << " >>> G4CascadeRecoilMaker::collide" << G4endl;
00086 
00087   // Available energy needed for "goodNucleus()" test at end
00088   inputEkin = bullet ? bullet->getKineticEnergy() : 0.;
00089 
00090   balance->setVerboseLevel(verboseLevel);
00091   balance->collide(bullet, target, output);
00092   fillRecoil();
00093 }
00094 
00095 // This is for use with G4IntraNucleiCascader
00096 
00097 void G4CascadeRecoilMaker::collide(G4InuclParticle* bullet,
00098                                    G4InuclParticle* target,
00099                                    G4CollisionOutput& output,
00100              const std::vector<G4CascadParticle>& cparticles) {
00101   if (verboseLevel > 1)
00102     G4cout << " >>> G4CascadeRecoilMaker::collide(<EP>,<CP>)" << G4endl;
00103 
00104   // Available energy needed for "goodNucleus()" test at end
00105   inputEkin = bullet ? bullet->getKineticEnergy() : 0.;
00106 
00107   balance->setVerboseLevel(verboseLevel);
00108   balance->collide(bullet, target, output, cparticles);
00109   fillRecoil();
00110 }
00111 
00112 
00113 // Used current event configuration to construct recoil nucleus
00114 // NOTE:  CheckBalance uses "final-initial", we want "initial-final"
00115 
00116 void G4CascadeRecoilMaker::fillRecoil() {
00117   recoilZ = -(balance->deltaQ());               // Charge "non-conservation"
00118   recoilA = -(balance->deltaB());               // Baryon "non-conservation"
00119   recoilMomentum = -(balance->deltaLV());
00120 
00121   theExcitons.clear();          // Discard previous exciton configuraiton
00122 
00123   // Bertini uses MeV for excitation energy
00124   if (!goodFragment()) excitationEnergy = 0.;
00125   else excitationEnergy = deltaM() * GeV;
00126 
00127   // Allow for very small negative mass difference, and round to zero
00128   if (std::abs(excitationEnergy) < excTolerance) excitationEnergy = 0.;
00129 
00130   if (verboseLevel > 2) {
00131     G4cout << "  recoil px " << recoilMomentum.px()
00132            << " py " << recoilMomentum.py() << " pz " << recoilMomentum.pz()
00133            << " E " << recoilMomentum.e() << " baryon " << recoilA
00134            << " charge " << recoilZ
00135            << "\n  recoil mass " << recoilMomentum.m()
00136            << " 'excitation' energy " << excitationEnergy << G4endl;
00137   }
00138 }
00139 
00140 
00141 // Construct physical nucleus from recoil parameters, if reasonable
00142 
00143 G4InuclNuclei* 
00144 G4CascadeRecoilMaker::makeRecoilNuclei(G4InuclParticle::Model model) {
00145   if (verboseLevel > 1) 
00146     G4cout << " >>> G4CascadeRecoilMaker::makeRecoilNuclei" << G4endl;
00147 
00148   if (!goodRecoil()) {
00149     if (verboseLevel > 2 && !wholeEvent())
00150       G4cout << theName << ": event recoil is not a physical nucleus" << G4endl;
00151 
00152     return 0;           // Null pointer means no fragment
00153   }
00154 
00155   theRecoilNuclei.fill(recoilMomentum, recoilA, recoilZ,
00156                        excitationEnergy, model);
00157   theRecoilNuclei.setExitonConfiguration(theExcitons);
00158 
00159   return &theRecoilNuclei;
00160 }
00161 
00162 
00163 // Construct pre-compound nuclear fragment from recoil parameters
00164 
00165 G4Fragment* G4CascadeRecoilMaker::makeRecoilFragment() {
00166   if (verboseLevel > 1) 
00167     G4cout << " >>> G4CascadeRecoilMaker::makeRecoilFragment" << G4endl;
00168 
00169   if (!goodRecoil()) {
00170     if (verboseLevel > 2 && !wholeEvent())
00171       G4cout << theName << ": event recoil is not a physical nucleus" << G4endl;
00172 
00173     return 0;           // Null pointer means no fragment
00174   }
00175 
00176   theRecoilFragment.SetZandA_asInt(recoilZ, recoilA);   // Note convention!
00177 
00178   // User may have overridden excitation energy; force four-momentum to match
00179   G4double fragMass = 
00180     G4InuclNuclei::getNucleiMass(recoilA,recoilZ) + excitationEnergy/GeV;
00181 
00182   G4LorentzVector fragMom; fragMom.setVectM(recoilMomentum.vect(), fragMass);
00183   theRecoilFragment.SetMomentum(fragMom*GeV);   // Bertini uses GeV!
00184 
00185   // Note:  exciton configuration has to be set piece by piece
00186   //            (arguments are Ntotal,Nproton in both cases)
00187   G4int nholes = theExcitons.protonHoles+theExcitons.neutronHoles;
00188   theRecoilFragment.SetNumberOfHoles(nholes, theExcitons.protonHoles);
00189 
00190   G4int nexcit = (theExcitons.protonQuasiParticles
00191                   + theExcitons.neutronQuasiParticles);
00192   theRecoilFragment.SetNumberOfExcitedParticle(nexcit,
00193                                     theExcitons.protonQuasiParticles);
00194 
00195   return &theRecoilFragment;
00196 }
00197 
00198 
00199 // Compute raw mass difference from recoil parameters
00200 
00201 G4double G4CascadeRecoilMaker::deltaM() const {
00202   G4double nucMass = G4InuclNuclei::getNucleiMass(recoilA,recoilZ);
00203   return (recoilMomentum.m() - nucMass);
00204 }
00205 
00206 
00207 // Data quality checks
00208 
00209 G4bool G4CascadeRecoilMaker::goodFragment() const {
00210   return (recoilA>0 && recoilZ>=0 && recoilA >= recoilZ);
00211 }
00212 
00213 G4bool G4CascadeRecoilMaker::goodRecoil() const {
00214   return (goodFragment() && excitationEnergy > -excTolerance);
00215 }
00216 
00217 G4bool G4CascadeRecoilMaker::wholeEvent() const {
00218   if (verboseLevel > 2) {
00219     G4cout << " >>> G4CascadeRecoilMaker::wholeEvent:"
00220            << " A " << recoilA << " Z " << recoilZ
00221            << " P " << recoilMomentum.rho() << " E " << recoilMomentum.e()
00222            << "\n wholeEvent returns "
00223            << (recoilA==0 && recoilZ==0 && 
00224                recoilMomentum.rho() < excTolerance/GeV &&
00225                std::abs(recoilMomentum.e()) < excTolerance/GeV) << G4endl;
00226   }
00227 
00228   return (recoilA==0 && recoilZ==0 && 
00229           recoilMomentum.rho() < excTolerance/GeV &&
00230           std::abs(recoilMomentum.e()) < excTolerance/GeV);
00231 }
00232 
00233 // Determine whether desired nuclear fragment is constructable outcome
00234 
00235 G4bool G4CascadeRecoilMaker::goodNucleus() const {
00236   if (verboseLevel > 2) {
00237     G4cout << " >>> G4CascadeRecoilMaker::goodNucleus" << G4endl;
00238   }
00239 
00240   const G4double minExcitation = 0.1*keV;
00241   const G4double reasonableExcitation = 7.0;    // Multiple of binding energy
00242   const G4double fractionalExcitation = 0.2;    // Fraction of input to excite
00243 
00244   if (!goodRecoil()) {
00245     if (verboseLevel>2) {
00246       if (!goodFragment()) G4cerr << " goodNucleus: invalid A/Z" << G4endl;
00247       else if (excitationEnergy < -excTolerance) 
00248         G4cerr << " goodNucleus: negative excitation" << G4endl;
00249     }
00250     return false;                               // Not a sensible nucleus
00251   }
00252 
00253   if (excitationEnergy <= minExcitation) return true;   // Effectively zero
00254 
00255   // Maximum possible excitation energy determined by initial energy
00256   G4double dm = bindingEnergy(recoilA,recoilZ);
00257   G4double exc_max0z = fractionalExcitation * inputEkin*GeV;
00258   G4double exc_dm    = reasonableExcitation * dm;
00259   G4double exc_max = (exc_max0z > exc_dm) ? exc_max0z : exc_dm;
00260   
00261   if (verboseLevel > 3) {
00262     G4cout << " eexs " << excitationEnergy << " max " << exc_max
00263            << " dm " << dm << G4endl;
00264   }
00265 
00266   if (verboseLevel > 2 && excitationEnergy >= exc_max)
00267     G4cerr << " goodNucleus: too much excitation" << G4endl;
00268 
00269   return (excitationEnergy < exc_max);          // Below maximum possible
00270 }

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