Jpp test-rotations-old
the software that should make you happy
Loading...
Searching...
No Matches
JRadiationFunction.hh
Go to the documentation of this file.
1#ifndef __JPHYSICS__JRADIATIONFUNCTION__
2#define __JPHYSICS__JRADIATIONFUNCTION__
3
7#include "JTools/JGrid.hh"
8
10
11
12/**
13 * \author mdejong
14 */
15
16namespace JPHYSICS {}
17namespace JPP { using namespace JPHYSICS; }
18
19namespace JPHYSICS {
20
21 /**
22 * Fast implementation of class JRadiation.
23 *
24 * In this, the methods
25 * - JRadiation::TotalCrossSectionEErad,
26 * - JRadiation::TotalCrossSectionGNrad,
27 * - JRadiation::CalculateACoeff and
28 * - JRadiation::IntegralofG
29 *
30 * are reimplemented using lookup tables.
31 */
33 public JRadiation
34 {
38
39 public:
40 /**
41 * Constructor.
42 *
43 * \param radiation JRadiation object
44 * \param number_of_bins number of bins
45 * \param Emin minimal muon energy [GeV]
46 * \param Emax maximal muon energy [GeV]
47 */
49 const unsigned int number_of_bins,
50 const double Emin,
51 const double Emax) :
52 JRadiation(radiation)
53 {
54 using namespace JTOOLS;
55
56 const double xmin = log(Emin);
57 const double xmax = log(Emax);
58
59
60 integral.configure(make_grid(number_of_bins, xmin, xmax));
61
62 for (JFunction2D_t::iterator i = integral.begin(); i != integral.end(); ++i) {
63
64 const double x = i->getX();
65 const double E = exp(x);
66
67 const double ymin = log(EminEErad/E);
68 const double ymax = 0.0;
69
70 if (ymax > ymin) {
71
72 i->getY().configure(make_grid(number_of_bins, ymin, ymax));
73
74 for (JFunction1D_t::iterator j = i->getY().begin(); j != i->getY().end(); ++j) {
75
76 const double y = j->getX();
77 const double eps = exp(y) * E;
78 const double z = JRadiation::IntegralofG(E, eps);
79
80 j->getY() = z;
81 }
82 }
83 }
84
86 integral.setExceptionHandler(new JFunction1D_t::JDefaultResult(0.0));
87
88
89 sigmaEE.configure(make_grid(number_of_bins, xmin, xmax));
90
91 for (JFunction1D_t::iterator i = sigmaEE.begin(); i != sigmaEE.end(); ++i) {
92
93 const double E = exp(i->getX());
94 const double y = JRadiation::TotalCrossSectionEErad(E);
95
96 i->getY() = y;
97 }
98
99 sigmaEE.compile();
100
101
102 sigmaGN.configure(make_grid(number_of_bins, xmin, xmax));
103
104 for (JFunction1D_t::iterator i = sigmaGN.begin(); i != sigmaGN.end(); ++i) {
105
106 const double E = exp(i->getX());
107 const double y = JRadiation::TotalCrossSectionGNrad(E);
108
109 i->getY() = y;
110 }
111
112 sigmaGN.compile();
113
114
115 Acoeff.configure(make_grid(number_of_bins, xmin, xmax));
116
117 for (JFunction1D_t::iterator i = Acoeff.begin(); i != Acoeff.end(); ++i) {
118
119 const double E = exp(i->getX());
120 const double y = JRadiation::CalculateACoeff(E);
121
122 i->getY() = y;
123 }
124
125 Acoeff.compile();
126 }
127
128
129 /**
130 * Pair production cross section.
131 *
132 * \param E muon energy [GeV]
133 * \return cross section [m^2/g]
134 */
135 virtual double TotalCrossSectionEErad(const double E) const override
136 {
137 const double x = log(E);
138
139 if (x >= sigmaEE. begin()->getX() &&
140 x <= sigmaEE.rbegin()->getX()) {
141
142 try {
143 return sigmaEE(x);
144 }
145 catch(std::exception& error) {}
146 }
147
149 }
150
151
152 /**
153 * Photo-nuclear cross section.
154 *
155 * \param E muon energy [GeV]
156 * \return cross section [m^2/g]
157 */
158 virtual double TotalCrossSectionGNrad(const double E) const override
159 {
160 const double x = log(E);
161
162 if (x >= sigmaGN. begin()->getX() &&
163 x <= sigmaGN.rbegin()->getX()) {
164
165 try {
166 return sigmaGN(x);
167 }
168 catch(std::exception& error) {}
169 }
170
172 }
173
174
175 /**
176 * Ionization a parameter.
177 *
178 * \param E muon energy [GeV]
179 * \return ionization coefficient [GeV*m^2/g]
180 */
181 virtual double CalculateACoeff(const double E) const override
182 {
183 const double x = log(E);
184
185 if (x >= Acoeff. begin()->getX() &&
186 x <= Acoeff.rbegin()->getX()) {
187
188 try {
189 return Acoeff(x);
190 }
191 catch(std::exception& error) {}
192 }
193
195 }
196
197 protected:
198 virtual double IntegralofG(const double E,
199 const double eps) const override
200 {
201 try {
202
203 const double x = log(E);
204 const double y = log(eps/E);
205 const double z = integral(x, y);
206
207 return z;
208 }
209 catch(std::exception& error) {}
210
211 return JRadiation::IntegralofG(E, eps);
212 }
213
218 };
219}
220
221#endif
Various implementations of functional maps.
Muon radiative cross sections.
Fast implementation of class JRadiation.
virtual double TotalCrossSectionGNrad(const double E) const override
Photo-nuclear cross section.
JTOOLS::JGridPolint1Function1D_t JFunction1D_t
virtual double CalculateACoeff(const double E) const override
Ionization a parameter.
virtual double TotalCrossSectionEErad(const double E) const override
Pair production cross section.
JRadiationFunction(const JRadiation &radiation, const unsigned int number_of_bins, const double Emin, const double Emax)
Constructor.
virtual double IntegralofG(const double E, const double eps) const override
JTOOLS::JMultiFunction< JFunction1D_t, JTOOLS::JMapList< JTOOLS::JPolint1FunctionalGridMap > > JFunction2D_t
Auxiliary class for the calculation of the muon radiative cross sections.
Definition JRadiation.hh:36
virtual double IntegralofG(const double E, const double eps) const
virtual double TotalCrossSectionEErad(const double E) const
Pair production cross section.
virtual double CalculateACoeff(double E) const
Ionization a parameter.
virtual double TotalCrossSectionGNrad(const double E) const
Photo-nuclear cross section.
const double EminEErad
Multidimensional interpolation method.
void setExceptionHandler(const typename function_type::supervisor_type &supervisor)
Set the supervisor for handling of exceptions.
void compile()
Compilation.
Auxiliary methods for light properties of deep-sea water.
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Auxiliary classes and methods for multi-dimensional interpolations and histograms.
Type definition of a 1st degree polynomial interpolation based on a JGridCollection with result type ...
Map list.
Definition JMapList.hh:25