#include "PropagationTimeModel.h"

Public Member Functions
	PropagationTimeModel (const fhicl::ParameterSet &VUVTimingParams, const fhicl::ParameterSet &VISTimingParams, CLHEP::HepRandomEngine &scintTimeEngine, const bool doReflectedLight=false, const bool GeoPropTimeOnly=false)

void	propagationTime (std::vector< double > &arrivalTimes, const geo::Point_t &x0, const size_t OpChannel, const bool Reflected=false)

Private Member Functions
void	generateVUVParams (const fhicl::ParameterSet &VUVTimingParams, CLHEP::HepRandomEngine &scintTimeEngine)

void	getVUVTimes (std::vector< double > &arrivalTimes, const double distance_in_cm, const size_t angle_bin)

void	getVUVTimesGeo (std::vector< double > &arrivalTimes, const double distance_in_cm) const

void	getVISTimes (std::vector< double > &arrivalTimes, const geo::Point_t &scintPoint, const geo::Point_t &opDetPoint)

geo::Point_t	cathodeCentre () const

std::vector< geo::Point_t >	opDetCenters () const

std::vector< int >	opDetOrientations () const

Private Attributes
const bool	fGeoPropTimeOnly

CLHEP::RandFlat	fUniformGen

geo::GeometryCore const &	fGeom

const double	fplane_depth

const geo::Point_t	fcathode_centre

const std::vector< geo::Point_t >	fOpDetCenter

const std::vector< int >	fOpDetOrientation

double	fstep_size

double	fvuv_vgroup_mean

double	fvuv_vgroup_max

double	fmin_d

double	finflexion_point_distance

double	fangle_bin_timing_vuv

std::vector< std::vector< CLHEP::RandGeneral > >	fVUVTimingGen

std::vector< std::vector< double > >	fVUV_max

std::vector< std::vector< double > >	fVUV_min

double	fvis_vmean

double	fangle_bin_timing_vis

std::vector< double >	fdistances_refl

std::vector< double >	fradial_distances_refl

std::vector< std::vector< std::vector< double > > >	fcut_off_pars

std::vector< std::vector< std::vector< double > > >	ftau_pars

Detailed Description

Definition at line 38 of file PropagationTimeModel.h.

Constructor & Destructor Documentation

phot::PropagationTimeModel::PropagationTimeModel	(	const fhicl::ParameterSet &	VUVTimingParams,
		const fhicl::ParameterSet &	VISTimingParams,
		CLHEP::HepRandomEngine &	scintTimeEngine,
		const bool	doReflectedLight = `false`,
		const bool	GeoPropTimeOnly = `false`
	)

Definition at line 61 of file PropagationTimeModel.cxx.

