Jpp  master_rocky
the software that should make you happy
JK40Rates.hh
Go to the documentation of this file.
1 #ifndef __JPHYSICS__JK40RATES___
2 #define __JPHYSICS__JK40RATES___
3 
4 #include <vector>
5 
7 
8 
9 /**
10  * \author mdejong
11  */
12 
13 namespace JPHYSICS {}
14 namespace JPP { using namespace JPHYSICS; }
15 
16 namespace JPHYSICS {
17 
18  /**
19  * Type definition of singles rate [Hz].
20  */
21  typedef double JRateL0_t;
22 
23  /**
24  * Type definition of count rate as a function of multiplicty [Hz]
25  * The multiples rate start counting at two-fold coincidences.
26  */
28 
29 
30  /**
31  * Type definition of multiplicity.
32  */
33  typedef size_t multiplicity_type;
34 
35 
36  /**
37  * Auxiliary class for K40 rates.
38  * The singles rate refers to the counting rate of a PMT and
39  * the multiples rate to the rate of genuine coincidences due to K40 decays.
40  */
41  struct JK40Rates {
42  /**
43  * Default constructor.
44  */
46  rateL0(0.0),
47  rateL1()
48  {}
49 
50 
51  /**
52  * Constructor.
53  *
54  * The multiples rates start counting at two-fold coincidences.
55  *
56  * \param rateL0_Hz singles rate [Hz]
57  * \param rateL1_Hz multiples rates [Hz]
58  */
59  JK40Rates(const JRateL0_t rateL0_Hz,
60  const JRateL1_t& rateL1_Hz = JRateL1_t()) :
61  rateL0(rateL0_Hz),
62  rateL1(rateL1_Hz)
63  {}
64 
65 
66  /**
67  * Get singles rate.
68  *
69  * \return rate [Hz]
70  */
71  double getSinglesRate() const
72  {
73  return rateL0;
74  }
75 
76 
77  /**
78  * Get multiples rate.
79  *
80  * \return rate [Hz]
81  */
83  {
84  return rateL1;
85  }
86 
87 
88  /**
89  * Get multiples rate at given multiplicity.
90  *
91  * \param M multiplicity (M >= JK40Rates::LOWER_L1_MULTIPLICITY)
92  * \return rate [Hz]
93  */
94  double getMultiplesRate(const multiplicity_type M) const
95  {
97  return rateL1[M - LOWER_L1_MULTIPLICITY];
98  else
99  return 0.0;
100  }
101 
102 
103  /**
104  * Get lower multiplicty.
105  *
106  * \return lower multiplicity
107  */
109  {
110  return LOWER_L1_MULTIPLICITY;
111  }
112 
113 
114  /**
115  * Get upper multiplicty.
116  *
117  * \return upper multiplicity
118  */
120  {
121  return rateL1.size() + 1;
122  }
123 
124 
125  /**
126  * Correct rates for global efficiency,
127  *
128  * \param QE global efficiency
129  */
130  void correct(const double QE)
131  {
132  if (QE > 0.0) {
133 
134  rateL0 /= QE;
135 
136  JRateL1_t buffer = rateL1;
137 
139 
140  // determine contribution from higher multiplicities
141 
142  double R = 0.0;
143 
144  for (multiplicity_type i = M + 1; i <= getUpperL1Multiplicity(); ++i) {
145  R += buffer[i - LOWER_L1_MULTIPLICITY] * JMATH::binomial(i, M) * pow(QE, M) * pow(1.0 - QE, i - M);
146  }
147 
148  if (getMultiplesRate(M) > R)
149  buffer[M - LOWER_L1_MULTIPLICITY] = (getMultiplesRate(M) - R) / pow(QE, M);
150  else
151  buffer[M - LOWER_L1_MULTIPLICITY] = 0.0;
152  }
153 
154  rateL1 = buffer;
155 
156  } else {
157 
158  rateL0 = 0.0;
159 
160  for (JRateL1_t::iterator i = rateL1.begin(); i != rateL1.end(); ++i) {
161  *i = 0.0;
162  }
163  }
164  }
165 
166 
167  /**
168  * Read K40 rates from input.
169  *
170  * \param in input stream
171  * \param object K40 rates
172  * \return input stream
173  */
174  friend inline std::istream& operator>>(std::istream& in, JK40Rates& object)
175  {
176  const double rateL0 = object.rateL0;
177 
178  if (in >> object.rateL0) {
179 
180  object.rateL1.clear();
181 
182  for (double x; in >> x; ) {
183  object.rateL1.push_back(x);
184  }
185 
186  } else {
187 
188  object.rateL0 = rateL0;
189  }
190 
191  return in;
192  }
193 
194 
195  /**
196  * Write K40 rates to output.
197  *
198  * \param out output stream
199  * \param object K40 rates
200  * \return output stream
201  */
202  friend inline std::ostream& operator<<(std::ostream& out, const JK40Rates& object)
203  {
204  out << object.rateL0;
205 
206  for (JRateL1_t::const_iterator i = object.rateL1.begin(); i != object.rateL1.end(); ++i) {
207  out << ' ' << *i;
208  }
209 
210  return out;
211  }
212 
213 
214  /**
215  * Lower L1 multiplicity
216  */
218 
219 
220  protected:
221  JRateL0_t rateL0; //!< singles rate [Hz]
222  JRateL1_t rateL1; //!< multiples rates [Hz]
223  };
224 }
225 
226 #endif
Auxiliary methods for mathematics.
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:97
double binomial(const size_t n, const size_t k)
Binomial function.
Auxiliary methods for light properties of deep-sea water.
size_t multiplicity_type
Type definition of multiplicity.
Definition: JK40Rates.hh:33
std::vector< double > JRateL1_t
Type definition of count rate as a function of multiplicty [Hz] The multiples rate start counting at ...
Definition: JK40Rates.hh:27
double JRateL0_t
Type definition of singles rate [Hz].
Definition: JK40Rates.hh:21
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Auxiliary class for K40 rates.
Definition: JK40Rates.hh:41
JRateL1_t rateL1
multiples rates [Hz]
Definition: JK40Rates.hh:222
double getSinglesRate() const
Get singles rate.
Definition: JK40Rates.hh:71
friend std::istream & operator>>(std::istream &in, JK40Rates &object)
Read K40 rates from input.
Definition: JK40Rates.hh:174
const JRateL1_t & getMultiplesRates() const
Get multiples rate.
Definition: JK40Rates.hh:82
JK40Rates(const JRateL0_t rateL0_Hz, const JRateL1_t &rateL1_Hz=JRateL1_t())
Constructor.
Definition: JK40Rates.hh:59
JK40Rates()
Default constructor.
Definition: JK40Rates.hh:45
double getMultiplesRate(const multiplicity_type M) const
Get multiples rate at given multiplicity.
Definition: JK40Rates.hh:94
friend std::ostream & operator<<(std::ostream &out, const JK40Rates &object)
Write K40 rates to output.
Definition: JK40Rates.hh:202
multiplicity_type getLowerL1Multiplicity() const
Get lower multiplicty.
Definition: JK40Rates.hh:108
multiplicity_type getUpperL1Multiplicity() const
Get upper multiplicty.
Definition: JK40Rates.hh:119
static const multiplicity_type LOWER_L1_MULTIPLICITY
Lower L1 multiplicity.
Definition: JK40Rates.hh:217
JRateL0_t rateL0
singles rate [Hz]
Definition: JK40Rates.hh:221
void correct(const double QE)
Correct rates for global efficiency,.
Definition: JK40Rates.hh:130