LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
phot::PropagationTimeModel Class Reference

#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_topDetCenters () const
 
std::vector< int > opDetOrientations () const
 

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

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_tfOpDetCenter
 
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 37 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 22 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(), and fhicl::ParameterSet::get().

27  : fGeoPropTimeOnly(GeoPropTimeOnly)
28  , fUniformGen(scintTimeEngine)
29  , fGeom(*(lar::providerFrom<geo::Geometry>()))
33  {
34  mf::LogInfo("PropagationTimeModel")
35  << "Initializing Photon propagation time model." << std::endl;
36  if (!fGeoPropTimeOnly) {
37  // Direct / VUV
38  mf::LogInfo("PropagationTimeModel") << "Using VUV timing parameterization";
39  fstep_size = VUVTimingParams.get<double>("step_size");
40  fmin_d = VUVTimingParams.get<double>("min_d");
41  fvuv_vgroup_mean = VUVTimingParams.get<double>("vuv_vgroup_mean");
42  fvuv_vgroup_max = VUVTimingParams.get<double>("vuv_vgroup_max");
43  finflexion_point_distance = VUVTimingParams.get<double>("inflexion_point_distance");
44  fangle_bin_timing_vuv = VUVTimingParams.get<double>("angle_bin_timing_vuv");
45  generateVUVParams(VUVTimingParams, scintTimeEngine);
46 
47  // Reflected / Visible
48  if (doReflectedLight) {
49  mf::LogInfo("PropagationTimeModel") << "Using VIS (reflected) timing parameterization";
50  fdistances_refl = VISTimingParams.get<std::vector<double>>("Distances_refl");
51  fradial_distances_refl = VISTimingParams.get<std::vector<double>>("Distances_radial_refl");
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");
55  fvis_vmean = VISTimingParams.get<double>("vis_vmean");
56  fangle_bin_timing_vis = VISTimingParams.get<double>("angle_bin_timing_vis");
57  }
58  }
59 
60  if (fGeoPropTimeOnly) {
61  mf::LogInfo("PropagationTimeModel") << "Using geometric VUV time propagation";
62  fvuv_vgroup_mean = VUVTimingParams.get<double>("vuv_vgroup_mean");
63  }
64 
65  mf::LogInfo("PropagationTimeModel") << "Photon propagation time model initalized." << std::endl;
66  }
Point_t GetCathodeCenter() const
Definition: TPCGeo.h:254
std::vector< geo::Point_t > opDetCenters() const
const std::vector< geo::Point_t > fOpDetCenter
const std::vector< int > fOpDetOrientation
std::vector< std::vector< std::vector< double > > > fcut_off_pars
geo::GeometryCore const & fGeom
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void generateVUVParams(const fhicl::ParameterSet &VUVTimingParams, CLHEP::HepRandomEngine &scintTimeEngine)
constexpr auto abs(T v)
Returns the absolute value of the argument.
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
T get(std::string const &key) const
Definition: ParameterSet.h:314
std::vector< double > fradial_distances_refl
std::vector< std::vector< std::vector< double > > > ftau_pars
std::vector< double > fdistances_refl
std::vector< int > opDetOrientations() const

Member Function Documentation

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

Definition at line 391 of file PropagationTimeModel.cxx.

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

392  {
393  larg4::ISTPC is_tpc = larg4::ISTPC{fGeom};
394  std::vector<geo::BoxBoundedGeo> activeVolumes = is_tpc.extractActiveLArVolume(fGeom);
395  geo::Point_t cathode_centre = {
396  fGeom.TPC().GetCathodeCenter().X(), activeVolumes[0].CenterY(), activeVolumes[0].CenterZ()};
397  return cathode_centre;
398  }
Point_t GetCathodeCenter() const
Definition: TPCGeo.h:254
static std::vector< geo::BoxBoundedGeo > extractActiveLArVolume(geo::GeometryCore const &geom)
Definition: ISTPC.cxx:45
geo::GeometryCore const & fGeom
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
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:180
double phot::PropagationTimeModel::finter_d ( const double *  x,
const double *  par 
)
staticprivate

Definition at line 435 of file PropagationTimeModel.cxx.

References y1, and y2.

Referenced by generateVUVParams().