References fangle_bin_timing_vis, fangle_bin_timing_vuv, fcut_off_pars, fdistances_refl, fGeoPropTimeOnly, finflexion_point_distance, fmin_d, fradial_distances_refl, fstep_size, ftau_pars, fvis_vmean, fvuv_vgroup_max, fvuv_vgroup_mean, generateVUVParams(), fhicl::ParameterSet::get(), and X.

     : fGeoPropTimeOnly(GeoPropTimeOnly)
     , fUniformGen(scintTimeEngine)
     , fGeom(*(lar::providerFrom<geo::Geometry>()))
     , fplane_depth(std::abs(fGeom.TPC().GetCathodeCenter().X()))
     , fOpDetCenter(opDetCenters())
     , fOpDetOrientation(opDetOrientations())
   {
     mf::LogInfo("PropagationTimeModel")
       << "Initializing Photon propagation time model." << std::endl;
     if (!fGeoPropTimeOnly) {
       // Direct / VUV
       mf::LogInfo("PropagationTimeModel") << "Using VUV timing parameterization";
       fstep_size = VUVTimingParams.get<double>("step_size");
       fmin_d = VUVTimingParams.get<double>("min_d");
       fvuv_vgroup_mean = VUVTimingParams.get<double>("vuv_vgroup_mean");
       fvuv_vgroup_max = VUVTimingParams.get<double>("vuv_vgroup_max");
       finflexion_point_distance = VUVTimingParams.get<double>("inflexion_point_distance");
       fangle_bin_timing_vuv = VUVTimingParams.get<double>("angle_bin_timing_vuv");
       generateVUVParams(VUVTimingParams, scintTimeEngine);
 
       // Reflected / Visible
       if (doReflectedLight) {
         mf::LogInfo("PropagationTimeModel") << "Using VIS (reflected) timing parameterization";
         fdistances_refl = VISTimingParams.get<std::vector<double>>("Distances_refl");
         fradial_distances_refl = VISTimingParams.get<std::vector<double>>("Distances_radial_refl");
         fcut_off_pars =
           VISTimingParams.get<std::vector<std::vector<std::vector<double>>>>("Cut_off");
         ftau_pars = VISTimingParams.get<std::vector<std::vector<std::vector<double>>>>("Tau");
         fvis_vmean = VISTimingParams.get<double>("vis_vmean");
         fangle_bin_timing_vis = VISTimingParams.get<double>("angle_bin_timing_vis");
       }
     }
 
     if (fGeoPropTimeOnly) {
       mf::LogInfo("PropagationTimeModel") << "Using geometric VUV time propagation";
       fvuv_vgroup_mean = VUVTimingParams.get<double>("vuv_vgroup_mean");
     }
 
     mf::LogInfo("PropagationTimeModel") << "Photon propagation time model initalized." << std::endl;
   }

Member Function Documentation

geo::Point_t phot::PropagationTimeModel::cathodeCentre ( ) const

private

Definition at line 430 of file PropagationTimeModel.cxx.

References larg4::ISTPC::extractActiveLArVolume(), fGeom, geo::TPCGeo::GetCathodeCenter(), and geo::GeometryCore::TPC().

   {
     larg4::ISTPC is_tpc = larg4::ISTPC{fGeom};
     std::vector<geo::BoxBoundedGeo> activeVolumes = is_tpc.extractActiveLArVolume(fGeom);
     geo::Point_t cathode_centre = {
       fGeom.TPC().GetCathodeCenter().X(), activeVolumes[0].CenterY(), activeVolumes[0].CenterZ()};
     return cathode_centre;
   }

void phot::PropagationTimeModel::generateVUVParams	(	const fhicl::ParameterSet &	VUVTimingParams,
		CLHEP::HepRandomEngine &	scintTimeEngine
	)

private

Definition at line 186 of file PropagationTimeModel.cxx.

References fangle_bin_timing_vuv, finflexion_point_distance, fmin_d, fstep_size, fVUV_max, fVUV_min, fvuv_vgroup_max, fvuv_vgroup_mean, fVUVTimingGen, fhicl::ParameterSet::get(), hh, phot::interpolate(), phot::interpolate3(), and lar::dump::vector().

