G4BigBanger Class Reference

#include <G4BigBanger.hh>

Inheritance diagram for G4BigBanger:

Collaboration diagram for G4BigBanger:

Public Member Functions
	G4BigBanger ()
virtual	~G4BigBanger ()
virtual void	deExcite (const G4Fragment &target, G4CollisionOutput &output)
Public Member Functions inherited from G4CascadeDeexciteBase
	G4CascadeDeexciteBase (const char *name)
virtual	~G4CascadeDeexciteBase ()
virtual void	setVerboseLevel (G4int verbose=0)
Public Member Functions inherited from G4VCascadeDeexcitation
	G4VCascadeDeexcitation (const G4String &name)
virtual	~G4VCascadeDeexcitation ()
virtual void	collide (G4InuclParticle *bullet, G4InuclParticle *target, G4CollisionOutput &globalOutput)
Public Member Functions inherited from G4VCascadeCollider
	G4VCascadeCollider (const G4String &name, G4int verbose=0)
virtual	~G4VCascadeCollider ()

Private Member Functions
void	generateBangInSCM (G4double etot, G4int a, G4int z)
void	generateMomentumModules (G4double etot, G4int a, G4int z)
G4double	xProbability (G4double x, G4int a) const
G4double	maxProbability (G4int a) const
G4double	generateX (G4int ia, G4double promax) const
	G4BigBanger (const G4BigBanger &)
G4BigBanger &	operator= (const G4BigBanger &)

Private Attributes
std::vector< G4InuclElementaryParticle >	particles
std::vector< G4double >	momModules
std::vector< G4LorentzVector >	scm_momentums

Additional Inherited Members
Protected Member Functions inherited from G4CascadeDeexciteBase
virtual G4bool	explosion (const G4Fragment &target) const
virtual G4bool	explosion (G4int A, G4int Z, G4double excitation) const
virtual G4bool	validateOutput (const G4Fragment &target, G4CollisionOutput &output)
virtual G4bool	validateOutput (const G4Fragment &target, const std::vector< G4InuclElementaryParticle > &particles)
virtual G4bool	validateOutput (const G4Fragment &target, const std::vector< G4InuclNuclei > &fragments)
void	getTargetData (const G4Fragment &target)
const G4Fragment &	makeFragment (G4LorentzVector mom, G4int A, G4int Z, G4double EX=0.)
const G4Fragment &	makeFragment (G4int A, G4int Z, G4double EX=0.)
Protected Member Functions inherited from G4VCascadeCollider
virtual void	setName (const G4String &name)
Protected Attributes inherited from G4CascadeDeexciteBase
G4CascadeCheckBalance *	balance
G4int	A
G4int	Z
G4LorentzVector	PEX
G4double	EEXS
G4Fragment	aFragment
Protected Attributes inherited from G4VCascadeCollider
G4String	theName
G4int	verboseLevel

Detailed Description

Definition at line 48 of file G4BigBanger.hh.

Constructor & Destructor Documentation

G4BigBanger() [1/2]

G4BigBanger::G4BigBanger

(

)

Definition at line 67 of file G4BigBanger.cc.

: G4CascadeDeexciteBase("G4BigBanger") {}

~G4BigBanger()

virtual G4BigBanger::~G4BigBanger ( )

inlinevirtual

Definition at line 51 of file G4BigBanger.hh.

{};

Here is the call graph for this function:

G4BigBanger() [2/2]

G4BigBanger::G4BigBanger ( const G4BigBanger &  )

private

Member Function Documentation

deExcite()

void G4BigBanger::deExcite ( const G4Fragment &  target, G4CollisionOutput &  output  )

virtual

Implements G4VCascadeDeexcitation.

