G4MuonVDNuclearModel Class Reference

#include <G4MuonVDNuclearModel.hh>

Inheritance diagram for G4MuonVDNuclearModel:

Collaboration diagram for G4MuonVDNuclearModel:

Public Member Functions
	G4MuonVDNuclearModel ()
	~G4MuonVDNuclearModel ()
G4HadFinalState *	ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
Public Member Functions inherited from G4HadronicInteraction
	G4HadronicInteraction (const G4String &modelName="HadronicModel")
virtual	~G4HadronicInteraction ()
virtual G4double	SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
virtual G4bool	IsApplicable (const G4HadProjectile &, G4Nucleus &)
G4double	GetMinEnergy () const
G4double	GetMinEnergy (const G4Material *aMaterial, const G4Element *anElement) const
void	SetMinEnergy (G4double anEnergy)
void	SetMinEnergy (G4double anEnergy, const G4Element *anElement)
void	SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)
G4double	GetMaxEnergy () const
G4double	GetMaxEnergy (const G4Material *aMaterial, const G4Element *anElement) const
void	SetMaxEnergy (const G4double anEnergy)
void	SetMaxEnergy (G4double anEnergy, const G4Element *anElement)
void	SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)
const G4HadronicInteraction *	GetMyPointer () const
virtual G4int	GetVerboseLevel () const
virtual void	SetVerboseLevel (G4int value)
const G4String &	GetModelName () const
void	DeActivateFor (const G4Material *aMaterial)
void	ActivateFor (const G4Material *aMaterial)
void	DeActivateFor (const G4Element *anElement)
void	ActivateFor (const G4Element *anElement)
G4bool	IsBlocked (const G4Material *aMaterial) const
G4bool	IsBlocked (const G4Element *anElement) const
void	SetRecoilEnergyThreshold (G4double val)
G4double	GetRecoilEnergyThreshold () const
G4bool	operator== (const G4HadronicInteraction &right) const
G4bool	operator!= (const G4HadronicInteraction &right) const
virtual const std::pair< G4double, G4double >	GetFatalEnergyCheckLevels () const
virtual std::pair< G4double, G4double >	GetEnergyMomentumCheckLevels () const
void	SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
virtual void	ModelDescription (std::ostream &outFile) const
virtual void	BuildPhysicsTable (const G4ParticleDefinition &)

Private Member Functions
G4DynamicParticle *	CalculateEMVertex (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
void	CalculateHadronicVertex (G4DynamicParticle *incident, G4Nucleus &target)
void	MakeSamplingTable ()

Private Attributes
G4double	CutFixed
G4int	NBIN
G4double	proba [5][8][1001]
G4double	ya [1001]
G4KokoulinMuonNuclearXS *	muNucXS
G4TheoFSGenerator *	ftfp
G4LundStringFragmentation *	theFragmentation
G4ExcitedStringDecay *	theStringDecay
G4CascadeInterface *	bert

Additional Inherited Members
Protected Member Functions inherited from G4HadronicInteraction
void	SetModelName (const G4String &nam)
G4bool	IsBlocked () const
void	Block ()
Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState	theParticleChange
G4int	verboseLevel
G4double	theMinEnergy
G4double	theMaxEnergy
G4bool	isBlocked

Detailed Description

Definition at line 49 of file G4MuonVDNuclearModel.hh.

Constructor & Destructor Documentation

G4MuonVDNuclearModel()

G4MuonVDNuclearModel::G4MuonVDNuclearModel

(

)

Definition at line 56 of file G4MuonVDNuclearModel.cc.

 : G4HadronicInteraction("G4MuonVDNuclearModel")
{
  muNucXS = (G4KokoulinMuonNuclearXS*)G4CrossSectionDataSetRegistry::Instance()->
    GetCrossSectionDataSet(G4KokoulinMuonNuclearXS::Default_Name());

  SetMinEnergy(0.0);
  SetMaxEnergy(1*PeV);
  CutFixed = 0.2*GeV;
  NBIN = 1000;

  for (G4int k = 0; k < 5; k++) {
    for (G4int j = 0; j < 8; j++) {
      for (G4int i = 0; i < 1001; i++) {
        proba[k][j][i] = 0.0;
        ya[i] = 0.0;
      }
    }
  }

  MakeSamplingTable();

  // reuse existing pre-compound model
  G4GeneratorPrecompoundInterface* precoInterface 
    = new G4GeneratorPrecompoundInterface();
  G4HadronicInteraction* p =
    G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
  G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p);
  if(!pre) { pre = new G4PreCompoundModel(); }
  precoInterface->SetDeExcitation(pre);

  // Build FTFP model
  ftfp = new G4TheoFSGenerator();
  ftfp->SetTransport(precoInterface);
  theFragmentation = new G4LundStringFragmentation();
  theStringDecay = new G4ExcitedStringDecay(theFragmentation);    
  G4FTFModel* theStringModel = new G4FTFModel;
  theStringModel->SetFragmentationModel(theStringDecay);
  ftfp->SetHighEnergyGenerator(theStringModel);

  // Build Bertini cascade
  bert = new G4CascadeInterface();
}