Referenced by PropagationTimeModel().

   {
     mf::LogInfo("PropagationTimeModel") << "Generating VUV parameters";
     double max_d = VUVTimingParams.get<double>("max_d");
     const size_t num_params =
       (max_d - fmin_d) / fstep_size; // for d < fmin_d, no parameterisaton, a
                                      // delta function is used instead
     const size_t num_angles = std::round(90. / fangle_bin_timing_vuv);
 
     // initialise vectors to contain range parameterisations
     double dummy[1] = {1.};
     fVUVTimingGen = std::vector(num_angles, std::vector(num_params, CLHEP::RandGeneral(dummy, 1)));
     fVUV_max = std::vector(num_angles, std::vector(num_params, 0.0));
     fVUV_min = std::vector(num_angles, std::vector(num_params, 0.0));
 
     std::vector<std::vector<double>> parameters[7];
     parameters[0] = std::vector(1, VUVTimingParams.get<std::vector<double>>("Distances_landau"));
     parameters[1] = VUVTimingParams.get<std::vector<std::vector<double>>>("Norm_over_entries");
     parameters[2] = VUVTimingParams.get<std::vector<std::vector<double>>>("Mpv");
     parameters[3] = VUVTimingParams.get<std::vector<std::vector<double>>>("Width");
     parameters[4] = std::vector(1, VUVTimingParams.get<std::vector<double>>("Distances_exp"));
     parameters[5] = VUVTimingParams.get<std::vector<std::vector<double>>>("Slope");
     parameters[6] = VUVTimingParams.get<std::vector<std::vector<double>>>("Expo_over_Landau_norm");
 
     // time range
     const double signal_t_range = 5000.;
 
     for (size_t index = 0; index < num_params; ++index) {
       double distance_in_cm = (index * fstep_size) + fmin_d;
 
       // direct path transport time
       double t_direct_mean = distance_in_cm / fvuv_vgroup_mean;
       double t_direct_min = distance_in_cm / fvuv_vgroup_max;
 
       // number of sampling points, for shorter distances, peak is
       // sharper so more sensitive sampling required
       int sampling;
       if (distance_in_cm < 2. * fmin_d)
         sampling = 10000;
       else if (distance_in_cm < 4. * fmin_d)
         sampling = 5000;
       else
         sampling = 1000;
 
       for (size_t angle_bin = 0; angle_bin < num_angles; ++angle_bin) {
         // Defining the model function(s) describing the photon
         // transportation timing vs distance.
         // Getting the landau parameters from the time parametrization
         std::array<double, 3> pars_landau;
         interpolate3(pars_landau,
                      parameters[0][0],
                      parameters[2][angle_bin],
                      parameters[3][angle_bin],
                      parameters[1][angle_bin],
                      distance_in_cm,
                      true);
         TF1 VUVTiming;
         // Deciding which time model to use (depends on the distance)
         // defining useful times for the VUV arrival time shapes
         if (distance_in_cm >= finflexion_point_distance) {
           // Set model: Landau
           double pars_far[4] = {t_direct_min, pars_landau[0], pars_landau[1], pars_landau[2]};
           VUVTiming = TF1("VUVTiming", model_far, 0., signal_t_range, 4);
           VUVTiming.SetParameters(pars_far);
         }
         else {
           // Set model: Landau + Exponential
           double pars_expo[2];
           // Getting the exponential parameters from the time parametrization
           pars_expo[1] =
             interpolate(parameters[4][0], parameters[5][angle_bin], distance_in_cm, true);
           pars_expo[0] =
             interpolate(parameters[4][0], parameters[6][angle_bin], distance_in_cm, true);
           pars_expo[0] *= pars_landau[2];
           pars_expo[0] = std::log(pars_expo[0]);
           // this is to find the intersection point between the two functions:
           TF1 fint = TF1("fint", finter_d, pars_landau[0], 4. * t_direct_mean, 5);
           double parsInt[5] = {
             pars_landau[0], pars_landau[1], pars_landau[2], pars_expo[0], pars_expo[1]};
           fint.SetParameters(parsInt);
           double t_int = fint.GetMinimumX();
           double minVal = fint.Eval(t_int);
           // the functions must intersect - output warning if they don't
           if (minVal > 0.015) {
             mf::LogWarning("PropagationTimeModel")
               << "WARNING: Parametrization of VUV light discontinuous for "
                  "distance = "
               << distance_in_cm << "\n"
               << "This shouldn't be happening " << std::endl;
           }
           double parsfinal[7] = {t_int,
                                  pars_landau[0],
                                  pars_landau[1],
                                  pars_landau[2],
                                  pars_expo[0],
                                  pars_expo[1],
                                  t_direct_min};
           VUVTiming = TF1("VUVTiming", model_close, 0., signal_t_range, 7);
           VUVTiming.SetParameters(parsfinal);
         }
 
         // store min/max, necessary to transform to the domain of the
         // original distribution
         // const size_t nq_max = 1;
         double xq_max[1] = {0.975}; // include 97.5% of tail
         double yq_max[1];
         VUVTiming.SetNpx(sampling);
         VUVTiming.GetQuantiles(1, yq_max, xq_max);
         double max = yq_max[0];
         double min = t_direct_min;
         VUVTiming.SetRange(min, max);
         fVUV_max[angle_bin][index] = max;
         fVUV_min[angle_bin][index] = min;
 
         // create the distributions that represent the parametrised timing,
         // and use those to create a RNG that samples said distributions
         auto hh = (TH1D*)VUVTiming.GetHistogram();
         std::vector<double> vuv_timings(sampling, 0.);
         for (int i = 0; i < sampling; ++i)
           vuv_timings[i] = hh->GetBinContent(i + 1);
         fVUVTimingGen[angle_bin][index] =
           CLHEP::RandGeneral(scintTimeEngine, vuv_timings.data(), vuv_timings.size());
       } // index < num_params
     }   // angle_bin < num_angles
   }

