Jpp  19.1.0
the software that should make you happy
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 
12 #include "JDetector/JTimeRange.hh"
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 
24 namespace JDETECTOR {}
25 namespace JPP { using namespace JDETECTOR; }
26 
27 namespace 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  */
95  {}
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.
Definition: JException.hh:712
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.
Definition: JException.hh:270
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
const double sigma[]
Definition: JQuadrature.cc:74
JDirection3D getDirection(const Vec &dir)
Get direction.
double getDot(const JNeutrinoDirection &first, const JNeutrinoDirection &second)
Dot product.
Definition: JAstronomy.hh:676
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).
Definition: JSTDTypes.hh:14
Data structure for PMT analogue signal.