Jpp  18.1.0
the software that should make you happy
 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) {
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 get_sigma();
281  }
282 
283 
284  /**
285  * Set intrinsic time smearing of K40 coincidences.
286  *
287  * \param sigma sigma [ns]
288  */
289  static void setSigma(const double sigma)
290  {
291  get_sigma() = sigma;
292  }
293 
294  private:
295  /**
296  * Get intrinsic time smearing of K40 coincidences.
297  *
298  * \return sigma [ns]
299  */
300  static double& get_sigma()
301  {
302  static double sigma = 0.5;
303 
304  return sigma;
305  }
306 
307 
308  /**
309  * This correlation matrix is a two-dimensional array in which element [i][j]
310  * corresponds to the probability of a genuine coincidence due to K40 decays,
311  * where i and j refer to the indices of the PMTs in the optical module.
312  */
314 
315  /**
316  * This probability vector is a one-dimensional array in which element [i]
317  * corresponds to the probability of a genuine coincidence due to K40 decays,
318  * where i refers to the index of the PMT in the optical module.
319  */
321 
322  /**
323  * This probability vector is a one-dimensional array in which element [i]
324  * corresponds to the probability of an additional hit,
325  * where i refers to the index of the PMT in the optical module.
326  */
328 
329  /**
330  * Multiples rate as a function of the multiplicity.
331  * The index i corresponds to multiplicity M = i + 2.
332  */
334  };
335 }
336 
337 #endif
static double getSigma()
Get intrinsic time smearing of K40 coincidences.
Exceptions.
Data structure for PMT analogue signal.
T getLowerLimit() const
Get lower limit.
Definition: JRange.hh:202
Auxiliary methods for geometrical methods.
Interface for simulation of K40 background.
Data structure for a composite optical module.
Definition: JModule.hh:68
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:712
JBuffer1D_t probabilityND
This probability vector is a one-dimensional array in which element [i] corresponds to the probabilit...
const double sigma[]
Definition: JQuadrature.cc:74
static unsigned int getRandomIndex(const JBuffer1D_t &buffer, const double random)
Get index based on random value.
Compiler version dependent expressions, macros, etc.
static double & get_sigma()
Get intrinsic time smearing of K40 coincidences.
virtual double getSinglesRate(const JPMTIdentifier &pmt) const =0
Get singles rate as a function of PMT.
int getID() const
Get identifier.
Definition: JObjectID.hh:50
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.
T getUpperLimit() const
Get upper limit.
Definition: JRange.hh:213
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:268
JBuffer1D_t probability1D
This probability vector is a one-dimensional array in which element [i] corresponds to the probabilit...
then usage $script< input file >[option[primary[working directory]]] nWhere option can be N
Definition: JMuonPostfit.sh:40
Auxiliary class for object identification.
Definition: JObjectID.hh:22
static void setSigma(const double sigma)
Set intrinsic time smearing of K40 coincidences.
int j
Definition: JPolint.hh:703
virtual double getProbability(const double ct) const =0
Get probability of coincidence.