436  {
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);
440  }
Float_t x
Definition: compare.C:6
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t y2[n_points_geant4]
Definition: compare.C:26
void phot::PropagationTimeModel::generateVUVParams ( const fhicl::ParameterSet VUVTimingParams,
CLHEP::HepRandomEngine &  scintTimeEngine 
)
private

Definition at line 147 of file PropagationTimeModel.cxx.

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

Referenced by PropagationTimeModel().

149  {
150  mf::LogInfo("PropagationTimeModel") << "Generating VUV parameters";
151  double max_d = VUVTimingParams.get<double>("max_d");
152  const size_t num_params =
153  (max_d - fmin_d) / fstep_size; // for d < fmin_d, no parameterisaton, a
154  // delta function is used instead
155  const size_t num_angles = std::round(90. / fangle_bin_timing_vuv);
156 
157  // initialise vectors to contain range parameterisations
158  double dummy[1] = {1.};
159  fVUVTimingGen = std::vector(num_angles, std::vector(num_params, CLHEP::RandGeneral(dummy, 1)));
160  fVUV_max = std::vector(num_angles, std::vector(num_params, 0.0));
161  fVUV_min = std::vector(num_angles, std::vector(num_params, 0.0));
162 
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");
171 
172  // time range
173  const double signal_t_range = 5000.;
174 
175  for (size_t index = 0; index < num_params; ++index) {
176  double distance_in_cm = (index * fstep_size) + fmin_d;
177 
178  // direct path transport time
179  double t_direct_mean = distance_in_cm / fvuv_vgroup_mean;
180  double t_direct_min = distance_in_cm / fvuv_vgroup_max;
181 
182  // number of sampling points, for shorter distances, peak is
183  // sharper so more sensitive sampling required
184  int sampling;
185  if (distance_in_cm < 2. * fmin_d)
186  sampling = 10000;
187  else if (distance_in_cm < 4. * fmin_d)
188  sampling = 5000;
189  else
190  sampling = 1000;
191 
192  for (size_t angle_bin = 0; angle_bin < num_angles; ++angle_bin) {
193  // Defining the model function(s) describing the photon
194  // transportation timing vs distance.
195  // Getting the landau parameters from the time parametrization
196  std::array<double, 3> pars_landau;
197  interpolate3(pars_landau,
198  parameters[0][0],
199  parameters[2][angle_bin],
200  parameters[3][angle_bin],
201  parameters[1][angle_bin],
202  distance_in_cm,
203  true);
204  TF1 VUVTiming;
205  // Deciding which time model to use (depends on the distance)
206  // defining useful times for the VUV arrival time shapes
207  if (distance_in_cm >= finflexion_point_distance) {
208  // Set model: Landau
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);
212  }
213  else {
214  // Set model: Landau + Exponential
215  double pars_expo[2];
216  // Getting the exponential parameters from the time parametrization
217  pars_expo[1] =
218  interpolate(parameters[4][0], parameters[5][angle_bin], distance_in_cm, true);
219  pars_expo[0] =
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]);
223  // this is to find the intersection point between the two functions:
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);
230  // the functions must intersect - output warning if they don't
231  if (minVal > 0.015) {
232  mf::LogWarning("PropagationTimeModel")
233  << "WARNING: Parametrization of VUV light discontinuous for "
234  "distance = "
235  << distance_in_cm << "\n"
236  << "This shouldn't be happening " << std::endl;
237  }
238  double parsfinal[7] = {t_int,
239  pars_landau[0],
240  pars_landau[1],
241  pars_landau[2],
242  pars_expo[0],
243  pars_expo[1],
244  t_direct_min};
245  VUVTiming = TF1("VUVTiming", model_close, 0., signal_t_range, 7);
246  VUVTiming.SetParameters(parsfinal);
247  }
248 
249  // store min/max, necessary to transform to the domain of the
250  // original distribution
251  // const size_t nq_max = 1;
252  double xq_max[1] = {0.975}; // include 97.5% of tail
253  double yq_max[1];
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);
259  fVUV_max[angle_bin][index] = max;
260  fVUV_min[angle_bin][index] = min;
261 
262  // create the distributions that represent the parametrised timing,
263  // and use those to create a RNG that samples said distributions
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);
268  fVUVTimingGen[angle_bin][index] =
269  CLHEP::RandGeneral(scintTimeEngine, vuv_timings.data(), vuv_timings.size());
270  } // index < num_params
271  } // angle_bin < num_angles
272  }
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)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
static double model_close(const double *x, const double *par)
std::vector< std::vector< CLHEP::RandGeneral > > fVUVTimingGen
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
std::vector< std::vector< double > > fVUV_min
T get(std::string const &key) const
Definition: ParameterSet.h:314
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
static double model_far(const double *x, const double *par)
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)
hh[ixs]
Definition: PlotSingle.C:6
std::vector< std::vector< double > > fVUV_max
void phot::PropagationTimeModel::getVISTimes ( std::vector< double > &  arrivalTimes,
const geo::Point_t scintPoint,
const geo::Point_t opDetPoint 
)
private

