Public Member Functions | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
PropagationTimeModel Class Reference

#include <PropagationTimeModel.h>

Public Member Functions

 PropagationTimeModel (fhicl::ParameterSet VUVTimingParams, fhicl::ParameterSet VISTimingParams, CLHEP::HepRandomEngine &ScintTimeEngine, bool doReflectedLight=false, bool GeoPropTimeOnly=false)
 
void propagationTime (std::vector< double > &arrival_time_dist, geo::Point_t const &x0, const size_t OpChannel, bool Reflected=false)
 

Private Member Functions

void Initialization ()
 
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)
 
void generateParam (const size_t index, const size_t angle_bin)
 
void getVISTimes (std::vector< double > &arrivalTimes, const TVector3 &ScintPoint, const TVector3 &OpDetPoint)
 
double fast_acos (double x) const
 
double interpolate (const std::vector< double > &xData, const std::vector< double > &yData, double x, bool extrapolate, size_t i=0) const
 
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, double x, bool extrapolate)
 

Static Private Member Functions

static double finter_d (const double *x, const double *par)
 
static double model_close (const double *x, const double *par)
 
static double model_far (const double *x, const double *par)
 

Private Attributes

fhicl::ParameterSet fVUVTimingParams
 
fhicl::ParameterSet fVISTimingParams
 
bool fdoReflectedLight
 
bool fGeoPropTimeOnly
 
larg4::ISTPC fISTPC
 
CLHEP::HepRandomEngine & fScintTimeEngine
 
double fplane_depth
 
TVector3 fcathode_centre
 
std::vector< geo::BoxBoundedGeofActiveVolumes
 
size_t nOpDets
 
std::vector< geo::Point_tfOpDetCenter
 
std::vector< int > fOpDetOrientation
 
double fstep_size
 
double fmax_d
 
double fmin_d
 
double fvuv_vgroup_mean
 
double fvuv_vgroup_max
 
double finflexion_point_distance
 
double fangle_bin_timing_vuv
 
std::vector< std::vector< double > > fparameters [7]
 
std::vector< std::vector< TF1 > > VUV_timing
 
std::vector< std::vector< double > > VUV_max
 
std::vector< std::vector< double > > VUV_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 28 of file PropagationTimeModel.h.

Constructor & Destructor Documentation

PropagationTimeModel::PropagationTimeModel ( fhicl::ParameterSet  VUVTimingParams,
fhicl::ParameterSet  VISTimingParams,
CLHEP::HepRandomEngine &  ScintTimeEngine,
bool  doReflectedLight = false,
bool  GeoPropTimeOnly = false 
)

Definition at line 18 of file PropagationTimeModel.cxx.

19  : fISTPC{*(lar::providerFrom<geo::Geometry>())}
20  , fScintTimeEngine(ScintTimeEngine)
21 {
22 
23  fVUVTimingParams = VUVTimingParams;
24  fdoReflectedLight = doReflectedLight;
25  fGeoPropTimeOnly = GeoPropTimeOnly;
26 
27  if (fdoReflectedLight) fVISTimingParams = VISTimingParams;
28 
29  // initialise parameters and geometry
30  std::cout << "Photon propagation time model initalized." << std::endl;
32 }
fhicl::ParameterSet fVISTimingParams
CLHEP::HepRandomEngine & fScintTimeEngine
fhicl::ParameterSet fVUVTimingParams
QTextStream & endl(QTextStream &s)

Member Function Documentation

double PropagationTimeModel::fast_acos ( double  x) const
private

Definition at line 416 of file PropagationTimeModel.cxx.

417 {
418  double negate = double(x < 0);
419  x = std::abs(x);
420  x -= double(x > 1.0) * (x - 1.0); // <- equivalent to min(1.0,x), but faster
421  double ret = -0.0187293;
422  ret = ret * x;
423  ret = ret + 0.0742610;
424  ret = ret * x;
425  ret = ret - 0.2121144;
426  ret = ret * x;
427  ret = ret + 1.5707288;
428  ret = ret * std::sqrt(1.0 - x);
429  ret = ret - 2. * negate * ret;
430  return negate * 3.14159265358979 + ret;
431 }
T abs(T value)
list x
Definition: train.py:276
double PropagationTimeModel::finter_d ( const double *  x,
const double *  par 
)
staticprivate

