Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JPMTAnalogueSignalProcessor.hh
Go to the documentation of this file.
1 #ifndef __JDETECTOR__JPMTANALOGUESIGNALPROCESSOR__
2 #define __JDETECTOR__JPMTANALOGUESIGNALPROCESSOR__
3 
4 #include <istream>
5 #include <cmath>
6 #include <limits>
7 
8 #include "JLang/JException.hh"
11 
12 /**
13  * \file
14  *
15  * PMT analogue signal processor.
16  * \author mdejong
17  */
18 namespace JDETECTOR {
19 
21 
22 
23  /**
24  * Specification for time-over-threshold corresponding to a one photo-electron pulse.
25  */
26  const double TIME_OVER_THRESHOLD_NS = 26.4; // [ns]
27 
28 
29  /**
30  * PMT analogue signal processor.
31  *
32  * This class provides for an implementation of the JDETECTOR::JPMTSignalProcessorInterface
33  * using a specific model for the analogue pulse of the PMT.\n
34  * In this, the leading edge of the analogue pulse from the PMT is assumed to be a Gaussian and the tail an exponential.\n
35  * The width of the Gaussian is referred to as the rise time and
36  * the inverse slope of the exponential to the decay time.\n
37  * The two functions are matched at a point where the values and first derivatives are identical.\n
38  *
39  * Note that the decay time is related to the rise time via the specification JDETECTOR::TIME_OVER_THRESHOLD_NS.
40  *
41  * The charge distribution is assumed to be a Gaussian which is centered at the specified gain
42  * and truncated by the specified threshold.
43  *
44  * The transition times are generated according the specified spread as follows.\n
45  * - If the specified transition time spread (TTS) is negative,
46  * the transition times are generated according to the measurements (see method JDETECTOR::getTransitionTime).\n
47  * - If the specified TTS is positive,
48  * the transition times are generated according a Gaussian with a sigma equals to the TTS.\n
49  * - If the TTS is zero, the transition times are generated without any spread.
50  */
53  public JPMTParameters
54  {
55  /**
56  * Constructor.
57  *
58  * \param parameters PMT parameters
59  */
62  JPMTParameters(parameters),
63  decayTime_ns(0.0),
64  t1(0.0),
65  y1(0.0),
66  x1(std::numeric_limits<double>::max())
67  {
68  configure();
69  }
70 
71 
72  /**
73  * Configure internal parameters.
74  *
75  * This method provides the implementations for
76  * - matching of the leading edge of the analogue pulse (Gaussian) and the tail (exponential); and
77  * - determination of number of photo-electrons above which the time-over-threshold
78  * linearly depends on the number of photo-electrons (apart from saturation).
79  *
80  * Note that this method will throw an error if the value of the rise time (i.e. width of the Gaussian)
81  * is too large with respect to the specification JDETECTOR::TIME_OVER_THRESHOLD_NS.
82  */
83  void configure()
84  {
85  static const int N = 100;
86  static const double precision = 1.0e-3;
87 
88  // check rise time
89 
91  THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::configure(): Invalid rise time [ns] " << riseTime_ns);
92  }
93 
94  // decay time
95 
96  const double y = -log(threshold);
97 
98  const double a = y;
99  const double b = riseTime_ns * sqrt(2.0*y) - TIME_OVER_THRESHOLD_NS;
100  const double c = 0.5*riseTime_ns*riseTime_ns;
101 
102  const double Q = b*b - 4.0*a*c;
103 
104  if (Q > 0.0)
105  decayTime_ns = (-b + sqrt(Q)) / (2.0*a);
106  else
107  decayTime_ns = -b / (2.0*a);
108 
109  // fix matching of Gaussian and exponential
110 
111  const double x = riseTime_ns / decayTime_ns;
112 
113  t1 = riseTime_ns*x;
114  y1 = exp(-0.5*x*x);
115 
116  // determine transition point to linear dependence of time-over-threshold as a function of number of photo-electrons
117 
118  const double v = slope / getDerivativeOfSaturation(getTimeOverThreshold(1.0));
119 
120  x1 = std::numeric_limits<double>::max(); // disable linearisation
121 
122  double xmin = 1.0;
123  double xmax = 1.0 / (getDerivative(1.0) * slope);
124 
125  for (int i = 0; i != N; ++i) {
126 
127  const double x = 0.5 * (xmin + xmax);
128  const double u = getDerivative(x) * v;
129 
130  if (fabs(1.0 - u) < precision) {
131  break;
132  }
133 
134  if (u < 1.0)
135  xmin = x;
136  else
137  xmax = x;
138  }
139 
140  x1 = 0.5 * (xmin + xmax);
141  }
142 
143 
144  /**
145  * Get decay time.
146  *
147  * \return decay time [ns]
148  */
149  double getDecayTime() const
150  {
151  return decayTime_ns;
152  }
153 
154 
155  /**
156  * Get time at transition point from Gaussian to exponential.
157  *
158  * \return time [ns]
159  */
160  double getT1() const
161  {
162  return t1;
163  }
164 
165 
166  /**
167  * Get amplitude at transition point from Gaussian to exponential.
168  *
169  * \return amplitude [npe]
170  */
171  double getY1() const
172  {
173  return y1;
174  }
175 
176 
177  /**
178  * Get transition point from a model dependent to linear relation between time-over-threshold and number of photo-electrons.
179  *
180  * \return number of photo-electrons [npe]
181  */
182  double getStartOfLinearisation() const
183  {
184  return x1;
185  }
186 
187 
188  /**
189  * Get amplitude at given time for a one photo-electron pulse.
190  *
191  * \param t1_ns time [ns]
192  * \return amplitude [npe]
193  */
194  double getAmplitude(const double t1_ns) const
195  {
196  if (t1_ns < t1) {
197 
198  const double x = t1_ns / riseTime_ns;
199 
200  return exp(-0.5*x*x); // Gaussian
201 
202  } else {
203 
204  const double x = t1_ns / decayTime_ns;
205 
206  return exp(-x) / y1; // exponential
207  }
208  }
209 
210 
211  /**
212  * Get time to pass from threshold to top of analogue pulse.\n
213  * In this, the leading edge of the analogue pulse is assumed to be Gaussian.
214  *
215  * \param npe number of photo-electrons
216  * \param th threshold [npe]
217  * \return time [ns]
218  */
219  double getRiseTime(const double npe, const double th) const
220  {
221  return riseTime_ns * sqrt(2.0*log(npe/th)); // Gaussian
222  }
223 
224 
225  /**
226  * Get time to pass from top of analogue pulse to threshold.\n
227  * In this, the trailing edge of the analogue pulse is assumed to be exponential.
228  *
229  * \param npe number of photo-electrons
230  * \param th threshold [npe]
231  * \return time [ns]
232  */
233  double getDecayTime(const double npe, const double th) const
234  {
235  if (npe*y1 > th)
236  return decayTime_ns * (log(npe/th) - log(y1)); // exponential
237  else
238  return riseTime_ns * sqrt(2.0*log(npe/th)); // Gaussian
239  }
240 
241 
242  /**
243  * Get maximal rise time for given threshold.
244  *
245  * Note that the rise time is entirely constraint by the specification JDETECTOR::TIME_OVER_THRESHOLD_NS.
246  *
247  * \param th threshold [npe]
248  * \return rise time [ns]
249  */
250  static double getMaximalRiseTime(const double th)
251  {
252  if (th > 0.0 && th < 1.0)
253  return 0.5 * TIME_OVER_THRESHOLD_NS / sqrt(-2.0*log(th));
254  else
255  THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::getMaximalRiseTime(): Invalid threshold " << th);
256  }
257 
258 
259  /**
260  * Get function value with saturation.
261  *
262  * \param y value
263  * \return value
264  */
265  double getSaturation(const double y) const
266  {
267  return saturation / sqrt(y*y + saturation*saturation);
268  }
269 
270 
271  /**
272  * Get derivative of function value with saturation.
273  *
274  * \param y value
275  * \return value
276  */
277  double getDerivativeOfSaturation(const double y) const
278  {
279  return (saturation*saturation + y*y) * sqrt(saturation*saturation + y*y) / (saturation * saturation * saturation);
280  }
281 
282 
283  /**
284  * Get gain spread for given number of photo-electrons.
285  *
286  * \param NPE number of photo-electrons
287  * \return gain spread
288  */
289  double getGainSpread(int NPE) const
290  {
291  return sqrt((double) NPE * gain) * gainSpread;
292  }
293 
294 
295  /**
296  * Get integral of probability.
297  *
298  * \param xmin minimum number of photo-electrons
299  * \param xmax maximum number of photo-electrons
300  * \param NPE true number of photo-electrons
301  * \return probability
302  */
303  double getIntegralOfProbability(const double xmin, const double xmax, const int NPE) const
304  {
305  double zmin = xmin;
306  double zmax = xmax;
307 
308  if (zmin < threshold) { zmin = threshold; }
309  if (zmax < threshold) { zmax = threshold; }
310 
311  const double mu = NPE * gain;
312  const double sigma = getGainSpread(NPE);
313 
314  return (0.5 * TMath::Erfc((zmin - mu) / sqrt(2.0) / sigma) -
315  0.5 * TMath::Erfc((zmax - mu) / sqrt(2.0) / sigma));
316  }
317 
318 
319  /**
320  * Set PMT parameters.
321  *
322  * \param parameters PMT parameters
323  */
324  void setPMTParameters(const JPMTParameters& parameters)
325  {
326  static_cast<JPMTParameters&>(*this).setPMTParameters(parameters);
327 
328  configure();
329  }
330 
331 
332  /**
333  * Read PMT signal from input.
334  *
335  * \param in input stream
336  * \param object PMT signal
337  * \return input stream
338  */
339  friend std::istream& operator>>(std::istream& in, JPMTAnalogueSignalProcessor& object)
340  {
341  in >> static_cast<JPMTParameters&>(object);
342 
343  object.configure();
344 
345  return in;
346  }
347 
348 
349  /**
350  * Apply relative QE.
351  *
352  * \return true if accepted; false if rejected
353  */
354  virtual bool applyQE() const
355  {
356  if (QE <= 0.0)
357  return false;
358  else if (QE < 1.0)
359  return gRandom->Rndm() < QE;
360  else
361  return true;
362  }
363 
364 
365  /**
366  * Get randomised time according transition time distribution.
367  *
368  * \param t_ns time [ns]
369  * \return time [ns]
370  */
371  virtual double getRandomTime(const double t_ns) const
372  {
373  if (TTS_ns < 0.0)
374  return t_ns + getTransitionTime(gRandom->Rndm());
375  else
376  return gRandom->Gaus(t_ns, TTS_ns);
377  }
378 
379 
380  /**
381  * Compare (arrival times of) photo-electrons.
382  *
383  * \param first first photo-electron
384  * \param second second photo-electron
385  * \return true if arrival times of photo-electrons are within rise time of analogue pulses; else false
386  */
387  virtual bool compare(const JPhotoElectron& first, const JPhotoElectron& second) const
388  {
389  return second.t_ns < first.t_ns + riseTime_ns;
390  }
391 
392 
393  /**
394  * Get randomised amplitude according gain and gain spread.
395  *
396  * \param NPE number of photo-electrons
397  * \return number of photo-electrons
398  */
399  virtual double getRandomAmplitude(const int NPE) const
400  {
401  double x;
402 
403  do {
404  x = gRandom->Gaus(NPE*gain, getGainSpread(NPE));
405  } while (x < 0.0);
406 
407  return x;
408  }
409 
410 
411  /**
412  * Get probability for given charge.
413  *
414  * \param npe observed number of photo-electrons
415  * \param NPE true number of photo-electrons
416  * \return probability [npe^-1]
417  */
418  virtual double getProbability(const double npe, const int NPE) const
419  {
420  if (npe >= threshold) {
421 
422  const double mu = NPE * gain;
423  const double sigma = getGainSpread(NPE);
424  const double V = 0.5 * TMath::Erfc((threshold - mu) / sqrt(2.0) / sigma);
425 
426  return TMath::Gaus(npe, mu, sigma, kTRUE) / V;
427  }
428 
429  return 0.0;
430  }
431 
432 
433  /**
434  * Apply threshold.
435  *
436  * \param npe number of photo-electrons
437  * \return true if pass; else false
438  */
439  virtual bool applyThreshold(const double npe) const
440  {
441  return npe >= threshold;
442  }
443 
444 
445  /**
446  * Get time to reach threshold.
447  *
448  * Note that the rise time is defined to be zero for a one photo-electron signal.
449  *
450  * \param npe number of photo-electrons
451  * \return time [ns]
452  */
453  virtual double getRiseTime(const double npe) const
454  {
455  if (slewing)
456  return ((getRiseTime(npe, getTH0()) - getRiseTime(npe, this->threshold)) -
457  (getRiseTime(1.0, getTH0()) - getRiseTime(1.0, this->threshold)));
458  else
459  return 0.0;
460  }
461 
462 
463  /**
464  * Get time-over-threshold (ToT).
465  *
466  * \param npe number of photo-electrons
467  * \return ToT [ns]
468  */
469  virtual double getTimeOverThreshold(const double npe) const
470  {
471  if (npe > threshold) {
472 
473  double tot = 0.0;
474 
475  if (npe*y1 > threshold) {
476 
477  if (npe <= getStartOfLinearisation()) {
478 
479  tot += getRiseTime (npe, threshold); // Gaussian
480  tot += getDecayTime(npe, threshold); // exponential
481 
482  } else {
483 
484  tot += getRiseTime (getStartOfLinearisation(), threshold); // Gaussian
485  tot += getDecayTime(getStartOfLinearisation(), threshold); // exponential
486 
487  tot += slope * (npe - getStartOfLinearisation()); // linear
488  }
489 
490  } else {
491 
492  tot += getRiseTime(npe, threshold); // Gaussian
493  tot += getRiseTime(npe, threshold); // Gaussian
494  }
495 
496  return tot * getSaturation(tot);
497 
498  } else {
499 
500  return 0.0;
501  }
502  }
503 
504 
505  /**
506  * Get derivative of number of photo-electrons to time-over-threshold.
507  *
508  * \param npe number of photo-electrons
509  * \return dnpe/dToT [ns^-1]
510  */
511  virtual double getDerivative(const double npe) const
512  {
513  if (npe >= threshold) {
514 
515  const double z = riseTime_ns / sqrt(2.0 * log(npe/threshold));
516 
517  double y = 0.0;
518 
519  if (npe*y1 > threshold) {
520 
521  if (npe <= getStartOfLinearisation())
522  y = npe / (z + decayTime_ns); // Gaussian + exponential
523  else
524  y = 1.0 / slope; // linear
525 
526  } else {
527 
528  y = npe / (2.0 * z); // Gaussian + Gaussian
529  }
530 
532 
533  } else {
534 
535  return 0.0;
536  }
537  }
538 
539 
540  /**
541  * Probability that a hit survives the simulation of the PMT.
542  *
543  * \param NPE number of photo-electrons
544  * \return probability
545  */
546  virtual double getSurvivalProbability(const int NPE) const
547  {
548  if (QE <= 0.0) {
549 
550  return 0.0;
551 
552  } else if (QE <= 1.0) {
553 
554  double P = 0.0;
555 
556  for (int i = 1; i <= NPE; ++i) {
557 
558  const double p = TMath::Binomial(NPE, i) * TMath::Power(QE, i) * TMath::Power(1.0 - QE, NPE - i);
559  const double mu = i * gain;
560  const double sigma = getGainSpread(i);
561  const double V = 0.5 * TMath::Erfc((threshold - mu) / sqrt(2.0) / sigma);
562  const double W = 0.5 * TMath::Erfc((0.0 - mu) / sqrt(2.0) / sigma);
563 
564  P += p*V/W;
565  }
566 
567  return P;
568 
569  } else {
570 
571  THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::getSurvivalProbability: Invalid QE " << QE);
572  }
573  }
574 
575 
576  /**
577  * Get lower threshold for rise time evaluation.
578  *
579  * \return threshold [npe]
580  */
581  static double getTH0()
582  {
583  return 0.1;
584  }
585 
586 
587  /**
588  * Get upper threshold for rise time evaluation.
589  *
590  * \return threshold [npe]
591  */
592  static double getTH1()
593  {
594  return 0.9;
595  }
596 
597  protected:
598  double decayTime_ns; //!< decay time [ns]
599  double t1; //!< time at match point [ns]
600  double y1; //!< amplitude at match point [npe]
601  /**
602  * Transition point from a logarithmic to a linear relation
603  * between time-over-threshold and number of photo-electrons.\n
604  * Measurements by B. Schermer and R. Bruijn at Nikhef.
605  */
606  double x1;
607  };
608 }
609 
610 #endif
double gain
gain [unit]
Exceptions.
double x1
Transition point from a logarithmic to a linear relation between time-over-threshold and number of ph...
void setPMTParameters(const JPMTParameters &parameters)
Set PMT parameters.
double getDerivativeOfSaturation(const double y) const
Get derivative of function value with saturation.
double gainSpread
gain spread [unit]
Data structure for single photo-electron.
const double TIME_OVER_THRESHOLD_NS
Specification for time-over-threshold corresponding to a one photo-electron pulse.
void setPMTParameters(const JPMTParameters &parameters)
Set PMT parameters.
virtual double getTimeOverThreshold(const double npe) const
Get time-over-threshold (ToT).
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:633
virtual bool applyQE() const
Apply relative QE.
double threshold
threshold [npe]
static double getTH0()
Get lower threshold for rise time evaluation.
double getY1() const
Get amplitude at transition point from Gaussian to exponential.
double riseTime_ns
rise time of analogue pulse [ns]
bool slewing
time slewing of analogue signal
static double getMaximalRiseTime(const double th)
Get maximal rise time for given threshold.
void configure()
Configure internal parameters.
JPMTAnalogueSignalProcessor(const JPMTParameters &parameters=JPMTParameters())
Constructor.
friend std::istream & operator>>(std::istream &in, JPMTAnalogueSignalProcessor &object)
Read PMT signal from input.
JPMTParameters()
Default constructor.
double getSaturation(const double y) const
Get function value with saturation.
virtual bool compare(const JPhotoElectron &first, const JPhotoElectron &second) const
Compare (arrival times of) photo-electrons.
virtual double getProbability(const double npe, const int NPE) const
Get probability for given charge.
virtual double getSurvivalProbability(const int NPE) const
Probability that a hit survives the simulation of the PMT.
static const JTransitionTimeGenerator getTransitionTime
Function object to generate transition time.
double getAmplitude(const double t1_ns) const
Get amplitude at given time for a one photo-electron pulse.
virtual bool applyThreshold(const double npe) const
Apply threshold.
virtual double getRiseTime(const double npe) const
Get time to reach threshold.
virtual double getRandomTime(const double t_ns) const
Get randomised time according transition time distribution.
double getGainSpread(int NPE) const
Get gain spread for given number of photo-electrons.
double TTS_ns
transition time spread [ns]
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:144
double getT1() const
Get time at transition point from Gaussian to exponential.
Data structure for PMT parameters.
double saturation
saturation [ns]
virtual double getRandomAmplitude(const int NPE) const
Get randomised amplitude according gain and gain spread.
virtual double getDerivative(const double npe) const
Get derivative of number of photo-electrons to time-over-threshold.
static double getTH1()
Get upper threshold for rise time evaluation.
double QE
relative quantum efficiency
double getDecayTime(const double npe, const double th) const
Get time to pass from top of analogue pulse to threshold.
double getIntegralOfProbability(const double xmin, const double xmax, const int NPE) const
Get integral of probability.
double getStartOfLinearisation() const
Get transition point from a model dependent to linear relation between time-over-threshold and number...
double y1
amplitude at match point [npe]
double slope
slope [ns/npe]
double getRiseTime(const double npe, const double th) const
Get time to pass from threshold to top of analogue pulse.