Jpp - the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JPMTSignalProcessorInterface.hh
Go to the documentation of this file.
1 
2 #ifndef __JDETECTOR__JPMTSIGNALPROCESSORINTERFACE__
3 #define __JDETECTOR__JPMTSIGNALPROCESSORINTERFACE__
4 
5 #include <algorithm>
6 #include <cmath>
7 #include <limits>
8 
9 #include "JLang/JCC.hh"
10 #include "JLang/JException.hh"
13 
14 
15 /**
16  * \author mdejong
17  */
18 
19 namespace JDETECTOR {}
20 namespace JPP { using namespace JDETECTOR; }
21 
22 namespace JDETECTOR {
23 
24 
26 
27 
28  /**
29  * PMT signal processor interface.
30  *
31  * This class supports the implementation of the PMT simulator interface using
32  * an alternative set of virtual methods.
33  * These methods constitute a user interface to the expected performance of a PMT.
34  *
35  * Each photon is converted to a photo-electron using member method JPMTSignalProcessorInterface::getRandomTime.
36  * For this, the data structure JPhotoElectron is used.
37  * Note that the quantum efficiency (QE) of the PMT actually is included
38  * in the simulation of the detector response.
39  * All photo-electrons are then time sorted.
40  * The photo-electrons which simultaneously arrive are merged.
41  * The corresponding condition is defined by member method JPMTSignalProcessorInterface::compare.
42  * The time of the combined signal is determined by the time of the first photo-electron and
43  * the rise time of the analogue pulse to the threshold of the discriminator.
44  * In this, the actual amplitude of the combined analogue signal and the calibration of the PMT are taken into account.
45  * The amplitude of the combined analogue signal is simulated using member method JPMTSignalProcessorInterface::getRandomCharge.
46  * For this, the data structure JPMTSignal is used.
47  *
48  * The analogue signals of the PMT are processed by a discriminator.
49  * This discriminator produces a time-over-threshold signal for each analogue signal that passes a preset threshold.
50  * The output signal is described by the time of the leading edge and the length of the time-over-threshold signal.
51  * For this, the data structure JPMTPulse is used.
52  * The determination of the time of the leading edge and the length of the time-over-threshold signal
53  * require a designated model.
54  * The class JPMTAnalogueSignalProcessor provides for an implementation of such a model.
55  *
56  * Overlapping time-over-threshold pulses are merged.
57  * The time of the leading edge is then set to the earliest leading edge and
58  * the time-over-threshold to the difference between
59  * the latest trailing edge and the earliest leading edge.
60  * The merging of hits is implemented in member method JPMTSignalProcessorInterface::merge.
61  *
62  * The default implementation of these virtual methods corresponds to no smearing
63  * and a time-over-threshold value equal to a fixed two photo-electron resolution.
64  * The width of the charge distribution and the two photo-electron resolution are implemented in methods
65  * JPMTSignalProcessorInterface::getQmin() and JPMTSignalProcessorInterface::getTmin(), respectively.
66  * For a realistic PMT simulation, the derived class should provide for
67  * an implementation of each of these virtual methods.
68  */
70  public:
71  /**
72  * Default constructor.
73  */
74  JPMTSignalProcessorInterface(const double max_charge = 100.0) :
75  MAX_CHARGE(max_charge)
76  {}
77 
78 
79  /**
80  * Virtual destructor.
81  */
83  {}
84 
85 
86  /**
87  * Threshold domain specifiers
88  */
90  BELOW_THRESHOLDBAND = -1, //!< below threshold band
91  IN_THRESHOLDBAND = 0, //!< inside threshold band
92  ABOVE_THRESHOLD = 2 //!< above threshold band
93  };
94 
95 
96  /**
97  * Process hits.
98  *
99  * \param calibration PMT calibration
100  * \param input PMT signals
101  * \param output PMT hits
102  */
103  void operator()(const JCalibration& calibration,
104  const JPMTData<JPMTSignal>& input,
105  JPMTData<JPMTPulse>& output) const
106  {
107  // apply transition time distribution to each photo-electron.
108 
109  buffer.clear();
110 
111  for (JPMTData<JPMTSignal>::const_iterator hit = input.begin(); hit != input.end(); ++hit) {
112  for (int i = 0; i != hit->npe; ++i) {
113  if (applyQE()) {
114  buffer.push_back(JPhotoElectron(getRandomTime(hit->t_ns)));
115  }
116  }
117  }
118 
119  if (!buffer.empty()) {
120 
122 
123  buffer.sort();
124 
125 
126  // generate PMT hits from time sequence of photo-electrons.
127 
128  for (JPMTData<JPhotoElectron>::const_iterator q = buffer.begin(), p = q++; q != buffer.end(); ++q) {
129 
130  while (compare(*p,*q)) {
131  ++q;
132  }
133 
134  const double npe = getRandomCharge(distance(p,q));
135  if (applyThreshold(npe) > BELOW_THRESHOLDBAND) {
136  output.push_back(JPMTPulse(putTime(p->t_ns + getRiseTime(npe), calibration), getTimeOverThreshold(npe)));
137  }
138 
139  p = q;
140  }
141 
142  // merge overlapping PMT hits.
143 
144  merge(output);
145  }
146  }
147 
148 
149  /**
150  * Apply relative QE.
151  * The default implementation returns always <tt>true</tt>.
152  *
153  * \return true
154  */
155  virtual bool applyQE() const
156  {
157  return true;
158  }
159 
160 
161  /**
162  * Get randomised time according transition time distribution.
163  * The default implementation returns the given value.
164  *
165  * \param t_ns time [ns]
166  * \return time [ns]
167  */
168  virtual double getRandomTime(const double t_ns) const
169  {
170  return t_ns;
171  }
172 
173 
174  /**
175  * Compare arrival times of photo-electrons.
176  *
177  * Two (or more) photo-electrons are merged if they are comparable.
178  * The default implementation returns uses JPMTSignalProcessorInterface::getTmin() as time window.
179  *
180  * \param first first photo-electron
181  * \param second second photo-electron
182  * \return false
183  */
184  virtual bool compare(const JPhotoElectron& first, const JPhotoElectron& second) const
185  {
186  return second.t_ns - first.t_ns <= getTmin();
187  }
188 
189 
190  /**
191  * Get randomised charge in accordance with gain and gain spread.
192  * The default implementation returns the given value.
193  *
194  * \param NPE number of photo-electrons
195  * \return number of photo-electrons
196  */
197  virtual double getRandomCharge(const int NPE) const
198  {
199  return (double) NPE;
200  }
201 
202 
203  /**
204  * Get probability for given charge.
205  *
206  * \param npe observed number of photo-electrons
207  * \param NPE true number of photo-electrons
208  * \return probability [npe^-1]
209  */
210  virtual double getChargeProbability(const double npe, const int NPE) const
211  {
212  return (fabs(npe - NPE) < 0.5*getQmin() ? 1.0/getQmin() : 0.0);
213  }
214 
215 
216  /**
217  * Apply threshold.
218  * The default implementation returns always <tt>ABOVE_THRESHOLD</tt>.
219  *
220  * \param npe number of photo-electrons
221  * \return threshold domain;
222  */
223  virtual JThresholdDomains applyThreshold(const double npe) const
224  {
225  return ABOVE_THRESHOLD;
226  }
227 
228 
229  /**
230  * Get time to pass threshold.
231  * The default implementation returns <tt>0</tt>.
232  *
233  * \param npe number of photo-electrons
234  * \return time [ns]
235  */
236  virtual double getRiseTime(const double npe) const
237  {
238  return 0.0;
239  }
240 
241 
242  /**
243  * Get time over threshold (ToT).
244  * The default implementation returns JPMTSignalProcessorInterface::getTmin().
245  *
246  * \param npe number of photo-electrons
247  * \return ToT [ns]
248  */
249  virtual double getTimeOverThreshold(const double npe) const
250  {
251  return getTmin();
252  }
253 
254 
255  /**
256  * Get derivative of number of photo-electrons to time over threshold.
257  * The default implementation returns infinity.
258  *
259  * \param npe number of photo-electrons
260  * \return dnpe/dToT [ns^-1]
261  */
262  virtual double getDerivative(const double npe) const
263  {
264  return std::numeric_limits<double>::max();
265  }
266 
267 
268  /**
269  * Probability that a hit survives the simulation of the PMT.
270  * The default implementation returns <tt>1</tt> if given value larger than <tt>0</tt>.
271  *
272  * \param NPE number of photo-electrons
273  * \return probability
274  */
275  virtual double getSurvivalProbability(const int NPE) const
276  {
277  if (NPE > 0)
278  return 1.0;
279  else
280  return 0.0;
281  }
282 
283 
284  /**
285  * Merging of PMT hits.
286  *
287  * Hits with overlapping time-over-threshold signals should -de facto- be combined.
288  * In this, the leading edge is maintained and the time-over-threshold is
289  * set to the difference between the overall trailing and leading edges.
290  * As a result, the number of PMT hits may be reduced.
291  *
292  * \param data PMT hits
293  */
294  virtual void merge(JPMTData<JPMTPulse>& data) const
295  {
296  using namespace std;
297 
298  JPMTData<JPMTPulse>::iterator out = data.begin();
299 
300  for (JPMTData<JPMTPulse>::iterator i = data.begin(); i != data.end(); ) {
301 
302  double t1 = i->t_ns;
303  double t2 = i->t_ns + i->tot_ns;
304 
305  while (++i != data.end() && i->t_ns < t2 + getTmin()) {
306  t2 = max(t2, i->t_ns + i->tot_ns);
307  }
308 
309  out->t_ns = t1;
310  out->tot_ns = t2 - t1;
311 
312  ++out;
313  }
314 
315  data.resize(distance(data.begin(), out));
316  }
317 
318 
319  /**
320  * Get number of photo-electrons.
321  * The default implementation finds the time-over-threshold iteratively.
322  *
323  * \param tot_ns time over threshold [ns]
324  * \param eps precision
325  * \return number of photo-electrons
326  */
327  virtual double getNPE(const double tot_ns,
328  const double eps = 1.0e-3) const
329  {
330  static const int N = 100;
331 
332  if (tot_ns > 0.0) {
333 
334  double xmin = 0.0;
335  double xmax = MAX_CHARGE;
336 
337  for (int i = 0; i != N; ++i) {
338 
339  const double x = 0.5 * (xmin + xmax);
340  const double y = this->getTimeOverThreshold(x);
341 
342  if (fabs(y - tot_ns) < eps) {
343  return x;
344  }
345 
346  if (y < tot_ns)
347  xmin = x;
348  else
349  xmax = x;
350  }
351 
352  return 0.5 * (xmin + xmax);
353 
354  } else {
355 
356  double xmin = 0.0;
357  double xmax = 2.0;
358 
359  for (int i = 0; i != N; ++i) {
360 
361  const double x = 0.5 * (xmin + xmax);
362  const double y = this->getTimeOverThreshold(x);
363 
364  if (y == 0.0)
365  xmin = x;
366  else
367  xmax = x;
368  }
369 
370  return 0.5 * (xmin + xmax);
371  }
372  }
373 
374 
375  /**
376  * Get two photo-electron resolution for time-over-threshold
377  *
378  * \return minimal time [ns]
379  */
380  static double getTmin()
381  {
382  return 1.0;
383  }
384 
385 
386  /**
387  * Get minimal width of charge distribution.
388  *
389  * \return minimal charge [npe]
390  */
391  static double getQmin()
392  {
393  return 1.0e-3;
394  }
395 
396 
397  protected:
398  double MAX_CHARGE; //!< Maximum charge [npe]
399 
400  private:
402 
403  };
404 
405 
406  /**
407  * Get charge probability.
408  *
409  * \param pmt PMT signal processor
410  * \param npe measured number of photo-electrons
411  * \param NPE expected number of photo-electrons
412  * \param precision precision
413  * \return probability
414  */
416  const double npe,
417  const double NPE,
418  const double precision = 1.0e-4)
419  {
420  int i = (int) (NPE - 5.0 * sqrt(NPE));
421 
422  if (i < 1) {
423  i = 1;
424  }
425 
426  double p = NPE * exp(-NPE) / (double) 1;
427 
428  for (int __i = 1; __i != i; ++__i) {
429  p *= NPE / __i;
430  }
431 
432  double P = 0.0;
433 
434  for (double p0 = 0.0; (p >= p0 || p > precision); ++i, p0 = p, p *= NPE / (double) i) {
435  P += pmt.getChargeProbability(npe, i) * p;
436  }
437 
438  return P;
439  }
440 }
441 
442 #endif
Exceptions.
Data structure for PMT digital pulse.
static double getQmin()
Get minimal width of charge distribution.
virtual double getRandomCharge(const int NPE) const
Get randomised charge in accordance with gain and gain spread.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Data structure for single photo-electron.
JPMTSignalProcessorInterface(const double max_charge=100.0)
Default constructor.
virtual void merge(JPMTData< JPMTPulse > &data) const
Merging of PMT hits.
virtual double getNPE(const double tot_ns, const double eps=1.0e-3) const
Get number of photo-electrons.
Data structure for time calibration.
virtual bool applyQE() const
Apply relative QE.
static double getTmin()
Get two photo-electron resolution for time-over-threshold.
then set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))*exp(-0.5 *(y-[1])*(y-[1])/([2]*[2]))" JF2 -o $WORKDIR/f2.root -F "$FORMULA" -@ "p0
void operator()(const JCalibration &calibration, const JPMTData< JPMTSignal > &input, JPMTData< JPMTPulse > &output) const
Process hits.
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
virtual double getRiseTime(const double npe) const
Get time to pass threshold.
virtual double getSurvivalProbability(const int NPE) const
Probability that a hit survives the simulation of the PMT.
Compiler version dependent expressions, macros, etc.
double putTime(const T &t1, const JCalibration &cal)
Get de-calibrated time.
std::vector< JElement_t >::iterator iterator
virtual ~JPMTSignalProcessorInterface()
Virtual destructor.
virtual JThresholdDomains applyThreshold(const double npe) const
Apply threshold.
virtual bool compare(const JPhotoElectron &first, const JPhotoElectron &second) const
Compare arrival times of photo-electrons.
virtual double getChargeProbability(const double npe, const int NPE) const
Get probability for given charge.
static JPhotoElectron getEndMarker()
Get end marker.
virtual double getRandomTime(const double t_ns) const
Get randomised time according transition time distribution.
virtual double getTimeOverThreshold(const double npe) const
Get time over threshold (ToT).
Template data structure for PMT I/O.
virtual double getDerivative(const double npe) const
Get derivative of number of photo-electrons to time over threshold.
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:162
std::vector< JElement_t >::const_iterator const_iterator
then usage $script[input file[working directory[option]]] nWhere option can be N
Definition: JMuonPostfit.sh:37
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:60