Definition at line 517 of file PropagationTimeModel.cxx.

518 {
519  double y1 = par[2] * TMath::Landau(x[0], par[0], par[1]);
520  double y2 = TMath::Exp(par[3] + x[0] * par[4]);
521 
522  return TMath::Abs(y1 - y2);
523 }
list x
Definition: train.py:276
void PropagationTimeModel::generateParam ( const size_t  index,
const size_t  angle_bin 
)
private

Definition at line 205 of file PropagationTimeModel.cxx.

206 {
207  // get distance
208  double distance_in_cm = (index * fstep_size) + fmin_d;
209 
210  // time range
211  const double signal_t_range = 5000.;
212 
213  // parameterisation TF1
214  TF1 fVUVTiming;
215 
216  // direct path transport time
217  double t_direct_mean = distance_in_cm / fvuv_vgroup_mean;
218  double t_direct_min = distance_in_cm / fvuv_vgroup_max;
219 
220  // Defining the model function(s) describing the photon transportation timing vs distance
221  // Getting the landau parameters from the time parametrization
222  std::array<double, 3> pars_landau;
223  interpolate3(pars_landau,
224  fparameters[0][0],
225  fparameters[2][angle_bin],
226  fparameters[3][angle_bin],
227  fparameters[1][angle_bin],
228  distance_in_cm,
229  true);
230  // Deciding which time model to use (depends on the distance)
231  // defining useful times for the VUV arrival time shapes
232  if (distance_in_cm >= finflexion_point_distance) {
233  double pars_far[4] = {t_direct_min, pars_landau[0], pars_landau[1], pars_landau[2]};
234  // Set model: Landau
235  fVUVTiming = TF1("fVUVTiming", model_far, 0, signal_t_range, 4);
236  fVUVTiming.SetParameters(pars_far);
237  }
238  else {
239  // Set model: Landau + Exponential
240  fVUVTiming = TF1("fVUVTiming", model_close, 0, signal_t_range, 7);
241  // Exponential parameters
242  double pars_expo[2];
243  // Getting the exponential parameters from the time parametrization
244  pars_expo[1] = interpolate(fparameters[4][0], fparameters[5][angle_bin], distance_in_cm, true);
245  pars_expo[0] = interpolate(fparameters[4][0], fparameters[6][angle_bin], distance_in_cm, true);
246  pars_expo[0] *= pars_landau[2];
247  pars_expo[0] = std::log(pars_expo[0]);
248  // this is to find the intersection point between the two functions:
249  TF1 fint = TF1("fint", finter_d, pars_landau[0], 4 * t_direct_mean, 5);
250  double parsInt[5] = {
251  pars_landau[0], pars_landau[1], pars_landau[2], pars_expo[0], pars_expo[1]};
252  fint.SetParameters(parsInt);
253  double t_int = fint.GetMinimumX();
254  double minVal = fint.Eval(t_int);
255  // the functions must intersect - output warning if they don't
256  if (minVal > 0.015) {
257  std::cout << "WARNING: Parametrization of VUV light discontinuous for distance = "
258  << distance_in_cm << std::endl;
259  std::cout << "WARNING: This shouldn't be happening " << std::endl;
260  }
261  double parsfinal[7] = {t_int,
262  pars_landau[0],
263  pars_landau[1],
264  pars_landau[2],
265  pars_expo[0],
266  pars_expo[1],
267  t_direct_min};
268  fVUVTiming.SetParameters(parsfinal);
269  }
270 
271  // set the number of points used to sample parameterisation
272  // for shorter distances, peak is sharper so more sensitive sampling required
273  int fsampling;
274  if (distance_in_cm < 50) fsampling = 10000;
275  else if (distance_in_cm < 100) fsampling = 5000;
276  else fsampling = 1000;
277  fVUVTiming.SetNpx(fsampling);
278 
279  // calculate max and min distance relevant to sample parameterisation
280  // max
281  const size_t nq_max = 1;
282  double xq_max[nq_max];
283  double yq_max[nq_max];
284  xq_max[0] = 0.975; // include 97.5% of tail
285  fVUVTiming.GetQuantiles(nq_max, yq_max, xq_max);
286  double max = yq_max[0];
287  // min
288  double min = t_direct_min;
289 
290  // store TF1 and min/max, this allows identical TF1 to be used every time sampling
291  // the first call of GetRandom generates the timing sampling and stores it in the TF1 object, this is the slow part
292  // all subsequent calls check if it has been generated previously and are ~100+ times quicker
293  VUV_timing[angle_bin][index] = fVUVTiming;
294  VUV_max[angle_bin][index] = max;
295  VUV_min[angle_bin][index] = min;
296 }
std::vector< std::vector< double > > VUV_max
static double model_close(const double *x, const double *par)
static double finter_d(const double *x, const double *par)
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, double x, bool extrapolate)
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, double x, bool extrapolate, size_t i=0) const
static double model_far(const double *x, const double *par)
static int max(int a, int b)
std::vector< std::vector< TF1 > > VUV_timing
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
std::vector< std::vector< double > > fparameters[7]
std::vector< std::vector< double > > VUV_min
QTextStream & endl(QTextStream &s)
void PropagationTimeModel::getVISTimes ( std::vector< double > &  arrivalTimes,
const TVector3 &  ScintPoint,
const TVector3 &  OpDetPoint 
)
private

