G4INCLPauliStandard.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 // INCL++ intra-nuclear cascade model
00027 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
00028 // Davide Mancusi, CEA
00029 // Alain Boudard, CEA
00030 // Sylvie Leray, CEA
00031 // Joseph Cugnon, University of Liege
00032 //
00033 #define INCLXX_IN_GEANT4_MODE 1
00034 
00035 #include "globals.hh"
00036 
00037 #include "G4INCLPauliStandard.hh"
00038 #include "G4INCLPauliBlocking.hh"
00039 #include "G4INCLGlobals.hh"
00040 #include "G4INCLRandom.hh"
00041 
00042 namespace G4INCL {
00043 
00044 
00045   PauliStandard::PauliStandard() :
00046     cellSize(std::pow(2.38*4.5*Math::pi,1./6.)*std::sqrt(PhysicalConstants::hc))
00047   {
00048     DEBUG("Initialising PauliStandard. cellSize=" << cellSize << std::endl);
00049   }
00050 
00051   PauliStandard::~PauliStandard() {}
00052 
00053   G4bool PauliStandard::isBlocked(ParticleList const pL, Nucleus const * const n) const {
00054     for(ParticleIter p = pL.begin(); p != pL.end(); ++p) {
00055       if( !(*p)->isNucleon() ) continue;
00056       if(getBlockingProbability(*p, n) > Random::shoot()) return true;
00057     }
00058     return false;
00059   }
00060 
00061   G4double PauliStandard::getBlockingProbability(Particle const * const particle, Nucleus const * const nucleus) const {
00062     const G4double r0 = nucleus->getDensity()->getNuclearRadius();
00063     const G4double pFermi = nucleus->getPotential()->getFermiMomentum(particle);
00064 
00065     const G4double pbl = cellSize * std::sqrt(pFermi/r0);
00066     const G4double rbl = pbl * r0/pFermi;
00067     const G4double maxVolR = rbl;
00068     const G4double maxVolP = pbl;
00069     G4double vol = std::pow(4.*Math::pi/3.0, 2)
00070       * std::pow(maxVolR*maxVolP/(Math::twoPi*PhysicalConstants::hc), 3);
00071 
00072     const G4double rdeq = nucleus->getUniverseRadius();
00073     const G4double rs = particle->getPosition().mag();
00074 
00075     if(rs - maxVolR > rdeq) {
00076       return 0.0;
00077     }
00078 
00079     if(rs + maxVolR > rdeq) {
00080       vol = vol * 0.5 * (rdeq - rs + maxVolR) / maxVolR;
00081     }
00082 
00083     // Get the list of particles that are currently inside the
00084     // nucleus.
00085     ParticleList particles = nucleus->getStore()->getParticles();
00086 
00087     G4int nl = 0;
00088     for(ParticleIter it = particles.begin(); it != particles.end(); ++it) {
00089       // Skip comparing with the same particle
00090       if( (*it)->getID() == particle->getID() ) continue;
00091 
00092       if((*it)->getType() == particle->getType()) {
00093         const ThreeVector dx2v = particle->getPosition() - (*it)->getPosition();
00094         const G4double dx2 = dx2v.mag2();
00095         if(dx2 > maxVolR * maxVolR) continue;
00096 
00097         const ThreeVector dp2v = particle->getMomentum() - (*it)->getMomentum();
00098         const G4double dp2 = dp2v.mag2();
00099         if(dp2 > maxVolP * maxVolP) continue;
00100 
00101         nl++;
00102       }
00103     }
00104     const G4double blockingProbability = ((G4double) nl) / vol / 2.0;
00105 
00106     if(blockingProbability > 1.0) return 1.0;
00107     else if(blockingProbability < 0.0) return 0.0;
00108     else return blockingProbability;
00109   }
00110 }

Generated on Mon May 27 17:48:36 2013 for Geant4 by  doxygen 1.4.7