void phot::PropagationTimeModel::getVISTimes	(	std::vector< double > &	arrivalTimes,
		const geo::Point_t &	scintPoint,
		const geo::Point_t &	opDetPoint
	)

private

Definition at line 315 of file PropagationTimeModel.cxx.

References util::abs(), util::counter(), fangle_bin_timing_vis, phot::fast_acos(), fcathode_centre, fcut_off_pars, fdistances_refl, fmin_d, fplane_depth, fradial_distances_refl, fstep_size, ftau_pars, fUniformGen, fvis_vmean, fVUV_min, fvuv_vgroup_max, getVUVTimes(), phot::interpolate(), util::pi(), r, and util::size().

Referenced by propagationTime().

   {
     // ***************************************************************************
     //     Calculation of earliest arrival times and corresponding unsmeared
     //     distribution
     // ***************************************************************************
 
     // set plane_depth for correct TPC:
     double plane_depth;
     if (scintPoint.X() < 0.)
       plane_depth = -fplane_depth;
     else
       plane_depth = fplane_depth;
 
     // calculate point of reflection for shortest path
     geo::Point_t bounce_point(plane_depth, scintPoint.Y(), scintPoint.Z());
 
     // calculate distance travelled by VUV light and by vis light
     double VUVdist = std::sqrt((bounce_point - scintPoint).Mag2());
     double Visdist = std::sqrt((opDetPoint - bounce_point).Mag2());
 
     // calculate times taken by VUV part of path
     int angle_bin_vuv = 0; // on-axis by definition
     getVUVTimes(arrivalTimes, VUVdist, angle_bin_vuv);
 
     // sum parts to get total transport times times
     for (size_t i = 0; i < arrivalTimes.size(); ++i) {
       arrivalTimes[i] += Visdist / fvis_vmean;
     }
 
     // ***************************************************************************
     //      Smearing of arrival time distribution
     // ***************************************************************************
     // calculate fastest time possible
     // vis part
     double vis_time = Visdist / fvis_vmean;
     // vuv part
     double vuv_time;
     if (VUVdist < fmin_d) { vuv_time = VUVdist / fvuv_vgroup_max; }
     else {
       // find index of required parameterisation
       const size_t index = std::round((VUVdist - fmin_d) / fstep_size);
       // find shortest time
       vuv_time = fVUV_min[angle_bin_vuv][index];
     }
     double fastest_time = vis_time + vuv_time;
 
     // calculate angle theta between bound_point and optical detector
     double cosine_theta = std::abs(opDetPoint.X() - bounce_point.X()) / Visdist;
     double theta = fast_acos(cosine_theta) * 180. / CLHEP::pi;
 
     // determine smearing parameters using interpolation of generated points:
     // 1). tau = exponential smearing factor, varies with distance and angle
     // 2). cutoff = largest smeared time allowed, preventing excessively large
     //     times caused by exponential distance to cathode
     double distance_cathode_plane = std::abs(plane_depth - scintPoint.X());
     // angular bin
     size_t theta_bin = theta / fangle_bin_timing_vis;
     // radial distance from centre of TPC (y,z plane)
     double r =
       std::hypot(scintPoint.Y() - fcathode_centre.Y(), scintPoint.Z() - fcathode_centre.Z());
 
     // cut-off and tau
     // cut-off
     // interpolate in d_c for each r bin
     std::vector<double> interp_vals(fcut_off_pars[theta_bin].size(), 0.0);
     for (size_t i = 0; i < fcut_off_pars[theta_bin].size(); i++) {
       interp_vals[i] =
         interpolate(fdistances_refl, fcut_off_pars[theta_bin][i], distance_cathode_plane, true);
     }
     // interpolate in r
     double cutoff = interpolate(fradial_distances_refl, interp_vals, r, true);
 
     // tau
     // interpolate in x for each r bin
     std::vector<double> interp_vals_tau(ftau_pars[theta_bin].size(), 0.0);
     for (size_t i = 0; i < ftau_pars[theta_bin].size(); i++) {
       interp_vals_tau[i] =
         interpolate(fdistances_refl, ftau_pars[theta_bin][i], distance_cathode_plane, true);
     }
     // interpolate in r
     double tau = interpolate(fradial_distances_refl, interp_vals_tau, r, true);
 
     // apply smearing:
     for (size_t i = 0; i < arrivalTimes.size(); ++i) {
       double arrival_time_smeared;
       // if time is already greater than cutoff, do not apply smearing
       if (arrivalTimes[i] >= cutoff) { continue; }
       // otherwise smear
       else {
         unsigned int counter = 0;
         // loop until time generated is within cutoff limit
         // most are within single attempt, very few take more than two
         do {
           // don't attempt smearings too many times
           if (counter >= 10) {
             arrival_time_smeared = arrivalTimes[i]; // don't smear
             break;
           }
           else {
             // generate random number in appropriate range
             double x = fUniformGen.fire(0.5, 1.0);
             // apply the exponential smearing
             arrival_time_smeared =
               arrivalTimes[i] + (arrivalTimes[i] - fastest_time) * (std::pow(x, -tau) - 1);
           }
           counter++;
         } while (arrival_time_smeared > cutoff);
       }
       arrivalTimes[i] = arrival_time_smeared;
     }
   }