Definition at line 301 of file PropagationTimeModel.cxx.

304 {
305  // *************************************************************************************************
306  // Calculation of earliest arrival times and corresponding unsmeared
307  // distribution
308  // *************************************************************************************************
309 
310  // set plane_depth for correct TPC:
311  double plane_depth;
312  if (ScintPoint[0] < 0) plane_depth = -fplane_depth;
313  else plane_depth = fplane_depth;
314 
315  // calculate point of reflection for shortest path
316  TVector3 bounce_point(plane_depth,ScintPoint[1],ScintPoint[2]);
317 
318  // calculate distance travelled by VUV light and by vis light
319  double VUVdist = (bounce_point - ScintPoint).Mag();
320  double Visdist = (OpDetPoint - bounce_point).Mag();
321 
322  // calculate times taken by VUV part of path
323  int angle_bin_vuv = 0; // on-axis by definition
324  getVUVTimes(arrivalTimes, VUVdist, angle_bin_vuv);
325 
326  // sum parts to get total transport times times
327  for (size_t i = 0; i < arrivalTimes.size(); ++i) {
328  arrivalTimes[i] += Visdist / fvis_vmean;
329  }
330 
331  // *************************************************************************************************
332  // Smearing of arrival time distribution
333  // *************************************************************************************************
334  // calculate fastest time possible
335  // vis part
336  double vis_time = Visdist / fvis_vmean;
337  // vuv part
338  double vuv_time;
339  if (VUVdist < fmin_d) {
340  vuv_time = VUVdist / fvuv_vgroup_max;
341  }
342  else {
343  // find index of required parameterisation
344  const size_t index = std::round((VUVdist - fmin_d) / fstep_size);
345  // find shortest time
346  vuv_time = VUV_min[angle_bin_vuv][index];
347  }
348  // sum
349  double fastest_time = vis_time + vuv_time;
350 
351  // calculate angle theta between bound_point and optical detector
352  double cosine_theta = std::abs(OpDetPoint[0] - bounce_point[0]) / Visdist;
353  double theta = fast_acos(cosine_theta) * 180. / CLHEP::pi;
354 
355  // determine smearing parameters using interpolation of generated points:
356  // 1). tau = exponential smearing factor, varies with distance and angle
357  // 2). cutoff = largest smeared time allowed, preventing excessively large
358  // times caused by exponential distance to cathode
359  double distance_cathode_plane = std::abs(plane_depth - ScintPoint[0]);
360  // angular bin
361  size_t theta_bin = theta / fangle_bin_timing_vis;
362  // radial distance from centre of TPC (y,z plane)
363  double r = std::hypot(ScintPoint[1] - fcathode_centre[1], ScintPoint[2] - fcathode_centre[2]);
364 
365  // cut-off and tau
366  // cut-off
367  // interpolate in d_c for each r bin
368  std::vector<double> interp_vals(fcut_off_pars[theta_bin].size(), 0.0);
369  for (size_t i = 0; i < fcut_off_pars[theta_bin].size(); i++){
370  interp_vals[i] = interpolate(fdistances_refl, fcut_off_pars[theta_bin][i], distance_cathode_plane, true);
371  }
372  // interpolate in r
373  double cutoff = interpolate(fradial_distances_refl, interp_vals, r, true);
374 
375  // tau
376  // interpolate in x for each r bin
377  std::vector<double> interp_vals_tau(ftau_pars[theta_bin].size(), 0.0);
378  for (size_t i = 0; i < ftau_pars[theta_bin].size(); i++){
379  interp_vals_tau[i] = interpolate(fdistances_refl, ftau_pars[theta_bin][i], distance_cathode_plane, true);
380  }
381  // interpolate in r
382  double tau = interpolate(fradial_distances_refl, interp_vals_tau, r, true);
383 
384  // apply smearing:
385  for (size_t i = 0; i < arrivalTimes.size(); ++i) {
386  double arrival_time_smeared;
387  // if time is already greater than cutoff, do not apply smearing
388  if (arrivalTimes[i] >= cutoff) { continue; }
389  // otherwise smear
390  else {
391  unsigned int counter = 0;
392  // loop until time generated is within cutoff limit
393  // most are within single attempt, very few take more than two
394  do {
395  // don't attempt smearings too many times
396  if (counter >= 10) {
397  arrival_time_smeared = arrivalTimes[i]; // don't smear
398  break;
399  }
400  else {
401  // generate random number in appropriate range
402  double x = CLHEP::RandFlat::shoot(&fScintTimeEngine, 0.5, 1.0);
403  // apply the exponential smearing
404  arrival_time_smeared =
405  arrivalTimes[i] + (arrivalTimes[i] - fastest_time) * (std::pow(x, -tau) - 1);
406  }
407  counter++;
408  } while (arrival_time_smeared > cutoff);
409  }
410  arrivalTimes[i] = arrival_time_smeared;
411  }
412 }
CLHEP::HepRandomEngine & fScintTimeEngine
void getVUVTimes(std::vector< double > &arrivalTimes, const double distance_in_cm, const size_t angle_bin)
constexpr T pow(T x)
Definition: pow.h:72
std::enable_if_t< std::is_arithmetic_v< T >, T > hypot(T x, T y)
Definition: hypot.h:60
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, double x, bool extrapolate, size_t i=0) const
T abs(T value)
std::vector< std::vector< std::vector< double > > > ftau_pars
std::vector< double > fradial_distances_refl
double fast_acos(double x) const
std::vector< std::vector< double > > VUV_min
std::vector< std::vector< std::vector< double > > > fcut_off_pars
list x
Definition: train.py:276
float pi
Definition: units.py:11
std::vector< double > fdistances_refl
void PropagationTimeModel::getVUVTimes ( std::vector< double > &  arrivalTimes,
const double  distance_in_cm,
const size_t  angle_bin 
)
private