Definition at line 276 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, util::size(), and x.

Referenced by propagationTime().

279  {
280  // ***************************************************************************
281  // Calculation of earliest arrival times and corresponding unsmeared
282  // distribution
283  // ***************************************************************************
284 
285  // set plane_depth for correct TPC:
286  double plane_depth;
287  if (scintPoint.X() < 0.)
288  plane_depth = -fplane_depth;
289  else
290  plane_depth = fplane_depth;
291 
292  // calculate point of reflection for shortest path
293  geo::Point_t bounce_point(plane_depth, scintPoint.Y(), scintPoint.Z());
294 
295  // calculate distance travelled by VUV light and by vis light
296  double VUVdist = std::sqrt((bounce_point - scintPoint).Mag2());
297  double Visdist = std::sqrt((opDetPoint - bounce_point).Mag2());
298 
299  // calculate times taken by VUV part of path
300  int angle_bin_vuv = 0; // on-axis by definition
301  getVUVTimes(arrivalTimes, VUVdist, angle_bin_vuv);
302 
303  // sum parts to get total transport times times
304  for (size_t i = 0; i < arrivalTimes.size(); ++i) {
305  arrivalTimes[i] += Visdist / fvis_vmean;
306  }
307 
308  // ***************************************************************************
309  // Smearing of arrival time distribution
310  // ***************************************************************************
311  // calculate fastest time possible
312  // vis part
313  double vis_time = Visdist / fvis_vmean;
314  // vuv part
315  double vuv_time;
316  if (VUVdist < fmin_d) { vuv_time = VUVdist / fvuv_vgroup_max; }
317  else {
318  // find index of required parameterisation
319  const size_t index = std::round((VUVdist - fmin_d) / fstep_size);
320  // find shortest time
321  vuv_time = fVUV_min[angle_bin_vuv][index];
322  }
323  double fastest_time = vis_time + vuv_time;
324 
325  // calculate angle theta between bound_point and optical detector
326  double cosine_theta = std::abs(opDetPoint.X() - bounce_point.X()) / Visdist;
327  double theta = fast_acos(cosine_theta) * 180. / CLHEP::pi;
328 
329  // determine smearing parameters using interpolation of generated points:
330  // 1). tau = exponential smearing factor, varies with distance and angle
331  // 2). cutoff = largest smeared time allowed, preventing excessively large
332  // times caused by exponential distance to cathode
333  double distance_cathode_plane = std::abs(plane_depth - scintPoint.X());
334  // angular bin
335  size_t theta_bin = theta / fangle_bin_timing_vis;
336  // radial distance from centre of TPC (y,z plane)
337  double r =
338  std::hypot(scintPoint.Y() - fcathode_centre.Y(), scintPoint.Z() - fcathode_centre.Z());
339 
340  // cut-off and tau
341  // cut-off
342  // interpolate in d_c for each r bin
343  std::vector<double> interp_vals(fcut_off_pars[theta_bin].size(), 0.0);
344  for (size_t i = 0; i < fcut_off_pars[theta_bin].size(); i++) {
345  interp_vals[i] =
346  interpolate(fdistances_refl, fcut_off_pars[theta_bin][i], distance_cathode_plane, true);
347  }
348  // interpolate in r
349  double cutoff = interpolate(fradial_distances_refl, interp_vals, r, true);
350 
351  // tau
352  // interpolate in x for each r bin
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++) {
355  interp_vals_tau[i] =
356  interpolate(fdistances_refl, ftau_pars[theta_bin][i], distance_cathode_plane, true);
357  }
358  // interpolate in r
359  double tau = interpolate(fradial_distances_refl, interp_vals_tau, r, true);
360 
361  // apply smearing:
362  for (size_t i = 0; i < arrivalTimes.size(); ++i) {
363  double arrival_time_smeared;
364  // if time is already greater than cutoff, do not apply smearing
365  if (arrivalTimes[i] >= cutoff) { continue; }
366  // otherwise smear
367  else {
368  unsigned int counter = 0;
369  // loop until time generated is within cutoff limit
370  // most are within single attempt, very few take more than two
371  do {
372  // don't attempt smearings too many times
373  if (counter >= 10) {
374  arrival_time_smeared = arrivalTimes[i]; // don't smear
375  break;
376  }
377  else {
378  // generate random number in appropriate range
379  double x = fUniformGen.fire(0.5, 1.0);
380  // apply the exponential smearing
381  arrival_time_smeared =
382  arrivalTimes[i] + (arrivalTimes[i] - fastest_time) * (std::pow(x, -tau) - 1);
383  }
384  counter++;
385  } while (arrival_time_smeared > cutoff);
386  }
387  arrivalTimes[i] = arrival_time_smeared;
388  }
389  }
Float_t x
Definition: compare.C:6
TRandom r
Definition: spectrum.C:23
std::vector< std::vector< std::vector< double > > > fcut_off_pars
constexpr auto abs(T v)
Returns the absolute value of the argument.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
Definition: counter.h:295
const geo::Point_t fcathode_centre
std::vector< std::vector< double > > fVUV_min
double fast_acos(double xin)
std::vector< double > fradial_distances_refl
std::vector< std::vector< std::vector< double > > > ftau_pars
void getVUVTimes(std::vector< double > &arrivalTimes, const double distance_in_cm, const size_t angle_bin)
std::vector< double > fdistances_refl
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:180
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, const double x, const bool extrapolate, size_t i)
void phot::PropagationTimeModel::getVUVTimes ( std::vector< double > &  arrivalTimes,
const double  distance_in_cm,
const size_t  angle_bin 
)
private

