G4BigBanger.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 // 20100114  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
00029 // 20100301  M. Kelsey -- In generateBangInSCM(), restore old G4CascMom calcs.
00030 //              for (N-1)th outgoing nucleon.
00031 // 20100319  M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff
00032 // 20100407  M. Kelsey -- Replace std::vector<> returns with data members.
00033 // 20100413  M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
00034 // 20100517  M. Kelsey -- Inherit from common base class, clean up code
00035 // 20100628  M. Kelsey -- Use old "bindingEnergy" fn as wrapper, add balance
00036 //              checking after bang.
00037 // 20100630  M. Kelsey -- Just do simple boost for target, instead of using
00038 //              G4LorentzConverter with dummy bullet.
00039 // 20100701  M. Kelsey -- Re-throw momentum list, not just angles!
00040 // 20100714  M. Kelsey -- Move conservation checking to base class
00041 // 20100726  M. Kelsey -- Move std::vector<> buffer to .hh file
00042 // 20100923  M. Kelsey -- Migrate to integer A and Z
00043 // 20110214  M. Kelsey -- Follow G4InuclParticle::Model enumerator migration
00044 // 20110806  M. Kelsey -- Pre-allocate buffers to reduce memory churn
00045 // 20110922  M. Kelsey -- Follow G4InuclParticle::print(ostream&) migration
00046 // 20120608  M. Kelsey -- Fix variable-name "shadowing" compiler warnings.
00047 
00048 #include <algorithm>
00049 
00050 #include "G4BigBanger.hh"
00051 #include "G4SystemOfUnits.hh"
00052 #include "G4CollisionOutput.hh"
00053 #include "G4InuclNuclei.hh"
00054 #include "G4InuclElementaryParticle.hh"
00055 #include "G4InuclSpecialFunctions.hh"
00056 #include "G4ParticleLargerEkin.hh"
00057 
00058 using namespace G4InuclSpecialFunctions;
00059 
00060 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
00061 
00062 G4BigBanger::G4BigBanger() : G4CascadeColliderBase("G4BigBanger") {}
00063 
00064 void
00065 G4BigBanger::collide(G4InuclParticle* /*bullet*/, G4InuclParticle* target,
00066                      G4CollisionOutput& output) {
00067 
00068   if (verboseLevel) G4cout << " >>> G4BigBanger::collide" << G4endl;
00069 
00070   // primitive explosion model A -> nucleons to prevent too exotic evaporation
00071 
00072   G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target);
00073   if (!nuclei_target) {
00074     G4cerr << " BigBanger -> try to bang not nuclei " << G4endl;
00075     return;
00076   }
00077 
00078   G4int A = nuclei_target->getA();
00079   G4int Z = nuclei_target->getZ();
00080 
00081   G4LorentzVector PEX = nuclei_target->getMomentum();
00082   G4double EEXS = nuclei_target->getExitationEnergy();
00083 
00084   G4ThreeVector toTheLabFrame = PEX.boostVector();      // From rest to lab
00085 
00086   // This "should" be difference between E-target and sum of m(nucleons)
00087   G4double etot = (EEXS - bindingEnergy(A,Z)) * MeV/GeV;  // To Bertini units
00088   if (etot < 0.0) etot = 0.0;
00089   
00090   if (verboseLevel > 2) {
00091     G4cout << " BigBanger: target\n" << *nuclei_target
00092            << "\n etot " << etot << G4endl;
00093   }
00094 
00095   if (verboseLevel > 3) {
00096     G4LorentzVector PEXrest = PEX;
00097     PEXrest.boost(-toTheLabFrame);
00098     G4cout << " target rest frame: px " << PEXrest.px() << " py "
00099            << PEXrest.py() << " pz " << PEXrest.pz() << " E " << PEXrest.e()
00100            << G4endl;
00101   }
00102 
00103   generateBangInSCM(etot, A, Z);
00104   
00105   if (verboseLevel > 2) {
00106     G4cout << " particles " << particles.size() << G4endl;
00107     for(G4int i = 0; i < G4int(particles.size()); i++) 
00108       G4cout << particles[i] << G4endl;
00109   }
00110 
00111   if (particles.empty()) {      // No bang!  Don't know why...
00112     G4cerr << " >>> G4BigBanger unable to process fragment "
00113            << nuclei_target->getDefinition()->GetParticleName() << G4endl;
00114 
00115     // FIXME:  This will violate baryon number, momentum, energy, etc.
00116     return;
00117   }
00118 
00119   // convert back to Lab
00120   G4LorentzVector totscm;
00121   G4LorentzVector totlab;
00122 
00123   if (verboseLevel > 2) G4cout << " BigBanger: boosting to lab" << G4endl;
00124 
00125   particleIterator ipart;
00126   for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
00127     G4LorentzVector mom = ipart->getMomentum();
00128     if (verboseLevel > 2) totscm += mom;
00129 
00130     mom.boost(toTheLabFrame);
00131     if (verboseLevel > 2) totlab += mom;
00132 
00133     ipart->setMomentum(mom); 
00134     if (verboseLevel > 2) G4cout << *ipart << G4endl;
00135   }
00136   
00137   std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
00138 
00139   validateOutput(0, target, particles);         // Checks <vector> directly
00140   
00141   if (verboseLevel > 2) {
00142     G4cout << " In SCM: total outgoing momentum " << G4endl 
00143            << " E " << totscm.e() << " px " << totscm.x()
00144            << " py " << totscm.y() << " pz " << totscm.z() << G4endl; 
00145     G4cout << " In Lab: mom cons " << G4endl 
00146            << " E " << PEX.e() - totlab.e()     // PEX now includes EEXS
00147            << " px " << PEX.x() - totlab.x()
00148            << " py " << PEX.y() - totlab.y() 
00149            << " pz " << PEX.z() - totlab.z() << G4endl; 
00150   }
00151 
00152   output.addOutgoingParticles(particles);
00153 }                    
00154 
00155 void G4BigBanger::generateBangInSCM(G4double etot, G4int a, G4int z) {
00156   if (verboseLevel > 3) {
00157     G4cout << " >>> G4BigBanger::generateBangInSCM" << G4endl;
00158   }
00159 
00160   const G4double ang_cut = 0.9999;
00161   const G4int itry_max = 1000;
00162   
00163   if (verboseLevel > 2) {
00164     G4cout << " a " << a << " z " << z << G4endl;
00165   }
00166 
00167   particles.clear();    // Reset output vector before filling
00168   
00169   if (a == 1) {         // Special -- bare nucleon doesn't really "explode"
00170     G4int knd = (z>0) ? 1 : 2;
00171     particles.push_back(G4InuclElementaryParticle(knd)); // zero momentum
00172     return;
00173   }
00174      
00175   // NOTE:  If distribution fails, need to regenerate magnitudes and angles!
00176   //*** generateMomentumModules(etot, a, z);
00177 
00178   scm_momentums.reserve(a);
00179   G4LorentzVector tot_mom;
00180 
00181   G4bool bad = true;
00182   G4int itry = 0;
00183   while(bad && itry < itry_max) {
00184     itry++;
00185     scm_momentums.clear();
00186 
00187     generateMomentumModules(etot, a, z);
00188     if (a == 2) {
00189       // This is only a three-vector, not a four-vector
00190       G4LorentzVector mom = generateWithRandomAngles(momModules[0]);
00191       scm_momentums.push_back(mom);
00192       scm_momentums.push_back(-mom);    // Only safe since three-vector!
00193       bad = false;
00194     } else {
00195       tot_mom *= 0.;            // Easy way to reset accumulator
00196 
00197       for(G4int i = 0; i < a-2; i++) {          // All but last two are thrown
00198       // This is only a three-vector, not a four-vector
00199         G4LorentzVector mom = generateWithRandomAngles(momModules[i]);
00200         scm_momentums.push_back(mom);
00201         tot_mom += mom;          
00202       };
00203 
00204       //                handle last two
00205       G4double tot_mod = tot_mom.rho(); 
00206       G4double ct = -0.5*(tot_mod*tot_mod + momModules[a-2]*momModules[a-2]
00207                           - momModules[a-1]*momModules[a-1]) / tot_mod
00208         / momModules[a-2];
00209 
00210       if (verboseLevel > 2) G4cout << " ct last " << ct << G4endl;
00211   
00212       if(std::fabs(ct) < ang_cut) {
00213         // This is only a three-vector, not a four-vector
00214         G4LorentzVector mom2 = generateWithFixedTheta(ct, momModules[a - 2]);
00215 
00216         // rotate to the normal system
00217         G4LorentzVector apr = tot_mom/tot_mod;
00218         G4double a_tr = std::sqrt(apr.x()*apr.x() + apr.y()*apr.y());
00219         G4LorentzVector mom;
00220         mom.setX(mom2.z()*apr.x() + ( mom2.x()*apr.y() + mom2.y()*apr.z()*apr.x())/a_tr);
00221         mom.setY(mom2.z()*apr.y() + (-mom2.x()*apr.x() + mom2.y()*apr.z()*apr.y())/a_tr);
00222         mom.setZ(mom2.z()*apr.z() - mom2.y()*a_tr);
00223 
00224         scm_momentums.push_back(mom);
00225 
00226         // and the last one (again, not actually a four-vector!)
00227         G4LorentzVector mom1 = -mom - tot_mom;
00228 
00229         scm_momentums.push_back(mom1);  
00230         bad = false;
00231       } // if (abs(ct) < ang_cut)
00232     }   // (a > 2)
00233   }     // while (bad && itry<itry_max)
00234 
00235   if (!bad) {
00236     particles.resize(a);        // Use assignment to avoid temporaries
00237     for(G4int i = 0; i < a; i++) {
00238       G4int knd = i < z ? 1 : 2;
00239       particles[i].fill(scm_momentums[i], knd, G4InuclParticle::BigBanger);
00240     };
00241   };
00242 
00243   if (verboseLevel > 2) {
00244     if (itry == itry_max) G4cout << " BigBanger -> can not generate bang " << G4endl;
00245   }
00246 
00247   return;
00248 }
00249            
00250 void G4BigBanger::generateMomentumModules(G4double etot, G4int a, G4int z) {
00251   if (verboseLevel > 3) {
00252     G4cout << " >>> G4BigBanger::generateMomentumModules" << G4endl;
00253   }
00254 
00255   // Proton and neutron masses
00256   const G4double mp = G4InuclElementaryParticle::getParticleMass(1);
00257   const G4double mn = G4InuclElementaryParticle::getParticleMass(2);
00258 
00259   momModules.clear();           // Reset buffer for filling
00260 
00261   G4double xtot = 0.0;
00262 
00263   if (a > 2) {                  // For "large" nuclei, energy is distributed
00264     G4double promax = maxProbability(a);
00265 
00266     momModules.resize(a, 0.);   // Pre-allocate to avoid memory churn
00267     for(G4int i = 0; i < a; i++) { 
00268       momModules[i] = generateX(a, promax);
00269       xtot += momModules[i];
00270       
00271       if (verboseLevel > 2) {
00272         G4cout << " i " << i << " x " << momModules[i] << G4endl;
00273       }
00274     }
00275   } else {                      // Two-body case is special, must be 50%
00276     xtot = 1.;
00277     momModules.push_back(0.5);
00278     momModules.push_back(0.5);
00279   }
00280 
00281   for(G4int i = 0; i < a; i++) {
00282     G4double mass = i < z ? mp : mn;
00283 
00284     momModules[i] *= etot/xtot;
00285     momModules[i] = std::sqrt(momModules[i] * (momModules[i] + 2.0 * mass));
00286 
00287     if (verboseLevel > 2) {
00288       G4cout << " i " << i << " pmod " << momModules[i] << G4endl;
00289     }
00290   };
00291 
00292   return;  
00293 }
00294 
00295 G4double G4BigBanger::xProbability(G4double x, G4int a) const {
00296   if (verboseLevel > 3) G4cout << " >>> G4BigBanger::xProbability" << G4endl;
00297 
00298   G4double ekpr = 0.0;
00299 
00300   if(x < 1.0 || x > 0.0) {
00301     ekpr = x * x;
00302 
00303     if (a%2 == 0) { // even A
00304       ekpr *= std::sqrt(1.0 - x) * std::pow((1.0 - x), (3*a-6)/2); 
00305     }
00306     else {
00307       ekpr *= std::pow((1.0 - x), (3*a-5)/2);
00308     };
00309   }; 
00310   
00311   return ekpr;
00312 }
00313 
00314 G4double G4BigBanger::maxProbability(G4int a) const {
00315   if (verboseLevel > 3) {
00316     G4cout << " >>> G4BigBanger::maxProbability" << G4endl;
00317   }
00318 
00319   return xProbability(2./3./(a-1.0), a);
00320 }
00321 
00322 G4double G4BigBanger::generateX(G4int a, G4double promax) const {
00323   if (verboseLevel > 3) G4cout << " >>> G4BigBanger::generateX" << G4endl;
00324 
00325   const G4int itry_max = 1000;
00326   G4int itry = 0;
00327   G4double x;
00328   
00329   while(itry < itry_max) {
00330     itry++;
00331     x = inuclRndm();
00332 
00333     if(xProbability(x, a) >= promax * inuclRndm()) return x;
00334   };
00335   if (verboseLevel > 2) {
00336     G4cout << " BigBanger -> can not generate x " << G4endl;
00337   }
00338 
00339   return maxProbability(a);
00340 }

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