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