13 #include "cetlib_except/exception.h" 24 CLHEP::HepRandomEngine& scintTimeEngine,
25 const bool doReflectedLight,
26 const bool GeoPropTimeOnly)
27 : fGeoPropTimeOnly(GeoPropTimeOnly)
28 , fUniformGen(scintTimeEngine)
30 , fplane_depth(
std::
abs(fGeom.TPC().GetCathodeCenter().
X()))
31 , fOpDetCenter(opDetCenters())
32 , fOpDetOrientation(opDetOrientations())
35 <<
"Initializing Photon propagation time model." << std::endl;
38 mf::LogInfo(
"PropagationTimeModel") <<
"Using VUV timing parameterization";
40 fmin_d = VUVTimingParams.
get<
double>(
"min_d");
48 if (doReflectedLight) {
49 mf::LogInfo(
"PropagationTimeModel") <<
"Using VIS (reflected) timing parameterization";
53 VISTimingParams.
get<std::vector<std::vector<std::vector<double>>>>(
"Cut_off");
54 ftau_pars = VISTimingParams.
get<std::vector<std::vector<std::vector<double>>>>(
"Tau");
61 mf::LogInfo(
"PropagationTimeModel") <<
"Using geometric VUV time propagation";
65 mf::LogInfo(
"PropagationTimeModel") <<
"Photon propagation time model initalized." << std::endl;
72 const size_t OpChannel,
80 std::hypot(x0.X() - opDetCenter.X(), x0.Y() - opDetCenter.Y(), x0.Z() - opDetCenter.Z());
83 cosine =
std::abs(x0.Y() - opDetCenter.Y()) / distance;
85 cosine =
std::abs(x0.X() - opDetCenter.X()) / distance;
99 std::hypot(x0.X() - opDetCenter.X(), x0.Y() - opDetCenter.Y(), x0.Z() - opDetCenter.Z());
103 throw cet::exception(
"PropagationTimeModel") <<
"Propagation time model not found.";
110 const double distance,
111 const size_t angle_bin)
116 for (
size_t i = 0; i < arrivalTimes.size(); ++i) {
117 arrivalTimes[i] = t_prop_correction;
124 for (
size_t i = 0; i < arrivalTimes.size(); ++i) {
136 const double distance)
const 140 for (
size_t i = 0; i < arrivalTimes.size(); ++i) {
141 arrivalTimes[i] = t_prop_correction;
148 CLHEP::HepRandomEngine& scintTimeEngine)
150 mf::LogInfo(
"PropagationTimeModel") <<
"Generating VUV parameters";
151 double max_d = VUVTimingParams.
get<
double>(
"max_d");
152 const size_t num_params =
158 double dummy[1] = {1.};
163 std::vector<std::vector<double>> parameters[7];
164 parameters[0] =
std::vector(1, VUVTimingParams.
get<std::vector<double>>(
"Distances_landau"));
165 parameters[1] = VUVTimingParams.
get<std::vector<std::vector<double>>>(
"Norm_over_entries");
166 parameters[2] = VUVTimingParams.
get<std::vector<std::vector<double>>>(
"Mpv");
167 parameters[3] = VUVTimingParams.
get<std::vector<std::vector<double>>>(
"Width");
168 parameters[4] =
std::vector(1, VUVTimingParams.
get<std::vector<double>>(
"Distances_exp"));
169 parameters[5] = VUVTimingParams.
get<std::vector<std::vector<double>>>(
"Slope");
170 parameters[6] = VUVTimingParams.
get<std::vector<std::vector<double>>>(
"Expo_over_Landau_norm");
173 const double signal_t_range = 5000.;
175 for (
size_t index = 0; index < num_params; ++index) {
185 if (distance_in_cm < 2. *
fmin_d)
187 else if (distance_in_cm < 4. *
fmin_d)
192 for (
size_t angle_bin = 0; angle_bin < num_angles; ++angle_bin) {
196 std::array<double, 3> pars_landau;
199 parameters[2][angle_bin],
200 parameters[3][angle_bin],
201 parameters[1][angle_bin],
209 double pars_far[4] = {t_direct_min, pars_landau[0], pars_landau[1], pars_landau[2]};
210 VUVTiming = TF1(
"VUVTiming",
model_far, 0., signal_t_range, 4);
211 VUVTiming.SetParameters(pars_far);
218 interpolate(parameters[4][0], parameters[5][angle_bin], distance_in_cm,
true);
220 interpolate(parameters[4][0], parameters[6][angle_bin], distance_in_cm,
true);
221 pars_expo[0] *= pars_landau[2];
222 pars_expo[0] = std::log(pars_expo[0]);
224 TF1 fint = TF1(
"fint",
finter_d, pars_landau[0], 4. * t_direct_mean, 5);
225 double parsInt[5] = {
226 pars_landau[0], pars_landau[1], pars_landau[2], pars_expo[0], pars_expo[1]};
227 fint.SetParameters(parsInt);
228 double t_int = fint.GetMinimumX();
229 double minVal = fint.Eval(t_int);
231 if (minVal > 0.015) {
233 <<
"WARNING: Parametrization of VUV light discontinuous for " 235 << distance_in_cm <<
"\n" 236 <<
"This shouldn't be happening " << std::endl;
238 double parsfinal[7] = {t_int,
245 VUVTiming = TF1(
"VUVTiming",
model_close, 0., signal_t_range, 7);
246 VUVTiming.SetParameters(parsfinal);
252 double xq_max[1] = {0.975};
254 VUVTiming.SetNpx(sampling);
255 VUVTiming.GetQuantiles(1, yq_max, xq_max);
256 double max = yq_max[0];
257 double min = t_direct_min;
258 VUVTiming.SetRange(min, max);
264 auto hh = (TH1D*)VUVTiming.GetHistogram();
265 std::vector<double> vuv_timings(sampling, 0.);
266 for (
int i = 0; i < sampling; ++i)
267 vuv_timings[i] =
hh->GetBinContent(i + 1);
269 CLHEP::RandGeneral(scintTimeEngine, vuv_timings.data(), vuv_timings.size());
287 if (scintPoint.X() < 0.)
293 geo::Point_t bounce_point(plane_depth, scintPoint.Y(), scintPoint.Z());
296 double VUVdist = std::sqrt((bounce_point - scintPoint).Mag2());
297 double Visdist = std::sqrt((opDetPoint - bounce_point).Mag2());
300 int angle_bin_vuv = 0;
304 for (
size_t i = 0; i < arrivalTimes.size(); ++i) {
321 vuv_time =
fVUV_min[angle_bin_vuv][index];
323 double fastest_time = vis_time + vuv_time;
326 double cosine_theta =
std::abs(opDetPoint.X() - bounce_point.X()) / Visdist;
333 double distance_cathode_plane =
std::abs(plane_depth - scintPoint.X());
344 for (
size_t i = 0; i <
fcut_off_pars[theta_bin].size(); i++) {
353 std::vector<double> interp_vals_tau(
ftau_pars[theta_bin].
size(), 0.0);
354 for (
size_t i = 0; i <
ftau_pars[theta_bin].size(); i++) {
362 for (
size_t i = 0; i < arrivalTimes.size(); ++i) {
363 double arrival_time_smeared;
365 if (arrivalTimes[i] >= cutoff) {
continue; }
374 arrival_time_smeared = arrivalTimes[i];
381 arrival_time_smeared =
382 arrivalTimes[i] + (arrivalTimes[i] - fastest_time) * (std::pow(x, -tau) - 1);
385 }
while (arrival_time_smeared > cutoff);
387 arrivalTimes[i] = arrival_time_smeared;
397 return cathode_centre;
402 std::vector<geo::Point_t> opDetCenter;
405 opDetCenter.push_back(opDet.
GetCenter());
412 std::vector<int> opDetOrientation;
416 opDetOrientation.push_back(0);
418 else if (opDet.
isBar()) {
421 opDetOrientation.push_back(1);
424 opDetOrientation.push_back(0);
428 opDetOrientation.push_back(0);
431 return opDetOrientation;
437 double y1 = par[2] * TMath::Landau(x[0], par[0], par[1]);
438 double y2 = TMath::Exp(par[3] + x[0] * par[4]);
439 return TMath::Abs(y1 - y2);
452 double y1 = par[3] * TMath::Landau(x[0], par[1], par[2]);
453 double y2 = TMath::Exp(par[4] + x[0] * par[5]);
454 if (x[0] <= par[6] || x[0] > par[0]) y1 = 0.;
455 if (x[0] < par[0]) y2 = 0.;
466 if (x[0] <= par[0])
return 0.;
467 return par[3] * TMath::Landau(x[0], par[1], par[2]);
Point_t GetCathodeCenter() const
std::vector< geo::Point_t > opDetCenters() const
const std::vector< geo::Point_t > fOpDetCenter
OpDetGeo const & OpDetGeoFromOpDet(unsigned int OpDet) const
Returns the geo::OpDetGeo object for the given detector number.
static std::vector< geo::BoxBoundedGeo > extractActiveLArVolume(geo::GeometryCore const &geom)
Utilities related to art service access.
const std::vector< int > fOpDetOrientation
CLHEP::RandFlat fUniformGen
Encapsulate the construction of a single cyostat.
std::vector< std::vector< std::vector< double > > > fcut_off_pars
void interpolate3(std::array< double, 3 > &inter, const std::vector< double > &xData, const std::vector< double > &yData1, const std::vector< double > &yData2, const std::vector< double > &yData3, const double x, const bool extrapolate)
Float_t y1[n_points_granero]
geo::GeometryCore const & fGeom
T::provider_type const * providerFrom()
Returns a constant pointer to the provider of specified service.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
bool isBar() const
Returns whether the detector shape is a bar (TGeoBBox).
void generateVUVParams(const fhicl::ParameterSet &VUVTimingParams, CLHEP::HepRandomEngine &scintTimeEngine)
constexpr auto abs(T v)
Returns the absolute value of the argument.
static double model_close(const double *x, const double *par)
std::vector< std::vector< CLHEP::RandGeneral > > fVUVTimingGen
Float_t y2[n_points_geant4]
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
bool isSphere() const
Returns whether the detector shape is a hemisphere (TGeoSphere).
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Access the description of detector geometry.
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
const geo::Point_t fcathode_centre
std::vector< std::vector< double > > fVUV_min
double fast_acos(double xin)
T get(std::string const &key) const
const double fplane_depth
void propagationTime(std::vector< double > &arrivalTimes, const geo::Point_t &x0, const size_t OpChannel, const bool Reflected=false)
void getVUVTimesGeo(std::vector< double > &arrivalTimes, const double distance_in_cm) const
std::vector< double > fradial_distances_refl
std::vector< std::vector< std::vector< double > > > ftau_pars
Test of util::counter and support utilities.
void getVUVTimes(std::vector< double > &arrivalTimes, const double distance_in_cm, const size_t angle_bin)
geo::Point_t cathodeCentre() const
double fangle_bin_timing_vuv
unsigned int NOpDets() const
Number of OpDets in the whole detector.
std::vector< double > fdistances_refl
double fangle_bin_timing_vis
PropagationTimeModel(const fhicl::ParameterSet &VUVTimingParams, const fhicl::ParameterSet &VISTimingParams, CLHEP::HepRandomEngine &scintTimeEngine, const bool doReflectedLight=false, const bool GeoPropTimeOnly=false)
Encapsulate the geometry of an optical detector.
geo::Point_t const & GetCenter() const
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
General LArSoft Utilities.
double finflexion_point_distance
LArSoft-specific namespace.
void getVISTimes(std::vector< double > &arrivalTimes, const geo::Point_t &scintPoint, const geo::Point_t &opDetPoint)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
static double model_far(const double *x, const double *par)
const bool fGeoPropTimeOnly
std::vector< int > opDetOrientations() const
static double finter_d(const double *x, const double *par)
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, const double x, const bool extrapolate, size_t i)
Namespace collecting geometry-related classes utilities.
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
std::vector< std::vector< double > > fVUV_max