Definition at line 109 of file PropagationTimeModel.cxx.

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

Referenced by getVISTimes(), and propagationTime().

112  {
113  if (distance < fmin_d) {
114  // times are fixed shift i.e. direct path only
115  double t_prop_correction = distance / fvuv_vgroup_mean;
116  for (size_t i = 0; i < arrivalTimes.size(); ++i) {
117  arrivalTimes[i] = t_prop_correction;
118  }
119  }
120  else {
121  // determine nearest parameterisation in discretisation
122  int index = std::round((distance - fmin_d) / fstep_size);
123  // randomly sample parameterisation for each photon
124  for (size_t i = 0; i < arrivalTimes.size(); ++i) {
125  arrivalTimes[i] = fVUVTimingGen[angle_bin][index].fire() *
126  (fVUV_max[angle_bin][index] - fVUV_min[angle_bin][index]) +
127  fVUV_min[angle_bin][index];
128  }
129  }
130  }
std::vector< std::vector< CLHEP::RandGeneral > > fVUVTimingGen
std::vector< std::vector< double > > fVUV_min
std::vector< std::vector< double > > fVUV_max
void phot::PropagationTimeModel::getVUVTimesGeo ( std::vector< double > &  arrivalTimes,
const double  distance_in_cm 
) const
private

Definition at line 135 of file PropagationTimeModel.cxx.

References fvuv_vgroup_mean.

Referenced by propagationTime().

137  {
138  // times are fixed shift i.e. direct path only
139  double t_prop_correction = distance / fvuv_vgroup_mean;
140  for (size_t i = 0; i < arrivalTimes.size(); ++i) {
141  arrivalTimes[i] = t_prop_correction;
142  }
143  }
double phot::PropagationTimeModel::model_close ( const double *  x,
const double *  par 
)
staticprivate

Definition at line 443 of file PropagationTimeModel.cxx.

References y1, and y2.

Referenced by generateVUVParams().

444  {
445  // par0 = joining point
446  // par1 = Landau MPV
447  // par2 = Landau width
448  // par3 = normalization
449  // par4 = Expo cte
450  // par5 = Expo tau
451  // par6 = t_min
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.;
456  return (y1 + y2);
457  }
Float_t x
Definition: compare.C:6
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t y2[n_points_geant4]
Definition: compare.C:26
double phot::PropagationTimeModel::model_far ( const double *  x,
const double *  par 
)
staticprivate

Definition at line 460 of file PropagationTimeModel.cxx.

Referenced by generateVUVParams().

