Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JK40DefaultSimulatorInterface.hh
Go to the documentation of this file.
1 #ifndef __JDETECTOR__JK40DEFAULTSIMULATORINTERFACE__
2 #define __JDETECTOR__JK40DEFAULTSIMULATORINTERFACE__
3 
4 #include "TRandom3.h"
5 
10 
11 #include "JMath/JMathToolkit.hh"
12 #include "JLang/JException.hh"
13 #include "JLang/JCC.hh"
14 
15 
16 /**
17  * \author mdejong
18  */
19 
20 namespace JDETECTOR {}
21 namespace JPP { using namespace JDETECTOR; }
22 
23 namespace JDETECTOR {
24 
26 
27 
28  /**
29  * Default K40 simulator interface.
30  *
31  * This class provides for a default implementation of the JK40Simulator interface
32  * which is based on a set of virtual methods.
33  * These methods constitute a user interface to the K40 background simulation.
34  */
36  public JK40Simulator
37  {
38  protected:
39  /**
40  * 1D-linear buffer.
41  */
42  struct JBuffer1D_t :
43  public std::vector<double>
44  {};
45 
46 
47  /**
48  * 2D-square buffer.
49  */
50  struct JBuffer2D_t :
51  public std::vector<JBuffer1D_t>
52  {
53  /**
54  * Resize.
55  *
56  * \param size size
57  */
58  void resize(size_t size)
59  {
61 
62  for (iterator i = begin(); i != end(); ++i) {
63  i->resize(size);
64  }
65  }
66  };
67 
68 
69  /**
70  * Default constructor.
71  */
73  {}
74 
75 
76  public:
77  /**
78  * Get singles rate as a function of PMT.
79  *
80  * \param pmt PMT identifier
81  * \return rate [Hz]
82  */
83  virtual double getSinglesRate(const JPMTIdentifier& pmt) const = 0;
84 
85 
86  /**
87  * Get multiples rate as a function of optical module.
88  *
89  * \param module optical module identifier
90  * \param M multiplicity (M >= 2)
91  * \return rate [Hz]
92  */
93  virtual double getMultiplesRate(const JModuleIdentifier& module, const int M) const = 0;
94 
95 
96  /**
97  * Get probability of coincidence.
98  *
99  * \param ct cosine space angle between PMT axes
100  * \return probability
101  */
102  virtual double getProbability(const double ct) const = 0;
103 
104 
105  /**
106  * Generate hits.
107  *
108  * \param module module
109  * \param period time window [ns]
110  * \param output background data
111  */
112  virtual void generateHits(const JModule& module,
113  const JTimeRange& period,
114  JModuleData& output) const
115  {
116 
117  // resize internal buffers
118 
119  const int N = module.size();
120 
122  probability1D.resize(N);
123  probabilityND.resize(N);
124 
125  rateL1_Hz.resize(N);
126 
127 
128  // generate singles
129 
130  for (int pmt = 0; pmt != N; ++pmt) {
131 
132  const double rateL0_Hz = getSinglesRate(JPMTIdentifier(module.getID(), pmt));
133 
134  if (rateL0_Hz > 0.0) {
135 
136  const double t_ns = 1.0e9 / rateL0_Hz; // [ns]
137 
138  for (double t1 = period.getLowerLimit() + gRandom->Exp(t_ns); t1 < period.getUpperLimit(); t1 += gRandom->Exp(t_ns)) {
139  output[pmt].push_back(JPMTSignal(t1, 1));
140  }
141  }
142  }
143 
144 
145  // generate coincidences
146 
147  double totalRateL1_Hz = 0.0;
148 
149  for (int i = 0; i != N; ++i) {
150  totalRateL1_Hz += rateL1_Hz[i] = getMultiplesRate(module.getID(), i + 2);
151  }
152 
153  if (totalRateL1_Hz > 0.0) {
154 
155  const double t_ns = 1.0e9 / totalRateL1_Hz; // [ns]
156 
157  double t1 = period.getLowerLimit() + gRandom->Exp(t_ns);
158 
159  if (t1 < period.getUpperLimit()) {
160 
161 
162  // configure correlation matrix
163 
164  for (int i = 0; i != N; ++i) {
165 
166  probability2D[i][i] = 0.0;
167 
168  for (int j = i + 1; j != N; ++j) {
169 
170  const double ct = JMATH::getDot(module[i], module[j]);
171  const double p = getProbability(ct);
172 
173  probability2D[i][j] = p;
174  probability2D[j][i] = p;
175  }
176  }
177 
178 
179  // determine probability of a coincidence as a function of the PMT number
180 
181  for (int i = 0; i != N; ++i) {
182 
183  probability1D[i] = 0.0;
184 
185  for (int j = 0; j != N; ++j) {
186  probability1D[i] += probability2D[i][j];
187  }
188  }
189 
190 
191  for ( ; t1 < period.getUpperLimit(); t1 += gRandom->Exp(t_ns)) {
192 
193 
194  // generate two-fold coincidence
195 
196  const unsigned int pmt1 = getRandomIndex(probability1D, gRandom->Rndm());
197  const unsigned int pmt2 = getRandomIndex(probability2D[pmt1], gRandom->Rndm());
198 
199  output[pmt1].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
200  output[pmt2].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
201 
202 
203  try {
204 
205  // generate larger than two-fold coincidences, if any
206 
207  unsigned int M = getRandomIndex(rateL1_Hz, gRandom->Rndm());
208 
209  if (M != 0) {
210 
212 
213  for (unsigned int pmtN = pmt2; M != 0; --M) {
214 
215  for (int i = 0; i != N; ++i) {
216  probabilityND[i] *= probability2D[pmtN][i];
217  }
218 
219  probabilityND[pmtN] = 0.0;
220 
221  pmtN = getRandomIndex(probabilityND, gRandom->Rndm());
222 
223  output[pmtN].insert(JPMTSignal(gRandom->Gaus(t1, getSigma()), 1));
224  }
225  }
226  }
227  catch (const JNumericalPrecision&) {}
228  }
229  }
230  }
231  }
232 
233 
234  /**
235  * Get index based on random value.
236  *
237  * It is assumed that the values of the input buffer monotonously decrease or increase.
238  * This method throws an exception when the summed values in the input buffer is zero.
239  *
240  * \param buffer input values
241  * \param random random value <0,1]
242  * \return index
243  */
244  static inline unsigned int getRandomIndex(const JBuffer1D_t& buffer, const double random)
245  {
246  double x = 0.0;
247 
248  for (JBuffer1D_t::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
249  x += *i;
250  }
251 
252  if (x > 0.0) {
253 
254  x *= random;
255 
256  unsigned int index = 0;
257 
258  for (JBuffer1D_t::const_iterator i = buffer.begin(); i != buffer.end() && (x -= *i) > 0.0; ++i, ++index) {}
259 
260  if (index == buffer.size()) {
261  --index;
262  }
263 
264  return index;
265 
266  } else {
267 
268  THROW(JNumericalPrecision, "getRandomIndex(): zero or negative probability.");
269  }
270  }
271 
272 
273  /**
274  * Get intrinsic time smearing of K40 coincidences.
275  *
276  * \return sigma [ns]
277  */
278  static double getSigma()
279  {
280  return 0.5;
281  }
282 
283  private:
284  /**
285  * This correlation matrix is a two-dimensional array in which element [i][j]
286  * corresponds to the probability of a genuine coincidence due to K40 decays,
287  * where i and j refer to the indices of the PMTs in the optical module.
288  */
290 
291  /**
292  * This probability vector is a one-dimensional array in which element [i]
293  * corresponds to the probability of a genuine coincidence due to K40 decays,
294  * where i refers to the index of the PMT in the optical module.
295  */
297 
298  /**
299  * This probability vector is a one-dimensional array in which element [i]
300  * corresponds to the probability of an additional hit,
301  * where i refers to the index of the PMT in the optical module.
302  */
304 
305  /**
306  * Multiples rate as a function of the multiplicity.
307  * The index i corresponds to multiplicity M = i + 2.
308  */
310  };
311 }
312 
313 #endif
static double getSigma()
Get intrinsic time smearing of K40 coincidences.
Exceptions.
Data structure for PMT analogue signal.
Auxiliary methods for geometrical methods.
Interface for simulation of K40 background.
Data structure for a composite optical module.
Definition: JModule.hh:47
Data structure for PMT data corresponding to a detector module.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:633
JBuffer1D_t probabilityND
This probability vector is a one-dimensional array in which element [i] corresponds to the probabilit...
JRange< double > JTimeRange
Type definition for time range.
static unsigned int getRandomIndex(const JBuffer1D_t &buffer, const double random)
Get index based on random value.
Compiler version dependent expressions, macros, etc.
virtual double getSinglesRate(const JPMTIdentifier &pmt) const =0
Get singles rate as a function of PMT.
int getID() const
Get identifier.
Definition: JObjectID.hh:54
JBuffer2D_t probability2D
This correlation matrix is a two-dimensional array in which element [i][j] corresponds to the probabi...
JBuffer1D_t rateL1_Hz
Multiples rate as a function of the multiplicity.
virtual void generateHits(const JModule &module, const JTimeRange &period, JModuleData &output) const
Generate hits.
double getDot(const JFirst_t &first, const JSecond_t &second)
Get dot product of objects.
virtual double getMultiplesRate(const JModuleIdentifier &module, const int M) const =0
Get multiples rate as a function of optical module.
Exception for numerical precision error.
Definition: JException.hh:234
JBuffer1D_t probability1D
This probability vector is a one-dimensional array in which element [i] corresponds to the probabilit...
Auxiliary class for object identification.
Definition: JObjectID.hh:27
virtual double getProbability(const double ct) const =0
Get probability of coincidence.