Effective gas: All quantities necessary to calculate PAI ionisation for gas mixture. More...

#include <TRT_PAI_effectiveGas.h>

Inheritance diagram for TRT_PAI_effectiveGas:

Collaboration diagram for TRT_PAI_effectiveGas:

Public Member Functions
	TRT_PAI_effectiveGas (TRT_PAI_gasMixture *gm, double Emin, double Emax, double tempK, double eps)
void	GasTab (const std::vector< float > &gamvec, std::vector< float > &EArray, std::vector< std::vector< float > > &fnArray, std::vector< float > &dndx)
	Tabulate double differential distribution.
bool	msgLvl (const MSG::Level lvl) const
	Test the output level.
MsgStream &	msg () const
	The standard message stream.
MsgStream &	msg (const MSG::Level lvl) const
	The standard message stream.
void	setLevel (MSG::Level lvl)
	Change the current logging level.

Private Member Functions
double	XSigma (double lnE, double dummy)
double	XFReal (double lnD, double lnE)
double	XGInt (double(TRT_PAI_effectiveGas::*pt2Func)(double, double), double lnLo, double lnHi, const double eps, double extraParameter)
double	dndedx (double lgE, double gamma)
void	initMessaging () const
	Initialize our message level and MessageSvc.

Private Attributes
std::vector< float >	m_lnELvls
std::vector< float >	m_lnFosc
std::vector< float >	m_lnIntegratedSigmas
std::vector< float >	m_lnEpsR
std::vector< float >	m_lnEpsI
const double	m_lnEmin {}
const double	m_lnEmax {}
const double	m_eps {}
double	m_S1 {}
double	m_Wp2 {}
double	m_S2 {}
double	m_ne {}
std::string	m_nm
	Message source name.
boost::thread_specific_ptr< MsgStream >	m_msg_tls
	MsgStream instance (a std::cout like with print-out levels)
std::atomic< IMessageSvc * >	m_imsg { nullptr }
	MessageSvc pointer.
std::atomic< MSG::Level >	m_lvl { MSG::NIL }
	Current logging level.
std::atomic_flag m_initialized	ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT
	Messaging initialized (initMessaging)

Detailed Description

Effective gas: All quantities necessary to calculate PAI ionisation for gas mixture.

Definition at line 18 of file TRT_PAI_effectiveGas.h.

Constructor & Destructor Documentation

◆ TRT_PAI_effectiveGas()

TRT_PAI_effectiveGas::TRT_PAI_effectiveGas	(	TRT_PAI_gasMixture *	gm,
		double	Emin,
		double	Emax,
		double	tempK,
		double	eps )

Parameters

gm	gas mixture
Emin	minimum energy transfer considered in calculations [eV]
Emax	maximum energy transfer considered in calculations [ev]
tempK	gas temperature [Kelvin]
eps	epsilon for numeric integration
mlog	Message stream

