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