Here is the call graph for this function:

~G4MuonVDNuclearModel()

G4MuonVDNuclearModel::~G4MuonVDNuclearModel

(

)

Definition at line 101 of file G4MuonVDNuclearModel.cc.

{
  delete theFragmentation;
  delete theStringDecay;
}

Member Function Documentation

ApplyYourself()

G4HadFinalState * G4MuonVDNuclearModel::ApplyYourself ( const G4HadProjectile &  aTrack, G4Nucleus &  targetNucleus  )

virtual

Implements G4HadronicInteraction.

Definition at line 109 of file G4MuonVDNuclearModel.cc.

{
  theParticleChange.Clear();

  // For very low energy, return initial track
  G4double epmax = aTrack.GetTotalEnergy() - 0.5*proton_mass_c2;
  if (epmax <= CutFixed) {
    theParticleChange.SetStatusChange(isAlive);
    theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
    theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
    return &theParticleChange;
  }

  // Produce recoil muon and transferred photon
  G4DynamicParticle* transferredPhoton = CalculateEMVertex(aTrack, targetNucleus);

  // Interact the gamma with the nucleus
  CalculateHadronicVertex(transferredPhoton, targetNucleus);
  return &theParticleChange;
}

Here is the call graph for this function:

CalculateEMVertex()

G4DynamicParticle * G4MuonVDNuclearModel::CalculateEMVertex ( const G4HadProjectile &  aTrack, G4Nucleus &  targetNucleus  )

private

Definition at line 133 of file G4MuonVDNuclearModel.cc.

