1 #ifndef __JDETECTOR__JPMTANALOGUESIGNALPROCESSOR__
2 #define __JDETECTOR__JPMTANALOGUESIGNALPROCESSOR__
63 x1(std::numeric_limits<double>::max())
82 static const int N = 100;
83 static const double precision = 1.0e-3;
105 const double Q = b*b - 4.0*a*c;
123 x1 = std::numeric_limits<double>::max();
128 for (
int i = 0; i !=
N; ++i) {
130 const double x = 0.5 * (xmin + xmax);
133 if (fabs(1.0 - u) < precision) {
143 x1 = 0.5 * (xmin + xmax);
203 return exp(-0.5*x*x);
255 if (th > 0.0 && th < 1.0)
258 THROW(
JValueOutOfRange,
"JPMTAnalogueSignalProcessor::getMaximalRiseTime(): Invalid threshold " << th);
329 const double mu = ((NPE-
k) * f *
gain) + (
k *
gain);
333 prob += weight * TMath::Gaus(npe, mu, sigma, kTRUE);
335 weight *= ( (NPE -
k) / ((
double) (
k+1)) *
362 if (zmin < th) { zmin = th; }
363 if (zmax < th) { zmax = th; }
371 const double mu = ((NPE-
k) * f *
gain) + (
k *
gain);
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));
378 weight *= ( (NPE -
k) / ((
double) (
k+1)) *
382 return cumulP / norm;
448 return gRandom->Rndm() <
QE;
465 return gRandom->Gaus(t_ns,
TTS_ns);
497 const double X = gRandom->Uniform();
500 for (
double sum_p = 0.0; k > 0; --
k) {
503 if (sum_p > X) {
break; }
505 weight *= ( k / ((double) (NPE - (k-1))) *
512 const double mu = ((NPE-
k) * f *
gain) + (k *
gain);
513 const double sigma = sqrt((NPE-k) * f + k) *
getGainSpread(1);
515 q = gRandom->Gaus(mu,sigma);
530 virtual double getNPE(
const double tot_ns,
531 const double eps = 1.0e-3)
const
547 const double z = (-B + sqrt(B*B - 4*A*C)) / (2*A);
660 THROW(
JValueOutOfRange,
"JPMTAnalogueSignalProcessor::getTimeOverThreshold: Invalid charge " << npe);
716 const double npe =
getNPE(tot_ns);
720 const double RthBand = PthBand * TMath::Gaus(tot_ns,
mean_ns,
sigma_ns, kTRUE);
721 const double RaboveTh = y *
v;
723 return RthBand + RaboveTh;
740 const double IthBand = PthBand * (0.5 * TMath::Erfc((Tmin -
mean_ns) / sqrt(2.0) /
sigma_ns) -
744 return IthBand + IaboveTh;
760 }
else if (
QE <= 1.0) {
766 for (
int i = 1; i <= NPE; ++i) {
768 const double p = (TMath::Binomial(NPE, i) *
769 TMath::Power(
QE, i) * TMath::Power(1.0 -
QE, NPE - i));
777 const double mu = ((i-
k) * f *
gain) + (
k *
gain);
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));
783 weight *= ( (i -
k) / ((
double) (
k+1)) *
787 P += p*Pabove/Ptotal;
*fatal Wrong number of arguments esac JCookie sh JRuns D $DETECTOR d sort n k
virtual double getRandomCharge(const int NPE) const
Get randomised charge according to gain and gain spread.
double MAX_CHARGE
Maximum charge [npe].
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.
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 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.
virtual double getTimeOverThreshold(const double npe) const
Get time-over-threshold (ToT).
JThresholdDomains
Threshold domain specifiers.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
esac print_variable DETECTOR INPUT_FILE OUTPUT_FILE CDF for TYPE in
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
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.
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 ¶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
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.
PMT signal processor interface.
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.
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.
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.
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
double getDecayTime(const double npe, const double th) const
Get time to pass from top of analogue pulse to threshold.
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
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.
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" -