ATLAS Offline Software
Loading...
Searching...
No Matches
TRT_PAI_effectiveGas Class Reference

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
gmgas mixture
Eminminimum energy transfer considered in calculations [eV]
Emaxmaximum energy transfer considered in calculations [ev]
tempKgas temperature [Kelvin]
epsepsilon for numeric integration
mlogMessage stream

Definition at line 23 of file TRT_PAI_effectiveGas.cxx.

28 : AthMessaging("TRT_PAI_effectiveGas"),
29 m_lnEmin(std::log(Emin)),
30 m_lnEmax(std::log(Emax)),
31 m_eps(eps)
32{
33 using namespace TRT_PAI_physicsConstants;
34
35 TRT_PAI_element* pe;
36 double rho = 0.;
37 double Zeff = 0.;
38 double Aeff = 0.;
39 double w = 0.;
40 double wtot = 0.;
41 int Nelem = gm->getNElements();
42
43 for ( int k=0; k<Nelem; k++ ) {
44 pe = gm->getElement(k);
45 w = gm->getElemWeight(k);
46 wtot += w;
47 rho += w * pe->getDensity(tempK);
48 Aeff += w * pe->getAtomicA();
49 Zeff += w * pe->getAtomicZ();
50 }
51
52 Aeff /= wtot;
53 Zeff /= wtot;
54
55 m_ne = Nav*rho*Zeff/Aeff; // Electron density
56 m_Wp2 = 4.*M_PI*r0*m_ne*he*he; // plasma freq**2 {ev}
57 m_S1 = mb/(2.*M_PI*M_PI*r0*he*Zeff); // x-section to F.osc
58 m_S2 = 2.*M_PI*r0*m_ne*Qe*Qe/erg; // dN/dx scale
59
60 ATH_MSG_DEBUG ( "effectiveGas: Electron density " << m_ne );
61 ATH_MSG_DEBUG ( "effectiveGas: Plasma freq**2 {eV} " << m_Wp2 );
62 ATH_MSG_DEBUG ( "effectiveGas: x-section to F.osc " << m_S1 );
63 ATH_MSG_DEBUG ( "effectiveGas: dN/dx scale " << m_S2 );
64
65 // Merge all energy levels into ELvls
66 std::vector<float> tempv;
67 for ( int k=0; k<Nelem; k++ ) {
68 pe = gm->getElement(k);
69 tempv = pe->getLnELvls();
70 m_lnELvls.insert( m_lnELvls.end(), tempv.begin(), tempv.end() );
71 }
72
73 // Insert also Emin and Emax into ELvls, then sort
74 m_lnELvls.push_back(m_lnEmin);
75 m_lnELvls.push_back(m_lnEmax);
76 sort(m_lnELvls.begin(),m_lnELvls.end());
77
78 std::vector<float>::iterator vit;
79
80 // Remove duplicate energies
81 vit = unique(m_lnELvls.begin(), m_lnELvls.end());
82 m_lnELvls.erase(vit, m_lnELvls.end());
83
84 // Remove everything below Emin
85 vit = lower_bound(m_lnELvls.begin(), m_lnELvls.end(), m_lnEmin);
86 m_lnELvls.erase(m_lnELvls.begin(), vit);
87
88 // Remove everything above Emax
89 vit = upper_bound(m_lnELvls.begin(), m_lnELvls.end(), m_lnEmax)-1;
90 m_lnELvls.erase(vit, m_lnELvls.end());
91
92 int NLvls = m_lnELvls.size();
93
94
95 // Expand vectors to correct dimension
96 m_lnFosc.resize(NLvls);
97 m_lnIntegratedSigmas.resize(NLvls);
98 m_lnEpsI.resize(NLvls);
99 m_lnEpsR.resize(NLvls);
100
101 // create array of effective cross sections (Fosc).
102
103 for ( int i=0; i<NLvls; i++ ) {
104 double fosc = 0.;
105 // all atoms with an absorbtion energylevel low enough contribute.
106 for ( int k=0; k<Nelem; k++ ) {
107 pe = gm->getElement(k);
108 if ( m_lnELvls[i] >= pe->getLnELvls()[0] ) {
109 double lnSig = TRT_PAI_utils::Interpolate( m_lnELvls[i],
110 pe->getLnELvls(),
111 pe->getLnSigmas() );
112 fosc += std::exp(lnSig) * gm->getElemWeight(k);
113 }
114 }
115 //clamp, just in case
116 fosc = std::clamp(fosc, 1e-300,1e+300);
117 m_lnFosc[i] = std::log(fosc);
118 }
119
120 // Create array of integrated cross-sections
121 double sigma = 0.;
122 double lnSigma = 0.;
123 double lnEHigh;
124 double lnELow = m_lnEmax;
125 double dsigma;
126
127 for ( int i=NLvls-1; i>=0; --i) {
128 lnEHigh = lnELow;
129 lnELow = m_lnELvls[i];
131 lnELow,
132 lnEHigh,
133 m_eps,
134 0.);
135 sigma += dsigma;
136 lnSigma = std::log(sigma);
137 m_lnIntegratedSigmas[i] = lnSigma;
138 }
139
140 // Normalize
141 for ( int i=0; i<NLvls; i++) {
142 m_lnFosc[i] -= lnSigma;
143 m_lnIntegratedSigmas[i] -= lnSigma;
144 }
145
146 double cnst = std::log(M_PI_2*m_Wp2);
147 for (int i=0; i<NLvls; i++) {
148 m_lnEpsI[i] = cnst - m_lnELvls[i] + m_lnFosc[i];
150 0.,
151 m_lnEmax,
152 m_eps,
153 m_lnELvls[i] );
154 m_lnEpsR[i] = m_Wp2 * std::exp(m_lnELvls[i]) * xint;
155 }
156 return;
157}
#define M_PI
#define ATH_MSG_DEBUG(x)
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
AthMessaging()
Default constructor:
std::vector< float > m_lnEpsI
double XFReal(double lnD, double lnE)
double XGInt(double(TRT_PAI_effectiveGas::*pt2Func)(double, double), double lnLo, double lnHi, const double eps, double extraParameter)
std::vector< float > m_lnELvls
double XSigma(double lnE, double dummy)
std::vector< float > m_lnIntegratedSigmas
std::vector< float > m_lnEpsR
std::vector< float > m_lnFosc
TRT_PAI_element * getElement(unsigned int n)
Get element no.
double getElemWeight(unsigned int n)
Get weight of element no.
int getNElements()
Get number of different element in this gas mixture.
const double erg
1 ev to erg {erg}
const double Nav
Avagadro constant.
const double Qe
Electron charge{ESU}.
const double r0
electron radius{cm}
float Interpolate(const float &xval, const std::vector< float > &xtabulated, const std::vector< float > &ytabulated)
Interpolation function.
DataModel_detail::iterator< DVL > unique(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of unique for DataVector/List.

Member Function Documentation

◆ dndedx()

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

Definition at line 237 of file TRT_PAI_effectiveGas.cxx.

237 {
238
239 using namespace TRT_PAI_physicsConstants;
240 using namespace TRT_PAI_utils;
241
242 double E = std::exp(lnE);
243 double gammaSq = gamma*gamma;
244 double betaSq = 1.-1./gammaSq;
245
246 double er = Interpolate(lnE, m_lnELvls, m_lnEpsR);
247 double ei = Interpolate(lnE, m_lnELvls, m_lnEpsI);
248
249 std::complex<double> Ceps1(er/(E*E), std::exp(ei));
250 std::complex<double> C1 = 1./gammaSq - Ceps1*betaSq;
251 std::complex<double> C2 = C1/(1.+Ceps1) * std::log(2.*betaSq*MeeV/(E*C1));
252
253 double x;
255 x = m_S2/betaSq * E * ( -2.*imag(C2)/(m_Wp2*M_PI) + (1.-std::exp(x))/(E*E) );
256
257 return x;
258}
#define 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
gamvecvector of gamma values for which to tabulate
Earrayvector of energy transfer values
fnArraydouble differential distribution
dndxThis array effectively contains the inverse of the mean free path as function of binned energy variable

Definition at line 262 of file TRT_PAI_effectiveGas.cxx.

266{
267 // For a given gamma factor, return dn/dx
268 // (by integrating dn^2/dEdx with respect to E)
269 // and tabulated values of dn^2/dEdx(?)....
270
271 double Rener = 0.05;
272
273 int NLvls = m_lnELvls.size();
274 int nGamVals = gamvec.size();
275
276 double lnEi = m_lnEmax;
277 double lnEo;
278 double lnEs = -99999999.;
279
280 fnArray.resize(nGamVals);
281 int nEVals = 0;
282
283 // Initialize
284 for ( int ig=0; ig<nGamVals; ++ig ) {
285 fnArray[ig].push_back(0.);
286 }
287
288 // Calculate integral from above and store the cumulative values in fnArray
289
290 for ( int ie=NLvls-1; ie>=0; --ie ) {
291 lnEo = lnEi;
292 lnEi = m_lnELvls[ie];
293 for ( int ig=0; ig<nGamVals; ++ig ) {
294 double ds =
295 XGInt( &TRT_PAI_effectiveGas::dndedx, lnEi, lnEo, m_eps, gamvec[ig]);
296 fnArray[ig][nEVals] += ds;
297 }
298 if ( std::abs(lnEs-lnEi) > Rener || ie==0 ) {
299 lnEs = lnEi;
300 lnEarray.push_back( lnEs );
301 if ( ie>0 ) {
302 for ( int ig=0; ig<nGamVals; ++ig ) {
303 fnArray[ig].push_back( fnArray[ig][nEVals] );
304 }
305 }
306 nEVals++;
307 }
308 }
309
310 // Copy the total integral into auxillary vector
311 dndx.resize(nGamVals);
312 for ( int ig = 0; ig < nGamVals; ++ig ) {
313 dndx[ig] = fnArray[ig][nEVals-1];
314 }
315
316 return;
317}
double dndedx(double lgE, double gamma)

◆ 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.

40{
42 // If user did not set an explicit level, set a default
43 if (m_lvl == MSG::NIL) {
44 m_lvl = m_imsg ?
45 static_cast<MSG::Level>( m_imsg.load()->outputLevel(m_nm) ) :
46 MSG::INFO;
47 }
48}
std::string m_nm
Message source name.
std::atomic< IMessageSvc * > m_imsg
MessageSvc pointer.
std::atomic< MSG::Level > m_lvl
Current logging level.
IMessageSvc * getMessageSvc(bool quiet=false)

◆ 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.

164{
165 MsgStream* ms = m_msg_tls.get();
166 if (!ms) {
167 if (!m_initialized.test_and_set()) initMessaging();
168 ms = new MsgStream(m_imsg,m_nm);
169 m_msg_tls.reset( ms );
170 }
171
172 ms->setLevel (m_lvl);
173 return *ms;
174}
boost::thread_specific_ptr< MsgStream > m_msg_tls
MsgStream instance (a std::cout like with print-out levels)
void initMessaging() const
Initialize our message level and MessageSvc.

◆ 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; }
MsgStream & msg() const
The standard message stream.

