diff --git a/MuonSpectrometer/MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h b/MuonSpectrometer/MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h index 186509f54164c41c66954bb669a994d6ba16a162..5db0f2ace53d2377e2959f7a8cf3f635145254c6 100644 --- a/MuonSpectrometer/MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h +++ b/MuonSpectrometer/MuonDigitization/sTGC_Digitization/sTGC_Digitization/sTgcDigitMaker.h @@ -83,7 +83,6 @@ class sTgcDigitMaker : public AthMessaging { double lowEdge; // low side of the interval in ns double kParameter; double thetaParameter; - double mostProbableTime; }; /** @@ -144,6 +143,8 @@ class sTgcDigitMaker : public AthMessaging { /** Find the gamma pdf parameters of a given distance */ GammaParameter getGammaParameter(double distance) const; + /** Get the most probable time of arrival */ + double getMostProbableArrivalTime(double distance) const; /** Energy threshold value for each chamber */ double m_energyThreshold[N_STATIONNAME][N_STATIONETA][N_STATIONPHI][N_MULTIPLET][N_GASGAP][N_CHANNELTYPE]{}; @@ -166,6 +167,8 @@ class sTgcDigitMaker : public AthMessaging { // Parameters of the gamma pdf required for determining digit time std::vector<GammaParameter> m_gammaParameter; + // 4th-order polymonial describing the most probable time as function of the distance of closest approach + std::vector<double> m_mostProbableArrivalTime; // Time offset to add to Strip timing std::vector<double> m_timeOffsetStrip; diff --git a/MuonSpectrometer/MuonDigitization/sTGC_Digitization/share/sTGC_Digitization_timeArrival.dat b/MuonSpectrometer/MuonDigitization/sTGC_Digitization/share/sTGC_Digitization_timeArrival.dat index 21db68524e9ef95c4ffafa7f82e016001eba4485..598864e505bc87ffa336a77201f4eeae641e74b7 100644 --- a/MuonSpectrometer/MuonDigitization/sTGC_Digitization/share/sTGC_Digitization_timeArrival.dat +++ b/MuonSpectrometer/MuonDigitization/sTGC_Digitization/share/sTGC_Digitization_timeArrival.dat @@ -3,37 +3,47 @@ # Time of arrival is obtained from Garfield simulation. It is parametrized as # a function of the distance of closest approach, i.e. the shortest distance # from a particle's trajectory to the nearest wire. -# The approach is the following: dividing the Garfield data into 9 intervals -# of 0.1 mm (distance of closest approach), then fit the time distribution in -# each interval. +# The approach is the following: dividing the Garfield data into 11 intervals, +# then fit the time distribution in each interval. # # In each interval, the time of arrival distribution can be described by a # gamma distribution with three parameters: # - the usual k and theta parameters, and -# - t0--the shift from zero. +# - t0, which is the earliest time of arrival for the interval and also the +# only parameter that varies with the distance of closest approach. # The gamma probability distribution function is: # P(t) = 1 / (\Gamma(k) * {\theta}^k) * (t - t0)^{k-1} * exp(-(t - t0) / {\theta}) -# where t is the time of arrival to be selected and t0 is the earliest time of -# arrival allowed. +# where t is the time of arrival to be selected. # -# The current implementation in Athena use the parameters k and theta, and the -# most probable value (peak of the gamma pdf) instead of t0. First, a gamma -# random number generator is launched to choose a time, say a value val, then -# it is shifted by the best-fitted most probable value minus (k-1)*theta. The -# best-fitted most probable value (mpv) already includes t0, while the product -# (k-1)*theta gives the most probable value of the gamma pdf begining at origin. -# In other word, t0 = mpv - (k-1)*theta and the resulting time of arrival is -# given by +# The current implementation in Athena uses the parameters k and theta, and the +# most probable value of time instead of t0. First, a gamma random number +# generator is launched to choose a time, say a value val, then it is shifted +# by the best-fitted most probable value minus (k-1)*theta. The best-fitted +# most probable value (mpv) already includes t0, while the product (k-1)*theta +# gives the most probable value of a gamma pdf begining at origin. +# In other word, t0 = mpv - (k-1)*theta for the intervals considered, except for +# the 1st interval (distance less than 0.050mm) where both t0 and mpv are +# set to zero nanosecond. Also, mpv is parametrized as a 4th-order polymonial +# in order to get a continuous time of arrival as a function of distance. +# So the resulting time of arrival is given by # time = val + (mpv - (k-1) * theta). # # Below are the parameters to be read. Line starting with '#' is not read. -# keyword low_edge[mm] k theta mpv[ns] -bin 0.0 1.000 0.833 0.0 -bin 0.1 2.342 0.455 0.899 -bin 0.2 3.474 0.391 2.220 -bin 0.3 4.153 0.376 3.733 -bin 0.4 5.778 0.321 5.441 -bin 0.5 5.046 0.364 7.186 -bin 0.6 5.981 0.392 9.300 -bin 0.7 5.771 0.533 11.997 -bin 0.8 3.068 1.488 16.095 +# IMPORTANT: the intervals must be sorted in ascending order of low_edge. +# keyword low_edge[mm] k theta +bin 0.00 0.500 1.460 +bin 0.05 3.236 0.335 +bin 0.15 7.270 0.219 +bin 0.25 8.406 0.236 +bin 0.35 10.809 0.223 +bin 0.45 8.017 0.267 +bin 0.55 5.859 0.366 +bin 0.65 8.973 0.290 +bin 0.75 9.319 0.290 +bin 0.80 6.282 0.517 +bin 0.85 4.680 1.044 + +# Most probable time is described using a 4th-order polymonial +# mpv(d) = p0 + p1 * d + p2 * d^2 + p3 * d^3 + p4 * d^4 +# keyword p0 p1 p2 p3 p4 +mpv 0.106 -2.031 71.36 -128.2 87.16 diff --git a/MuonSpectrometer/MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx b/MuonSpectrometer/MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx index 0aabb85c363424da0b89d71a781d13d2a2571702..0313ba7796c7cd07c5567884f9633fc1ea27e4fc 100644 --- a/MuonSpectrometer/MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx +++ b/MuonSpectrometer/MuonDigitization/sTGC_Digitization/src/sTgcDigitMaker.cxx @@ -231,33 +231,55 @@ std::unique_ptr<sTgcDigitCollection> sTgcDigitMaker::executeDigi(const sTGCSimHi // Distance should be in the range [0, 0.9] mm, unless particle passes through // the wire plane near the edges double wire_pitch = detEl->wirePitch(); - if ((dist_wire > -9.) && (std::abs(dist_wire) > (wire_pitch / 2))) { - ATH_MSG_DEBUG("Distance to the nearest wire (" << std::abs(dist_wire) << ") is greater than expected."); + // Absolute value of the distance + double abs_dist_wire = std::abs(dist_wire); + if ((dist_wire > -9.) && (abs_dist_wire > (wire_pitch / 2))) { + ATH_MSG_DEBUG("Distance to the nearest wire (" << abs_dist_wire << ") is greater than expected."); + } + + // Do not digitize hits that are too far from the nearest wire + if (abs_dist_wire > wire_pitch) { + return nullptr; } // Get the gamma pdf parameters associated with the distance of closest approach. - GammaParameter gamma_par = getGammaParameter(std::abs(dist_wire)); + //GammaParameter gamma_par = getGammaParameter(abs_dist_wire); + double par_kappa = (getGammaParameter(abs_dist_wire)).kParameter; + double par_theta = (getGammaParameter(abs_dist_wire)).thetaParameter; + double most_prob_time = getMostProbableArrivalTime(abs_dist_wire); // Compute the most probable value of the gamma pdf - double gamma_mpv = (gamma_par.kParameter - 1) * gamma_par.thetaParameter; - // If the most probable value is near zero, then ensure it is zero - if ((gamma_par.mostProbableTime) < 0.1) {gamma_mpv = 0.;} - double t0_par = gamma_par.mostProbableTime - gamma_mpv; + double gamma_mpv = (par_kappa - 1) * par_theta; + // If the most probable value is less than zero, then set it to zero + if (gamma_mpv < 0.) {gamma_mpv = 0.;} + double t0_par = most_prob_time - gamma_mpv; // Digit time follows a gamma distribution, so a value val is - // chosen using a gamma random generator then shifted by t0 + // chosen using a gamma random generator then is shifted by t0 // to account for drift time. // Note: CLHEP::RandGamma takes the parameters k and lambda, // where lambda = 1 / theta. - double digit_time = t0_par + CLHEP::RandGamma::shoot(rndmEngine, gamma_par.kParameter, 1/gamma_par.thetaParameter); - if (digit_time < 0.0) { - // Ensure the digit time is positive - digit_time = -1.0 * digit_time; + double digit_time = t0_par + CLHEP::RandGamma::shoot(rndmEngine, par_kappa, 1/par_theta); + + // Sometimes, digit_time is negative because t0_par can be negative. + // In such case, discard the negative value and shoot RandGamma for another value. + // However, if that has already been done many times then set digit_time to zero + // in order to avoid runaway loop. + const int shoot_limit = 4; + int shoot_counter = 0; + while (digit_time < 0.) { + if (shoot_counter > shoot_limit) { + digit_time = 0.; + break; + } + digit_time = t0_par + CLHEP::RandGamma::shoot(rndmEngine, par_kappa, 1/par_theta); + ++shoot_counter; } + ATH_MSG_DEBUG("sTgcDigitMaker distance = " << dist_wire << ", time = " << digit_time - << ", k parameter = " << gamma_par.kParameter - << ", theta parameter = " << gamma_par.thetaParameter - << ", most probable time = " << gamma_par.mostProbableTime); + << ", k parameter = " << par_kappa + << ", theta parameter = " << par_theta + << ", most probable time = " << most_prob_time); //// HV efficiency correction if (m_doEfficiencyCorrection){ @@ -1167,9 +1189,12 @@ StatusCode sTgcDigitMaker::readFileOfTimeArrival() { std::istringstream iss(line); iss >> key; if (key.compare("bin") == 0) { - iss >> param.lowEdge >> param.kParameter >> param.thetaParameter >> param.mostProbableTime; + iss >> param.lowEdge >> param.kParameter >> param.thetaParameter; m_gammaParameter.push_back(param); - } + } else if (key.compare("mpv") == 0) { + double mpt; + while (iss >> mpt) {m_mostProbableArrivalTime.push_back(mpt);} + } } // Close the file @@ -1182,12 +1207,12 @@ sTgcDigitMaker::GammaParameter sTgcDigitMaker::getGammaParameter(double distance double d = distance; if (d < 0.) { - ATH_MSG_WARNING("getGammaParameter: expecting a positive distance, but got negative value: " << d - << ". Proceed to the calculation with the absolute value."); + ATH_MSG_WARNING("getGammaParameter: expecting a positive distance, but got a negative value: " << d + << ". Proceed to the calculation using its absolute value."); d = -1.0 * d; } - // Find the parameters assuming the container is sorted + // Find the parameters assuming the container is sorted in ascending order of 'lowEdge' int index{-1}; for (auto& par: m_gammaParameter) { if (distance < par.lowEdge) { @@ -1198,6 +1223,25 @@ sTgcDigitMaker::GammaParameter sTgcDigitMaker::getGammaParameter(double distance return m_gammaParameter.at(index); } +//+++++++++++++++++++++++++++++++++++++++++++++++ +double sTgcDigitMaker::getMostProbableArrivalTime(double distance) const { + + double d = distance; + if (d < 0.) { + ATH_MSG_WARNING("getMostProbableArrivalTime: expecting a positive distance, but got a negative value: " << d + << ". Proceed to the calculation using its absolute value."); + d = -1.0 * d; + } + + double mpt = m_mostProbableArrivalTime.at(0) + + m_mostProbableArrivalTime.at(1) * d + + m_mostProbableArrivalTime.at(2) * d * d + + m_mostProbableArrivalTime.at(3) * d * d * d + + m_mostProbableArrivalTime.at(4) * d * d * d * d; + return mpt; +} + +//+++++++++++++++++++++++++++++++++++++++++++++++ StatusCode sTgcDigitMaker::readFileOfTimeOffsetStrip() { // Verify the file sTGC_Digitization_timeOffsetStrip.dat exists const std::string file_name = "sTGC_Digitization_timeOffsetStrip.dat"; @@ -1242,6 +1286,7 @@ StatusCode sTgcDigitMaker::readFileOfTimeOffsetStrip() { return StatusCode::SUCCESS; } +//+++++++++++++++++++++++++++++++++++++++++++++++ double sTgcDigitMaker::getTimeOffsetStrip(int neighbor_index) const { if ((!m_timeOffsetStrip.empty()) && (neighbor_index >= 0)) { // Return the last element if out of range