{
  // Select sampling table
  G4double KineticEnergy = aTrack.GetKineticEnergy();
  G4double TotalEnergy = aTrack.GetTotalEnergy();
  G4double Mass = G4MuonMinus::MuonMinus()->GetPDGMass();
  G4double lnZ = G4Log(G4double(targetNucleus.GetZ_asInt() ) );

  G4double epmin = CutFixed;
  G4double epmax = TotalEnergy - 0.5*proton_mass_c2;
  G4double m0 = 0.2*GeV;

  G4double delmin = 1.e10;
  G4double del;
  G4int izz = 0;
  G4int itt = 0;
  G4int NBINminus1 = NBIN - 1;

  G4int nzdat = 5;
  G4double zdat[] = {1.,4.,13.,29.,92.};
  for (G4int iz = 0; iz < nzdat; iz++) {
    del = std::abs(lnZ-G4Log(zdat[iz]));
    if (del < delmin) {
      delmin = del;
      izz = iz;
    }
  }
 
  G4int ntdat = 8;
  G4double tdat[] = {1.e3,1.e4,1.e5,1.e6,1.e7,1.e8,1.e9,1.e10};
  delmin = 1.e10;
  for (G4int it = 0; it < ntdat; it++) {
    del = std::abs(G4Log(KineticEnergy)-G4Log(tdat[it]) );
    if (del < delmin) {
      delmin = del;
      itt = it;
    }
  }

  // Sample the energy transfer according to the probability table
  G4double r = G4UniformRand();

  G4int iy = -1;
  G4ExceptionDescription ed;
  ed << " While count exceeded " << G4endl;
 
  do {
       iy += 1 ;
       if (iy > 10000) {
         G4Exception("G4RPGAntiProtonInelastic::Cascade()", "HAD_RPG_100", JustWarning, ed);
         break;
       }
     } while (((proba[izz][itt][iy]) < r)&&(iy < NBINminus1)) ;  /* Loop checking, 01.09.2015, D.Wright */

  // Sampling is done uniformly in y in the bin

  G4double y; 
  if (iy < NBIN)
    y = ya[iy] + G4UniformRand() * (ya[iy+1] - ya[iy]);
  else
    y = ya[iy];

  G4double x = G4Exp(y);
  G4double ep = epmin*G4Exp(x*G4Log(epmax/epmin) );

  // Sample scattering angle of mu, but first t should be sampled.
  G4double yy = ep/TotalEnergy;
  G4double tmin = Mass*Mass*yy*yy/(1.-yy);
  G4double tmax = 2.*proton_mass_c2*ep;
  G4double t1;
  G4double t2;
  if (m0 < ep) {
    t1 = m0*m0;
    t2 = ep*ep;
  } else {
    t1 = ep*ep;
    t2 = m0*m0;
  }

  G4double w1 = tmax*t1;
  G4double w2 = tmax+t1;
  G4double w3 = tmax*(tmin+t1)/(tmin*w2);
  G4double y1 = 1.-yy;
  G4double y2 = 0.5*yy*yy;
  G4double y3 = y1+y2;

  G4double t;
  G4double rej;

  // Now sample t
  G4int ntry = 0;
  G4ExceptionDescription eda;
  eda << " While count exceeded " << G4endl;
 
  do
  {
    ntry += 1;
    if (ntry > 10000) {
      G4Exception("G4MuonVDNuclearModel::CalculateEMVertex()", "HAD_RPG_100", JustWarning, eda);
      break;
    }
  
    t = w1/(w2*G4Exp(G4UniformRand()*G4Log(w3))-tmax);
    rej = (1.-t/tmax)*(y1*(1.-tmin/t)+y2)/(y3*(1.-t/t2)); 
  } while (G4UniformRand() > rej) ;   /* Loop checking, 01.09.2015, D.Wright */

  // compute angle from t
  G4double sinth2 =
             0.5*(t-tmin)/(2.*(TotalEnergy*(TotalEnergy-ep)-Mass*Mass)-tmin);
  G4double theta = std::acos(1. - 2.*sinth2);

  G4double phi = twopi*G4UniformRand();
  G4double sinth = std::sin(theta);
  G4double dirx = sinth*std::cos(phi);
  G4double diry = sinth*std::sin(phi);
  G4double dirz = std::cos(theta);
  G4ThreeVector finalDirection(dirx,diry,dirz);
  G4ThreeVector ParticleDirection(aTrack.Get4Momentum().vect().unit() );
  finalDirection.rotateUz(ParticleDirection);

  G4double NewKinEnergy = KineticEnergy - ep;
  G4double finalMomentum = std::sqrt(NewKinEnergy*(NewKinEnergy+2.*Mass) );
  G4double Ef = NewKinEnergy + Mass;
  G4double initMomentum = std::sqrt(KineticEnergy*(TotalEnergy+Mass) );

  // Set energy and direction of scattered primary in theParticleChange
  theParticleChange.SetStatusChange(isAlive);
  theParticleChange.SetEnergyChange(NewKinEnergy);
  theParticleChange.SetMomentumChange(finalDirection);

  // Now create the emitted gamma 
  G4LorentzVector primaryMomentum(initMomentum*ParticleDirection, TotalEnergy);
  G4LorentzVector fsMomentum(finalMomentum*finalDirection, Ef);
  G4LorentzVector momentumTransfer = primaryMomentum - fsMomentum;

  G4DynamicParticle* gamma = 
           new G4DynamicParticle(G4Gamma::Gamma(), momentumTransfer);
 
  return gamma;
}

Here is the call graph for this function:

Here is the caller graph for this function:

CalculateHadronicVertex()

void G4MuonVDNuclearModel::CalculateHadronicVertex ( G4DynamicParticle *  incident, G4Nucleus &  target  )

private

Definition at line 277 of file G4MuonVDNuclearModel.cc.

{
  G4HadFinalState* hfs = 0;
  G4double gammaE = incident->GetTotalEnergy();

  if (gammaE < 10*GeV) {
    G4HadProjectile projectile(*incident);
    hfs = bert->ApplyYourself(projectile, target);
  } else {
    // convert incident gamma to a pi0
    G4double piMass = G4PionZero::PionZero()->GetPDGMass();
    G4double piKE = incident->GetTotalEnergy() - piMass;
    G4double piMom = std::sqrt(piKE*(piKE + 2*piMass) );
    G4ThreeVector piMomentum(incident->GetMomentumDirection() );
    piMomentum *= piMom;
    G4DynamicParticle theHadron(G4PionZero::PionZero(), piMomentum);
    G4HadProjectile projectile(theHadron);
    hfs = ftfp->ApplyYourself(projectile, target);
  }

  delete incident;

  // Copy secondaries from sub-model to model
  theParticleChange.AddSecondaries(hfs);
} 