Definition at line 24 of file TRT_PAI_effectiveGas.cxx.

  : AthMessaging("TRT_PAI_effectiveGas"),
    m_lnEmin(std::log(Emin)),
    m_lnEmax(std::log(Emax)),
    m_eps(eps)
{
  // Tell clang to optimize assuming that FP may trap.
  CXXUTILS_TRAPPING_FP;
 
  using namespace TRT_PAI_physicsConstants;
 
  TRT_PAI_element* pe;
  double rho  = 0.;
  double Zeff = 0.;
  double Aeff = 0.;
  double w    = 0.;
  double wtot = 0.;
  int Nelem   = gm->getNElements();
 
  for ( int k=0; k<Nelem; k++ ) {
    pe    = gm->getElement(k);
    w     = gm->getElemWeight(k);
    wtot += w;
    rho  += w * pe->getDensity(tempK);
    Aeff += w * pe->getAtomicA();
    Zeff += w * pe->getAtomicZ();
  }
  if (wtot == 0.)[[unlikely]]{
    ATH_MSG_ERROR("TRT_PAI_effectiveGas::TRT_PAI_effectiveGas: wtot is zero.");
    return;
  }
  Aeff /= wtot;
  Zeff /= wtot;
 
  m_ne  = Nav*rho*Zeff/Aeff;         // Electron density
  m_Wp2 = 4.*M_PI*r0*m_ne*he*he;       // plasma freq**2 {ev}
  m_S1  = mb/(2.*M_PI*M_PI*r0*he*Zeff);  // x-section to F.osc
  m_S2  = 2.*M_PI*r0*m_ne*Qe*Qe/erg;   // dN/dx scale
 
  ATH_MSG_DEBUG ( "effectiveGas: Electron density     " << m_ne );
  ATH_MSG_DEBUG ( "effectiveGas: Plasma freq**2 {eV}  " << m_Wp2 );
  ATH_MSG_DEBUG ( "effectiveGas: x-section to F.osc   " << m_S1 );
  ATH_MSG_DEBUG ( "effectiveGas: dN/dx scale          " << m_S2 );
 
  // Merge all energy levels into ELvls
  std::vector<float> tempv;
  for ( int k=0; k<Nelem; k++ ) {
    pe    = gm->getElement(k);
    tempv = pe->getLnELvls();
    m_lnELvls.insert( m_lnELvls.end(), tempv.begin(), tempv.end() );
  }
 
  // Insert also Emin and Emax into ELvls, then sort
  m_lnELvls.push_back(m_lnEmin);
  m_lnELvls.push_back(m_lnEmax);
  sort(m_lnELvls.begin(),m_lnELvls.end());
 
  std::vector<float>::iterator vit;
 
  // Remove duplicate energies
  vit = unique(m_lnELvls.begin(), m_lnELvls.end());
  m_lnELvls.erase(vit, m_lnELvls.end());
 
  // Remove everything below Emin
  vit = lower_bound(m_lnELvls.begin(), m_lnELvls.end(), m_lnEmin);
  m_lnELvls.erase(m_lnELvls.begin(), vit);
 
  // Remove everything above Emax
  vit = upper_bound(m_lnELvls.begin(), m_lnELvls.end(), m_lnEmax)-1;
  m_lnELvls.erase(vit, m_lnELvls.end());
 
  int NLvls = m_lnELvls.size();
 
 
  // Expand vectors to correct dimension
  m_lnFosc.resize(NLvls);
  m_lnIntegratedSigmas.resize(NLvls);
  m_lnEpsI.resize(NLvls);
  m_lnEpsR.resize(NLvls);
 
  // create array of effective cross sections (Fosc).
 
  for ( int i=0; i<NLvls; i++ ) {
    double fosc = 0.;
    // all atoms with an absorbtion energylevel low enough contribute.
    for ( int k=0; k<Nelem; k++ ) {
      pe = gm->getElement(k);
      if ( m_lnELvls[i] >= pe->getLnELvls()[0] ) {
    double lnSig = TRT_PAI_utils::Interpolate( m_lnELvls[i],
                           pe->getLnELvls(),
                           pe->getLnSigmas() );
    fosc += std::exp(lnSig) * gm->getElemWeight(k);
      }
    }
    //clamp, just in case
    fosc = std::clamp(fosc, 1e-300,1e+300);
    m_lnFosc[i] = std::log(fosc);
  }
 
  // Create array of integrated cross-sections
  double sigma = 0.;
  double lnSigma = 0.;
  double lnEHigh;
  double lnELow = m_lnEmax;
  double dsigma;
 
  for ( int i=NLvls-1; i>=0; --i) {
    lnEHigh = lnELow;
    lnELow  = m_lnELvls[i];
    dsigma  = XGInt( &TRT_PAI_effectiveGas::XSigma,
             lnELow,
             lnEHigh,
             m_eps,
             0.);
    sigma  += dsigma;
    lnSigma = std::log(sigma);
    m_lnIntegratedSigmas[i] = lnSigma;
  }
 
  // Normalize
  for ( int i=0; i<NLvls; i++) {
    m_lnFosc[i]             -= lnSigma;
    m_lnIntegratedSigmas[i] -= lnSigma;
  }
 
  double cnst = std::log(M_PI_2*m_Wp2);
  for (int i=0; i<NLvls; i++) {
    m_lnEpsI[i]  =  cnst - m_lnELvls[i] + m_lnFosc[i];
    float xint = XGInt(&TRT_PAI_effectiveGas::XFReal,
               0.,
               m_lnEmax,
               m_eps,
               m_lnELvls[i] );
    m_lnEpsR[i]  = m_Wp2 * std::exp(m_lnELvls[i]) * xint;
  }
  return;
}