Definition at line 169 of file PropagationTimeModel.cxx.

170 {
171  if (distance < fmin_d) {
172  // times are fixed shift i.e. direct path only
173  double t_prop_correction = distance / fvuv_vgroup_mean;
174  for (size_t i = 0; i < arrivalTimes.size(); ++i) {
175  arrivalTimes[i] = t_prop_correction;
176  }
177  }
178  else {
179  // determine nearest parameterisation in discretisation
180  int index = std::round((distance - fmin_d) / fstep_size);
181  // check whether required parameterisation has been generated, generating if not
182  if (VUV_timing[angle_bin][index].GetNdim() == 0) generateParam(index, angle_bin);
183  // randomly sample parameterisation for each photon
184  for (size_t i = 0; i < arrivalTimes.size(); ++i) {
185  arrivalTimes[i] = VUV_timing[angle_bin][index].GetRandom(VUV_min[angle_bin][index], VUV_max[angle_bin][index]);
186  }
187  }
188 }
std::vector< std::vector< double > > VUV_max
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
std::vector< std::vector< TF1 > > VUV_timing
std::vector< std::vector< double > > VUV_min
void generateParam(const size_t index, const size_t angle_bin)
void PropagationTimeModel::getVUVTimesGeo ( std::vector< double > &  arrivalTimes,
const double  distance_in_cm 
)
private

Definition at line 193 of file PropagationTimeModel.cxx.

