G4Fancy3DNucleus Class Reference

#include <G4Fancy3DNucleus.hh>

Inheritance diagram for G4Fancy3DNucleus:

Collaboration diagram for G4Fancy3DNucleus:

Public Member Functions
	G4Fancy3DNucleus ()
	~G4Fancy3DNucleus ()
void	Init (G4int theA, G4int theZ)
G4bool	StartLoop ()
G4Nucleon *	GetNextNucleon ()
const std::vector< G4Nucleon > &	GetNucleons ()
G4int	GetMassNumber ()
G4double	GetMass ()
G4int	GetCharge ()
G4double	GetNuclearRadius ()
G4double	GetNuclearRadius (const G4double maxRelativeDensity)
G4double	GetOuterRadius ()
G4double	AddExcitationEnergy (G4double)
G4double	GetExcitationEnergy ()
G4double	CoulombBarrier ()
void	DoLorentzBoost (const G4LorentzVector &theBoost)
void	DoLorentzBoost (const G4ThreeVector &theBeta)
void	DoLorentzContraction (const G4LorentzVector &theBoost)
void	DoLorentzContraction (const G4ThreeVector &theBeta)
void	CenterNucleons ()
void	DoTranslation (const G4ThreeVector &theShift)
const G4VNuclearDensity *	GetNuclearDensity () const
void	SortNucleonsIncZ ()
void	SortNucleonsDecZ ()
Public Member Functions inherited from G4V3DNucleus
	G4V3DNucleus ()
virtual	~G4V3DNucleus ()
std::pair< G4double, G4double >	ChooseImpactXandY (G4double maxImpact)
std::pair< G4double, G4double >	RefetchImpactXandY ()

Private Member Functions
	G4Fancy3DNucleus (const G4Fancy3DNucleus &right)
const G4Fancy3DNucleus &	operator= (const G4Fancy3DNucleus &right)
int	operator== (const G4Fancy3DNucleus &right) const
int	operator!= (const G4Fancy3DNucleus &right) const
void	ChooseNucleons ()
void	ChoosePositions ()
void	ChooseFermiMomenta ()
G4double	BindingEnergy ()
G4bool	ReduceSum ()

Private Attributes
G4int	myA
G4int	myZ
std::vector< G4Nucleon >	theNucleons
G4int	currentNucleon
G4VNuclearDensity *	theDensity
G4FermiMomentum	theFermi
const G4double	nucleondistance
G4double	excitationEnergy
std::vector< G4ThreeVector >	places
std::vector< G4ThreeVector >	momentum
std::vector< G4double >	fermiM
std::vector< G4Fancy3DNucleusHelper >	testSums

Detailed Description

Definition at line 54 of file G4Fancy3DNucleus.hh.

Constructor & Destructor Documentation

G4Fancy3DNucleus() [1/2]

G4Fancy3DNucleus::G4Fancy3DNucleus

(

)

Definition at line 53 of file G4Fancy3DNucleus.cc.

  : myA(0), myZ(0), theNucleons(250), currentNucleon(-1), theDensity(0), 
    nucleondistance(0.8*fermi),excitationEnergy(0.),
    places(250), momentum(250), fermiM(250), testSums(250)
{
//G4cout <<"G4Fancy3DNucleus::G4Fancy3DNucleus()"<<G4endl;
}

~G4Fancy3DNucleus()

G4Fancy3DNucleus::~G4Fancy3DNucleus

(

)

Definition at line 61 of file G4Fancy3DNucleus.cc.

{
  if(theDensity) delete theDensity;
}

Here is the call graph for this function:

G4Fancy3DNucleus() [2/2]

G4Fancy3DNucleus::G4Fancy3DNucleus ( const G4Fancy3DNucleus &  right )

private

Member Function Documentation

AddExcitationEnergy()

G4double G4Fancy3DNucleus::AddExcitationEnergy ( G4double  anE )

inline

Definition at line 130 of file G4Fancy3DNucleus.hh.

{
   excitationEnergy +=anE;
   return excitationEnergy;
}

BindingEnergy()

G4double G4Fancy3DNucleus::BindingEnergy ( )

private

Definition at line 161 of file G4Fancy3DNucleus.cc.

{
  return G4NucleiProperties::GetBindingEnergy(myA,myZ);
}