Member Function Documentation

◆ dndedx()

double TRT_PAI_effectiveGas::dndedx	(	double	lgE,
		double	gamma )

private

Definition at line 244 of file TRT_PAI_effectiveGas.cxx.

                                                            {
 
  using namespace TRT_PAI_physicsConstants;
  using namespace TRT_PAI_utils;
 
  double E       = std::exp(lnE);
  double gammaSq = gamma*gamma;
  double betaSq  = 1.-1./gammaSq;
 
  double er = Interpolate(lnE, m_lnELvls, m_lnEpsR);
  double ei = Interpolate(lnE, m_lnELvls, m_lnEpsI);
 
  std::complex<double> Ceps1(er/(E*E), std::exp(ei));
  std::complex<double> C1 = 1./gammaSq - Ceps1*betaSq;
  std::complex<double> C2 = C1/(1.+Ceps1) * std::log(2.*betaSq*MeeV/(E*C1));
 
  double x;
  x = Interpolate(lnE, m_lnELvls, m_lnIntegratedSigmas);
  x = m_S2/betaSq * E * ( -2.*imag(C2)/(m_Wp2*M_PI) + (1.-std::exp(x))/(E*E) );
 
  return x;
}

◆ GasTab()

void TRT_PAI_effectiveGas::GasTab	(	const std::vector< float > &	gamvec,
		std::vector< float > &	EArray,
		std::vector< std::vector< float > > &	fnArray,
		std::vector< float > &	dndx )

Tabulate double differential distribution.

Parameters

gamvec	vector of gamma values for which to tabulate
Earray	vector of energy transfer values
fnArray	double differential distribution
dndx	This array effectively contains the inverse of the mean free path as function of binned energy variable

Definition at line 269 of file TRT_PAI_effectiveGas.cxx.

{
  // For a given gamma factor, return dn/dx
  //     (by integrating dn^2/dEdx with respect to E)
  // and tabulated values of dn^2/dEdx(?)....
 
  double Rener = 0.05;
 
  int NLvls    = m_lnELvls.size();
  int nGamVals = gamvec.size();
 
  double lnEi = m_lnEmax;
  double lnEo;
  double lnEs = -99999999.;
 
  fnArray.resize(nGamVals);
  int nEVals = 0;
 
  // Initialize
  for ( int ig=0; ig<nGamVals; ++ig ) {
    fnArray[ig].push_back(0.);
  }
 
  // Calculate integral from above and store the cumulative values in fnArray
 
  for ( int ie=NLvls-1; ie>=0; --ie ) {
    lnEo = lnEi;
    lnEi = m_lnELvls[ie];
    for ( int ig=0; ig<nGamVals; ++ig ) {
      double ds =
    XGInt( &TRT_PAI_effectiveGas::dndedx, lnEi, lnEo, m_eps, gamvec[ig]);
      fnArray[ig][nEVals] += ds;
    }
    if ( std::abs(lnEs-lnEi) > Rener || ie==0 ) {
      lnEs = lnEi;
      lnEarray.push_back( lnEs );
      if ( ie>0 ) {
    for ( int ig=0; ig<nGamVals; ++ig ) {
      fnArray[ig].push_back( fnArray[ig][nEVals] );
    }
      }
      nEVals++;
    }
  }
 
  // Copy the total integral into auxillary vector
  dndx.resize(nGamVals);
  for ( int ig = 0; ig < nGamVals; ++ig ) {
    dndx[ig] = fnArray[ig][nEVals-1];
  }
 
  return;
}

◆ initMessaging()

void AthMessaging::initMessaging ( ) const

privateinherited

Initialize our message level and MessageSvc.

This method should only be called once.

Definition at line 39 of file AthMessaging.cxx.

{
  m_imsg = Athena::getMessageSvc();
  // If user did not set an explicit level, set a default
  if (m_lvl == MSG::NIL) {
    m_lvl = m_imsg ?
      static_cast<MSG::Level>( m_imsg.load()->outputLevel(m_nm) ) :
      MSG::INFO;
  }
}

◆ msg() [1/2]

MsgStream & AthMessaging::msg ( ) const

inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 163 of file AthMessaging.h.