194 {
195  // times are fixed shift i.e. direct path only
196  double t_prop_correction = distance / fvuv_vgroup_mean;
197  for (size_t i = 0; i < arrivalTimes.size(); ++i) {
198  arrivalTimes[i] = t_prop_correction;
199  }
200 }
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
void PropagationTimeModel::Initialization ( )
private

Definition at line 36 of file PropagationTimeModel.cxx.

37 {
38  // load photon propagation time model parameters
39  if (!fGeoPropTimeOnly) {
40 
41  // Direct / VUV
42  mf::LogInfo("PropagationTimeModel") << "Using VUV timing parameterization";
43 
44  fparameters[0] = std::vector(1, fVUVTimingParams.get<std::vector<double>>("Distances_landau"));
45  fparameters[1] = fVUVTimingParams.get<std::vector<std::vector<double>>>("Norm_over_entries");
46  fparameters[2] = fVUVTimingParams.get<std::vector<std::vector<double>>>("Mpv");
47  fparameters[3] = fVUVTimingParams.get<std::vector<std::vector<double>>>("Width");
48  fparameters[4] = std::vector(1, fVUVTimingParams.get<std::vector<double>>("Distances_exp"));
49  fparameters[5] = fVUVTimingParams.get<std::vector<std::vector<double>>>("Slope");
50  fparameters[6] = fVUVTimingParams.get<std::vector<std::vector<double>>>("Expo_over_Landau_norm");
51 
52  fstep_size = fVUVTimingParams.get<double>("step_size");
53  fmax_d = fVUVTimingParams.get<double>("max_d");
54  fmin_d = fVUVTimingParams.get<double>("min_d");
55  fvuv_vgroup_mean = fVUVTimingParams.get<double>("vuv_vgroup_mean");
56  fvuv_vgroup_max = fVUVTimingParams.get<double>("vuv_vgroup_max");
57  finflexion_point_distance = fVUVTimingParams.get<double>("inflexion_point_distance");
58  fangle_bin_timing_vuv = fVUVTimingParams.get<double>("angle_bin_timing_vuv");
59 
60  // create vector of empty TF1s that will be replaced with the sampled
61  // parameterisations that are generated as they are required
62  const size_t num_params = (fmax_d - fmin_d) / fstep_size; // for d < fmin_d, no parameterisaton, a delta function is used instead
63  size_t num_angles = std::round(90/fangle_bin_timing_vuv);
64  VUV_timing = std::vector(num_angles, std::vector(num_params, TF1()));
65 
66  // initialise vectors to contain range parameterisations sampled to in each case
67  // when using TF1->GetRandom(xmin,xmax), must be in same range otherwise sampling
68  // is regenerated, this is the slow part!
69  VUV_max = std::vector(num_angles, std::vector(num_params, 0.0));
70  VUV_min = std::vector(num_angles, std::vector(num_params, 0.0));
71 
72 
73  // Reflected / Visible
74  if (fdoReflectedLight) {
75  mf::LogInfo("PropagationTimeModel") << "Using VIS (reflected) timing parameterization";
76 
77  fdistances_refl = fVISTimingParams.get<std::vector<double>>("Distances_refl");
78  fradial_distances_refl = fVISTimingParams.get<std::vector<double>>("Distances_radial_refl");
79  fcut_off_pars = fVISTimingParams.get<std::vector<std::vector<std::vector<double>>>>("Cut_off");
80  ftau_pars = fVISTimingParams.get<std::vector<std::vector<std::vector<double>>>>("Tau");
81  fvis_vmean = fVISTimingParams.get<double>("vis_vmean");
82  fangle_bin_timing_vis = fVISTimingParams.get<double>("angle_bin_timing_vis");
83  }
84  }
85 
86  if (fGeoPropTimeOnly) {
87  mf::LogInfo("PropagationTimeModel") << "Using geometric VUV time propagation";
88  fvuv_vgroup_mean = fVUVTimingParams.get<double>("vuv_vgroup_mean");
89  }
90 
91  // access information from the geometry service
92  geo::GeometryCore const& geom = *(lar::providerFrom<geo::Geometry>());
93 
94  // get TPC information
95  fplane_depth = std::abs(geom.TPC(0, 0).GetCathodeCenter().X());
97  fcathode_centre = {geom.TPC(0, 0).GetCathodeCenter().X(),
98  fActiveVolumes[0].CenterY(),
99  fActiveVolumes[0].CenterZ()};
100 
101  // get PDS information
102  nOpDets = geom.NOpDets();
103 
104  fOpDetOrientation.reserve(nOpDets); fOpDetCenter.reserve(nOpDets);
105  for (size_t const i : util::counter(nOpDets)) {
106 
107  geo::OpDetGeo const& opDet = geom.OpDetGeoFromOpDet(i);
108  fOpDetCenter.push_back(opDet.GetCenter());
109 
110  if (opDet.isSphere()) { // dome PMTs
111  fOpDetOrientation.push_back(0);
112  }
113  else if (opDet.isBar()) {
114  // determine orientation to get correction OpDet dimensions
115  if (opDet.Width() > opDet.Height()) { // laterals, Y dimension smallest
116  fOpDetOrientation.push_back(1);
117  }
118  else { // anode/cathode (default), X dimension smallest
119  fOpDetOrientation.push_back(0);
120  }
121  }
122  else {
123  fOpDetOrientation.push_back(0);
124  }
125  }
126 }
static std::vector< geo::BoxBoundedGeo > extractActiveLArVolume(geo::GeometryCore const &geom)
Definition: ISTPC.cxx:49
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
fhicl::ParameterSet fVISTimingParams
bool isBar() const
Returns whether the detector shape is a bar (TGeoBBox).
Definition: OpDetGeo.h:195
struct vector vector
std::vector< std::vector< double > > VUV_max
void GetCenter(double *xyz, double localz=0.0) const
Definition: OpDetGeo.cxx:40
std::vector< geo::Point_t > fOpDetCenter
bool isSphere() const
Returns whether the detector shape is a hemisphere (TGeoSphere).
Definition: OpDetGeo.h:198
T abs(T value)
std::vector< std::vector< std::vector< double > > > ftau_pars
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
Definition: counter.h:285
std::vector< double > fradial_distances_refl
std::vector< geo::BoxBoundedGeo > fActiveVolumes
T get(std::string const &key) const
Definition: ParameterSet.h:271
std::vector< int > fOpDetOrientation
Description of geometry of one entire detector.
std::vector< std::vector< TF1 > > VUV_timing
fhicl::ParameterSet fVUVTimingParams
std::vector< std::vector< double > > fparameters[7]
std::vector< std::vector< double > > VUV_min
std::vector< std::vector< std::vector< double > > > fcut_off_pars
std::vector< double > fdistances_refl
double Width() const
Definition: OpDetGeo.h:82
double Height() const
Definition: OpDetGeo.h:83
double PropagationTimeModel::interpolate ( const std::vector< double > &  xData,
const std::vector< double > &  yData,
double  x,
bool  extrapolate,
size_t  i = 0 
) const
private