Here is the call graph for this function:

Here is the caller graph for this function:

CenterNucleons()

void G4Fancy3DNucleus::CenterNucleons

(

)

Definition at line 238 of file G4Fancy3DNucleus.cc.

{
    G4ThreeVector   center;
    
    for (G4int i=0; i<myA; i++ )
    {
       center+=theNucleons[i].GetPosition();
    }   
    center /= -myA;
    DoTranslation(center);
}

Here is the call graph for this function:

ChooseFermiMomenta()

void G4Fancy3DNucleus::ChooseFermiMomenta ( )

private

Definition at line 357 of file G4Fancy3DNucleus.cc.

{
    G4int i;
    G4double density;

    // Pre-allocate buffers for filling by index
    momentum.resize(myA, G4ThreeVector(0.,0.,0.));
    fermiM.resize(myA, 0.*GeV);

    for (G4int ntry=0; ntry<1 ; ntry ++ )
    {
    for (i=0; i < myA; i++ )    // momenta for all, including last, in case we swap nucleons
    {
       density = theDensity->GetDensity(theNucleons[i].GetPosition());
       fermiM[i] = theFermi.GetFermiMomentum(density);
       G4ThreeVector mom=theFermi.GetMomentum(density);
       if (theNucleons[i].GetDefinition() == G4Proton::Proton())
       {
          G4double eMax = std::sqrt(sqr(fermiM[i]) +sqr(theNucleons[i].GetDefinition()->GetPDGMass()) )
                          - CoulombBarrier();
          if ( eMax > theNucleons[i].GetDefinition()->GetPDGMass() )
          {
              G4double pmax2= sqr(eMax) - sqr(theNucleons[i].GetDefinition()->GetPDGMass());
          fermiM[i] = std::sqrt(pmax2);
          while ( mom.mag2() > pmax2 )  /* Loop checking, 30-Oct-2015, G.Folger */
          {
              mom=theFermi.GetMomentum(density, fermiM[i]);
          }
          }  else
          {
              G4cerr << "G4Fancy3DNucleus: difficulty finding proton momentum" << G4endl;
          mom=G4ThreeVector(0,0,0);
          }

       }
       momentum[i]= mom;
    }

    if ( ReduceSum() ) break;
//       G4cout <<" G4FancyNucleus: iterating to find momenta: "<< ntry<< G4endl;
    }

//     G4ThreeVector sum;
//     for (G4int index=0; index<myA;sum+=momentum[index++])
//     ;
//     G4cout << "final sum / mag() " << sum << " / " << sum.mag() << G4endl;

    G4double energy;
    for ( i=0; i< myA ; i++ )
    {
       energy = theNucleons[i].GetParticleType()->GetPDGMass()
            - BindingEnergy()/myA;
       G4LorentzVector tempV(momentum[i],energy);
       theNucleons[i].SetMomentum(tempV);
       // GF 11-05-2011: set BindingEnergy to be T of Nucleon with p , ~ p**2/2m
       //theNucleons[i].SetBindingEnergy(
       //     0.5*sqr(fermiM[i])/theNucleons[i].GetParticleType()->GetPDGMass());
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

ChooseNucleons()

void G4Fancy3DNucleus::ChooseNucleons ( )

private

Definition at line 267 of file G4Fancy3DNucleus.cc.

{
    G4int protons=0,nucleons=0;
    
    while (nucleons < myA )  /* Loop checking, 30-Oct-2015, G.Folger */
    {
      if ( protons < myZ && G4UniformRand() < (G4double)(myZ-protons)/(G4double)(myA-nucleons) )
      {
         protons++;
         theNucleons[nucleons++].SetParticleType(G4Proton::Proton());
      }
      else if ( (nucleons-protons) < (myA-myZ) )
      {
           theNucleons[nucleons++].SetParticleType(G4Neutron::Neutron());
      }
      else G4cout << "G4Fancy3DNucleus::ChooseNucleons not efficient" << G4endl;
    }
    return;
}

Here is the call graph for this function:

Here is the caller graph for this function:

ChoosePositions()

void G4Fancy3DNucleus::ChoosePositions ( )

private

Definition at line 287 of file G4Fancy3DNucleus.cc.

{
    G4int i=0;
    G4ThreeVector   aPos, delta;
    G4bool      freeplace;
    const G4double nd2=sqr(nucleondistance);
    G4double maxR=GetNuclearRadius(0.001);   //  there are no nucleons at a
                                            //  relative Density of 0.01
    G4int jr=0;
    G4int jx,jy;
    G4double arand[600];
    G4double *prand=arand;

    places.clear();             // Reset data buffer
  G4int interationsLeft=1000*myA;
  while ( (i < myA) && (--interationsLeft>0))  /* Loop checking, 30-Oct-2015, G.Folger */
    {
       do
       {    
        if ( jr < 3 ) 
        {
            jr=std::min(600,9*(myA - i));
            G4RandFlat::shootArray(jr,prand);
            //CLHEP::RandFlat::shootArray(jr, prand );
        }
        jx=--jr;
        jy=--jr;
        aPos.set((2*arand[jx]-1.), (2*arand[jy]-1.), (2*arand[--jr]-1.));
       } while (aPos.mag2() > 1. );  /* Loop checking, 30-Oct-2015, G.Folger */
       aPos *=maxR;
       G4double density=theDensity->GetRelativeDensity(aPos);
       if (G4UniformRand() < density)
       {
          freeplace= true;
          std::vector<G4ThreeVector>::iterator iplace;
          for( iplace=places.begin(); iplace!=places.end() && freeplace;++iplace)
          {
            delta = *iplace - aPos;
        freeplace= delta.mag2() > nd2;
          }
          
          if ( freeplace )
          {
          G4double pFermi=theFermi.GetFermiMomentum(theDensity->GetDensity(aPos));
            // protons must at least have binding energy of CoulombBarrier, so
            //  assuming the Fermi energy corresponds to a potential, we must place these such
            //  that the Fermi Energy > CoulombBarrier
          if (theNucleons[i].GetDefinition() == G4Proton::Proton())
          {
            G4double nucMass = theNucleons[i].GetDefinition()->GetPDGMass();
            G4double eFermi= std::sqrt( sqr(pFermi) + sqr(nucMass) )
                              - nucMass;
            if (eFermi <= CoulombBarrier() ) freeplace=false;
          }
          }
          if ( freeplace )
          {
          theNucleons[i].SetPosition(aPos);
          places.push_back(aPos);
          ++i;
          }
       }
    }
    if (interationsLeft<=0) {
       G4Exception("model/util/G4Fancy3DNucleus.cc", "mod_util001", FatalException,
             "Problem to place nucleons");
    }

}

Here is the call graph for this function:

Here is the caller graph for this function:

CoulombBarrier()

G4double G4Fancy3DNucleus::CoulombBarrier ( )

virtual

Implements G4V3DNucleus.

Definition at line 501 of file G4Fancy3DNucleus.cc.

{
  static const G4double cfactor = (1.44/1.14) * MeV;
  return cfactor*myZ/(1.0 + G4Pow::GetInstance()->Z13(myA));
}

Here is the call graph for this function:

Here is the caller graph for this function:

DoLorentzBoost() [1/2]

void G4Fancy3DNucleus::DoLorentzBoost ( const G4LorentzVector &  theBoost )

virtual

Implements G4V3DNucleus.

Definition at line 200 of file G4Fancy3DNucleus.cc.

{
    for (G4int i=0; i<myA; i++){
        theNucleons[i].Boost(theBoost);
    }
}

DoLorentzBoost() [2/2]

void G4Fancy3DNucleus::DoLorentzBoost ( const G4ThreeVector &  theBeta )

virtual

Implements G4V3DNucleus.

Definition at line 207 of file G4Fancy3DNucleus.cc.

{
    for (G4int i=0; i<myA; i++){
        theNucleons[i].Boost(theBeta);
    }
}

DoLorentzContraction() [1/2]

void G4Fancy3DNucleus::DoLorentzContraction ( const G4LorentzVector &  theBoost )

virtual

Implements G4V3DNucleus.

Definition at line 228 of file G4Fancy3DNucleus.cc.

{
    if (theBoost.e() !=0 ) {
       G4ThreeVector beta = theBoost.vect()/theBoost.e();
       DoLorentzContraction(beta);
    }
}

Here is the call graph for this function:

DoLorentzContraction() [2/2]

void G4Fancy3DNucleus::DoLorentzContraction ( const G4ThreeVector &  theBeta )

virtual

Implements G4V3DNucleus.

Definition at line 214 of file G4Fancy3DNucleus.cc.

{
    G4double beta2=theBeta.mag2();
    if (beta2 > 0) {
       G4double factor=(1-std::sqrt(1-beta2))/beta2; // (gamma-1)/gamma/beta**2
       G4ThreeVector rprime;
       for (G4int i=0; i< myA; i++) {
         rprime = theNucleons[i].GetPosition() - 
           factor * (theBeta*theNucleons[i].GetPosition()) * theBeta;  
         theNucleons[i].SetPosition(rprime);
       }
    }    
}

Here is the call graph for this function:

DoTranslation()

void G4Fancy3DNucleus::DoTranslation ( const G4ThreeVector &  theShift )

virtual

Implements G4V3DNucleus.

Definition at line 250 of file G4Fancy3DNucleus.cc.

{
  G4ThreeVector tempV;
  for (G4int i=0; i<myA; i++ )
    {
      tempV = theNucleons[i].GetPosition() + theShift;
      theNucleons[i].SetPosition(tempV);
    }   
}

Here is the caller graph for this function:

GetCharge()

G4int G4Fancy3DNucleus::GetCharge ( )

inlinevirtual

Implements G4V3DNucleus.

Definition at line 121 of file G4Fancy3DNucleus.hh.

{
    return myZ;
}

GetExcitationEnergy()

G4double G4Fancy3DNucleus::GetExcitationEnergy ( void  )

inline

Definition at line 136 of file G4Fancy3DNucleus.hh.

{
   return excitationEnergy;
}

GetMass()

G4double G4Fancy3DNucleus::GetMass ( )

virtual

Implements G4V3DNucleus.

Definition at line 191 of file G4Fancy3DNucleus.cc.

{
        return   myZ*G4Proton::Proton()->GetPDGMass() + 
                 (myA-myZ)*G4Neutron::Neutron()->GetPDGMass() -
                 BindingEnergy();
}

Here is the call graph for this function:

GetMassNumber()

G4int G4Fancy3DNucleus::GetMassNumber ( )

inlinevirtual

Implements G4V3DNucleus.

Definition at line 126 of file G4Fancy3DNucleus.hh.

{
    return myA;
}

GetNextNucleon()

G4Nucleon * G4Fancy3DNucleus::GetNextNucleon ( )

virtual

Implements G4V3DNucleus.

Definition at line 126 of file G4Fancy3DNucleus.cc.

{
  return ( (currentNucleon>=0 && currentNucleon<myA) ? 
       &theNucleons[currentNucleon++] : 0 );
}

Here is the caller graph for this function:

GetNuclearDensity()

const G4VNuclearDensity * G4Fancy3DNucleus::GetNuclearDensity ( ) const

virtual

Implements G4V3DNucleus.

Definition at line 260 of file G4Fancy3DNucleus.cc.

{
    return theDensity;
}

GetNuclearRadius() [1/2]

G4double G4Fancy3DNucleus::GetNuclearRadius ( )

virtual

Implements G4V3DNucleus.

Definition at line 167 of file G4Fancy3DNucleus.cc.

{
    return GetNuclearRadius(0.5);
}

Here is the caller graph for this function:

GetNuclearRadius() [2/2]

G4double G4Fancy3DNucleus::GetNuclearRadius ( const G4double  maxRelativeDensity )

virtual

Implements G4V3DNucleus.

Definition at line 172 of file G4Fancy3DNucleus.cc.

{   
    return theDensity->GetRadius(maxRelativeDensity);
}

Here is the call graph for this function:

GetNucleons()

const std::vector< G4Nucleon > & G4Fancy3DNucleus::GetNucleons ( )

virtual

Implements G4V3DNucleus.

Definition at line 132 of file G4Fancy3DNucleus.cc.

{
  return theNucleons;
}

Here is the caller graph for this function:

GetOuterRadius()

G4double G4Fancy3DNucleus::GetOuterRadius ( )

virtual

Implements G4V3DNucleus.

Definition at line 177 of file G4Fancy3DNucleus.cc.

{
    G4double maxradius2=0;
    
    for (int i=0; i<myA; i++)
    {
       if ( theNucleons[i].GetPosition().mag2() > maxradius2 )
       {
        maxradius2=theNucleons[i].GetPosition().mag2();
       }
    }
    return std::sqrt(maxradius2)+nucleondistance;
}

Here is the caller graph for this function:

Init()

void G4Fancy3DNucleus::Init ( G4int  theA, G4int  theZ  )

virtual

Implements G4V3DNucleus.

Definition at line 77 of file G4Fancy3DNucleus.cc.

{
//  G4cout << "G4Fancy3DNucleus::Init(theA, theZ) called"<<G4endl;
  currentNucleon=-1;
  theNucleons.clear();

  myZ = theZ;
  myA= theA;
  excitationEnergy=0;

  theNucleons.resize(myA);  // Pre-loads vector with empty elements
  
//  G4cout << "myA, myZ" << myA << ", " << myZ << G4endl;

  if(theDensity) delete theDensity;
  if ( myA < 17 ) {
     theDensity = new G4NuclearShellModelDensity(myA, myZ);
  } else {
     theDensity = new G4NuclearFermiDensity(myA, myZ);
  }

  theFermi.Init(myA, myZ);
  
  ChooseNucleons();
  
  ChoosePositions();
  
//  CenterNucleons();   // This would introduce a bias 

  ChooseFermiMomenta();
  
  G4double Ebinding= BindingEnergy()/myA;
  
  for (G4int aNucleon=0; aNucleon < myA; aNucleon++)
  {
    theNucleons[aNucleon].SetBindingEnergy(Ebinding);
  }
  
  
  return;
}

Here is the call graph for this function:

Here is the caller graph for this function:

operator!=()

int G4Fancy3DNucleus::operator!= ( const G4Fancy3DNucleus &  right ) const

private

operator=()

const G4Fancy3DNucleus& G4Fancy3DNucleus::operator= ( const G4Fancy3DNucleus &  right )

private

operator==()

int G4Fancy3DNucleus::operator== ( const G4Fancy3DNucleus &  right ) const

private

ReduceSum()

G4bool G4Fancy3DNucleus::ReduceSum ( )

private

Definition at line 418 of file G4Fancy3DNucleus.cc.

{
    G4ThreeVector sum;
    G4double PFermi=fermiM[myA-1];

    for (G4int i=0; i < myA-1 ; i++ )
         { sum+=momentum[i]; }

// check if have to do anything at all..
    if ( sum.mag() <= PFermi )
    {
        momentum[myA-1]=-sum;
        return true;
    }

// find all possible changes in momentum, changing only the component parallel to sum
    G4ThreeVector testDir=sum.unit();
    testSums.clear();
    testSums.resize(myA-1);     // Allocate block for filling below

    G4ThreeVector delta;
    for (G4int aNucleon=0; aNucleon < myA-1; aNucleon++) {
      delta = 2.*((momentum[aNucleon]*testDir)*testDir);

      testSums[aNucleon].Fill(delta, delta.mag(), aNucleon);
    }

    std::sort(testSums.begin(), testSums.end());

//    reduce Momentum Sum until the next would be allowed.
    G4int index=testSums.size();
    while ( (sum-testSums[--index].Vector).mag()>PFermi && index>0)  /* Loop checking, 30-Oct-2015, G.Folger */
    {
      // Only take one which improve, ie. don't change sign and overshoot...
      if ( sum.mag() > (sum-testSums[index].Vector).mag() ) {
        momentum[testSums[index].Index]-=testSums[index].Vector;
        sum-=testSums[index].Vector;
      }
    }

    if ( (sum-testSums[index].Vector).mag() <= PFermi )
    {
        G4int best=-1;
        G4double pBest=2*PFermi; // anything larger than PFermi
        for ( G4int aNucleon=0; aNucleon<=index; aNucleon++)
        {
            // find the momentum closest to choosen momentum for last Nucleon.
            G4double pTry=(testSums[aNucleon].Vector-sum).mag();
            if ( pTry < PFermi
             &&  std::abs(momentum[myA-1].mag() - pTry ) < pBest )
                {
               pBest=std::abs(momentum[myA-1].mag() - pTry );
               best=aNucleon;
            }
        }
        if ( best < 0 )  
        {
          G4String text = "G4Fancy3DNucleus.cc: Logic error in ReduceSum()";
              throw G4HadronicException(__FILE__, __LINE__, text);
        }
        momentum[testSums[best].Index]-=testSums[best].Vector;
        momentum[myA-1]=testSums[best].Vector-sum;

        return true;

    }

    // try to compensate momentum using another Nucleon....
    G4int swapit=-1;
    while (swapit<myA-1)  /* Loop checking, 30-Oct-2015, G.Folger */
    {
      if ( fermiM[++swapit] > PFermi ) break;
    }
    if (swapit == myA-1 ) return false;

    // Now we have a nucleon with a bigger Fermi Momentum.
    // Exchange with last nucleon.. and iterate.
    std::swap(theNucleons[swapit], theNucleons[myA-1]);
    std::swap(momentum[swapit], momentum[myA-1]);
    std::swap(fermiM[swapit], fermiM[myA-1]);
    return ReduceSum();
}

Here is the call graph for this function:

Here is the caller graph for this function:

SortNucleonsDecZ()

void G4Fancy3DNucleus::SortNucleonsDecZ ( )

virtual

Implements G4V3DNucleus.

Definition at line 152 of file G4Fancy3DNucleus.cc.

{
  if (theNucleons.size() < 2 ) return;      // Avoid unnecessary work
  SortNucleonsIncZ();

  std::reverse(theNucleons.begin(), theNucleons.end());
}

Here is the call graph for this function:

SortNucleonsIncZ()

void G4Fancy3DNucleus::SortNucleonsIncZ ( )

virtual

Implements G4V3DNucleus.

Definition at line 144 of file G4Fancy3DNucleus.cc.

{
  if (theNucleons.size() < 2 ) return;      // Avoid unnecesary work

  std::sort(theNucleons.begin(), theNucleons.end(),
        G4Fancy3DNucleusHelperForSortInZ); 
}

Here is the call graph for this function:

Here is the caller graph for this function:

StartLoop()

G4bool G4Fancy3DNucleus::StartLoop ( )

virtual

Implements G4V3DNucleus.

Definition at line 119 of file G4Fancy3DNucleus.cc.

{
    currentNucleon=0;
    return (theNucleons.size()>0);
}   

Here is the caller graph for this function:

Member Data Documentation

currentNucleon

G4int G4Fancy3DNucleus::currentNucleon

private

Definition at line 108 of file G4Fancy3DNucleus.hh.

excitationEnergy

G4double G4Fancy3DNucleus::excitationEnergy

private

Definition at line 112 of file G4Fancy3DNucleus.hh.

fermiM

std::vector<G4double> G4Fancy3DNucleus::fermiM

private

Definition at line 116 of file G4Fancy3DNucleus.hh.

momentum

std::vector<G4ThreeVector> G4Fancy3DNucleus::momentum

private

Definition at line 115 of file G4Fancy3DNucleus.hh.

myA

G4int G4Fancy3DNucleus::myA

private

Definition at line 104 of file G4Fancy3DNucleus.hh.

myZ

G4int G4Fancy3DNucleus::myZ

private

Definition at line 105 of file G4Fancy3DNucleus.hh.

nucleondistance

const G4double G4Fancy3DNucleus::nucleondistance

private

Definition at line 111 of file G4Fancy3DNucleus.hh.

places

std::vector<G4ThreeVector> G4Fancy3DNucleus::places

private

Definition at line 114 of file G4Fancy3DNucleus.hh.

testSums

std::vector<G4Fancy3DNucleusHelper> G4Fancy3DNucleus::testSums

private

Definition at line 117 of file G4Fancy3DNucleus.hh.

theDensity

G4VNuclearDensity* G4Fancy3DNucleus::theDensity

private

Definition at line 109 of file G4Fancy3DNucleus.hh.

theFermi

G4FermiMomentum G4Fancy3DNucleus::theFermi

private

Definition at line 110 of file G4Fancy3DNucleus.hh.

theNucleons

std::vector<G4Nucleon> G4Fancy3DNucleus::theNucleons

private

Definition at line 106 of file G4Fancy3DNucleus.hh.

---

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

• G4Fancy3DNucleus.hh
• G4Fancy3DNucleus.cc