461  {
462  // par1 = Landau MPV
463  // par2 = Landau width
464  // par3 = normalization
465  // par0 = t_min
466  if (x[0] <= par[0]) return 0.;
467  return par[3] * TMath::Landau(x[0], par[1], par[2]);
468  }
Float_t x
Definition: compare.C:6
std::vector< geo::Point_t > phot::PropagationTimeModel::opDetCenters ( ) const
private

Definition at line 400 of file PropagationTimeModel.cxx.

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

401  {
402  std::vector<geo::Point_t> opDetCenter;
403  for (size_t const i : util::counter(fGeom.NOpDets())) {
404  geo::OpDetGeo const& opDet = fGeom.OpDetGeoFromOpDet(i);
405  opDetCenter.push_back(opDet.GetCenter());
406  }
407  return opDetCenter;
408  }
OpDetGeo const & OpDetGeoFromOpDet(unsigned int OpDet) const
Returns the geo::OpDetGeo object for the given detector number.
geo::GeometryCore const & fGeom
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
Definition: counter.h:295
unsigned int NOpDets() const
Number of OpDets in the whole detector.
geo::Point_t const & GetCenter() const
Definition: OpDetGeo.h:72
std::vector< int > phot::PropagationTimeModel::opDetOrientations ( ) const
private

Definition at line 410 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().

411  {
412  std::vector<int> opDetOrientation;
413  for (size_t const i : util::counter(fGeom.NOpDets())) {
414  geo::OpDetGeo const& opDet = fGeom.OpDetGeoFromOpDet(i);
415  if (opDet.isSphere()) { // dome PMTs
416  opDetOrientation.push_back(0);
417  }
418  else if (opDet.isBar()) {
419  // determine orientation to get correction OpDet dimensions
420  if (opDet.Width() > opDet.Height()) { // laterals, Y dimension smallest
421  opDetOrientation.push_back(1);
422  }
423  else { // anode/cathode (default), X dimension smallest
424  opDetOrientation.push_back(0);
425  }
426  }
427  else {
428  opDetOrientation.push_back(0);
429  }
430  }
431  return opDetOrientation;
432  }
OpDetGeo const & OpDetGeoFromOpDet(unsigned int OpDet) const
Returns the geo::OpDetGeo object for the given detector number.
geo::GeometryCore const & fGeom
bool isBar() const
Returns whether the detector shape is a bar (TGeoBBox).
Definition: OpDetGeo.h:180
bool isSphere() const
Returns whether the detector shape is a hemisphere (TGeoSphere).
Definition: OpDetGeo.h:183
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
Definition: counter.h:295
unsigned int NOpDets() const
Number of OpDets in the whole detector.
double Width() const
Definition: OpDetGeo.h:79
double Height() const
Definition: OpDetGeo.h:80
void phot::PropagationTimeModel::propagationTime ( std::vector< double > &  arrivalTimes,
const geo::Point_t x0,
const size_t  OpChannel,
const bool  Reflected = false 
)

Definition at line 70 of file PropagationTimeModel.cxx.

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

74  {
75  if (!fGeoPropTimeOnly) {
76  // Get VUV photons transport time distribution from the parametrization
77  geo::Point_t const& opDetCenter = fOpDetCenter[OpChannel];
78  if (!Reflected) {
79  double distance =
80  std::hypot(x0.X() - opDetCenter.X(), x0.Y() - opDetCenter.Y(), x0.Z() - opDetCenter.Z());
81  double cosine;
82  if (fOpDetOrientation[OpChannel] == 1)
83  cosine = std::abs(x0.Y() - opDetCenter.Y()) / distance;
84  else
85  cosine = std::abs(x0.X() - opDetCenter.X()) / distance;
86 
87  double theta = fast_acos(cosine) * 180. / CLHEP::pi;
88  int angle_bin = theta / fangle_bin_timing_vuv;
89  getVUVTimes(arrivalTimes, distance, angle_bin); // in ns
90  }
91  else {
92  getVISTimes(arrivalTimes, x0, opDetCenter); // in ns
93  }
94  }
95  else if (fGeoPropTimeOnly && !Reflected) {
96  // Get VUV photons arrival time geometrically
97  geo::Point_t const& opDetCenter = fOpDetCenter[OpChannel];
98  double distance =
99  std::hypot(x0.X() - opDetCenter.X(), x0.Y() - opDetCenter.Y(), x0.Z() - opDetCenter.Z());
100  getVUVTimesGeo(arrivalTimes, distance); // in ns
101  }
102  else {
103  throw cet::exception("PropagationTimeModel") << "Propagation time model not found.";
104  }
105  }
const std::vector< geo::Point_t > fOpDetCenter
const std::vector< int > fOpDetOrientation
constexpr auto abs(T v)
Returns the absolute value of the argument.
double fast_acos(double xin)
void getVUVTimesGeo(std::vector< double > &arrivalTimes, const double distance_in_cm) const
void getVUVTimes(std::vector< double > &arrivalTimes, const double distance_in_cm, const size_t angle_bin)
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:180
void getVISTimes(std::vector< double > &arrivalTimes, const geo::Point_t &scintPoint, const geo::Point_t &opDetPoint)
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33