Definition at line 69 of file G4BigBanger.cc.

                                                {
  if (verboseLevel) G4cout << " >>> G4BigBanger::deExcite" << G4endl;

  getTargetData(target);
  G4ThreeVector toTheLabFrame = PEX.boostVector();  // From rest to lab

  // This "should" be difference between E-target and sum of m(nucleons)
  G4double etot = (EEXS - bindingEnergy(A,Z)) * MeV/GeV;  // To Bertini units
  if (etot < 0.0) etot = 0.0;
  
  if (verboseLevel > 2) {
    G4cout << " BigBanger: target\n" << target
       << "\n etot " << etot << G4endl;
  }

  if (verboseLevel > 3) {
    G4LorentzVector PEXrest = PEX;
    PEXrest.boost(-toTheLabFrame);
    G4cout << " target rest frame: px " << PEXrest.px() << " py "
       << PEXrest.py() << " pz " << PEXrest.pz() << " E " << PEXrest.e()
       << G4endl;
  }

  generateBangInSCM(etot, A, Z);
  
  if (verboseLevel > 2) {
    G4cout << " particles " << particles.size() << G4endl;
    for(G4int i = 0; i < G4int(particles.size()); i++) 
      G4cout << particles[i] << G4endl;
  }

  if (particles.empty()) {  // No bang!  Don't know why...
    G4cerr << " >>> G4BigBanger unable to process fragment "
       << target << G4endl;

    // FIXME:  This will violate baryon number, momentum, energy, etc.
    return;
  }

  // convert back to Lab
  G4LorentzVector totscm;
  G4LorentzVector totlab;

  if (verboseLevel > 2) G4cout << " BigBanger: boosting to lab" << G4endl;

  particleIterator ipart;
  for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
    G4LorentzVector mom = ipart->getMomentum();
    if (verboseLevel > 2) totscm += mom;

    mom.boost(toTheLabFrame);
    if (verboseLevel > 2) totlab += mom;

    ipart->setMomentum(mom); 
    if (verboseLevel > 2) G4cout << *ipart << G4endl;
  }
  
  std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());

  validateOutput(target, particles);        // Checks <vector> directly
  
  if (verboseLevel > 2) {
    G4cout << " In SCM: total outgoing momentum " << G4endl 
       << " E " << totscm.e() << " px " << totscm.x()
       << " py " << totscm.y() << " pz " << totscm.z() << G4endl; 
    G4cout << " In Lab: mom cons " << G4endl 
       << " E " << PEX.e() - totlab.e() // PEX now includes EEXS
       << " px " << PEX.x() - totlab.x()
       << " py " << PEX.y() - totlab.y() 
       << " pz " << PEX.z() - totlab.z() << G4endl; 
  }

  globalOutput.addOutgoingParticles(particles);
}            

Here is the call graph for this function:

Here is the caller graph for this function:

generateBangInSCM()

void G4BigBanger::generateBangInSCM ( G4double  etot, G4int  a, G4int  z  )

private

Definition at line 145 of file G4BigBanger.cc.

                                                                   {
  if (verboseLevel > 3) {
    G4cout << " >>> G4BigBanger::generateBangInSCM" << G4endl;
  }

  const G4double ang_cut = 0.9999;
  const G4int itry_max = 1000;
  
  if (verboseLevel > 2) {
    G4cout << " a " << a << " z " << z << G4endl;
  }

  particles.clear();    // Reset output vector before filling
  
  if (a == 1) {     // Special -- bare nucleon doesn't really "explode"
    G4int knd = (z>0) ? 1 : 2;
    particles.push_back(G4InuclElementaryParticle(knd)); // zero momentum
    return;
  }
     
  // NOTE:  If distribution fails, need to regenerate magnitudes and angles!
  //*** generateMomentumModules(etot, a, z);

  scm_momentums.reserve(a);
  G4LorentzVector tot_mom;

  G4bool bad = true;
  G4int itry = 0;
  while(bad && itry < itry_max) {   /* Loop checking 08.06.2015 MHK */
    itry++;
    scm_momentums.clear();

    generateMomentumModules(etot, a, z);
    if (a == 2) {
      // This is only a three-vector, not a four-vector
      G4LorentzVector mom = generateWithRandomAngles(momModules[0]);
      scm_momentums.push_back(mom);
      scm_momentums.push_back(-mom);    // Only safe since three-vector!
      bad = false;
    } else {
      tot_mom *= 0.;        // Easy way to reset accumulator

      for(G4int i = 0; i < a-2; i++) {      // All but last two are thrown
      // This is only a three-vector, not a four-vector
    G4LorentzVector mom = generateWithRandomAngles(momModules[i]);
    scm_momentums.push_back(mom);
    tot_mom += mom;      
      };

      //                handle last two
      G4double tot_mod = tot_mom.rho(); 
      G4double ct = -0.5*(tot_mod*tot_mod + momModules[a-2]*momModules[a-2]
              - momModules[a-1]*momModules[a-1]) / tot_mod
    / momModules[a-2];

      if (verboseLevel > 2) G4cout << " ct last " << ct << G4endl;
  
      if(std::fabs(ct) < ang_cut) {
    // This is only a three-vector, not a four-vector
    G4LorentzVector mom2 = generateWithFixedTheta(ct, momModules[a - 2]);

    // rotate to the normal system
    G4LorentzVector apr = tot_mom/tot_mod;
    G4double a_tr = std::sqrt(apr.x()*apr.x() + apr.y()*apr.y());
    G4LorentzVector mom;
    mom.setX(mom2.z()*apr.x() + ( mom2.x()*apr.y() + mom2.y()*apr.z()*apr.x())/a_tr);
    mom.setY(mom2.z()*apr.y() + (-mom2.x()*apr.x() + mom2.y()*apr.z()*apr.y())/a_tr);
    mom.setZ(mom2.z()*apr.z() - mom2.y()*a_tr);

    scm_momentums.push_back(mom);

    // and the last one (again, not actually a four-vector!)
    G4LorentzVector mom1 = -mom - tot_mom;

    scm_momentums.push_back(mom1);  
    bad = false;
      } // if (abs(ct) < ang_cut)
    }   // (a > 2)
  } // while (bad && itry<itry_max)

  if (!bad) {
    particles.resize(a);    // Use assignment to avoid temporaries
    for(G4int i = 0; i < a; i++) {
      G4int knd = i < z ? 1 : 2;
      particles[i].fill(scm_momentums[i], knd, G4InuclParticle::BigBanger);
    };
  };

  if (verboseLevel > 2) {
    if (itry == itry_max) G4cout << " BigBanger -> can not generate bang " << G4endl;
  }

  return;
}