void phot::PropagationTimeModel::getVUVTimes	(	std::vector< double > &	arrivalTimes,
		const double	distance_in_cm,
		const size_t	angle_bin
	)

private

Definition at line 148 of file PropagationTimeModel.cxx.

References fmin_d, fstep_size, fVUV_max, fVUV_min, fvuv_vgroup_mean, and fVUVTimingGen.

Referenced by getVISTimes(), and propagationTime().

   {
     if (distance < fmin_d) {
       // times are fixed shift i.e. direct path only
       double t_prop_correction = distance / fvuv_vgroup_mean;
       for (size_t i = 0; i < arrivalTimes.size(); ++i) {
         arrivalTimes[i] = t_prop_correction;
       }
     }
     else {
       // determine nearest parameterisation in discretisation
       int index = std::round((distance - fmin_d) / fstep_size);
       // randomly sample parameterisation for each photon
       for (size_t i = 0; i < arrivalTimes.size(); ++i) {
         arrivalTimes[i] = fVUVTimingGen[angle_bin][index].fire() *
                             (fVUV_max[angle_bin][index] - fVUV_min[angle_bin][index]) +
                           fVUV_min[angle_bin][index];
       }
     }
   }

void phot::PropagationTimeModel::getVUVTimesGeo	(	std::vector< double > &	arrivalTimes,
		const double	distance_in_cm
	)		const

private

Definition at line 174 of file PropagationTimeModel.cxx.

References fvuv_vgroup_mean.

Referenced by propagationTime().

   {
     // times are fixed shift i.e. direct path only
     double t_prop_correction = distance / fvuv_vgroup_mean;
     for (size_t i = 0; i < arrivalTimes.size(); ++i) {
       arrivalTimes[i] = t_prop_correction;
     }
   }