{
  MsgStream* ms = m_msg_tls.get();
  if (!ms) {
    if (!m_initialized.test_and_set()) initMessaging();
    ms = new MsgStream(m_imsg,m_nm);
    m_msg_tls.reset( ms );
  }
 
  ms->setLevel (m_lvl);
  return *ms;
}

◆ msg() [2/2]

MsgStream & AthMessaging::msg ( const MSG::Level lvl ) const

inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 178 of file AthMessaging.h.

179{ return msg() << lvl; }

AthMessaging::msg

MsgStream & msg() const

The standard message stream.

Definition AthMessaging.h:163

◆ msgLvl()

bool AthMessaging::msgLvl ( const MSG::Level lvl ) const

inlineinherited

Test the output level.

Parameters

lvl	The message level to test against

Returns: boolean Indicating if messages at given level will be printed

Return values

true	Messages at level "lvl" will be printed

Definition at line 151 of file AthMessaging.h.

{
  if (m_lvl <= lvl) {
    msg() << lvl;
    return true;
  } else {
    return false;
  }
}

◆ setLevel()

void AthMessaging::setLevel ( MSG::Level lvl )

inherited

Change the current logging level.

Use this rather than msg().setLevel() for proper operation with MT.

Definition at line 28 of file AthMessaging.cxx.

{
  m_lvl = lvl;
}

◆ XFReal()

double TRT_PAI_effectiveGas::XFReal	(	double	lnD,
		double	lnE )

private

Definition at line 180 of file TRT_PAI_effectiveGas.cxx.

                                                           {
 
  double fp = 0.;
  if ( lnE+lnD > m_lnEmin ) {
    fp = TRT_PAI_utils::Interpolate(lnE+lnD, m_lnELvls, m_lnFosc);
    fp = std::exp(std::max(fp, -99.0));
  }
 
  double fm = 0.;
  if (lnE-lnD > m_lnEmin ) {
    fm = TRT_PAI_utils::Interpolate(lnE-lnD, m_lnELvls, m_lnFosc);
    fm = std::exp(std::max(fm, -99.0));
  }
 
  double x = std::exp(lnD);
  return x/(x*x-1.)*(fp-fm);
}

◆ XGInt()

double TRT_PAI_effectiveGas::XGInt	(	double(TRT_PAI_effectiveGas::*	pt2Func )(double, double),
		double	lnLo,
		double	lnHi,
		const double	eps,
		double	extraParameter )

private

Definition at line 200 of file TRT_PAI_effectiveGas.cxx.

                                             {
  //
  // Pavel Nevskis "famous integration procedure"
  //
  // Applies the Gauss-Legendre Quadrature method with n=4
 
  // These constants belong to the Gauss-Legendre quadrature method:
  const int M = 4;                                              // 4 points
  const double U[M]={-.8611363,-.3399810, .3399810 ,.8611363};  // abscissas
  const double W[M]={ .3478548, .6521452, .6521452, .3478548};  // weights
 
  int N=10;
  double OTB=0,Y,D,result;
 
  double Dt = 0.5*(lnHi-lnLo);
 
  do {
    Y   = OTB;
    OTB = 0.;
    D   = Dt/N;
    for (int i=1; i<=N; i++) {
      for (int k=0; k<M; k++) {
    double arg = lnLo + D*(2*i-1+U[k]);
    OTB += W[k] * (this->*pt2Func)(arg, extraParameter );
      }
    }
    OTB *= D;
    result = OTB;
    N *= 2;
    if ( N>100000 ) {
      ATH_MSG_WARNING( "effectiveGas::XGInt: integrate Divergence! " << std::abs(OTB-Y) << " " << std::abs(eps*OTB) );
      break;
    }
  } while ( (std::abs(OTB-Y) > std::abs(eps*OTB)) && (eps>0) );
 
  return result;
}

◆ XSigma()

double TRT_PAI_effectiveGas::XSigma	(	double	lnE,
		double	dummy )

private

Definition at line 168 of file TRT_PAI_effectiveGas.cxx.

                                                             {
  double xsig = dummy;
  xsig = TRT_PAI_utils::Interpolate(std::max(lnE,m_lnEmin),
                    m_lnELvls,
                    m_lnFosc);
  xsig = std::max(xsig,-99.);
  xsig = std::exp( xsig + lnE );
  return xsig;
}

