1 #ifndef __JDETECTOR__JPMTANALOGUESIGNALPROCESSOR__
2 #define __JDETECTOR__JPMTANALOGUESIGNALPROCESSOR__
67 x1(std::numeric_limits<double>::max())
86 static const int N = 100;
87 static const double precision = 1.0e-4;
109 const double Q = b*b - 4.0*a*c;
129 x1 = std::numeric_limits<double>::max();
134 for (
int i = 0; i !=
N; ++i) {
136 const double x = 0.5 * (xmin + xmax);
139 if (fabs(1.0 - u) < precision) {
149 x1 = 0.5 * (xmin + xmax);
211 return exp(-0.5*x*x);
263 if (th > 0.0 && th < 1.0)
266 THROW(
JValueOutOfRange,
"JPMTAnalogueSignalProcessor::getMaximalRiseTime(): Invalid threshold " << th);
337 const double mu = ((NPE-
k) * f *
gain) + (
k *
gain);
341 prob += weight * TMath::Gaus(npe, mu, sigma, kTRUE);
343 weight *= ( (NPE -
k) / ((
double) (
k+1)) *
370 if (zmin < th) { zmin = th; }
371 if (zmax < th) { zmax = th; }
379 const double mu = ((NPE-
k) * f *
gain) + (
k *
gain);
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));
386 weight *= ( (NPE -
k) / ((
double) (
k+1)) *
390 return cumulP / norm;
456 return gRandom->Rndm() <
QE;
473 return gRandom->Gaus(t_ns,
TTS_ns);
507 const double X = gRandom->Uniform();
510 for (
double sum_p = 0.0; k > 0; --
k) {
513 if (sum_p > X) {
break; }
515 weight *= ( k / ((double) (NPE - (k-1))) *
522 const double mu = ((NPE-
k) * f *
gain) + (k *
gain);
523 const double sigma = sqrt((NPE-k) * f + k) *
getGainSpread(1);
525 q = gRandom->Gaus(mu,sigma);
540 virtual double getNPE(
const double tot_ns,
541 const double eps = 1.0e-3)
const
557 const double z = (-B + sqrt(B*B - 4*A*C)) / (2*A);
670 THROW(
JValueOutOfRange,
"JPMTAnalogueSignalProcessor::getTimeOverThreshold: Invalid charge " << npe);
726 const double npe =
getNPE(tot_ns);
730 const double RthBand = PthBand * TMath::Gaus(tot_ns,
mean_ns,
sigma_ns, kTRUE);
731 const double RaboveTh = y *
v;
733 return RthBand + RaboveTh;
750 const double IthBand = PthBand * (0.5 * TMath::Erfc((Tmin -
mean_ns) / sqrt(2.0) /
sigma_ns) -
754 return IthBand + IaboveTh;
770 }
else if (
QE <= 1.0) {
776 for (
int i = 1; i <= NPE; ++i) {
778 const double p = (TMath::Binomial(NPE, i) *
779 TMath::Power(
QE, i) * TMath::Power(1.0 -
QE, NPE - i));
787 const double mu = ((i-
k) * f *
gain) + (
k *
gain);
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));
793 weight *= ( (i -
k) / ((
double) (
k+1)) *
797 P += p*Pabove/Ptotal;
856 const double precision = 1.0e-4)
858 int i = (int) (NPE - 5.0 * sqrt(NPE) - 0.5);
859 int M = (int) (NPE + 5.0 * sqrt(NPE) + 0.5);
861 if (i < 1) { i = 1; }
862 if (M < i) { M = i; }
864 double p = NPE *
exp(-NPE) / (double) 1;
866 for (
int __i = 1; __i != i; ++__i) {
872 for (
double p0 = 0.0; (p >= p0 || p > precision) && i != M; ++i, p0 = p, p *= NPE / (double) i) {
virtual double getRandomCharge(const int NPE) const override
Get randomised charge according to gain and gain spread.
double MAX_CHARGE
Maximum charge [npe].
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 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 ¶meters)
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 decayTime_ns
decay time [ns]
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 ¶meters)
Set PMT parameters.
JThresholdDomains
Threshold domain specifiers.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
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
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 ¶meters=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 .
then break fi done getCenter read X Y Z let X
PMT signal processor interface.
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.
double getDecayTime() const
Get decay time.
PMT analogue signal processor.
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]
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.
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.
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
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 t1
time at match point [ns]
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
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.