Jpp 19.3.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
JK40DefaultSimulatorInterface.hh
Go to the documentation of this file.
1#ifndef __JDETECTOR__JK40DEFAULTSIMULATORINTERFACE__
2#define __JDETECTOR__JK40DEFAULTSIMULATORINTERFACE__
3
4#include <vector>
5#include <set>
6
7#include "TRandom3.h"
8
13
14#include "JMath/JMathToolkit.hh"
15#include "JLang/JComparator.hh"
16#include "JLang/JException.hh"
17#include "JLang/JCC.hh"
18
19
20/**
21 * \author mdejong
22 */
23
24namespace JDETECTOR {}
25namespace JPP { using namespace JDETECTOR; }
26
27namespace JDETECTOR {
28
30
31
32 /**
33 * Default K40 simulator interface.
34 *
35 * This class provides for a default implementation of the JK40Simulator interface
36 * which is based on a set of virtual methods.
37 * These methods constitute a user interface to the K40 background simulation.
38 */
40 public JK40Simulator
41 {
42 /**
43 * PMT pair.
44 */
45 struct pair_type {
46 size_t pmt1; //!< first PMT
47 size_t pmt2; //!< second PMT
48 double P; //!< probability of coincidence
49 };
50
51
52 /**
53 * Auxiliary data structure for the probabilities of genuine coincidences due to K40 decays as a function of the PMT pair.\n
54 * The first element in this list should correspond to an arbitray PMT pair with zero probability before using the function operator.
55 */
56 mutable struct :
58 {
59 /**
60 * Get random PMT pair according probabilities of genuine coincidences due to K40 decays.
61 *
62 * \param random random value <0,1]
63 * \return PMT pair
64 */
65 const pair_type& operator()(const double random) const
66 {
67 double P = this->rbegin()->P;
68
69 if (P <= 0.0) {
70 THROW(JNumericalPrecision, "Probability " << P);
71 }
72
73 P *= random;
74
75 int imin = 0;
76 int imax = this->size() - 1;
77
78 for (int i = (imax + imin) / 2; imax - imin != 1; i = (imax + imin) / 2) {
79 if ((*this)[i].P < P)
80 imin = i;
81 else
82 imax = i;
83 }
84
85 return (*this)[imax];
86 }
88
89
90 protected:
91 /**
92 * Default constructor.
93 */
96
97
98 public:
99 /**
100 * Get singles rate as a function of PMT.
101 *
102 * \param pmt PMT identifier
103 * \return rate [Hz]
104 */
105 virtual double getSinglesRate(const JPMTIdentifier& pmt) const = 0;
106
107
108 /**
109 * Get multiples rate as a function of optical module.
110 *
111 * \param module optical module identifier
112 * \param M multiplicity (M >= 2)
113 * \return rate [Hz]
114 */
115 virtual double getMultiplesRate(const JModuleIdentifier& module, const int M) const = 0;
116
117
118 /**
119 * Get probability of coincidence.
120 *
121 * \param ct cosine space angle between PMT axes
122 * \return probability
123 */
124 virtual double getProbability(const double ct) const = 0;
125
126
127 /**
128 * Generate hits.
129 *
130 * \param module module
131 * \param period time window [ns]
132 * \param output background data
133 */
134 virtual void generateHits(const JModule& module,
135 const JTimeRange& period,
136 JModuleData& output) const
137 {
138 using namespace std;
139 using namespace JPP;
140
141
142 // resize internal buffers
143
144 const size_t N = module.size();
145 const size_t M = (N * (N - 1)) / 2;
146
147 rateL1_Hz .resize(N);
148 probabilityL1.resize(M + 1);
149
150
151 // generate singles
152
153 for (size_t pmt = 0; pmt != N; ++pmt) {
154
155 const double rateL0_Hz = getSinglesRate(JPMTIdentifier(module.getID(), pmt));
156
157 if (rateL0_Hz > 0.0) {
158
159 const double t_ns = 1.0e9 / rateL0_Hz; // [ns]
160
161 for (double t1 = period.getLowerLimit() + gRandom->Exp(t_ns); t1 < period.getUpperLimit(); t1 += gRandom->Exp(t_ns)) {
162 output[pmt].push_back(JPMTSignal(t1, 1));
163 }
164 }
165 }
166
167
168 // generate coincidences
169
170 double totalRateL1_Hz = 0.0;
171
172 for (size_t i = 0; i != N; ++i) {
173 totalRateL1_Hz += rateL1_Hz[i] = getMultiplesRate(module.getID(), i + 2);
174 }
175
176 if (totalRateL1_Hz > 0.0) {
177
178 const double t_ns = 1.0e9 / totalRateL1_Hz; // [ns]
179
180 double t1 = period.getLowerLimit() + gRandom->Exp(t_ns);
181
182 if (t1 < period.getUpperLimit()) {
183
184 // configure pair-wise propabilities
185
186 probabilityL1[0] = { N, N, 0.0};
187
188 size_t i = 0;
189 double P = 0.0;
190
191 for (size_t pmt1 = 0; pmt1 != N; ++pmt1) {
192 for (size_t pmt2 = 0; pmt2 != pmt1; ++pmt2) {
193
194 const double ct = getDot(module[pmt1].getDirection(), module[pmt2].getDirection());
195 const double p = getProbability(ct);
196
197 i += 1;
198 P += p;
199
200 probabilityL1[i] = { pmt1, pmt2, P};
201 }
202 }
203
204 for ( ; t1 < period.getUpperLimit(); t1 += gRandom->Exp(t_ns)) {
205
206 try {
207
208 // generate two-fold coincidence
209
210 const pair_type& pair = probabilityL1(gRandom->Rndm());
211
212 output[pair.pmt1].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
213 output[pair.pmt2].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
214
215 // generate larger than two-fold coincidences, if any
216
217 size_t M = 0;
218
219 for (double R = totalRateL1_Hz * gRandom->Rndm(); M != N && (R -= rateL1_Hz[M]) > 0.0; ++M) {}
220
221 if (M != 0) {
222
223 set<size_t> buffer = { pair.pmt1, pair.pmt2 }; // hit PMTs
224
225 for ( ; M != 0; --M) {
226
227 vector<double> probability1D(N, 1.0);
228
229 double P = 0.0;
230
231 for (size_t i = 0; i != N; ++i) {
232
233 if (buffer.count(i) == 0) {
234
235 for (set<size_t>::const_iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt) {
236
237 const double ct = getDot(module[i].getDirection(), module[*pmt].getDirection());
238
239 probability1D[i] *= getProbability(ct);
240 }
241
242 P += probability1D[i];
243
244 } else {
245
246 probability1D[i] = 0.0;
247 }
248 }
249
250 if (P > 0.0) {
251
252 size_t pmt = 0;
253
254 for (P *= gRandom->Rndm(); pmt != N && (P -= probability1D[pmt]) > 0.0; ++pmt) {}
255
256 if (pmt != N) {
257
258 output[pmt].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
259
260 buffer.insert(pmt);
261 }
262
263 } else {
264
265 break;
266 }
267 }
268 }
269 }
270 catch (const JNumericalPrecision&) {}
271 }
272 }
273 }
274 }
275
276
277 /**
278 * Get intrinsic time smearing of K40 coincidences.
279 *
280 * \return sigma [ns]
281 */
282 static double getSigma()
283 {
284 return get_sigma();
285 }
286
287
288 /**
289 * Set intrinsic time smearing of K40 coincidences.
290 *
291 * \param sigma sigma [ns]
292 */
293 static void setSigma(const double sigma)
294 {
295 get_sigma() = sigma;
296 }
297
298 private:
299 /**
300 * Get intrinsic time smearing of K40 coincidences.
301 *
302 * \return sigma [ns]
303 */
304 static double& get_sigma()
305 {
306 static double sigma = 0.5;
307
308 return sigma;
309 }
310
311
312 /**
313 * Multiples rate as a function of the multiplicity.
314 * The index i corresponds to multiplicity M = i + 2.
315 */
317 };
318}
319
320#endif
Compiler version dependent expressions, macros, etc.
Exceptions.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Auxiliary methods for geometrical methods.
virtual void generateHits(const JModule &module, const JTimeRange &period, JModuleData &output) const
Generate hits.
static double & get_sigma()
Get intrinsic time smearing of K40 coincidences.
static double getSigma()
Get intrinsic time smearing of K40 coincidences.
std::vector< double > rateL1_Hz
Multiples rate as a function of the multiplicity.
static void setSigma(const double sigma)
Set intrinsic time smearing of K40 coincidences.
virtual double getSinglesRate(const JPMTIdentifier &pmt) const =0
Get singles rate as a function of PMT.
virtual double getMultiplesRate(const JModuleIdentifier &module, const int M) const =0
Get multiples rate as a function of optical module.
virtual double getProbability(const double ct) const =0
Get probability of coincidence.
Interface for simulation of K40 background.
Data structure for PMT data corresponding to a detector module.
Data structure for a composite optical module.
Definition JModule.hh:75
Exception for numerical precision error.
Auxiliary class for object identification.
Definition JObjectID.hh:25
int getID() const
Get identifier.
Definition JObjectID.hh:50
T getLowerLimit() const
Get lower limit.
Definition JRange.hh:202
T getUpperLimit() const
Get upper limit.
Definition JRange.hh:213
file Auxiliary data structures and methods for detector calibration.
Definition JAnchor.hh:12
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Data structure for PMT analogue signal.