std::vector< geo::Point_t > phot::PropagationTimeModel::opDetCenters ( ) const

private

Definition at line 439 of file PropagationTimeModel.cxx.

References util::counter(), fGeom, geo::OpDetGeo::GetCenter(), geo::GeometryCore::NOpDets(), and geo::GeometryCore::OpDetGeoFromOpDet().

   {
     std::vector<geo::Point_t> opDetCenter;
     for (size_t const i : util::counter(fGeom.NOpDets())) {
       geo::OpDetGeo const& opDet = fGeom.OpDetGeoFromOpDet(i);
       opDetCenter.push_back(opDet.GetCenter());
     }
     return opDetCenter;
   }

std::vector< int > phot::PropagationTimeModel::opDetOrientations ( ) const

private

Definition at line 449 of file PropagationTimeModel.cxx.

References util::counter(), fGeom, geo::OpDetGeo::Height(), geo::OpDetGeo::isBar(), geo::OpDetGeo::isSphere(), geo::GeometryCore::NOpDets(), geo::GeometryCore::OpDetGeoFromOpDet(), and geo::OpDetGeo::Width().

   {
     std::vector<int> opDetOrientation;
     for (size_t const i : util::counter(fGeom.NOpDets())) {
       geo::OpDetGeo const& opDet = fGeom.OpDetGeoFromOpDet(i);
       if (opDet.isSphere()) { // dome PMTs
         opDetOrientation.push_back(0);
       }
       else if (opDet.isBar()) {
         // determine orientation to get correction OpDet dimensions
         if (opDet.Width() > opDet.Height()) { // laterals, Y dimension smallest
           opDetOrientation.push_back(1);
         }
         else { // anode/cathode (default), X dimension smallest
           opDetOrientation.push_back(0);
         }
       }
       else {
         opDetOrientation.push_back(0);
       }
     }
     return opDetOrientation;
   }

void phot::PropagationTimeModel::propagationTime	(	std::vector< double > &	arrivalTimes,
		const geo::Point_t &	x0,
		const size_t	OpChannel,
		const bool	Reflected = `false`
	)

Definition at line 109 of file PropagationTimeModel.cxx.

References util::abs(), fangle_bin_timing_vuv, phot::fast_acos(), fGeoPropTimeOnly, fOpDetCenter, fOpDetOrientation, getVISTimes(), getVUVTimes(), getVUVTimesGeo(), and util::pi().

   {
     if (!fGeoPropTimeOnly) {
       // Get VUV photons transport time distribution from the parametrization
       geo::Point_t const& opDetCenter = fOpDetCenter[OpChannel];
       if (!Reflected) {
         double distance =
           std::hypot(x0.X() - opDetCenter.X(), x0.Y() - opDetCenter.Y(), x0.Z() - opDetCenter.Z());
         double cosine;
         if (fOpDetOrientation[OpChannel] == 1)
           cosine = std::abs(x0.Y() - opDetCenter.Y()) / distance;
         else
           cosine = std::abs(x0.X() - opDetCenter.X()) / distance;
 
         double theta = fast_acos(cosine) * 180. / CLHEP::pi;
         int angle_bin = theta / fangle_bin_timing_vuv;
         getVUVTimes(arrivalTimes, distance, angle_bin); // in ns
       }
       else {
         getVISTimes(arrivalTimes, x0, opDetCenter); // in ns
       }
     }
     else if (fGeoPropTimeOnly && !Reflected) {
       // Get VUV photons arrival time geometrically
       geo::Point_t const& opDetCenter = fOpDetCenter[OpChannel];
       double distance =
         std::hypot(x0.X() - opDetCenter.X(), x0.Y() - opDetCenter.Y(), x0.Z() - opDetCenter.Z());
       getVUVTimesGeo(arrivalTimes, distance); // in ns
     }
     else {
       throw cet::exception("PropagationTimeModel") << "Propagation time model not found.";
     }
   }