Jpp  15.0.2
the software that should make you happy
 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 "TRandom3.h"
9 #include "TMath.h"
10 
11 #include "JLang/JException.hh"
16 
17 /**
18  * \file
19  *
20  * PMT analogue signal processor.
21  * \author mdejong
22  */
23 namespace JDETECTOR {
24 
26 
27 
28  /**
29  * PMT analogue signal processor.
30  *
31  * This class provides for an implementation of the JDETECTOR::JPMTSignalProcessorInterface
32  * using a specific model for the analogue pulse of the PMT.\n
33  * In this, the leading edge of the analogue pulse from the PMT is assumed to be a Gaussian and the tail an exponential.\n
34  * The width of the Gaussian is referred to as the rise time and
35  * the inverse slope of the exponential to the decay time.\n
36  * The two functions are matched at a point where the values and first derivatives are identical.\n
37  *
38  * Note that the decay time is related to the rise time via the specification JDETECTOR::TIME_OVER_THRESHOLD_NS.
39  *
40  * The charge distribution is assumed to be a Gaussian which is centered at the specified gain
41  * and truncated by the specified threshold.
42  *
43  * The transition times are generated according the specified spread as follows.
44  * - If the specified transition-time spread (TTS) is negative,
45  * the transition times are generated according to measurements (see method JDETECTOR::getTransitionTime).\n
46  * In this, the negated integral value of the TTS corresponds to the option
47  * which in turn corresponds to the detector identifier of the measurements.
48  * - If the specified TTS is positive,
49  * the transition times are generated according a Gaussian with a sigma equals to the given TTS.
50  * - If the TTS is zero, the transition times are generated without any spread.
51  */
54  public JPMTParameters
55  {
56  /**
57  * Constructor.
58  *
59  * \param parameters PMT parameters
60  */
64  decayTime_ns(0.0),
65  t1(0.0),
66  y1(0.0),
67  x1(std::numeric_limits<double>::max())
68  {
69  configure();
70  }
71 
72 
73  /**
74  * Configure internal parameters.
75  *
76  * This method provides the implementations for
77  * - matching of the leading edge of the analogue pulse (Gaussian) and the tail (exponential); and
78  * - determination of number of photo-electrons above which the time-over-threshold
79  * linearly depends on the number of photo-electrons (apart from saturation).
80  *
81  * Note that this method will throw an error if the value of the rise time (i.e.\ width of the Gaussian)
82  * is too large with respect to the specification JDETECTOR::TIME_OVER_THRESHOLD_NS.
83  */
84  void configure()
85  {
86  static const int N = 100;
87  static const double precision = 1.0e-4;
88 
89  // check thresholdband
90 
91  if (threshold - thresholdBand < getTH0()) {
92  THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::configure(): Invalid thresholdband [npe] " << thresholdBand);
93  }
94 
95  // check rise time
96 
98  THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::configure(): Invalid rise time [ns] " << riseTime_ns);
99  }
100 
101  // decay time
102 
103  const double y = -log(threshold);
104 
105  const double a = y;
106  const double b = riseTime_ns * sqrt(2.0*y) - TIME_OVER_THRESHOLD_NS;
107  const double c = 0.5*riseTime_ns*riseTime_ns;
108 
109  const double Q = b*b - 4.0*a*c;
110 
111  if (Q > 0.0)
112  decayTime_ns = (-b + sqrt(Q)) / (2.0*a);
113  else
114  decayTime_ns = -b / (2.0*a);
115 
116  // fix matching of Gaussian and exponential
117 
118  const double x = riseTime_ns / decayTime_ns;
119 
120  t1 = riseTime_ns*x;
121  y1 = exp(-0.5*x*x);
122 
123  // determine transition point to linear dependence of time-over-threshold as a function of number of photo-electrons
124 
125  const double xs = saturation; // disable saturation
126 
127  saturation = 1.0e50;
128 
129  x1 = std::numeric_limits<double>::max(); // disable linearisation
130 
131  double xmin = 1.0;
132  double xmax = 1.0 / (getDerivative(1.0) * slope);
133 
134  for (int i = 0; i != N; ++i) {
135 
136  const double x = 0.5 * (xmin + xmax);
137  const double u = getDerivative(x) * slope;
138 
139  if (fabs(1.0 - u) < precision) {
140  break;
141  }
142 
143  if (u < 1.0)
144  xmin = x;
145  else
146  xmax = x;
147  }
148 
149  x1 = 0.5 * (xmin + xmax);
150 
151  saturation = xs;
152  }
153 
154 
155  /**
156  * Get decay time.
157  *
158  * \return decay time [ns]
159  */
160  double getDecayTime() const
161  {
162  return decayTime_ns;
163  }
164 
165 
166  /**
167  * Get time at transition point from Gaussian to exponential.
168  *
169  * \return time [ns]
170  */
171  double getT1() const
172  {
173  return t1;
174  }
175 
176 
177  /**
178  * Get amplitude at transition point from Gaussian to exponential.
179  *
180  * \return amplitude [npe]
181  */
182  double getY1() const
183  {
184  return y1;
185  }
186 
187 
188  /**
189  * Get transition point from a model-dependent to linear relation between time-over-threshold and number of photo-electrons.
190  *
191  * \return number of photo-electrons [npe]
192  */
193  double getStartOfLinearisation() const
194  {
195  return x1;
196  }
197 
198 
199  /**
200  * Get amplitude at given time for a one photo-electron pulse.
201  *
202  * \param t1_ns time [ns]
203  * \return amplitude [npe]
204  */
205  double getAmplitude(const double t1_ns) const
206  {
207  if (t1_ns < t1) {
208 
209  const double x = t1_ns / riseTime_ns;
210 
211  return exp(-0.5*x*x); // Gaussian
212 
213  } else {
214 
215  const double x = t1_ns / decayTime_ns;
216 
217  return exp(-x) / y1; // exponential
218  }
219  }
220 
221 
222  /**
223  * Get time to pass from threshold to top of analogue pulse.\n
224  * In this, the leading edge of the analogue pulse is assumed to be Gaussian.
225  *
226  * \param npe number of photo-electrons
227  * \param th threshold [npe]
228  * \return time [ns]
229  */
230  double getRiseTime(const double npe, const double th) const
231  {
232  return riseTime_ns * sqrt(2.0*log(npe/th)); // Gaussian
233  }
234 
235 
236  /**
237  * Get time to pass from top of analogue pulse to threshold.\n
238  * In this, the trailing edge of the analogue pulse is assumed to be exponential.
239  *
240  * \param npe number of photo-electrons
241  * \param th threshold [npe]
242  * \return time [ns]
243  */
244  double getDecayTime(const double npe, const double th) const
245  {
246  if (npe*y1 > th)
247  return decayTime_ns * (log(npe/th) - log(y1)); // exponential
248  else
249  return riseTime_ns * sqrt(2.0*log(npe/th)); // Gaussian
250  }
251 
252 
253  /**
254  * Get maximal rise time for given threshold.
255  *
256  * Note that the rise time is entirely constrained by the specification JDETECTOR::TIME_OVER_THRESHOLD_NS.
257  *
258  * \param th threshold [npe]
259  * \return rise time [ns]
260  */
261  static double getMaximalRiseTime(const double th)
262  {
263  if (th > 0.0 && th < 1.0)
264  return 0.5 * TIME_OVER_THRESHOLD_NS / sqrt(-2.0*log(th));
265  else
266  THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::getMaximalRiseTime(): Invalid threshold " << th);
267  }
268 
269 
270  /**
271  * Get time-over-threshold with saturation.
272  *
273  * \param T time-over-threshold without saturation
274  * \return time-over-threshold with saturation
275  */
276  double applySaturation(const double T) const
277  {
278  return saturation / sqrt(T*T + saturation*saturation) * T;
279  }
280 
281 
282  /**
283  * Get time-over-threshold without saturation.
284  *
285  * \param tot_ns time-over-threshold with saturation
286  * \return time-over-threshold without saturation
287  */
288  double removeSaturation(const double tot_ns) const
289  {
290  return saturation / sqrt(saturation*saturation - tot_ns*tot_ns) * tot_ns;
291  }
292 
293 
294  /**
295  * Get derivative of saturation factor.
296  *
297  * \param T time-over-threshold without saturation
298  * \return derivative of saturation factor
299  */
300  double getDerivativeOfSaturation(const double T) const
301  {
302  return saturation * saturation * saturation / ((saturation*saturation + T*T) * sqrt(saturation*saturation + T*T));
303  }
304 
305 
306  /**
307  * Get gain spread for given number of photo-electrons.
308  *
309  * \param NPE number of photo-electrons
310  * \return gain spread
311  */
312  double getGainSpread(int NPE) const
313  {
314  return sqrt((double) NPE * gain) * gainSpread;
315  }
316 
317 
318  /**
319  * Get probability density for given charge.
320  *
321  * \param npe observed number of photo-electrons
322  * \param NPE true number of photo-electrons
323  * \return probability [npe^-1]
324  */
325  virtual double getChargeProbability(const double npe, const int NPE) const override
326  {
327  if (applyThreshold(npe) > BELOW_THRESHOLDBAND) {
328 
329  const double f = gainSpread * gainSpread;
330 
331  double norm = 0.0;
332  double prob = 0.0;
333  double weight = (PunderAmplified > 0.0 ? pow(PunderAmplified, NPE) : 1.0);
334 
335  for (int k = (PunderAmplified > 0.0 ? 0 : NPE); k <= NPE; ++k) {
336 
337  const double mu = ((NPE-k) * f * gain) + (k * gain);
338  const double sigma = sqrt((NPE-k) * f + k) * getGainSpread(1);
339 
340  norm += weight * (0.5 * TMath::Erfc(((threshold-thresholdBand) - mu) / sqrt(2.0) / sigma));
341  prob += weight * TMath::Gaus(npe, mu, sigma, kTRUE);
342 
343  weight *= ( (NPE - k) / ((double) (k+1)) *
344  (1.0 - PunderAmplified) / PunderAmplified );
345  }
346 
347  return prob / norm;
348  }
349 
350  return 0.0;
351  }
352 
353 
354  /**
355  * Get integral of probability.
356  *
357  * \param xmin minimum number of photo-electrons
358  * \param xmax maximum number of photo-electrons
359  * \param NPE true number of photo-electrons
360  * \return probability
361  */
362  double getIntegralOfChargeProbability(const double xmin, const double xmax, const int NPE) const
363  {
364  double zmin = xmin;
365  double zmax = xmax;
366 
367  const double th = threshold - thresholdBand;
368  const double f = gainSpread * gainSpread;
369 
370  if (zmin < th) { zmin = th; }
371  if (zmax < th) { zmax = th; }
372 
373  double norm = 0.0;
374  double cumulP = 0.0;
375  double weight = (PunderAmplified > 0.0 ? pow(PunderAmplified, NPE) : 1.0);
376 
377  for (int k = (PunderAmplified > 0.0 ? 0 : NPE); k <= NPE; ++k) {
378 
379  const double mu = ((NPE-k) * f * gain) + (k * gain);
380  const double sigma = sqrt((NPE-k) * f + k) * getGainSpread(1);
381 
382  norm += weight * (0.5 * TMath::Erfc((th - mu) / sqrt(2.0) / sigma));
383  cumulP += weight * (0.5 * TMath::Erfc((zmin - mu) / sqrt(2.0) / sigma) -
384  0.5 * TMath::Erfc((zmax - mu) / sqrt(2.0) / sigma));
385 
386  weight *= ( (NPE - k) / ((double) (k+1)) *
387  (1.0 - PunderAmplified) / PunderAmplified);
388  }
389 
390  return cumulP / norm;
391  }
392 
393 
394  /**
395  * Get integral of probability in specific threshold domain
396  *
397  * \param domain threshold domain
398  * \param NPE true number of photo-electrons
399  * \return probability
400  */
401  double getIntegralOfChargeProbability(JThresholdDomains domain, const int NPE) const
402  {
403  switch (domain) {
404  case ABOVE_THRESHOLD:
406 
407  case IN_THRESHOLDBAND:
409 
410  default:
411  return 0.0;
412  }
413  }
414 
415 
416  /**
417  * Set PMT parameters.
418  *
419  * \param parameters PMT parameters
420  */
422  {
423  static_cast<JPMTParameters&>(*this).setPMTParameters(parameters);
424 
425  configure();
426  }
427 
428 
429  /**
430  * Read PMT signal from input.
431  *
432  * \param in input stream
433  * \param object PMT signal
434  * \return input stream
435  */
436  friend std::istream& operator>>(std::istream& in, JPMTAnalogueSignalProcessor& object)
437  {
438  in >> static_cast<JPMTParameters&>(object);
439 
440  object.configure();
441 
442  return in;
443  }
444 
445 
446  /**
447  * Apply relative QE.
448  *
449  * \return true if accepted; false if rejected
450  */
451  virtual bool applyQE() const override
452  {
453  if (QE <= 0.0)
454  return false;
455  else if (QE < 1.0)
456  return gRandom->Rndm() < QE;
457  else
458  return true;
459  }
460 
461 
462  /**
463  * Get randomised time according transition time distribution.
464  *
465  * \param t_ns time [ns]
466  * \return time [ns]
467  */
468  virtual double getRandomTime(const double t_ns) const override
469  {
470  if (TTS_ns < 0.0)
471  return t_ns + getTransitionTime(gRandom->Rndm(), -lrint(TTS_ns));
472  else if (TTS_ns > 0.0)
473  return gRandom->Gaus(t_ns, TTS_ns);
474  else
475  return t_ns;
476  }
477 
478 
479  /**
480  * Compare (arrival times of) photo-electrons.
481  *
482  * \param first first photo-electron
483  * \param second second photo-electron
484  * \return true if arrival times of photo-electrons are within rise time of analogue pulses; else false
485  */
486  virtual bool compare(const JPhotoElectron& first, const JPhotoElectron& second) const override
487  {
488  return second.t_ns < first.t_ns + riseTime_ns;
489  }
490 
491 
492  /**
493  * Get randomised charge according to gain and gain spread.
494  *
495  * \param NPE number of photo-electrons
496  * \return number of photo-electrons
497  */
498  virtual double getRandomCharge(const int NPE) const override
499  {
500  double q;
501  do {
502 
503  // Determine which contribution to sample from
504  int k = NPE;
505  if (PunderAmplified > 0.0) {
506 
507  const double X = gRandom->Uniform();
508  double weight = pow(1.0 - PunderAmplified, NPE);
509 
510  for (double sum_p = 0.0; k > 0; --k) {
511 
512  sum_p += weight;
513  if (sum_p > X) { break; }
514 
515  weight *= ( k / ((double) (NPE - (k-1))) *
516  PunderAmplified / (1.0 - PunderAmplified) );
517  }
518  }
519 
520  // Sample from chosen contribution
521  const double f = gainSpread * gainSpread;
522  const double mu = ((NPE-k) * f * gain) + (k * gain);
523  const double sigma = sqrt((NPE-k) * f + k) * getGainSpread(1);
524 
525  q = gRandom->Gaus(mu,sigma);
526 
527  } while (q < 0.0);
528 
529  return q;
530  }
531 
532 
533  /**
534  * Get number of photo-electrons.
535  *
536  * \param tot_ns time over threshold [ns]
537  * \param eps precision
538  * \return number of photo-electrons
539  */
540  virtual double getNPE(const double tot_ns,
541  const double eps = 1.0e-3) const
542  {
543  const double t = removeSaturation(tot_ns);
544 
545  double npe;
546  if (t <= 2*getRiseTime(threshold/y1,threshold)) { // Gaussian + Gaussian
547 
548  npe = threshold * exp(t*t/riseTime_ns/riseTime_ns/8.0);
549 
550  } else if (t <= (getRiseTime(getStartOfLinearisation(), threshold) +
552  // Gaussian + Exponential
553 
554  const double A = decayTime_ns;
555  const double B = sqrt(2.0) * riseTime_ns;
556  const double C = -(decayTime_ns*log(y1) + t);
557  const double z = (-B + sqrt(B*B - 4*A*C)) / (2*A);
558 
559  npe = threshold * exp(z*z);
560 
561  } else { // linear
562 
563  const double T = (getRiseTime(getStartOfLinearisation(), threshold) +
565 
566  npe = getStartOfLinearisation() + (t - T) / slope;
567  }
568 
569  return (npe <= MAX_CHARGE ? npe : MAX_CHARGE);
570  }
571 
572 
573  /**
574  * Apply threshold.
575  *
576  * \param npe number of photo-electrons
577  * \return true if pass; else false
578  */
579  virtual JThresholdDomains applyThreshold(const double npe) const override
580  {
581  if (npe > threshold) {
582 
583  return ABOVE_THRESHOLD;
584 
585  } else if (npe > threshold - thresholdBand) {
586 
587  return IN_THRESHOLDBAND;
588 
589  } else {
590 
591  return BELOW_THRESHOLDBAND;
592  }
593  }
594 
595 
596  /**
597  * Get time to reach threshold.
598  *
599  * Note that the rise time is defined to be zero for a one photo-electron signal.
600  *
601  * \param npe number of photo-electrons
602  * \return time [ns]
603  */
604  virtual double getRiseTime(const double npe) const override
605  {
606  if (slewing) {
607 
608  switch (applyThreshold(npe)) {
609  case IN_THRESHOLDBAND:
610  return ((getRiseTime(npe, getTH0()) - getRiseTime(npe, threshold-thresholdBand)) -
611  (getRiseTime(1.0, getTH0()) - getRiseTime(1.0, threshold-thresholdBand))) + this->mean_ns;
612 
613  case ABOVE_THRESHOLD:
614  return ((getRiseTime(npe, getTH0()) - getRiseTime(npe, threshold)) -
615  (getRiseTime(1.0, getTH0()) - getRiseTime(1.0, threshold)));
616 
617  default:
618  THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::getRiseTime: Invalid charge " << npe);
619  }
620 
621  } else {
622 
623  return 0.0;
624  }
625  }
626 
627 
628  /**
629  * Get time-over-threshold (ToT).
630  *
631  * \param npe number of photo-electrons
632  * \return ToT [ns]
633  */
634  virtual double getTimeOverThreshold(const double npe) const override
635  {
636 
637  switch (applyThreshold(npe)) {
638  case IN_THRESHOLDBAND: {
639 
640  return gRandom->Gaus(mean_ns, sigma_ns);
641  }
642 
643  case ABOVE_THRESHOLD: {
644 
645  double tot = 0.0;
646 
647  if (npe*y1 <= threshold) {
648 
649  tot += getRiseTime(npe, threshold); // Gaussian
650  tot += getRiseTime(npe, threshold); // Gaussian
651 
652  } else if (npe <= getStartOfLinearisation()) {
653 
654  tot += getRiseTime (npe, threshold); // Gaussian
655  tot += getDecayTime(npe, threshold); // exponential
656 
657  } else {
658 
659  tot += getRiseTime (getStartOfLinearisation(), threshold); // Gaussian
660  tot += getDecayTime(getStartOfLinearisation(), threshold); // exponential
661 
662  tot += slope * (npe - getStartOfLinearisation()); // linear
663  }
664 
665  return applySaturation(tot);
666  }
667 
668  default: {
669 
670  THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::getTimeOverThreshold: Invalid charge " << npe);
671  }}
672  }
673 
674 
675  /**
676  * Get derivative of number of photo-electrons to time-over-threshold.
677  *
678  * \param npe number of photo-electrons
679  * \return dnpe/dToT [ns^-1]
680  */
681  virtual double getDerivative(const double npe) const override
682  {
683  switch(applyThreshold(npe)) {
684  case ABOVE_THRESHOLD: {
685 
686  const double z = riseTime_ns / sqrt(2.0 * log(npe/threshold));
687 
688  double y = 0.0;
689 
690  if (npe*y1 > threshold) {
691 
692  if (npe <= getStartOfLinearisation())
693  y = npe / (z + decayTime_ns); // Gaussian + exponential
694  else
695  y = 1.0 / slope; // linear
696 
697  } else {
698 
699  y = npe / (2.0 * z); // Gaussian + Gaussian
700  }
701 
702  const double tot_ns = getTimeOverThreshold(npe);
703 
704  return y / getDerivativeOfSaturation(removeSaturation(tot_ns));
705 
706  }
707 
708  default:
709  return 0.0;
710  }
711  }
712 
713 
714  /**
715  * Get probability of having a pulse with specific time-over-threshold
716  *
717  * \param tot_ns time-over-threshold with saturation [ns]
718  * \param NPE true number of photo-electrons
719  * \return probability [ns^-1]
720  */
721  double getTimeOverThresholdProbability(const double tot_ns, const int NPE) const
722  {
723 
724  const double PthBand = getIntegralOfChargeProbability(IN_THRESHOLDBAND, NPE);
725 
726  const double npe = getNPE(tot_ns);
727  const double y = getChargeProbability(npe, NPE);
728  const double v = getDerivative(npe);
729 
730  const double RthBand = PthBand * TMath::Gaus(tot_ns, mean_ns, sigma_ns, kTRUE);
731  const double RaboveTh = y * v;
732 
733  return RthBand + RaboveTh;
734  }
735 
736 
737  /**
738  * Get cumulative probability of time-over-threshold distribution
739  *
740  * \param Tmin minimum time-over-threshold (with saturation) [ns]
741  * \param Tmax maximum time-over-threshold (with saturation) [ns]
742  * \param NPE true number of photo-electrons
743  * \return probability [ns^-1]
744  */
745  double getIntegralOfTimeOverThresholdProbability(const double Tmin, const double Tmax, const int NPE) const
746  {
747 
748  const double PthBand = getIntegralOfChargeProbability(IN_THRESHOLDBAND, NPE);
749 
750  const double IthBand = PthBand * (0.5 * TMath::Erfc((Tmin - mean_ns) / sqrt(2.0) / sigma_ns) -
751  0.5 * TMath::Erfc((Tmax - mean_ns) / sqrt(2.0) / sigma_ns));
752  const double IaboveTh = getIntegralOfChargeProbability(getNPE(Tmin), getNPE(Tmax), NPE);
753 
754  return IthBand + IaboveTh;
755  }
756 
757 
758  /**
759  * Probability that a hit survives the simulation of the PMT.
760  *
761  * \param NPE number of photo-electrons
762  * \return probability
763  */
764  virtual double getSurvivalProbability(const int NPE) const override
765  {
766  if (QE <= 0.0) {
767 
768  return 0.0;
769 
770  } else if (QE <= 1.0) {
771 
772  double P = 0.0;
773 
774  const double f = gainSpread * gainSpread;
775 
776  for (int i = 1; i <= NPE; ++i) {
777 
778  const double p = (TMath::Binomial(NPE, i) *
779  TMath::Power(QE, i) * TMath::Power(1.0 - QE, NPE - i));
780 
781  double Ptotal = 0.0;
782  double Pabove = 0.0;
783  double weight = (PunderAmplified > 0.0 ? pow(PunderAmplified, i) : 1.0);
784 
785  for (int k = (PunderAmplified > 0.0 ? 0 : NPE); k <= NPE; ++k) {
786 
787  const double mu = ((i-k) * f * gain) + (k * gain);
788  const double sigma = sqrt((i-k) * f + k) * getGainSpread(1);
789 
790  Ptotal += weight * (0.5 * TMath::Erfc(( 0.0 - mu) / sqrt(2.0) / sigma));
791  Pabove += weight * (0.5 * TMath::Erfc((threshold - mu) / sqrt(2.0) / sigma));
792 
793  weight *= ( (i - k) / ((double) (k+1)) *
794  (1.0 - PunderAmplified) / PunderAmplified );
795  }
796 
797  P += p*Pabove/Ptotal;
798  }
799 
800  return P;
801 
802  } else {
803 
804  THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::getSurvivalProbability: Invalid QE " << QE);
805  }
806  }
807 
808 
809  /**
810  * Get lower threshold for rise time evaluation.
811  *
812  * \return threshold [npe]
813  */
814  static double getTH0()
815  {
816  return 0.1;
817  }
818 
819 
820  /**
821  * Get upper threshold for rise time evaluation.
822  *
823  * \return threshold [npe]
824  */
825  static double getTH1()
826  {
827  return 0.9;
828  }
829 
830  protected:
831 
832  double decayTime_ns; //!< decay time [ns]
833  double t1; //!< time at match point [ns]
834  double y1; //!< amplitude at match point [npe]
835  /**
836  * Transition point from a logarithmic to a linear relation
837  * between time-over-threshold and number of photo-electrons.\n
838  * Measurements by B. Schermer and R. Bruijn at Nikhef.
839  */
840  double x1;
841  };
842 
843 
844  /**
845  * Get time-over-threshold probability.
846  *
847  * \param pmt PMT signal processor
848  * \param tot_ns time-over-threshold [ns]
849  * \param NPE expected number of photo-electrons
850  * \param precision precision
851  * \return probability
852  */
854  const double tot_ns,
855  const double NPE,
856  const double precision = 1.0e-4)
857  {
858  int i = (int) (NPE - 5.0 * sqrt(NPE) - 0.5);
859  int M = (int) (NPE + 5.0 * sqrt(NPE) + 0.5);
860 
861  if (i < 1) { i = 1; }
862  if (M < i) { M = i; }
863 
864  double p = NPE * exp(-NPE) / (double) 1;
865 
866  for (int __i = 1; __i != i; ++__i) {
867  p *= NPE / __i;
868  }
869 
870  double P = 0.0;
871 
872  for (double p0 = 0.0; (p >= p0 || p > precision) && i != M; ++i, p0 = p, p *= NPE / (double) i) {
873  P += pmt.getTimeOverThresholdProbability(tot_ns, i) * p;
874  }
875 
876  return P;
877  }
878 }
879 
880 #endif
virtual double getRandomCharge(const int NPE) const override
Get randomised charge according to gain and gain spread.
double gain
gain [unit]
Exceptions.
then fatal No hydrophone data file $HYDROPHONE_TXT fi sort gr k
Q(UTCMax_s-UTCMin_s)-livetime_s
virtual double getRiseTime(const double npe) const override
Get time to reach threshold.
do $JPP JMEstimator M
Definition: JMEstimator.sh:37
double x1
Transition point from a logarithmic to a linear relation between time-over-threshold and number of ph...
double getDerivativeOfSaturation(const double T) const
Get derivative of saturation factor.
void setPMTParameters(const JPMTParameters &parameters)
Set PMT parameters.
double getTimeOverThresholdProbability(const JPMTAnalogueSignalProcessor &pmt, const double tot_ns, const double NPE, const double precision=1.0e-4)
Get time-over-threshold probability.
Time calibration (including definition of sign of time offset).
then JMuonPostfit f
double getIntegralOfChargeProbability(const double xmin, const double xmax, const int NPE) const
Get integral of probability.
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.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
*fatal Wrong number of arguments esac JCookie sh typeset Z DETECTOR typeset Z SOURCE_RUN typeset Z TARGET_RUN set_variable PARAMETERS_FILE $WORKDIR parameters
Definition: diff-Tuna.sh:38
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 getIntegralOfChargeProbability(JThresholdDomains domain, const int NPE) const
Get integral of probability in specific threshold domain.
double PunderAmplified
probability of underamplified hit
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable STRING $argv[2] set_array QUANTILES set_variable FORMULA *[0] exp(-0.5 *(x-[1])*(x-[1])/([2]*[2]))" set_variable MODULE `getModule -a $DETECTOR -L "$STRING 0"` typeset -Z 4 STRING JOpera1D -f hydrophone.root
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
double thresholdBand
threshold-band [npe]
double riseTime_ns
rise time of analogue pulse [ns]
static const double C
Physics constants.
virtual double getDerivative(const double npe) const override
Get derivative of number of photo-electrons to time-over-threshold.
bool slewing
time slewing of analogue signal
double mean_ns
mean time-over-threshold of threshold-band hits [ns]
virtual double getChargeProbability(const double npe, const int NPE) const override
Get probability density for given charge.
static double getMaximalRiseTime(const double th)
Get maximal rise time for given threshold.
void configure()
Configure internal parameters.
JPMTAnalogueSignalProcessor(const JPMTParameters &parameters=JPMTParameters())
Constructor.
double sigma_ns
time-over-threshold standard deviation of threshold-band hits [ns]
friend std::istream & operator>>(std::istream &in, JPMTAnalogueSignalProcessor &object)
Read PMT signal from input.
JPMTParameters()
Default constructor.
do set_variable OUTPUT_DIRECTORY $WORKDIR T
getTransitionTime
Function object to generate transition time.
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:98
then break fi done getCenter read X Y Z let X
double getAmplitude(const double t1_ns) const
Get amplitude at given time for a one photo-electron pulse.
double getGainSpread(int NPE) const
Get gain spread for given number of photo-electrons.
virtual double getNPE(const double tot_ns, const double eps=1.0e-3) const
Get number of photo-electrons.
virtual double getTimeOverThreshold(const double npe) const override
Get time-over-threshold (ToT).
virtual double getRandomTime(const double t_ns) const override
Get randomised time according transition time distribution.
double TTS_ns
transition time spread [ns]
then JCalibrateToT a
Definition: JTuneHV.sh:116
double getIntegralOfTimeOverThresholdProbability(const double Tmin, const double Tmax, const int NPE) const
Get cumulative probability of time-over-threshold distribution.
double getTimeOverThresholdProbability(const double tot_ns, const int NPE) const
Get probability of having a pulse with specific time-over-threshold.
double applySaturation(const double T) const
Get time-over-threshold with saturation.
Exception for accessing a value in a collection that is outside of its range.
Definition: JException.hh:162
double getT1() const
Get time at transition point from Gaussian to exponential.
Data structure for PMT parameters.
double saturation
saturation [ns]
double removeSaturation(const double tot_ns) const
Get time-over-threshold without saturation.
data_type v[N+1][M+1]
Definition: JPolint.hh:740
double u[N+1]
Definition: JPolint.hh:739
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:41
static double getTH1()
Get upper threshold for rise time evaluation.
double QE
relative quantum efficiency
virtual bool compare(const JPhotoElectron &first, const JPhotoElectron &second) const override
Compare (arrival times of) photo-electrons.
double getDecayTime(const double npe, const double th) const
Get time to pass from top of analogue pulse to threshold.
virtual bool applyQE() const override
Apply relative QE.
double getStartOfLinearisation() const
Get transition point from a model-dependent to linear relation between time-over-threshold and number...
source $JPP_DIR setenv csh $JPP_DIR eval JShellParser o a A
virtual double getSurvivalProbability(const int NPE) const override
Probability that a hit survives the simulation of the PMT.
virtual JThresholdDomains applyThreshold(const double npe) const override
Apply threshold.
std::vector< double > weight
Definition: JAlgorithm.hh:417
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:62
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.