Member Data Documentation

◆ ATLAS_THREAD_SAFE

std::atomic_flag m_initialized AthMessaging::ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT

mutableprivateinherited

Messaging initialized (initMessaging)

Definition at line 141 of file AthMessaging.h.

◆ m_eps

const double TRT_PAI_effectiveGas::m_eps {}

private

Definition at line 62 of file TRT_PAI_effectiveGas.h.

62{}; // epsilon for numerical integration.

◆ m_imsg

std::atomic<IMessageSvc*> AthMessaging::m_imsg { nullptr }

mutableprivateinherited

MessageSvc pointer.

Definition at line 135 of file AthMessaging.h.

135{ nullptr };

◆ m_lnELvls

std::vector<float> TRT_PAI_effectiveGas::m_lnELvls

private

Definition at line 55 of file TRT_PAI_effectiveGas.h.

◆ m_lnEmax

const double TRT_PAI_effectiveGas::m_lnEmax {}

private

Definition at line 61 of file TRT_PAI_effectiveGas.h.

61{};

◆ m_lnEmin

const double TRT_PAI_effectiveGas::m_lnEmin {}

private

Definition at line 60 of file TRT_PAI_effectiveGas.h.

60{};

◆ m_lnEpsI

std::vector<float> TRT_PAI_effectiveGas::m_lnEpsI

private

Definition at line 59 of file TRT_PAI_effectiveGas.h.

◆ m_lnEpsR

std::vector<float> TRT_PAI_effectiveGas::m_lnEpsR

private

Definition at line 58 of file TRT_PAI_effectiveGas.h.

◆ m_lnFosc

std::vector<float> TRT_PAI_effectiveGas::m_lnFosc

private

Definition at line 56 of file TRT_PAI_effectiveGas.h.

◆ m_lnIntegratedSigmas

std::vector<float> TRT_PAI_effectiveGas::m_lnIntegratedSigmas

private

Definition at line 57 of file TRT_PAI_effectiveGas.h.

◆ m_lvl

std::atomic<MSG::Level> AthMessaging::m_lvl { MSG::NIL }

mutableprivateinherited

Current logging level.

Definition at line 138 of file AthMessaging.h.

138{ MSG::NIL };

◆ m_msg_tls

boost::thread_specific_ptr<MsgStream> AthMessaging::m_msg_tls

mutableprivateinherited

MsgStream instance (a std::cout like with print-out levels)

Definition at line 132 of file AthMessaging.h.

◆ m_ne

double TRT_PAI_effectiveGas::m_ne {}

private

Definition at line 66 of file TRT_PAI_effectiveGas.h.

66{} ; // Electron density

◆ m_nm

std::string AthMessaging::m_nm

privateinherited

Message source name.

Definition at line 129 of file AthMessaging.h.

◆ m_S1

double TRT_PAI_effectiveGas::m_S1 {}

private

Definition at line 63 of file TRT_PAI_effectiveGas.h.

63{} ; // x-section to F.osc

◆ m_S2

double TRT_PAI_effectiveGas::m_S2 {}

private

Definition at line 65 of file TRT_PAI_effectiveGas.h.

65{} ; // dN/dx scale

◆ m_Wp2

double TRT_PAI_effectiveGas::m_Wp2 {}

private

Definition at line 64 of file TRT_PAI_effectiveGas.h.

64{}; // plasma freq**2 {ev}

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

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ TRT_PAI_effectiveGas()

Member Function Documentation

◆ dndedx()

◆ GasTab()

◆ initMessaging()

◆ msg() [1/2]

◆ msg() [2/2]

◆ msgLvl()

◆ setLevel()

◆ XFReal()

◆ XGInt()

◆ XSigma()

Member Data Documentation

◆ ATLAS_THREAD_SAFE

◆ m_eps

◆ m_imsg

◆ m_lnELvls

◆ m_lnEmax

◆ m_lnEmin

◆ m_lnEpsI

◆ m_lnEpsR

◆ m_lnFosc

◆ m_lnIntegratedSigmas

◆ m_lvl

◆ m_msg_tls

◆ m_ne

◆ m_nm

◆ m_S1

◆ m_S2

◆ m_Wp2