1 #ifndef __JDETECTOR__JPMTANALOGUESIGNALPROCESSOR__
2 #define __JDETECTOR__JPMTANALOGUESIGNALPROCESSOR__
77 x1(std::numeric_limits<double>::max())
96 static const int N = 100;
97 static const double precision = 1.0e-4;
118 const double Q = b*b - 4.0*a*
c;
138 x1 = std::numeric_limits<double>::max();
143 for (
int i = 0; i !=
N; ++i) {
145 const double x = 0.5 * (xmin + xmax);
148 if (fabs(1.0 - u) < precision) {
158 x1 = 0.5 * (xmin + xmax);
220 return exp(-0.5*x*x);
272 if (th > 0.0 && th < 1.0)
275 THROW(
JValueOutOfRange,
"JPMTAnalogueSignalProcessor::getMaximalRiseTime(): Invalid threshold " << th);
302 return std::numeric_limits<double>::max();
347 if (zmin < th) { zmin = th; }
348 if (zmax < th) { zmax = th; }
356 const double mu = ((NPE-
k) * f *
gain) + (
k *
gain);
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));
363 weight *= ( (NPE -
k) / ((
double) (
k+1)) *
367 return cumulP / norm;
457 return gRandom->Rndm() <
QE;
474 return gRandom->Gaus(t_ns,
TTS_ns);
513 const double X = gRandom->Uniform();
516 for (
double sum_p = 0.0; k > 0; --
k) {
519 if (sum_p > X) {
break; }
521 weight *= ( k / ((double) (NPE - (k-1))) *
528 const double mu = ((NPE-
k) * f *
gain) + (k *
gain);
529 const double sigma = sqrt((NPE-k) * f + k) *
getGainSpread(1);
531 q = gRandom->Gaus(mu,sigma);
558 const double mu = ((NPE-
k) * f *
gain) + (
k *
gain);
562 prob += weight * TMath::Gaus(npe, mu, sigma, kTRUE);
564 weight *= ( (NPE -
k) / ((
double) (
k+1)) *
662 THROW(
JValueOutOfRange,
"JPMTAnalogueSignalProcessor::getTimeOverThreshold: Invalid charge " << npe);
718 }
else if (
QE <= 1.0) {
724 for (
int i = 1; i <= NPE; ++i) {
726 const double p = (TMath::Binomial(NPE, i) *
727 TMath::Power(
QE, i) * TMath::Power(1.0 -
QE, NPE - i));
735 const double mu = ((i-
k) * f *
gain) + (
k *
gain);
738 Ptotal += weight * (0.5 * TMath::Erfc((0.0 - mu) / sqrt(2.0) / sigma));
741 weight *= ( (i -
k) / ((
double) (
k+1)) *
745 P += p*Pabove/Ptotal;
763 virtual double getNPE(
const double tot_ns)
const override
765 double tot = std::numeric_limits<double>::max();
778 }
else if (tot <= TOT) {
783 const double z = (-b + sqrt(b*b - 4*a*c)) / (2*a);
805 const double npe =
getNPE(tot_ns);
809 const double RthBand = PthBand * TMath::Gaus(tot_ns,
mean_ns,
sigma_ns, kTRUE);
810 const double RaboveTh = y *
v;
812 return RthBand + RaboveTh;
829 const double IthBand = PthBand * (0.5 * TMath::Erfc((Tmin -
mean_ns) / sqrt(2.0) /
sigma_ns) -
833 return IthBand + IaboveTh;
884 const double precision = 1.0e-4)
886 int i = (int) (NPE - 5.0 * sqrt(NPE) - 0.5);
887 int M = (int) (NPE + 5.0 * sqrt(NPE) + 0.5);
889 if (i < 1) { i = 1; }
890 if (M < i) { M = i; }
892 double p = NPE *
exp(-NPE) / (double) 1;
894 for (
int __i = 1; __i != i; ++__i) {
900 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 getDerivativeOfSaturation(const double tot_ns) const
Get derivative of saturation factor.
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.
double x1
Transition point from a logarithmic to a linear relation between time-over-threshold and number of ph...
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.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
o $QUALITY_ROOT d $DEBUG!JPlot1D f
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 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 ¶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.
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
virtual double getNPE(const double tot_ns) const override
Get number of photo-electrons.
PMT signal processor interface.
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.
double getDecayTime() const
Get decay time.
PMT analogue signal processor.
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.
JThresholdDomain
Threshold domain specifiers.
$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.
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...
virtual double getSurvivalProbability(const int NPE) const override
Probability that a hit survives the simulation of the PMT.
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.