Here is the call graph for this function:

Here is the caller graph for this function:

MakeSamplingTable()

void G4MuonVDNuclearModel::MakeSamplingTable ( )

private

Definition at line 305 of file G4MuonVDNuclearModel.cc.

{
  G4double adat[] = {1.01,9.01,26.98,63.55,238.03};
  G4double zdat[] = {1.,4.,13.,29.,92.};
  G4int nzdat = 5;

  G4double tdat[] = {1.e3,1.e4,1.e5,1.e6,1.e7,1.e8,1.e9,1.e10};
  G4int ntdat = 8;

  G4int nbin;
  G4double KineticEnergy;
  G4double TotalEnergy;
  G4double Maxep;
  G4double CrossSection;

  G4double c;
  G4double y;
  G4double ymin,ymax;
  G4double dy,yy;
  G4double dx,x;
  G4double ep;

  G4double AtomicNumber;
  G4double AtomicWeight;

  for (G4int iz = 0; iz < nzdat; iz++) {
    AtomicNumber = zdat[iz];
    AtomicWeight = adat[iz]*(g/mole);  

    for (G4int it = 0; it < ntdat; it++) {
      KineticEnergy = tdat[it];
      TotalEnergy = KineticEnergy + G4MuonMinus::MuonMinus()->GetPDGMass();
      Maxep = TotalEnergy - 0.5*proton_mass_c2;

      CrossSection = 0.0;

      // Calculate the differential cross section
      // numerical integration in log .........
      c = G4Log(Maxep/CutFixed);
      ymin = -5.0;
      ymax = 0.0;
      dy = (ymax-ymin)/NBIN; 

      nbin=-1;              

      y = ymin - 0.5*dy;
      yy = ymin - dy;
      for (G4int i = 0; i < NBIN; i++) {
        y += dy;
        x = G4Exp(y);
        yy += dy;
        dx = G4Exp(yy+dy)-G4Exp(yy);
      
        ep = CutFixed*G4Exp(c*x);

        CrossSection +=
           ep*dx*muNucXS->ComputeDDMicroscopicCrossSection(KineticEnergy,
                               AtomicNumber,
                               AtomicWeight, ep);
        if (nbin < NBIN) {
          nbin += 1;
          ya[nbin] = y;
          proba[iz][it][nbin] = CrossSection;
        }
      }
      ya[NBIN] = 0.;
 
      if (CrossSection > 0.0) {
        for (G4int ib = 0; ib <= nbin; ib++) proba[iz][it][ib] /= CrossSection;
      }
    } // loop on it
  } // loop on iz

  // G4cout << " Kokoulin XS = "
  //       << muNucXS->ComputeDDMicroscopicCrossSection(1*GeV, 20.0, 
  //         40.0*g/mole, 0.3*GeV)/millibarn
  //       << G4endl; 
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

bert

G4CascadeInterface* G4MuonVDNuclearModel::bert

private

Definition at line 80 of file G4MuonVDNuclearModel.hh.

CutFixed

G4double G4MuonVDNuclearModel::CutFixed

private

Definition at line 69 of file G4MuonVDNuclearModel.hh.

ftfp

G4TheoFSGenerator* G4MuonVDNuclearModel::ftfp

private

Definition at line 77 of file G4MuonVDNuclearModel.hh.

muNucXS

G4KokoulinMuonNuclearXS* G4MuonVDNuclearModel::muNucXS

private

Definition at line 75 of file G4MuonVDNuclearModel.hh.

NBIN

G4int G4MuonVDNuclearModel::NBIN

private

Definition at line 70 of file G4MuonVDNuclearModel.hh.

proba

G4double G4MuonVDNuclearModel::proba[5][8][1001]

private

Definition at line 72 of file G4MuonVDNuclearModel.hh.

theFragmentation

G4LundStringFragmentation* G4MuonVDNuclearModel::theFragmentation

private

Definition at line 78 of file G4MuonVDNuclearModel.hh.

theStringDecay

G4ExcitedStringDecay* G4MuonVDNuclearModel::theStringDecay

private

Definition at line 79 of file G4MuonVDNuclearModel.hh.

ya

G4double G4MuonVDNuclearModel::ya[1001]

private

Definition at line 73 of file G4MuonVDNuclearModel.hh.

---

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

• G4MuonVDNuclearModel.hh
• G4MuonVDNuclearModel.cc