Here is the call graph for this function:

Here is the caller graph for this function:

generateMomentumModules()

void G4BigBanger::generateMomentumModules ( G4double  etot, G4int  a, G4int  z  )

private

Definition at line 240 of file G4BigBanger.cc.

                                                                         {
  if (verboseLevel > 3) {
    G4cout << " >>> G4BigBanger::generateMomentumModules" << G4endl;
  }

  // Proton and neutron masses
  const G4double mp = G4InuclElementaryParticle::getParticleMass(1);
  const G4double mn = G4InuclElementaryParticle::getParticleMass(2);

  momModules.clear();       // Reset buffer for filling

  G4double xtot = 0.0;

  if (a > 2) {          // For "large" nuclei, energy is distributed
    G4double promax = maxProbability(a);

    momModules.resize(a, 0.);   // Pre-allocate to avoid memory churn
    for(G4int i = 0; i < a; i++) { 
      momModules[i] = generateX(a, promax);
      xtot += momModules[i];
      
      if (verboseLevel > 2) {
    G4cout << " i " << i << " x " << momModules[i] << G4endl;
      }
    }
  } else {          // Two-body case is special, must be 50%
    xtot = 1.;
    momModules.push_back(0.5);
    momModules.push_back(0.5);
  }

  for(G4int i = 0; i < a; i++) {
    G4double mass = i < z ? mp : mn;

    momModules[i] *= etot/xtot;
    momModules[i] = std::sqrt(momModules[i] * (momModules[i] + 2.0 * mass));

    if (verboseLevel > 2) {
      G4cout << " i " << i << " pmod " << momModules[i] << G4endl;
    }
  };

  return;  
}

Here is the call graph for this function:

Here is the caller graph for this function:

generateX()

G4double G4BigBanger::generateX ( G4int  ia, G4double  promax  ) const

private

Definition at line 313 of file G4BigBanger.cc.

                                                              {
  if (verboseLevel > 3) G4cout << " >>> G4BigBanger::generateX" << G4endl;

  const G4int itry_max = 1000;
  G4int itry = 0;
  G4double x;
  
  while(itry < itry_max) {  /* Loop checking 08.06.2015 MHK */
    itry++;
    x = inuclRndm();

    if(xProbability(x, a) >= promax * inuclRndm()) return x;
  };
  if (verboseLevel > 2) {
    G4cout << " BigBanger -> can not generate x " << G4endl;
  }

  return maxProbability(a);
}

Here is the call graph for this function:

Here is the caller graph for this function:

maxProbability()

G4double G4BigBanger::maxProbability ( G4int  a ) const

private

Definition at line 305 of file G4BigBanger.cc.

                                                  {
  if (verboseLevel > 3) {
    G4cout << " >>> G4BigBanger::maxProbability" << G4endl;
  }

  return xProbability(2./3./(a-1.0), a);
}

Here is the call graph for this function:

Here is the caller graph for this function:

operator=()

G4BigBanger& G4BigBanger::operator= ( const G4BigBanger &  )

private

xProbability()

G4double G4BigBanger::xProbability ( G4double  x, G4int  a  ) const

private

Definition at line 285 of file G4BigBanger.cc.

                                                            {
  if (verboseLevel > 3) G4cout << " >>> G4BigBanger::xProbability" << G4endl;

  G4Pow* theG4Pow = G4Pow::GetInstance();   // For convenience

  G4double ekpr = 0.0;
  if(x < 1.0 || x > 0.0) {
    ekpr = x * x;

    if (a%2 == 0) { // even A
      ekpr *= std::sqrt(1.0 - x) * theG4Pow->powN((1.0 - x), (3*a-6)/2); 
    }
    else {
      ekpr *= theG4Pow->powN((1.0 - x), (3*a-5)/2);
    };
  }; 
  
  return ekpr;
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

momModules

std::vector<G4double> G4BigBanger::momModules

private

Definition at line 64 of file G4BigBanger.hh.

particles

std::vector<G4InuclElementaryParticle> G4BigBanger::particles

private

Definition at line 63 of file G4BigBanger.hh.

scm_momentums

std::vector<G4LorentzVector> G4BigBanger::scm_momentums

private

Definition at line 65 of file G4BigBanger.hh.

---

The documentation for this class was generated from the following files:

• G4BigBanger.hh
• G4BigBanger.cc