Definition at line 439 of file PropagationTimeModel.cxx.

444 {
445  if (i == 0) {
446  size_t size = xData.size();
447  if (x >= xData[size - 2]) { // special case: beyond right end
448  i = size - 2;
449  }
450  else {
451  while (x > xData[i + 1])
452  i++;
453  }
454  }
455  double xL = xData[i];
456  double xR = xData[i + 1];
457  double yL = yData[i];
458  double yR = yData[i + 1]; // points on either side (unless beyond ends)
459  if (!extrapolate) { // if beyond ends of array and not extrapolating
460  if (x < xL) return yL;
461  if (x > xR) return yL;
462  }
463  const double dydx = (yR - yL) / (xR - xL); // gradient
464  return yL + dydx * (x - xL); // linear interpolation
465 }
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
list x
Definition: train.py:276
void PropagationTimeModel::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,
double  x,
bool  extrapolate 
)
private

Definition at line 469 of file PropagationTimeModel.cxx.

476 {
477  size_t size = xData.size();
478  size_t i = 0; // find left end of interval for interpolation
479  if (x >= xData[size - 2]) { // special case: beyond right end
480  i = size - 2;
481  }
482  else {
483  while (x > xData[i + 1])
484  i++;
485  }
486  double xL = xData[i];
487  double xR = xData[i + 1]; // points on either side (unless beyond ends)
488  double yL1 = yData1[i];
489  double yR1 = yData1[i + 1];
490  double yL2 = yData2[i];
491  double yR2 = yData2[i + 1];
492  double yL3 = yData3[i];
493  double yR3 = yData3[i + 1];
494 
495  if (!extrapolate) { // if beyond ends of array and not extrapolating
496  if (x < xL) {
497  inter[0] = yL1;
498  inter[1] = yL2;
499  inter[2] = yL3;
500  return;
501  }
502  if (x > xR) {
503  inter[0] = yL1;
504  inter[1] = yL2;
505  inter[2] = yL3;
506  return;
507  }
508  }
509  const double m = (x - xL) / (xR - xL);
510  inter[0] = m * (yR1 - yL1) + yL1;
511  inter[1] = m * (yR2 - yL2) + yL2;
512  inter[2] = m * (yR3 - yL3) + yL3;
513 }
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
list x
Definition: train.py:276
double PropagationTimeModel::model_close ( const double *  x,
const double *  par 
)
staticprivate