Member Data Documentation

double phot::PropagationTimeModel::fangle_bin_timing_vis
private

Definition at line 106 of file PropagationTimeModel.h.

Referenced by getVISTimes(), and PropagationTimeModel().

double phot::PropagationTimeModel::fangle_bin_timing_vuv
private

Definition at line 96 of file PropagationTimeModel.h.

Referenced by generateVUVParams(), propagationTime(), and PropagationTimeModel().

const geo::Point_t phot::PropagationTimeModel::fcathode_centre
private

Definition at line 89 of file PropagationTimeModel.h.

Referenced by getVISTimes().

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

Definition at line 109 of file PropagationTimeModel.h.

Referenced by getVISTimes(), and PropagationTimeModel().

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

Definition at line 107 of file PropagationTimeModel.h.

Referenced by getVISTimes(), and PropagationTimeModel().

geo::GeometryCore const& phot::PropagationTimeModel::fGeom
private

Definition at line 87 of file PropagationTimeModel.h.

Referenced by cathodeCentre(), opDetCenters(), and opDetOrientations().

const bool phot::PropagationTimeModel::fGeoPropTimeOnly
private

Definition at line 81 of file PropagationTimeModel.h.

Referenced by propagationTime(), and PropagationTimeModel().

double phot::PropagationTimeModel::finflexion_point_distance
private

Definition at line 96 of file PropagationTimeModel.h.

Referenced by generateVUVParams(), and PropagationTimeModel().

double phot::PropagationTimeModel::fmin_d
private
const std::vector<geo::Point_t> phot::PropagationTimeModel::fOpDetCenter
private

Definition at line 92 of file PropagationTimeModel.h.

Referenced by propagationTime().

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

Definition at line 93 of file PropagationTimeModel.h.

Referenced by propagationTime().

const double phot::PropagationTimeModel::fplane_depth
private

Definition at line 88 of file PropagationTimeModel.h.

Referenced by getVISTimes().

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

Definition at line 108 of file PropagationTimeModel.h.

Referenced by getVISTimes(), and PropagationTimeModel().

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

Definition at line 110 of file PropagationTimeModel.h.

Referenced by getVISTimes(), and PropagationTimeModel().

CLHEP::RandFlat phot::PropagationTimeModel::fUniformGen
private

Definition at line 84 of file PropagationTimeModel.h.

Referenced by getVISTimes().

double phot::PropagationTimeModel::fvis_vmean
private

Definition at line 106 of file PropagationTimeModel.h.

Referenced by getVISTimes(), and PropagationTimeModel().

std::vector<std::vector<double> > phot::PropagationTimeModel::fVUV_max
private

Definition at line 102 of file PropagationTimeModel.h.

Referenced by generateVUVParams(), and getVUVTimes().

std::vector<std::vector<double> > phot::PropagationTimeModel::fVUV_min
private

Definition at line 103 of file PropagationTimeModel.h.

Referenced by generateVUVParams(), getVISTimes(), and getVUVTimes().

double phot::PropagationTimeModel::fvuv_vgroup_max
private

Definition at line 96 of file PropagationTimeModel.h.

Referenced by generateVUVParams(), getVISTimes(), and PropagationTimeModel().

double phot::PropagationTimeModel::fvuv_vgroup_mean
private
std::vector<std::vector<CLHEP::RandGeneral> > phot::PropagationTimeModel::fVUVTimingGen
private

Definition at line 99 of file PropagationTimeModel.h.

Referenced by generateVUVParams(), and getVUVTimes().


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