◆ msgLvl()

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

Test the output level.

Parameters
lvlThe message level to test against
Returns
boolean Indicating if messages at given level will be printed
Return values
trueMessages at level "lvl" will be printed

Definition at line 151 of file AthMessaging.h.

152{
153 if (m_lvl <= lvl) {
154 msg() << lvl;
155 return true;
156 } else {
157 return false;
158 }
159}

◆ 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.

29{
30 m_lvl = lvl;
31}

◆ XFReal()

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

Definition at line 173 of file TRT_PAI_effectiveGas.cxx.

173 {
174
175 double fp = 0.;
176 if ( lnE+lnD > m_lnEmin ) {
178 fp = std::exp(std::max(fp, -99.0));
179 }
180
181 double fm = 0.;
182 if (lnE-lnD > m_lnEmin ) {
184 fm = std::exp(std::max(fm, -99.0));
185 }
186
187 double x = std::exp(lnD);
188 return x/(x*x-1.)*(fp-fm);
189}

◆ 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 193 of file TRT_PAI_effectiveGas.cxx.

197 {
198 //
199 // Pavel Nevskis "famous integration procedure"
200 //
201 // Applies the Gauss-Legendre Quadrature method with n=4
202
203 // These constants belong to the Gauss-Legendre quadrature method:
204 const int M = 4; // 4 points
205 const double U[M]={-.8611363,-.3399810, .3399810 ,.8611363}; // abscissas
206 const double W[M]={ .3478548, .6521452, .6521452, .3478548}; // weights
207
208 int N=10;
209 double OTB=0,Y,D,result;
210
211 double Dt = 0.5*(lnHi-lnLo);
212
213 do {
214 Y = OTB;
215 OTB = 0.;
216 D = Dt/N;
217 for (int i=1; i<=N; i++) {
218 for (int k=0; k<M; k++) {
219 double arg = lnLo + D*(2*i-1+U[k]);
220 OTB += W[k] * (this->*pt2Func)(arg, extraParameter );
221 }
222 }
223 OTB *= D;
224 result = OTB;
225 N *= 2;
226 if ( N>100000 ) {
227 ATH_MSG_WARNING( "effectiveGas::XGInt: integrate Divergence! " << std::abs(OTB-Y) << " " << std::abs(eps*OTB) );
228 break;
229 }
230 } while ( (std::abs(OTB-Y) > std::abs(eps*OTB)) && (eps>0) );
231
232 return result;
233}
#define ATH_MSG_WARNING(x)

◆ XSigma()

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

Definition at line 161 of file TRT_PAI_effectiveGas.cxx.

161 {
162 double xsig = dummy;
163 xsig = TRT_PAI_utils::Interpolate(std::max(lnE,m_lnEmin),
164 m_lnELvls,
165 m_lnFosc);
166 xsig = std::max(xsig,-99.);
167 xsig = std::exp( xsig + lnE );
168 return xsig;
169}

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.

◆ 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.

◆ m_lnEmin

const double TRT_PAI_effectiveGas::m_lnEmin
private

Definition at line 60 of file TRT_PAI_effectiveGas.h.

◆ 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.

◆ 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.

◆ m_S2

double TRT_PAI_effectiveGas::m_S2
private

Definition at line 65 of file TRT_PAI_effectiveGas.h.

◆ m_Wp2

double TRT_PAI_effectiveGas::m_Wp2
private

Definition at line 64 of file TRT_PAI_effectiveGas.h.


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