Definition at line 527 of file PropagationTimeModel.cxx.

528 {
529  // par0 = joining point
530  // par1 = Landau MPV
531  // par2 = Landau width
532  // par3 = normalization
533  // par4 = Expo cte
534  // par5 = Expo tau
535  // par6 = t_min
536 
537  double y1 = par[3] * TMath::Landau(x[0], par[1], par[2]);
538  double y2 = TMath::Exp(par[4] + x[0] * par[5]);
539  if (x[0] <= par[6] || x[0] > par[0]) y1 = 0.;
540  if (x[0] < par[0]) y2 = 0.;
541 
542  return (y1 + y2);
543 }
list x
Definition: train.py:276
double PropagationTimeModel::model_far ( const double *  x,
const double *  par 
)
staticprivate

Definition at line 547 of file PropagationTimeModel.cxx.

548 {
549  // par1 = Landau MPV
550  // par2 = Landau width
551  // par3 = normalization
552  // par0 = t_min
553 
554  double y = par[3] * TMath::Landau(x[0], par[1], par[2]);
555  if (x[0] <= par[0]) y = 0.;
556 
557  return y;
558 }
list x
Definition: train.py:276
void PropagationTimeModel::propagationTime ( std::vector< double > &  arrival_time_dist,
geo::Point_t const &  x0,
const size_t  OpChannel,
bool  Reflected = false 
)

Definition at line 131 of file PropagationTimeModel.cxx.

135 {
136  if (!fGeoPropTimeOnly) {
137  // Get VUV photons transport time distribution from the parametrization
138  geo::Point_t const& opDetCenter = fOpDetCenter[OpChannel];
139  if (!Reflected) {
140  double distance = std::hypot(x0.X() - opDetCenter.X(), x0.Y() - opDetCenter.Y(), x0.Z() - opDetCenter.Z());
141  double cosine;
142  if (fOpDetOrientation[OpChannel] == 1) cosine = std::abs(x0.Y() - opDetCenter.Y()) / distance;
143  else cosine = std::abs(x0.X() - opDetCenter.X()) / distance;
144 
145  double theta = fast_acos(cosine)*180./CLHEP::pi;
146  int angle_bin = theta/fangle_bin_timing_vuv;
147  getVUVTimes(arrival_time_dist, distance, angle_bin); // in ns
148  }
149  else {
150  getVISTimes(arrival_time_dist, geo::vect::toTVector3(x0),
151  geo::vect::toTVector3(opDetCenter)); // in ns
152  }
153  }
154  else if (fGeoPropTimeOnly && !Reflected) {
155  // Get VUV photons arrival time geometrically
156  geo::Point_t const& opDetCenter = fOpDetCenter[OpChannel];
157  double distance = std::hypot(x0.X() - opDetCenter.X(), x0.Y() - opDetCenter.Y(), x0.Z() - opDetCenter.Z());
158  getVUVTimesGeo(arrival_time_dist, distance); // in ns
159  }
160  else {
161  throw cet::exception("PropagationTimeModel")
162  << "Propagation time model not found.";
163  }
164 }
Index OpChannel(Index detNum, Index channel)
void getVUVTimes(std::vector< double > &arrivalTimes, const double distance_in_cm, const size_t angle_bin)
void getVISTimes(std::vector< double > &arrivalTimes, const TVector3 &ScintPoint, const TVector3 &OpDetPoint)
std::vector< geo::Point_t > fOpDetCenter
std::enable_if_t< std::is_arithmetic_v< T >, T > hypot(T x, T y)
Definition: hypot.h:60
T abs(T value)
double fast_acos(double x) const
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
std::vector< int > fOpDetOrientation
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:184
void getVUVTimesGeo(std::vector< double > &arrivalTimes, const double distance_in_cm)
TVector3 toTVector3(Vector const &v)
Converts a vector into a TVector3.
float pi
Definition: units.py:11
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

