#ifndef __JDETECTOR__JPMTANALOGUESIGNALPROCESSOR__ #define __JDETECTOR__JPMTANALOGUESIGNALPROCESSOR__ #include #include #include #include "TRandom3.h" #include "TMath.h" #include "JLang/JException.hh" #include "JDetector/JCalibration.hh" #include "JDetector/JPMTParameters.hh" #include "JDetector/JPMTSignalProcessorInterface.hh" #include "JDetector/JPMTTransitTimeGenerator.hh" /** * \file * * PMT analogue signal processor. * \author mdejong */ namespace JDETECTOR { using JLANG::JValueOutOfRange; /** * PMT analogue signal processor. * * This class provides for an implementation of the JDETECTOR::JPMTSignalProcessorInterface * using a specific model for the analogue pulse of the PMT.\n * In this, the leading edge of the analogue pulse from the PMT is assumed to be a Gaussian and the tail an exponential.\n * The width of the Gaussian is referred to as the rise time and * the inverse slope of the exponential to the decay time.\n * The two functions are matched at a point where the values and first derivatives are identical.\n * * Note that the decay time is related to the rise time via the specification JDETECTOR::TIME_OVER_THRESHOLD_NS. * * The charge distribution is assumed to be a Gaussian which is centered at the specified gain * and truncated by the specified threshold. * * The transition times are generated according the specified spread as follows. * - If the specified transition-time spread (TTS) is negative, * the transition times are generated according to measurements (see method JDETECTOR::getTransitionTime).\n * In this, the negated integral value of the TTS corresponds to the option * which in turn corresponds to the detector identifier of the measurements. * - If the specified TTS is positive, * the transition times are generated according a Gaussian with a sigma equals to the given TTS. * - If the TTS is zero, the transition times are generated without any spread. */ struct JPMTAnalogueSignalProcessor : public JPMTSignalProcessorInterface, public JPMTParameters { /** * Threshold domain specifiers */ enum JThresholdDomain { BELOW_THRESHOLD = -1, //!< below threshold THRESHOLDBAND = 0, //!< inside threshold band ABOVE_THRESHOLD = 2 //!< above threshold }; /** * Constructor. * * \param parameters PMT parameters */ JPMTAnalogueSignalProcessor(const JPMTParameters& parameters = JPMTParameters()) : JPMTSignalProcessorInterface(), JPMTParameters(parameters), decayTime_ns(0.0), t1(0.0), y1(0.0), x1(std::numeric_limits::max()) { configure(); } /** * Configure internal parameters. * * This method provides the implementations for * - matching of the leading edge of the analogue pulse (Gaussian) and the tail (exponential); and * - determination of number of photo-electrons above which the time-over-threshold * linearly depends on the number of photo-electrons (apart from saturation). * * Note that this method will throw an error if the value of the rise time (i.e.\ width of the Gaussian) * is too large with respect to the specification JDETECTOR::TIME_OVER_THRESHOLD_NS. */ void configure() { static const int N = 100; static const double precision = 1.0e-4; // check thresholdband if (threshold - thresholdBand < getTH0()) { THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::configure(): Invalid thresholdband [npe] " << thresholdBand); } // check rise time if (riseTime_ns > getMaximalRiseTime(threshold)) { THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::configure(): Invalid rise time [ns] " << riseTime_ns); } // decay time const double y = -log(threshold); const double a = y; const double b = riseTime_ns * sqrt(2.0*y) - TIME_OVER_THRESHOLD_NS; const double c = 0.5*riseTime_ns*riseTime_ns; const double Q = b*b - 4.0*a*c; if (Q > 0.0) decayTime_ns = (-b + sqrt(Q)) / (2.0*a); else decayTime_ns = -b / (2.0*a); // fix matching of Gaussian and exponential const double x = riseTime_ns / decayTime_ns; t1 = riseTime_ns*x; y1 = exp(-0.5*x*x); // determine transition point to linear dependence of time-over-threshold as a function of number of photo-electrons const double xs = saturation; // disable saturation saturation = 1.0e50; x1 = std::numeric_limits::max(); // disable linearisation double xmin = 1.0; double xmax = 1.0 / (getDerivative(1.0) * slope); for (int i = 0; i != N; ++i) { const double x = 0.5 * (xmin + xmax); const double u = getDerivative(x) * slope; if (fabs(1.0 - u) < precision) { break; } if (u < 1.0) xmin = x; else xmax = x; } x1 = 0.5 * (xmin + xmax); saturation = xs; // restore saturation } /** * Get decay time. * * \return decay time [ns] */ double getDecayTime() const { return decayTime_ns; } /** * Get time at transition point from Gaussian to exponential. * * \return time [ns] */ double getT1() const { return t1; } /** * Get amplitude at transition point from Gaussian to exponential. * * \return amplitude [npe] */ double getY1() const { return y1; } /** * Get transition point from a model-dependent to linear relation between time-over-threshold and number of photo-electrons. * * \return number of photo-electrons [npe] */ double getStartOfLinearisation() const { return x1; } /** * Get amplitude at given time for a one photo-electron pulse. * * \param t1_ns time [ns] * \return amplitude [npe] */ double getAmplitude(const double t1_ns) const { if (t1_ns < t1) { const double x = t1_ns / riseTime_ns; return exp(-0.5*x*x); // Gaussian } else { const double x = t1_ns / decayTime_ns; return exp(-x) / y1; // exponential } } /** * Get time to pass from threshold to top of analogue pulse.\n * In this, the leading edge of the analogue pulse is assumed to be Gaussian. * * \param npe number of photo-electrons * \param th threshold [npe] * \return time [ns] */ double getRiseTime(const double npe, const double th) const { return riseTime_ns * sqrt(2.0*log(npe/th)); // Gaussian } /** * Get time to pass from top of analogue pulse to threshold.\n * In this, the trailing edge of the analogue pulse is assumed to be exponential. * * \param npe number of photo-electrons * \param th threshold [npe] * \return time [ns] */ double getDecayTime(const double npe, const double th) const { if (npe*y1 > th) return decayTime_ns * (log(npe/th) - log(y1)); // exponential else return riseTime_ns * sqrt(2.0*log(npe/th)); // Gaussian } /** * Get maximal rise time for given threshold. * * Note that the rise time is entirely constrained by the specification JDETECTOR::TIME_OVER_THRESHOLD_NS. * * \param th threshold [npe] * \return rise time [ns] */ static double getMaximalRiseTime(const double th) { if (th > 0.0 && th < 1.0) return 0.5 * TIME_OVER_THRESHOLD_NS / sqrt(-2.0*log(th)); else THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::getMaximalRiseTime(): Invalid threshold " << th); } /** * Get time-over-threshold with saturation. * * \param tot_ns time-over-threshold without saturation * \return time-over-threshold with saturation */ double applySaturation(const double tot_ns) const { return saturation / sqrt(tot_ns*tot_ns + saturation*saturation) * tot_ns; } /** * Get time-over-threshold without saturation. * * \param tot_ns time-over-threshold with saturation * \return time-over-threshold without saturation */ double removeSaturation(const double tot_ns) const { if (tot_ns < saturation) return saturation / sqrt(saturation*saturation - tot_ns*tot_ns) * tot_ns; else return std::numeric_limits::max(); //THROW(JValueOutOfRange, "Time-over-threshold exceeds saturation " << tot_ns << " >= " << saturation); } /** * Get derivative of saturation factor. * * \param tot_ns time-over-threshold without saturation * \return derivative of saturation factor */ double getDerivativeOfSaturation(const double tot_ns) const { return saturation * saturation * saturation / ((saturation*saturation + tot_ns*tot_ns) * sqrt(saturation*saturation + tot_ns*tot_ns)); } /** * Get gain spread for given number of photo-electrons. * * \param NPE number of photo-electrons * \return gain spread */ double getGainSpread(int NPE) const { return sqrt((double) NPE * gain) * gainSpread; } /** * Get integral of probability. * * \param xmin minimum number of photo-electrons * \param xmax maximum number of photo-electrons * \param NPE true number of photo-electrons * \return probability */ double getIntegralOfChargeProbability(const double xmin, const double xmax, const int NPE) const { double zmin = xmin; double zmax = xmax; const double th = threshold - thresholdBand; const double f = gainSpread * gainSpread; if (zmin < th) { zmin = th; } if (zmax < th) { zmax = th; } double norm = 0.0; double cumulP = 0.0; double weight = (PunderAmplified > 0.0 ? pow(PunderAmplified, NPE) : 1.0); for (int k = (PunderAmplified > 0.0 ? 0 : NPE); k <= NPE; ++k) { const double mu = ((NPE-k) * f * gain) + (k * gain); const double sigma = sqrt((NPE-k) * f + k) * getGainSpread(1); norm += weight * (0.5 * TMath::Erfc((th - mu) / sqrt(2.0) / sigma)); cumulP += weight * (0.5 * TMath::Erfc((zmin - mu) / sqrt(2.0) / sigma) - 0.5 * TMath::Erfc((zmax - mu) / sqrt(2.0) / sigma)); weight *= ( (NPE - k) / ((double) (k+1)) * (1.0 - PunderAmplified) / PunderAmplified); } return cumulP / norm; } /** * Get integral of probability in specific threshold domain * * \param domain threshold domain * \param NPE true number of photo-electrons * \return probability */ double getIntegralOfChargeProbability(const JThresholdDomain domain, const int NPE) const { switch (domain) { case ABOVE_THRESHOLD: return 1.0 - getIntegralOfChargeProbability(threshold-thresholdBand, threshold, NPE); case THRESHOLDBAND: return getIntegralOfChargeProbability(threshold-thresholdBand, threshold, NPE); default: return 0.0; } } /** * Set PMT parameters. * * \param parameters PMT parameters */ void setPMTParameters(const JPMTParameters& parameters) { static_cast(*this).setPMTParameters(parameters); configure(); } /** * Read PMT signal from input. * * \param in input stream * \param object PMT signal * \return input stream */ friend std::istream& operator>>(std::istream& in, JPMTAnalogueSignalProcessor& object) { in >> static_cast(object); object.configure(); return in; } /** * Get threshold domain. * * \param npe number of photo-electrons * \return threshold domain */ JThresholdDomain getThresholdDomain(const double npe) const { if (npe > threshold) { return ABOVE_THRESHOLD; } else if (npe > threshold - thresholdBand) { return THRESHOLDBAND; } else { return BELOW_THRESHOLD; } } /** * Apply relative QE. * * \return true if accepted; false if rejected */ virtual bool applyQE() const override { if (QE <= 0.0) return false; else if (QE < 1.0) return gRandom->Rndm() < QE; else return true; } /** * Get randomised time according transition time distribution. * * \param t_ns time [ns] * \return time [ns] */ virtual double getRandomTime(const double t_ns) const override { if (TTS_ns < 0.0) return t_ns + getTransitionTime(gRandom->Rndm(), -lrint(TTS_ns)); else if (TTS_ns > 0.0) return gRandom->Gaus(t_ns, TTS_ns); else return t_ns; } /** * Compare arrival times of photo-electrons. * This implementation uses the internal rise time as two photo-electron resolution. * * Two (or more) photo-electrons are merged if they are comparable. * * \param first first photo-electron * \param second second photo-electron * \return true if arrival times of photo-electrons are within two photo-electron resolution; else false */ virtual bool compare(const JPhotoElectron& first, const JPhotoElectron& second) const override { return second.t_ns < first.t_ns + riseTime_ns; } /** * Get randomised charge according to gain and gain spread. * * \param NPE number of photo-electrons * \return number of photo-electrons */ virtual double getRandomCharge(const int NPE) const override { double q; do { // Determine which contribution to sample from int k = NPE; if (PunderAmplified > 0.0) { const double X = gRandom->Uniform(); double weight = pow(1.0 - PunderAmplified, NPE); for (double sum_p = 0.0; k > 0; --k) { sum_p += weight; if (sum_p > X) { break; } weight *= ( k / ((double) (NPE - (k-1))) * PunderAmplified / (1.0 - PunderAmplified) ); } } // Sample from chosen contribution const double f = gainSpread * gainSpread; const double mu = ((NPE-k) * f * gain) + (k * gain); const double sigma = sqrt((NPE-k) * f + k) * getGainSpread(1); q = gRandom->Gaus(mu,sigma); } while (q < 0.0); return q; } /** * Get probability density for given charge. * * \param npe observed number of photo-electrons * \param NPE true number of photo-electrons * \return probability [npe^-1] */ virtual double getChargeProbability(const double npe, const int NPE) const override { if (getThresholdDomain(npe) > BELOW_THRESHOLD) { const double f = gainSpread * gainSpread; double norm = 0.0; double prob = 0.0; double weight = (PunderAmplified > 0.0 ? pow(PunderAmplified, NPE) : 1.0); for (int k = (PunderAmplified > 0.0 ? 0 : NPE); k <= NPE; ++k) { const double mu = ((NPE-k) * f * gain) + (k * gain); const double sigma = sqrt((NPE-k) * f + k) * getGainSpread(1); norm += weight * (0.5 * TMath::Erfc(((threshold-thresholdBand) - mu) / sqrt(2.0) / sigma)); prob += weight * TMath::Gaus(npe, mu, sigma, kTRUE); weight *= ( (NPE - k) / ((double) (k+1)) * (1.0 - PunderAmplified) / PunderAmplified ); } return prob / norm; } return 0.0; } /** * Apply threshold. * * \param npe number of photo-electrons * \return true if pass; else false */ virtual bool applyThreshold(const double npe) const override { return getThresholdDomain(npe) > BELOW_THRESHOLD; } /** * Get time to reach threshold. * * Note that the rise time is defined to be zero for a one photo-electron signal. * * \param npe number of photo-electrons * \return time [ns] */ virtual double getRiseTime(const double npe) const override { if (slewing) { switch (getThresholdDomain(npe)) { case THRESHOLDBAND: return ((getRiseTime(npe, getTH0()) - getRiseTime(npe, threshold-thresholdBand)) - (getRiseTime(1.0, getTH0()) - getRiseTime(1.0, threshold-thresholdBand))) + this->mean_ns; case ABOVE_THRESHOLD: return ((getRiseTime(npe, getTH0()) - getRiseTime(npe, threshold)) - (getRiseTime(1.0, getTH0()) - getRiseTime(1.0, threshold))); default: THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::getRiseTime: Invalid charge " << npe); } } else { return 0.0; } } /** * Get time-over-threshold (ToT). * * \param npe number of photo-electrons * \return ToT [ns] */ virtual double getTimeOverThreshold(const double npe) const override { switch (getThresholdDomain(npe)) { case THRESHOLDBAND: { return gRandom->Gaus(mean_ns, sigma_ns); } case ABOVE_THRESHOLD: { double tot = 0.0; if (npe*y1 <= threshold) { tot += getRiseTime(npe, threshold); // Gaussian tot += getRiseTime(npe, threshold); // Gaussian } else if (npe <= getStartOfLinearisation()) { tot += getRiseTime (npe, threshold); // Gaussian tot += getDecayTime(npe, threshold); // exponential } else { tot += getRiseTime (getStartOfLinearisation(), threshold); // Gaussian tot += getDecayTime(getStartOfLinearisation(), threshold); // exponential tot += slope * (npe - getStartOfLinearisation()); // linear } return applySaturation(tot); } default: { THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::getTimeOverThreshold: Invalid charge " << npe); }} } /** * Get derivative of number of photo-electrons to time-over-threshold. * * \param npe number of photo-electrons * \return dnpe/dToT [ns^-1] */ virtual double getDerivative(const double npe) const override { switch (getThresholdDomain(npe)) { case ABOVE_THRESHOLD: { const double z = riseTime_ns / sqrt(2.0 * log(npe/threshold)); double y = 0.0; if (npe*y1 > threshold) { if (npe <= getStartOfLinearisation()) y = npe / (z + decayTime_ns); // Gaussian + exponential else y = 1.0 / slope; // linear } else { y = npe / (2.0 * z); // Gaussian + Gaussian } const double tot_ns = getTimeOverThreshold(npe); return y / getDerivativeOfSaturation(removeSaturation(tot_ns)); } default: return 0.0; } } /** * Probability that a hit survives the simulation of the PMT. * * \param NPE number of photo-electrons * \return probability */ virtual double getSurvivalProbability(const int NPE) const override { if (QE <= 0.0) { return 0.0; } else if (QE <= 1.0) { double P = 0.0; const double f = gainSpread * gainSpread; for (int i = 1; i <= NPE; ++i) { const double p = (TMath::Binomial(NPE, i) * TMath::Power(QE, i) * TMath::Power(1.0 - QE, NPE - i)); double Ptotal = 0.0; double Pabove = 0.0; double weight = (PunderAmplified > 0.0 ? pow(PunderAmplified, i) : 1.0); for (int k = (PunderAmplified > 0.0 ? 0 : NPE); k <= NPE; ++k) { const double mu = ((i-k) * f * gain) + (k * gain); const double sigma = sqrt((i-k) * f + k) * getGainSpread(1); Ptotal += weight * (0.5 * TMath::Erfc((0.0 - mu) / sqrt(2.0) / sigma)); Pabove += weight * (0.5 * TMath::Erfc((threshold - thresholdBand - mu) / sqrt(2.0) / sigma)); weight *= ( (i - k) / ((double) (k+1)) * (1.0 - PunderAmplified) / PunderAmplified ); } P += p*Pabove/Ptotal; } return P; } else { THROW(JValueOutOfRange, "JPMTAnalogueSignalProcessor::getSurvivalProbability: Invalid QE " << QE); } } /** * Get number of photo-electrons. * * \param tot_ns time over threshold [ns] * \return number of photo-electrons */ virtual double getNPE(const double tot_ns) const override { if (tot_ns >= saturation) { return std::numeric_limits::max(); } const double tot = removeSaturation(tot_ns); const double TOT = (getRiseTime (getStartOfLinearisation(), threshold) + getDecayTime(getStartOfLinearisation(), threshold)); if (tot <= 2*getRiseTime(threshold/y1,threshold)) { // Gaussian + Gaussian return threshold * exp(tot*tot/riseTime_ns/riseTime_ns/8.0); } else if (tot <= TOT) { // Gaussian + Exponential const double a = decayTime_ns; const double b = sqrt(2.0) * riseTime_ns; const double c = -(decayTime_ns*log(y1) + tot); const double z = (-b + sqrt(b*b - 4*a*c)) / (2*a); return threshold * exp(z*z); } else { // linear return getStartOfLinearisation() + (tot - TOT) / slope; } } /** * Get probability of having a pulse with specific time-over-threshold * * \param tot_ns time-over-threshold with saturation [ns] * \param NPE true number of photo-electrons * \return probability [ns^-1] */ double getTimeOverThresholdProbability(const double tot_ns, const int NPE) const { const double PthBand = getIntegralOfChargeProbability(THRESHOLDBAND, NPE); const double npe = getNPE(tot_ns); const double y = getChargeProbability(npe, NPE); const double v = getDerivative(npe); const double RthBand = PthBand * TMath::Gaus(tot_ns, mean_ns, sigma_ns, kTRUE); const double RaboveTh = y * v; return RthBand + RaboveTh; } /** * Get cumulative probability of time-over-threshold distribution * * \param Tmin minimum time-over-threshold (with saturation) [ns] * \param Tmax maximum time-over-threshold (with saturation) [ns] * \param NPE true number of photo-electrons * \return probability [ns^-1] */ double getIntegralOfTimeOverThresholdProbability(const double Tmin, const double Tmax, const int NPE) const { const double PthBand = getIntegralOfChargeProbability(THRESHOLDBAND, NPE); const double IthBand = PthBand * (0.5 * TMath::Erfc((Tmin - mean_ns) / sqrt(2.0) / sigma_ns) - 0.5 * TMath::Erfc((Tmax - mean_ns) / sqrt(2.0) / sigma_ns)); const double IaboveTh = getIntegralOfChargeProbability(getNPE(Tmin), getNPE(Tmax), NPE); return IthBand + IaboveTh; } /** * Get lower threshold for rise time evaluation. * * \return threshold [npe] */ static double getTH0() { return 0.1; } /** * Get upper threshold for rise time evaluation. * * \return threshold [npe] */ static double getTH1() { return 0.9; } protected: double decayTime_ns; //!< decay time [ns] double t1; //!< time at match point [ns] double y1; //!< amplitude at match point [npe] /** * Transition point from a logarithmic to a linear relation * between time-over-threshold and number of photo-electrons.\n * Measurements by B. Schermer and R. Bruijn at Nikhef. */ double x1; }; /** * Get time-over-threshold probability. * * \param pmt PMT signal processor * \param tot_ns time-over-threshold [ns] * \param NPE expected number of photo-electrons * \param precision precision * \return probability */ inline double getTimeOverThresholdProbability(const JPMTAnalogueSignalProcessor& pmt, const double tot_ns, const double NPE, const double precision = 1.0e-4) { int i = (int) (NPE - 5.0 * sqrt(NPE) - 0.5); int M = (int) (NPE + 5.0 * sqrt(NPE) + 0.5); if (i < 1) { i = 1; } if (M < i) { M = i; } double p = NPE * exp(-NPE) / (double) 1; for (int __i = 1; __i != i; ++__i) { p *= NPE / __i; } double P = 0.0; for (double p0 = 0.0; (p >= p0 || p > precision) && i != M; ++i, p0 = p, p *= NPE / (double) i) { P += pmt.getTimeOverThresholdProbability(tot_ns, i) * p; } return P; } } #endif