Member Data Documentation

std::vector<geo::BoxBoundedGeo> PropagationTimeModel::fActiveVolumes
private

Definition at line 106 of file PropagationTimeModel.h.

double PropagationTimeModel::fangle_bin_timing_vis
private

Definition at line 123 of file PropagationTimeModel.h.

double PropagationTimeModel::fangle_bin_timing_vuv
private

Definition at line 114 of file PropagationTimeModel.h.

TVector3 PropagationTimeModel::fcathode_centre
private

Definition at line 105 of file PropagationTimeModel.h.

std::vector<std::vector<std::vector<double> > > PropagationTimeModel::fcut_off_pars
private

Definition at line 126 of file PropagationTimeModel.h.

std::vector<double> PropagationTimeModel::fdistances_refl
private

Definition at line 124 of file PropagationTimeModel.h.

bool PropagationTimeModel::fdoReflectedLight
private

Definition at line 94 of file PropagationTimeModel.h.

bool PropagationTimeModel::fGeoPropTimeOnly
private

Definition at line 95 of file PropagationTimeModel.h.

double PropagationTimeModel::finflexion_point_distance
private

Definition at line 114 of file PropagationTimeModel.h.

larg4::ISTPC PropagationTimeModel::fISTPC
private

Definition at line 98 of file PropagationTimeModel.h.

double PropagationTimeModel::fmax_d
private

Definition at line 114 of file PropagationTimeModel.h.

double PropagationTimeModel::fmin_d
private

Definition at line 114 of file PropagationTimeModel.h.

std::vector<geo::Point_t> PropagationTimeModel::fOpDetCenter
private

Definition at line 110 of file PropagationTimeModel.h.

std::vector<int> PropagationTimeModel::fOpDetOrientation
private

Definition at line 111 of file PropagationTimeModel.h.

std::vector<std::vector<double> > PropagationTimeModel::fparameters[7]
private

Definition at line 115 of file PropagationTimeModel.h.

double PropagationTimeModel::fplane_depth
private

Definition at line 104 of file PropagationTimeModel.h.

std::vector<double> PropagationTimeModel::fradial_distances_refl
private

Definition at line 125 of file PropagationTimeModel.h.

CLHEP::HepRandomEngine& PropagationTimeModel::fScintTimeEngine
private

Definition at line 101 of file PropagationTimeModel.h.

double PropagationTimeModel::fstep_size
private

Definition at line 114 of file PropagationTimeModel.h.

std::vector<std::vector<std::vector<double> > > PropagationTimeModel::ftau_pars
private

Definition at line 127 of file PropagationTimeModel.h.

double PropagationTimeModel::fvis_vmean
private

Definition at line 123 of file PropagationTimeModel.h.

fhicl::ParameterSet PropagationTimeModel::fVISTimingParams
private

Definition at line 91 of file PropagationTimeModel.h.

double PropagationTimeModel::fvuv_vgroup_max
private

Definition at line 114 of file PropagationTimeModel.h.

double PropagationTimeModel::fvuv_vgroup_mean
private

Definition at line 114 of file PropagationTimeModel.h.

fhicl::ParameterSet PropagationTimeModel::fVUVTimingParams
private

Definition at line 90 of file PropagationTimeModel.h.

size_t PropagationTimeModel::nOpDets
private

Definition at line 109 of file PropagationTimeModel.h.

std::vector<std::vector<double> > PropagationTimeModel::VUV_max
private

Definition at line 119 of file PropagationTimeModel.h.

std::vector<std::vector<double> > PropagationTimeModel::VUV_min
private

Definition at line 120 of file PropagationTimeModel.h.

std::vector<std::vector<TF1> > PropagationTimeModel::VUV_timing
private

Definition at line 117 of file PropagationTimeModel.h.


The documentation for this class was generated from the following files: