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

#include "SemiAnalyticalModel.h"

Classes

struct  Dims
 
struct  OpticalDetector
 

Public Member Functions

 SemiAnalyticalModel (const fhicl::ParameterSet &VUVHits, const fhicl::ParameterSet &VISHits, const bool doReflectedLight=false, const bool includeAnodeReflections=false, const bool useXeAbsorption=false)
 
void detectedDirectVisibilities (std::vector< double > &DetectedVisibilities, geo::Point_t const &ScintPoint) const
 
void detectedReflectedVisibilities (std::vector< double > &ReflDetectedVisibilities, geo::Point_t const &ScintPoint, bool AnodeMode=false) const
 

Private Member Functions

double VUVAbsorptionLength () const
 
double VUVVisibility (geo::Point_t const &ScintPoint, OpticalDetector const &opDet) const
 
double VISVisibility (geo::Point_t const &ScintPoint, OpticalDetector const &opDet, const double cathode_visibility, geo::Point_t const &hotspot, bool AnodeMode=false) const
 
double Gaisser_Hillas (const double x, const double *par) const
 
double Rectangle_SolidAngle (const double a, const double b, const double d) const
 
double Rectangle_SolidAngle (Dims const &o, geo::Vector_t const &v, const int OpDetOrientation) const
 
double Disk_SolidAngle (const double d, const double h, const double b) const
 
double Omega_Dome_Model (const double distance, const double theta) const
 
bool isOpDetInSameTPC (geo::Point_t const &ScintPoint, geo::Point_t const &OpDetPoint) const
 
std::vector< OpticalDetectoropticalDetectors () const
 

Private Attributes

geo::GeometryCore const & fGeom
 
const larg4::ISTPC fISTPC
 
const int fNTPC
 
const std::vector< geo::BoxBoundedGeofActiveVolumes
 
const TVector3 fcathode_centre
 
const TVector3 fanode_centre
 
double fplane_depth
 
double fanode_plane_depth
 
const size_t fNOpDets
 
const std::vector< OpticalDetectorfOpDetector
 
double fradius
 
Dims fcathode_plane
 
Dims fanode_plane
 
double fdelta_angulo_vuv
 
bool fIsFlatPDCorr
 
std::vector< std::vector< double > > fGHvuvpars_flat
 
std::vector< double > fborder_corr_angulo_flat
 
std::vector< std::vector< double > > fborder_corr_flat
 
bool fIsFlatPDCorrLat
 
std::vector< std::vector< double > > fGHvuvpars_flat_lateral
 
std::vector< double > fborder_corr_angulo_flat_lateral
 
std::vector< std::vector< double > > fborder_corr_flat_lateral
 
bool fIsDomePDCorr
 
std::vector< std::vector< double > > fGHvuvpars_dome
 
std::vector< double > fborder_corr_angulo_dome
 
std::vector< std::vector< double > > fborder_corr_dome
 
bool fApplyFieldCageTransparency
 
double fFieldCageTransparencyLateral
 
double fFieldCageTransparencyCathode
 
const bool fDoReflectedLight
 
const bool fIncludeAnodeReflections
 
double fdelta_angulo_vis
 
double fAnodeReflectivity
 
std::vector< double > fvis_distances_x_flat
 
std::vector< double > fvis_distances_r_flat
 
std::vector< std::vector< std::vector< double > > > fvispars_flat
 
std::vector< double > fvis_distances_x_flat_lateral
 
std::vector< double > fvis_distances_r_flat_lateral
 
std::vector< std::vector< std::vector< double > > > fvispars_flat_lateral
 
std::vector< double > fvis_distances_x_dome
 
std::vector< double > fvis_distances_r_dome
 
std::vector< std::vector< std::vector< double > > > fvispars_dome
 
const bool fUseXeAbsorption
 
double fvuv_absorption_length
 

Detailed Description

Definition at line 33 of file SemiAnalyticalModel.h.

Constructor & Destructor Documentation

phot::SemiAnalyticalModel::SemiAnalyticalModel ( const fhicl::ParameterSet VUVHits,
const fhicl::ParameterSet VISHits,
const bool  doReflectedLight = false,
const bool  includeAnodeReflections = false,
const bool  useXeAbsorption = false 
)

Definition at line 28 of file SemiAnalyticalModel.cxx.

References util::abs(), larg4::ISTPC::extractActiveLArVolume(), fActiveVolumes, fanode_centre, fanode_plane, fanode_plane_depth, fAnodeReflectivity, fApplyFieldCageTransparency, fborder_corr_angulo_dome, fborder_corr_angulo_flat, fborder_corr_angulo_flat_lateral, fborder_corr_dome, fborder_corr_flat, fborder_corr_flat_lateral, fcathode_centre, fcathode_plane, fdelta_angulo_vis, fdelta_angulo_vuv, fDoReflectedLight, fFieldCageTransparencyCathode, fFieldCageTransparencyLateral, fGeom, fGHvuvpars_dome, fGHvuvpars_flat, fGHvuvpars_flat_lateral, fIncludeAnodeReflections, geo::TPCGeo::FirstPlane(), fIsDomePDCorr, fIsFlatPDCorr, fIsFlatPDCorrLat, fISTPC, fNOpDets, fNTPC, fOpDetector, fplane_depth, fradius, fUseXeAbsorption, fvis_distances_r_dome, fvis_distances_r_flat, fvis_distances_r_flat_lateral, fvis_distances_x_dome, fvis_distances_x_flat, fvis_distances_x_flat_lateral, fvispars_dome, fvispars_flat, fvispars_flat_lateral, fvuv_absorption_length, geo::TPCGeo::GetCathodeCenter(), geo::PlaneGeo::GetCenter(), phot::SemiAnalyticalModel::Dims::h, geo::GeometryCore::NOpDets(), geo::GeometryCore::NTPC(), opticalDetectors(), geo::GeometryCore::TPC(), VUVAbsorptionLength(), and phot::SemiAnalyticalModel::Dims::w.

33  : fGeom{*(lar::providerFrom<geo::Geometry>())}
34  , fISTPC{fGeom}
35  , fNTPC(fGeom.NTPC())
38  fActiveVolumes[0].CenterY(),
39  fActiveVolumes[0].CenterZ()}
41  fActiveVolumes[0].CenterY(),
42  fActiveVolumes[0].CenterZ()}
45  , fDoReflectedLight(doReflectedLight)
46  , fIncludeAnodeReflections(includeAnodeReflections)
47  , fUseXeAbsorption(useXeAbsorption)
48  {
49  mf::LogInfo("SemiAnalyticalModel") << "Initializing Semi-analytical model." << std::endl;
50 
51  // Load Gaisser-Hillas corrections for VUV semi-analytic hits
52  mf::LogInfo("SemiAnalyticalModel") << "Using VUV visibility parameterization";
53 
54  fIsFlatPDCorr = VUVHitsParams.get<bool>("FlatPDCorr", false);
55  fIsFlatPDCorrLat = VUVHitsParams.get<bool>("FlatPDCorrLat", false);
56  fIsDomePDCorr = VUVHitsParams.get<bool>("DomePDCorr", false);
57  fdelta_angulo_vuv = VUVHitsParams.get<double>("delta_angulo_vuv", 10);
58  fradius = VUVHitsParams.get<double>("PMT_radius", 10.16);
59  fApplyFieldCageTransparency = VUVHitsParams.get<bool>("ApplyFieldCageTransparency", false);
60  fFieldCageTransparencyLateral = VUVHitsParams.get<double>("FieldCageTransparencyLateral", 1.0);
61  fFieldCageTransparencyCathode = VUVHitsParams.get<double>("FieldCageTransparencyCathode", 1.0);
62 
64  throw cet::exception("SemiAnalyticalModel")
65  << "Both isFlatPDCorr/isFlatPDCorrLat and isDomePDCorr parameters are "
66  "false, at least one type of parameterisation is required for the "
67  "semi-analytic light simulation."
68  << "\n";
69  }
70  if (fIsFlatPDCorr) {
71  fGHvuvpars_flat = VUVHitsParams.get<std::vector<std::vector<double>>>("GH_PARS_flat");
72  fborder_corr_angulo_flat = VUVHitsParams.get<std::vector<double>>("GH_border_angulo_flat");
73  fborder_corr_flat = VUVHitsParams.get<std::vector<std::vector<double>>>("GH_border_flat");
74  }
75  if (fIsFlatPDCorrLat) {
77  VUVHitsParams.get<std::vector<std::vector<double>>>("GH_PARS_flat_lateral");
79  VUVHitsParams.get<std::vector<double>>("GH_border_angulo_flat_lateral");
81  VUVHitsParams.get<std::vector<std::vector<double>>>("GH_border_flat_lateral");
82  }
83  if (fIsDomePDCorr) {
84  fGHvuvpars_dome = VUVHitsParams.get<std::vector<std::vector<double>>>("GH_PARS_dome");
85  fborder_corr_angulo_dome = VUVHitsParams.get<std::vector<double>>("GH_border_angulo_dome");
86  fborder_corr_dome = VUVHitsParams.get<std::vector<std::vector<double>>>("GH_border_dome");
87  }
88 
89  // Load corrections for VIS semi-analytic hits
90  if (fDoReflectedLight) {
91  mf::LogInfo("SemiAnalyticalModel") << "Using VIS (reflected) visibility parameterization";
92  fdelta_angulo_vis = VISHitsParams.get<double>("delta_angulo_vis");
93 
94  if (fIsFlatPDCorr) {
95  fvis_distances_x_flat = VISHitsParams.get<std::vector<double>>("VIS_distances_x_flat");
96  fvis_distances_r_flat = VISHitsParams.get<std::vector<double>>("VIS_distances_r_flat");
98  VISHitsParams.get<std::vector<std::vector<std::vector<double>>>>("VIS_correction_flat");
99  }
100  if (fIsDomePDCorr) {
101  fvis_distances_x_dome = VISHitsParams.get<std::vector<double>>("VIS_distances_x_dome");
102  fvis_distances_r_dome = VISHitsParams.get<std::vector<double>>("VIS_distances_r_dome");
103  fvispars_dome =
104  VISHitsParams.get<std::vector<std::vector<std::vector<double>>>>("VIS_correction_dome");
105  }
106 
107  // set cathode plane struct for solid angle function
108  fcathode_plane.h = fActiveVolumes[0].SizeY();
109  fcathode_plane.w = fActiveVolumes[0].SizeZ();
111  }
112 
113  // Load corrections for Anode reflections configuration
115  mf::LogInfo("SemiAnalyticalModel") << "Using anode reflections parameterization";
116  fdelta_angulo_vis = VISHitsParams.get<double>("delta_angulo_vis");
117  fAnodeReflectivity = VISHitsParams.get<double>("AnodeReflectivity");
118 
119  if (fIsFlatPDCorr) {
120  fvis_distances_x_flat = VISHitsParams.get<std::vector<double>>("VIS_distances_x_flat");
121  fvis_distances_r_flat = VISHitsParams.get<std::vector<double>>("VIS_distances_r_flat");
122  fvispars_flat =
123  VISHitsParams.get<std::vector<std::vector<std::vector<double>>>>("VIS_correction_flat");
124  }
125 
126  if (fIsFlatPDCorrLat) {
128  VISHitsParams.get<std::vector<double>>("VIS_distances_x_flat_lateral");
130  VISHitsParams.get<std::vector<double>>("VIS_distances_r_flat_lateral");
131  fvispars_flat_lateral = VISHitsParams.get<std::vector<std::vector<std::vector<double>>>>(
132  "VIS_correction_flat_lateral");
133  }
134 
135  // set anode plane struct for solid angle function
136  fanode_plane.h = fActiveVolumes[0].SizeY();
137  fanode_plane.w = fActiveVolumes[0].SizeZ();
139  }
140 
141  // set absorption length
143  mf::LogInfo("SemiAnalyticalModel")
144  << "Setting absorption length to: " << fvuv_absorption_length << std::endl;
145 
146  mf::LogInfo("SemiAnalyticalModel") << "Semi-analytical model initialized." << std::endl;
147  }
std::vector< double > fvis_distances_x_dome
std::vector< std::vector< double > > fborder_corr_flat_lateral
Point_t GetCathodeCenter() const
Definition: TPCGeo.h:254
unsigned int NTPC(CryostatID const &cryoid=cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:686
static std::vector< geo::BoxBoundedGeo > extractActiveLArVolume(geo::GeometryCore const &geom)
Definition: ISTPC.cxx:45
Point_t const & GetCenter() const
Returns the centre of the wire plane in world coordinates [cm].
Definition: PlaneGeo.h:435
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
geo::PlaneGeo const & FirstPlane() const
Returns the first wire plane (the closest to TPC center).
Definition: TPCGeo.h:212
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::vector< std::vector< std::vector< double > > > fvispars_dome
std::vector< std::vector< double > > fGHvuvpars_flat_lateral
std::vector< double > fvis_distances_r_flat
const std::vector< geo::BoxBoundedGeo > fActiveVolumes
std::vector< std::vector< double > > fborder_corr_dome
TPCGeo const & TPC(TPCID const &tpcid=tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:722
geo::GeometryCore const & fGeom
std::vector< double > fborder_corr_angulo_flat
std::vector< double > fvis_distances_r_dome
std::vector< double > fborder_corr_angulo_dome
std::vector< std::vector< std::vector< double > > > fvispars_flat_lateral
std::vector< double > fvis_distances_x_flat_lateral
unsigned int NOpDets() const
Number of OpDets in the whole detector.
std::vector< std::vector< double > > fGHvuvpars_flat
std::vector< std::vector< double > > fGHvuvpars_dome
std::vector< double > fvis_distances_x_flat
std::vector< std::vector< std::vector< double > > > fvispars_flat
std::vector< double > fvis_distances_r_flat_lateral
std::vector< OpticalDetector > opticalDetectors() const
const std::vector< OpticalDetector > fOpDetector
std::vector< std::vector< double > > fborder_corr_flat
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< double > fborder_corr_angulo_flat_lateral

Member Function Documentation

void phot::SemiAnalyticalModel::detectedDirectVisibilities ( std::vector< double > &  DetectedVisibilities,
geo::Point_t const &  ScintPoint 
) const

Definition at line 175 of file SemiAnalyticalModel.cxx.

References util::counter(), fNOpDets, fOpDetector, isOpDetInSameTPC(), and VUVVisibility().

177  {
178  DetectedVisibilities.resize(fNOpDets);
179  for (size_t const OpDet : util::counter(fNOpDets)) {
180  if (!isOpDetInSameTPC(ScintPoint, fOpDetector[OpDet].center)) {
181  DetectedVisibilities[OpDet] = 0.;
182  continue;
183  }
184 
185  DetectedVisibilities[OpDet] = VUVVisibility(ScintPoint, fOpDetector[OpDet]);
186  ;
187  }
188  }
bool isOpDetInSameTPC(geo::Point_t const &ScintPoint, geo::Point_t const &OpDetPoint) const
double VUVVisibility(geo::Point_t const &ScintPoint, OpticalDetector const &opDet) const
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 std::vector< OpticalDetector > fOpDetector
void phot::SemiAnalyticalModel::detectedReflectedVisibilities ( std::vector< double > &  ReflDetectedVisibilities,
geo::Point_t const &  ScintPoint,
bool  AnodeMode = false 
) const

Definition at line 317 of file SemiAnalyticalModel.cxx.

References util::abs(), util::counter(), fanode_centre, fanode_plane, fanode_plane_depth, fborder_corr_angulo_flat, fborder_corr_flat, fcathode_centre, fcathode_plane, fGHvuvpars_flat, fIsFlatPDCorr, fNOpDets, fOpDetector, fplane_depth, fvuv_absorption_length, Gaisser_Hillas(), phot::interpolate(), isOpDetInSameTPC(), util::pi(), r, Rectangle_SolidAngle(), and VISVisibility().

321  {
322  // 1). calculate visibility of VUV photons on
323  // reflective foils via solid angle + Gaisser-Hillas
324  // corrections:
325 
326  // get scintpoint coords relative to centre of cathode plane and set plane
327  // dimensions
328  geo::Vector_t ScintPoint_relative;
329  Dims plane_dimensions;
330  double plane_depth;
331  if (AnodeMode) {
332  plane_dimensions = fanode_plane;
333  plane_depth = fanode_plane_depth;
334  ScintPoint_relative.SetCoordinates(std::abs(ScintPoint.X() - fanode_plane_depth),
335  std::abs(ScintPoint.Y() - fanode_centre[1]),
336  std::abs(ScintPoint.Z() - fanode_centre[2]));
337  }
338  else {
339  plane_dimensions = fcathode_plane;
340  plane_depth = ScintPoint.X() < 0. ? -fplane_depth : fplane_depth;
341  ScintPoint_relative.SetCoordinates(std::abs(ScintPoint.X() - plane_depth),
342  std::abs(ScintPoint.Y() - fcathode_centre[1]),
343  std::abs(ScintPoint.Z() - fcathode_centre[2]));
344  }
345 
346  // calculate solid angle of anode/cathode from the scintillation point,
347  // orientation always = 0 (anode/cathode)
348  double solid_angle_cathode = Rectangle_SolidAngle(plane_dimensions, ScintPoint_relative, 0);
349 
350  // calculate distance and angle between ScintPoint and hotspot
351  // vast majority of hits in hotspot region directly infront of scintpoint,
352  // therefore consider attenuation for this distance and on axis GH instead of
353  // for the centre coordinate
354  double distance_cathode = std::abs(plane_depth - ScintPoint.X());
355  // calculate hits on cathode plane via geometric acceptance
356  double cathode_visibility_geo = std::exp(-1. * distance_cathode / fvuv_absorption_length) *
357  (solid_angle_cathode / (4. * CLHEP::pi));
358 
359  // determine Gaisser-Hillas correction including border effects
360  // use flat correction
361  double r = std::hypot(ScintPoint.Y() - fcathode_centre[1], ScintPoint.Z() - fcathode_centre[2]);
362  double pars_ini[4] = {0, 0, 0, 0};
363  double s1 = 0;
364  double s2 = 0;
365  double s3 = 0;
366  if (fIsFlatPDCorr) {
367  pars_ini[0] = fGHvuvpars_flat[0][0];
368  pars_ini[1] = fGHvuvpars_flat[1][0];
369  pars_ini[2] = fGHvuvpars_flat[2][0];
370  pars_ini[3] = fGHvuvpars_flat[3][0];
374  }
375  else {
376  throw cet::exception("SemiAnalyticalModel")
377  << "Error: flat optical detector VUV correction required for reflected "
378  "semi-analytic hits. - configuration error in semi-analytical model."
379  << "\n";
380  }
381 
382  // add border correction
383  pars_ini[0] = pars_ini[0] + s1 * r;
384  pars_ini[1] = pars_ini[1] + s2 * r;
385  pars_ini[2] = pars_ini[2] + s3 * r;
386  pars_ini[3] = pars_ini[3];
387 
388  // calculate corrected number of hits
389  double GH_correction = Gaisser_Hillas(distance_cathode, pars_ini);
390  const double cathode_visibility_rec = GH_correction * cathode_visibility_geo;
391 
392  // 2). detemine visibility of each PD
393  const geo::Point_t hotspot = {plane_depth, ScintPoint.Y(), ScintPoint.Z()};
394  ReflDetectedVisibilities.resize(fNOpDets);
395  for (size_t const OpDet : util::counter(fNOpDets)) {
396  if (!isOpDetInSameTPC(ScintPoint, fOpDetector[OpDet].center)) {
397  ReflDetectedVisibilities[OpDet] = 0.;
398  continue;
399  }
400 
401  ReflDetectedVisibilities[OpDet] =
402  VISVisibility(ScintPoint, fOpDetector[OpDet], cathode_visibility_rec, hotspot, AnodeMode);
403  }
404  }
TRandom r
Definition: spectrum.C:23
double VISVisibility(geo::Point_t const &ScintPoint, OpticalDetector const &opDet, const double cathode_visibility, geo::Point_t const &hotspot, bool AnodeMode=false) const
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:160
bool isOpDetInSameTPC(geo::Point_t const &ScintPoint, geo::Point_t const &OpDetPoint) const
double Rectangle_SolidAngle(const double a, const double b, const double d) const
constexpr auto abs(T v)
Returns the absolute value of the argument.
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
std::vector< double > fborder_corr_angulo_flat
double Gaisser_Hillas(const double x, const double *par) const
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
std::vector< std::vector< double > > fGHvuvpars_flat
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
const std::vector< OpticalDetector > fOpDetector
std::vector< std::vector< double > > fborder_corr_flat
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, const double x, const bool extrapolate, size_t i)
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
double phot::SemiAnalyticalModel::Disk_SolidAngle ( const double  d,
const double  h,
const double  b 
) const
private

Definition at line 542 of file SemiAnalyticalModel.cxx.

References d, e, phot::isApproximatelyEqual(), phot::isApproximatelyZero(), phot::isDefinitelyGreaterThan(), phot::isDefinitelyLessThan(), and util::pi().

Referenced by VISVisibility(), and VUVVisibility().

543  {
544  if (b <= 0. || d < 0. || h <= 0.) return 0.;
545  const double leg2 = (b + d) * (b + d);
546  const double aa = std::sqrt(h * h / (h * h + leg2));
547  if (isApproximatelyZero(d)) { return 2. * CLHEP::pi * (1. - aa); }
548  double bb = 2. * std::sqrt(b * d / (h * h + leg2));
549  double cc = 4. * b * d / leg2;
550 
551  if (isDefinitelyGreaterThan(d, b)) {
552  try {
553  return 2. * aa *
554  (std::sqrt(1. - cc) * boost::math::ellint_3(bb, cc, noLDoublePromote()) -
555  boost::math::ellint_1(bb, noLDoublePromote()));
556  }
557  catch (std::domain_error& e) {
558  if (isApproximatelyEqual(d, b, 1e-9)) {
559  mf::LogWarning("SemiAnalyticalModel")
560  << "Elliptic Integral in Disk_SolidAngle() given parameters "
561  "outside domain."
562  << "\nbb: " << bb << "\ncc: " << cc << "\nException message: " << e.what()
563  << "\nRelax condition and carry on.";
564  return CLHEP::pi - 2. * aa * boost::math::ellint_1(bb, noLDoublePromote());
565  }
566  else {
567  mf::LogError("SemiAnalyticalModel")
568  << "Elliptic Integral inside Disk_SolidAngle() given parameters "
569  "outside domain.\n"
570  << "\nbb: " << bb << "\ncc: " << cc << "Exception message: " << e.what();
571  return 0.;
572  }
573  }
574  }
575  if (isDefinitelyLessThan(d, b)) {
576  try {
577  return 2. * CLHEP::pi -
578  2. * aa *
579  (boost::math::ellint_1(bb, noLDoublePromote()) +
580  std::sqrt(1. - cc) * boost::math::ellint_3(bb, cc, noLDoublePromote()));
581  }
582  catch (std::domain_error& e) {
583  if (isApproximatelyEqual(d, b, 1e-9)) {
584  mf::LogWarning("SemiAnalyticalModel")
585  << "Elliptic Integral in Disk_SolidAngle() given parameters "
586  "outside domain."
587  << "\nbb: " << bb << "\ncc: " << cc << "\nException message: " << e.what()
588  << "\nRelax condition and carry on.";
589  return CLHEP::pi - 2. * aa * boost::math::ellint_1(bb, noLDoublePromote());
590  }
591  else {
592  mf::LogError("SemiAnalyticalModel")
593  << "Elliptic Integral inside Disk_SolidAngle() given parameters "
594  "outside domain.\n"
595  << "\nbb: " << bb << "\ncc: " << cc << "Exception message: " << e.what();
596  return 0.;
597  }
598  }
599  }
600  if (isApproximatelyEqual(d, b)) {
601  return CLHEP::pi - 2. * aa * boost::math::ellint_1(bb, noLDoublePromote());
602  }
603  return 0.;
604  }
static constexpr bool isDefinitelyGreaterThan(TReal a, TReal b, TReal tolerance=std::numeric_limits< TReal >::epsilon())
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
boost::math::policies::policy< boost::math::policies::promote_double< false > > noLDoublePromote
Float_t d
Definition: plot.C:235
static constexpr bool isApproximatelyEqual(TReal a, TReal b, TReal tolerance=std::numeric_limits< TReal >::epsilon())
static constexpr bool isDefinitelyLessThan(TReal a, TReal b, TReal tolerance=std::numeric_limits< TReal >::epsilon())
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
Float_t e
Definition: plot.C:35
static constexpr bool isApproximatelyZero(TReal a, TReal tolerance=std::numeric_limits< TReal >::epsilon())
double phot::SemiAnalyticalModel::Gaisser_Hillas ( const double  x,
const double *  par 
) const
private

Definition at line 529 of file SemiAnalyticalModel.cxx.

Referenced by detectedReflectedVisibilities(), and VUVVisibility().

530  {
531  double X_mu_0 = par[3];
532  double Normalization = par[0];
533  double Diff = par[1] - X_mu_0;
534  double Term = std::pow((x - X_mu_0) / Diff, Diff / par[2]);
535  double Exponential = std::exp((par[1] - x) / par[2]);
536 
537  return (Normalization * Term * Exponential);
538  }
Float_t x
Definition: compare.C:6
bool phot::SemiAnalyticalModel::isOpDetInSameTPC ( geo::Point_t const &  ScintPoint,
geo::Point_t const &  OpDetPoint 
) const
private

Definition at line 730 of file SemiAnalyticalModel.cxx.

References util::abs(), and fNTPC.

Referenced by detectedDirectVisibilities(), and detectedReflectedVisibilities().

732  {
733  // method working for SBND, uBooNE, DUNE-HD 1x2x6 and DUNE-VD 1x8x6
734  // will need to be replaced to work in full DUNE geometry, ICARUS geometry
735  // check x coordinate has same sign or is close to zero, otherwise return
736  // false
737  if (((ScintPoint.X() < 0.) != (OpDetPoint.X() < 0.)) && std::abs(OpDetPoint.X()) > 10. &&
738  fNTPC == 2) { // TODO: replace with geometry service
739  return false;
740  }
741  return true;
742  }
constexpr auto abs(T v)
Returns the absolute value of the argument.
double phot::SemiAnalyticalModel::Omega_Dome_Model ( const double  distance,
const double  theta 
) const
private

Definition at line 695 of file SemiAnalyticalModel.cxx.

References fradius, and util::pi().

Referenced by VISVisibility(), and VUVVisibility().

696  {
697  // this function calculates the solid angle of a semi-sphere of radius b,
698  // as a correction to the analytic formula of the on-axix solid angle,
699  // as we move off-axis an angle theta. We have used 9-angular bins
700  // with delta_theta width.
701 
702  // par0 = Radius correction close
703  // par1 = Radius correction far
704  // par2 = breaking distance betwween "close" and "far"
705 
706  constexpr double par0[9] = {0., 0., 0., 0., 0., 0.597542, 1.00872, 1.46993, 2.04221};
707  constexpr double par1[9] = {
708  0., 0., 0.19569, 0.300449, 0.555598, 0.854939, 1.39166, 2.19141, 2.57732};
709  const double delta_theta = 10.; // TODO: should this be fdelta_angulo_vuv?
710  int j = int(theta / delta_theta);
711  // PMT radius
712  const double b = fradius; // cm
713  // distance form which the model parameters break (empirical value)
714  const double d_break = 5 * b; // par2
715 
716  if (distance >= d_break) {
717  double R_apparent_far = b - par1[j];
718  double ratio_square = (R_apparent_far * R_apparent_far) / (distance * distance);
719  return (2 * CLHEP::pi * (1 - std::sqrt(1 - ratio_square)));
720  }
721  else {
722  double R_apparent_close = b - par0[j];
723  double ratio_square = (R_apparent_close * R_apparent_close) / (distance * distance);
724  return (2 * CLHEP::pi * (1 - std::sqrt(1 - ratio_square)));
725  }
726  }
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
std::vector< SemiAnalyticalModel::OpticalDetector > phot::SemiAnalyticalModel::opticalDetectors ( ) const
private

Definition at line 744 of file SemiAnalyticalModel.cxx.

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

Referenced by SemiAnalyticalModel().

745  {
746  std::vector<SemiAnalyticalModel::OpticalDetector> opticalDetector;
747  for (size_t const i : util::counter(fNOpDets)) {
748  geo::OpDetGeo const& opDet = fGeom.OpDetGeoFromOpDet(i);
749  geo::Point_t center = opDet.GetCenter();
750  double height;
751  double length;
752  int type;
753  int orientation;
754  if (opDet.isSphere()) { // dome PMTs
755  type = 1; // dome
756  orientation = 0; // anode/cathode (default)
757  length = -1;
758  height = -1;
759  }
760  else if (opDet.isBar()) {
761  type = 0; // (X)Arapucas/Bars
762  // determine orientation to get correction OpDet dimensions
763  length = opDet.Length();
764  height = opDet.Height();
765  if (opDet.Width() > opDet.Length()) { // laterals along x-y plane, Z dimension smallest
766  orientation = 2;
767  length = opDet.Width();
768  }
769  else if (opDet.Width() > opDet.Height()) { // laterals along x-z plane, Y dimension smallest
770  orientation = 1;
771  height = opDet.Width();
772  }
773  else { // anode/cathode (default), X dimension smallest
774  orientation = 0;
775  }
776  }
777  else {
778  type = 2; // disk PMTs
779  orientation = 0; // anode/cathode (default)
780  length = -1;
781  height = -1;
782  }
783  opticalDetector.emplace_back(
784  SemiAnalyticalModel::OpticalDetector{height, length, center, type, orientation});
785  }
786  return opticalDetector;
787  }
OpDetGeo const & OpDetGeoFromOpDet(unsigned int OpDet) const
Returns the geo::OpDetGeo object for the given detector number.
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
geo::GeometryCore const & fGeom
double Length() const
Definition: OpDetGeo.h:78
geo::Point_t const & GetCenter() const
Definition: OpDetGeo.h:72
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 Width() const
Definition: OpDetGeo.h:79
double Height() const
Definition: OpDetGeo.h:80
double phot::SemiAnalyticalModel::Rectangle_SolidAngle ( const double  a,
const double  b,
const double  d 
) const
private

Definition at line 608 of file SemiAnalyticalModel.cxx.

References d, and phot::fast_acos().

Referenced by detectedReflectedVisibilities(), Rectangle_SolidAngle(), VISVisibility(), and VUVVisibility().

611  {
612  double aa = a / (2. * d);
613  double bb = b / (2. * d);
614  double aux = (1. + aa * aa + bb * bb) / ((1. + aa * aa) * (1. + bb * bb));
615  return 4. * fast_acos(std::sqrt(aux));
616  }
double fast_acos(double xin)
Float_t d
Definition: plot.C:235
double phot::SemiAnalyticalModel::Rectangle_SolidAngle ( Dims const &  o,
geo::Vector_t const &  v,
const int  OpDetOrientation 
) const
private

Definition at line 618 of file SemiAnalyticalModel.cxx.

References util::abs(), phot::SemiAnalyticalModel::Dims::h, phot::isApproximatelyZero(), phot::isDefinitelyGreaterThan(), Rectangle_SolidAngle(), and phot::SemiAnalyticalModel::Dims::w.

621  {
622  // v is the position of the track segment with respect to
623  // the center position of the arapuca window
624 
625  // solid angle calculation depends on orientation of PD, set correct distances
626  // to use
627  double d1;
628  double d2;
629  if (OpDetOrientation == 2) {
630  // lateral PD, arapuca plane fixed in z direction
631  d1 = std::abs(v.X());
632  d2 = std::abs(v.Z());
633  }
634  else if (OpDetOrientation == 1) {
635  // lateral PD, arapuca plane fixed in y direction
636  d1 = std::abs(v.X());
637  d2 = std::abs(v.Y());
638  }
639  else {
640  // anode/cathode PD, arapuca plane fixed in x direction [default]
641  d1 = std::abs(v.Y());
642  d2 = std::abs(v.X());
643  }
644  // arapuca plane fixed in x direction
645  if (isApproximatelyZero(d1) && isApproximatelyZero(v.Z())) {
646  return Rectangle_SolidAngle(o.h, o.w, d2);
647  }
648  if (isDefinitelyGreaterThan(d1, o.h * .5) &&
649  isDefinitelyGreaterThan(std::abs(v.Z()), o.w * .5)) {
650  double A = d1 - o.h * .5;
651  double B = std::abs(v.Z()) - o.w * .5;
652  double to_return = (Rectangle_SolidAngle(2. * (A + o.h), 2. * (B + o.w), d2) -
653  Rectangle_SolidAngle(2. * A, 2. * (B + o.w), d2) -
654  Rectangle_SolidAngle(2. * (A + o.h), 2. * B, d2) +
655  Rectangle_SolidAngle(2. * A, 2. * B, d2)) *
656  .25;
657  return to_return;
658  }
659  if ((d1 <= o.h * .5) && (std::abs(v.Z()) <= o.w * .5)) {
660  double A = -d1 + o.h * .5;
661  double B = -std::abs(v.Z()) + o.w * .5;
662  double to_return = (Rectangle_SolidAngle(2. * (o.h - A), 2. * (o.w - B), d2) +
663  Rectangle_SolidAngle(2. * A, 2. * (o.w - B), d2) +
664  Rectangle_SolidAngle(2. * (o.h - A), 2. * B, d2) +
665  Rectangle_SolidAngle(2. * A, 2. * B, d2)) *
666  .25;
667  return to_return;
668  }
669  if (isDefinitelyGreaterThan(d1, o.h * .5) && (std::abs(v.Z()) <= o.w * .5)) {
670  double A = d1 - o.h * .5;
671  double B = -std::abs(v.Z()) + o.w * .5;
672  double to_return = (Rectangle_SolidAngle(2. * (A + o.h), 2. * (o.w - B), d2) -
673  Rectangle_SolidAngle(2. * A, 2. * (o.w - B), d2) +
674  Rectangle_SolidAngle(2. * (A + o.h), 2. * B, d2) -
675  Rectangle_SolidAngle(2. * A, 2. * B, d2)) *
676  .25;
677  return to_return;
678  }
679  if ((d1 <= o.h * .5) && isDefinitelyGreaterThan(std::abs(v.Z()), o.w * .5)) {
680  double A = -d1 + o.h * .5;
681  double B = std::abs(v.Z()) - o.w * .5;
682  double to_return = (Rectangle_SolidAngle(2. * (o.h - A), 2. * (B + o.w), d2) -
683  Rectangle_SolidAngle(2. * (o.h - A), 2. * B, d2) +
684  Rectangle_SolidAngle(2. * A, 2. * (B + o.w), d2) -
685  Rectangle_SolidAngle(2. * A, 2. * B, d2)) *
686  .25;
687  return to_return;
688  }
689 
690  return 0.;
691  }
double Rectangle_SolidAngle(const double a, const double b, const double d) const
constexpr auto abs(T v)
Returns the absolute value of the argument.
static constexpr bool isDefinitelyGreaterThan(TReal a, TReal b, TReal tolerance=std::numeric_limits< TReal >::epsilon())
static constexpr bool isApproximatelyZero(TReal a, TReal tolerance=std::numeric_limits< TReal >::epsilon())
double phot::SemiAnalyticalModel::VISVisibility ( geo::Point_t const &  ScintPoint,
OpticalDetector const &  opDet,
const double  cathode_visibility,
geo::Point_t const &  hotspot,
bool  AnodeMode = false 
) const
private

Definition at line 406 of file SemiAnalyticalModel.cxx.

References util::abs(), phot::SemiAnalyticalModel::OpticalDetector::center, Disk_SolidAngle(), fanode_plane_depth, fAnodeReflectivity, phot::fast_acos(), fcathode_centre, fdelta_angulo_vis, fIsDomePDCorr, fIsFlatPDCorr, fIsFlatPDCorrLat, fplane_depth, fradius, fvis_distances_r_dome, fvis_distances_r_flat, fvis_distances_r_flat_lateral, fvis_distances_x_dome, fvis_distances_x_flat, fvis_distances_x_flat_lateral, fvispars_dome, fvispars_flat, fvispars_flat_lateral, phot::SemiAnalyticalModel::OpticalDetector::h, phot::interpolate2(), Omega_Dome_Model(), phot::SemiAnalyticalModel::OpticalDetector::orientation, util::pi(), r, Rectangle_SolidAngle(), phot::SemiAnalyticalModel::OpticalDetector::type, and phot::SemiAnalyticalModel::OpticalDetector::w.

Referenced by detectedReflectedVisibilities().

411  {
412  // set correct plane_depth
413  double plane_depth;
414  if (AnodeMode)
415  plane_depth = fanode_plane_depth;
416  else
417  plane_depth = ScintPoint.X() < 0. ? -fplane_depth : fplane_depth;
418 
419  // calculate visibility of the optical
420  // detector from the hotspot using solid angle:
421 
422  geo::Vector_t const emission_relative = hotspot - opDet.center;
423 
424  // calculate distances and angles for application of corrections
425  // distance from hotspot to optical detector
426  const double distance_vis = emission_relative.R();
427  // angle between hotspot and optical detector
428  double cosine_vis;
429  if (opDet.orientation == 2) { // lateral fixed at z
430  cosine_vis = std::abs(emission_relative.Z()) / distance_vis;
431  }
432  else if (opDet.orientation == 1) { // lateral fixed at y
433  cosine_vis = std::abs(emission_relative.Y()) / distance_vis;
434  }
435  else { // anode/cathode (default)
436  cosine_vis = std::abs(emission_relative.X()) / distance_vis;
437  }
438  const double theta_vis = fast_acos(cosine_vis) * 180. / CLHEP::pi;
439 
440  // calculate solid angle of optical channel
441  double solid_angle_detector = 0.;
442  // ARAPUCAS/Bars (rectangle)
443  if (opDet.type == 0) {
444  // get hotspot coordinates relative to opDet
445  geo::Vector_t const abs_emission_relative{std::abs(emission_relative.X()),
446  std::abs(emission_relative.Y()),
447  std::abs(emission_relative.Z())};
448  solid_angle_detector =
449  Rectangle_SolidAngle(Dims{opDet.h, opDet.w}, abs_emission_relative, opDet.orientation);
450  }
451  // PMTS (dome)
452  else if (opDet.type == 1) {
453  solid_angle_detector = Omega_Dome_Model(distance_vis, theta_vis);
454  }
455  // PMTs (disk)
456  else if (opDet.type == 2) {
457  const double zy_offset = std::sqrt(emission_relative.Y() * emission_relative.Y() +
458  emission_relative.Z() * emission_relative.Z());
459  const double x_distance = std::abs(emission_relative.X());
460  solid_angle_detector = Disk_SolidAngle(zy_offset, x_distance, fradius);
461  }
462  else {
463  throw cet::exception("SemiAnalyticalModel")
464  << "Error: Invalid optical detector shape requested - configuration "
465  "error in semi-analytical model, only rectangular, dome or disk "
466  "optical detectors are supported."
467  << "\n";
468  }
469 
470  // calculate number of hits via geometeric acceptance
471  double visibility_geo = (solid_angle_detector / (2. * CLHEP::pi)) *
472  cathode_visibility; // 2*pi due to presence of reflective foils
473 
474  // determine correction factor, depending on PD type
475  const size_t k = (theta_vis / fdelta_angulo_vis); // off-set angle bin
476  double r;
477  if ((opDet.orientation == 1 && fIsFlatPDCorrLat) ||
478  (opDet.orientation == 2 && fIsFlatPDCorrLat))
479  r = std::hypot(ScintPoint.X() - fcathode_centre[0], ScintPoint.Z() - fcathode_centre[2]);
480  else
481  r = std::hypot(ScintPoint.Y() - fcathode_centre[1], ScintPoint.Z() - fcathode_centre[2]);
482  double d_c = std::abs(ScintPoint.X() - plane_depth); // distance to cathode
483  double border_correction = 0;
484  // flat PDs
485  if ((opDet.type == 0 || opDet.type == 2) && (fIsFlatPDCorr || fIsFlatPDCorrLat)) {
486  // cathode/anode case
487  if (opDet.orientation == 0 && fIsFlatPDCorr) {
488 
489  border_correction =
491  }
492  // laterals case
493  else if ((opDet.orientation == 1 && fIsFlatPDCorrLat) ||
494  (opDet.orientation == 2 && fIsFlatPDCorrLat)) {
495  border_correction = interpolate2(fvis_distances_x_flat_lateral,
498  d_c,
499  r,
500  k);
501  }
502  else {
503  throw cet::exception("SemiAnalyticalModel")
504  << "Error: Invalid optical detector shape requested or corrections "
505  "are missing - configuration error in semi-analytical model."
506  << "\n";
507  }
508  }
509  // dome PDs
510  else if (opDet.type == 1 && fIsDomePDCorr) {
511  border_correction =
513  }
514  else {
515  throw cet::exception("SemiAnalyticalModel")
516  << "Error: Invalid optical detector shape requested or corrections are "
517  "missing - configuration error in semi-analytical model."
518  << "\n";
519  }
520 
521  // apply anode reflectivity factor
522  if (AnodeMode) border_correction *= fAnodeReflectivity;
523 
524  return border_correction * visibility_geo / cosine_vis;
525  }
std::vector< double > fvis_distances_x_dome
TRandom r
Definition: spectrum.C:23
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:160
double Rectangle_SolidAngle(const double a, const double b, const double d) const
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::vector< std::vector< std::vector< double > > > fvispars_dome
std::vector< double > fvis_distances_r_flat
double fast_acos(double xin)
double interpolate2(const std::vector< double > &xDistances, const std::vector< double > &rDistances, const std::vector< std::vector< std::vector< double >>> &parameters, const double x, const double r, const size_t k)
std::vector< double > fvis_distances_r_dome
std::vector< std::vector< std::vector< double > > > fvispars_flat_lateral
std::vector< double > fvis_distances_x_flat_lateral
double Omega_Dome_Model(const double distance, const double theta) const
double Disk_SolidAngle(const double d, const double h, const double b) const
std::vector< double > fvis_distances_x_flat
std::vector< std::vector< std::vector< double > > > fvispars_flat
std::vector< double > fvis_distances_r_flat_lateral
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
double phot::SemiAnalyticalModel::VUVAbsorptionLength ( ) const
private

Definition at line 149 of file SemiAnalyticalModel.cxx.

References fUseXeAbsorption, and phot::interpolate().

Referenced by SemiAnalyticalModel().

150  {
151  // determine LAr absorption length in cm
152  std::map<double, double> abs_length_spectrum =
153  lar::providerFrom<detinfo::LArPropertiesService>()->AbsLengthSpectrum();
154  std::vector<double> x_v, y_v;
155  for (auto elem : abs_length_spectrum) {
156  x_v.push_back(elem.first);
157  y_v.push_back(elem.second);
158  }
159  double vuv_absorption_length;
160  if (fUseXeAbsorption)
161  vuv_absorption_length =
162  interpolate(x_v, y_v, 7.1, false); // 7.1 eV: peak of Xe VUV emission spectrum
163  else
164  vuv_absorption_length =
165  interpolate(x_v, y_v, 9.7, false); // 9.7 eV: peak of Ar VUV emission spectrum
166  if (vuv_absorption_length <= 0) {
167  throw cet::exception("SemiAnalyticalModel")
168  << "Error: VUV Absorption Length is 0 or negative.\n";
169  }
170  return vuv_absorption_length;
171  }
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, const double x, const bool extrapolate, size_t i)
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
double phot::SemiAnalyticalModel::VUVVisibility ( geo::Point_t const &  ScintPoint,
OpticalDetector const &  opDet 
) const
private

Definition at line 190 of file SemiAnalyticalModel.cxx.

References util::abs(), phot::SemiAnalyticalModel::OpticalDetector::center, Disk_SolidAngle(), phot::fast_acos(), fborder_corr_angulo_dome, fborder_corr_angulo_flat, fborder_corr_angulo_flat_lateral, fborder_corr_dome, fborder_corr_flat, fborder_corr_flat_lateral, fcathode_centre, fdelta_angulo_vuv, fGHvuvpars_dome, fGHvuvpars_flat, fGHvuvpars_flat_lateral, fIsDomePDCorr, fIsFlatPDCorr, fIsFlatPDCorrLat, fradius, fvuv_absorption_length, Gaisser_Hillas(), phot::SemiAnalyticalModel::OpticalDetector::h, phot::interpolate(), Omega_Dome_Model(), phot::SemiAnalyticalModel::OpticalDetector::orientation, util::pi(), r, Rectangle_SolidAngle(), phot::SemiAnalyticalModel::OpticalDetector::type, and phot::SemiAnalyticalModel::OpticalDetector::w.

Referenced by detectedDirectVisibilities().

192  {
193  // distance and angle between ScintPoint and OpDet center
194  geo::Vector_t const relative = ScintPoint - opDet.center;
195  const double distance = relative.R();
196  double cosine;
197  if (opDet.orientation == 2)
198  cosine = std::abs(relative.Z()) / distance;
199  else if (opDet.orientation == 1)
200  cosine = std::abs(relative.Y()) / distance;
201  else
202  cosine = std::abs(relative.X()) / distance;
203  const double theta = fast_acos(cosine) * 180. / CLHEP::pi;
204 
205  double solid_angle = 0.;
206  // ARAPUCAS/Bars (rectangle)
207  if (opDet.type == 0) {
208  // get scintillation point coordinates relative to arapuca window centre
209  geo::Vector_t const abs_relative{
210  std::abs(relative.X()), std::abs(relative.Y()), std::abs(relative.Z())};
211  solid_angle = Rectangle_SolidAngle(Dims{opDet.h, opDet.w}, abs_relative, opDet.orientation);
212  }
213  // PMTs (dome)
214  else if (opDet.type == 1) {
215  solid_angle = Omega_Dome_Model(distance, theta);
216  }
217  // PMTs (disk)
218  else if (opDet.type == 2) {
219  const double zy_offset = std::sqrt(relative.Y() * relative.Y() + relative.Z() * relative.Z());
220  const double x_distance = std::abs(relative.X());
221  solid_angle = Disk_SolidAngle(zy_offset, x_distance, fradius);
222  }
223  else {
224  throw cet::exception("SemiAnalyticalModel")
225  << "Error: Invalid optical detector shape requested - configuration "
226  "error in semi-analytical model, only rectangular, dome or disk "
227  "optical detectors are supported."
228  << "\n";
229  }
230 
231  // calculate visibility by geometric acceptance
232  // accounting for solid angle and LAr absorbtion length
233  double visibility_geo =
234  std::exp(-1. * distance / fvuv_absorption_length) * (solid_angle / (4 * CLHEP::pi));
235 
236  // apply Gaisser-Hillas correction for Rayleigh scattering distance
237  // and angular dependence offset angle bin
238  const size_t j = (theta / fdelta_angulo_vuv);
239 
240  // determine GH parameters, accounting for border effects
241  // radial distance from centre of detector (Y-Z standard / X-Z laterals)
242  double r = 0;
243  if (opDet.orientation == 2)
244  r = std::hypot(ScintPoint.X() - fcathode_centre[0], ScintPoint.Y() - fcathode_centre[1]);
245  else if (opDet.orientation == 1)
246  r = std::hypot(ScintPoint.X() - fcathode_centre[0], ScintPoint.Z() - fcathode_centre[2]);
247  else
248  r = std::hypot(ScintPoint.Y() - fcathode_centre[1], ScintPoint.Z() - fcathode_centre[2]);
249 
250  double pars_ini[4] = {0, 0, 0, 0};
251  double s1 = 0;
252  double s2 = 0;
253  double s3 = 0;
254  // flat PDs
255  if ((opDet.type == 0 || opDet.type == 2) && (fIsFlatPDCorr || fIsFlatPDCorrLat)) {
256  if ((opDet.orientation == 1 && fIsFlatPDCorrLat) ||
257  (opDet.orientation == 2 && fIsFlatPDCorrLat)) { // laterals
258  pars_ini[0] = fGHvuvpars_flat_lateral[0][j];
259  pars_ini[1] = fGHvuvpars_flat_lateral[1][j];
260  pars_ini[2] = fGHvuvpars_flat_lateral[2][j];
261  pars_ini[3] = fGHvuvpars_flat_lateral[3][j];
262  s1 =
264  s2 =
266  s3 =
268  }
269  else if (opDet.orientation == 0 && fIsFlatPDCorr) { // cathode/anode
270  pars_ini[0] = fGHvuvpars_flat[0][j];
271  pars_ini[1] = fGHvuvpars_flat[1][j];
272  pars_ini[2] = fGHvuvpars_flat[2][j];
273  pars_ini[3] = fGHvuvpars_flat[3][j];
277  }
278  else {
279  throw cet::exception("SemiAnalyticalModel")
280  << "Error: flat optical detectors are found, but parameters are "
281  "missing - configuration error in semi-analytical model."
282  << "\n";
283  }
284  }
285  // dome PDs
286  else if (opDet.type == 1 && fIsDomePDCorr) {
287  pars_ini[0] = fGHvuvpars_dome[0][j];
288  pars_ini[1] = fGHvuvpars_dome[1][j];
289  pars_ini[2] = fGHvuvpars_dome[2][j];
290  pars_ini[3] = fGHvuvpars_dome[3][j];
294  }
295  else {
296  throw cet::exception("SemiAnalyticalModel")
297  << "Error: Invalid optical detector shape requested or corrections are "
298  "missing - configuration error in semi-analytical model."
299  << "\n";
300  }
301 
302  // add border correction to parameters
303  pars_ini[0] = pars_ini[0] + s1 * r;
304  pars_ini[1] = pars_ini[1] + s2 * r;
305  pars_ini[2] = pars_ini[2] + s3 * r;
306  pars_ini[3] = pars_ini[3];
307 
308  // calculate correction
309  double GH_correction = Gaisser_Hillas(distance, pars_ini);
310 
311  // determine corrected visibility of photo-detector
312  return GH_correction * visibility_geo / cosine;
313  }
TRandom r
Definition: spectrum.C:23
std::vector< std::vector< double > > fborder_corr_flat_lateral
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:160
double Rectangle_SolidAngle(const double a, const double b, const double d) const
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::vector< std::vector< double > > fGHvuvpars_flat_lateral
std::vector< std::vector< double > > fborder_corr_dome
double fast_acos(double xin)
std::vector< double > fborder_corr_angulo_flat
std::vector< double > fborder_corr_angulo_dome
double Gaisser_Hillas(const double x, const double *par) const
std::vector< std::vector< double > > fGHvuvpars_flat
double Omega_Dome_Model(const double distance, const double theta) const
std::vector< std::vector< double > > fGHvuvpars_dome
double Disk_SolidAngle(const double d, const double h, const double b) const
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
std::vector< std::vector< double > > fborder_corr_flat
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, const double x, const bool extrapolate, size_t i)
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< double > fborder_corr_angulo_flat_lateral

Member Data Documentation

const std::vector<geo::BoxBoundedGeo> phot::SemiAnalyticalModel::fActiveVolumes
private

Definition at line 101 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel().

const TVector3 phot::SemiAnalyticalModel::fanode_centre
private

Definition at line 102 of file SemiAnalyticalModel.h.

Referenced by detectedReflectedVisibilities(), and SemiAnalyticalModel().

Dims phot::SemiAnalyticalModel::fanode_plane
private

Definition at line 110 of file SemiAnalyticalModel.h.

Referenced by detectedReflectedVisibilities(), and SemiAnalyticalModel().

double phot::SemiAnalyticalModel::fanode_plane_depth
private
double phot::SemiAnalyticalModel::fAnodeReflectivity
private

Definition at line 140 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VISVisibility().

bool phot::SemiAnalyticalModel::fApplyFieldCageTransparency
private

Definition at line 131 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel().

std::vector<double> phot::SemiAnalyticalModel::fborder_corr_angulo_dome
private

Definition at line 128 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VUVVisibility().

std::vector<double> phot::SemiAnalyticalModel::fborder_corr_angulo_flat
private
std::vector<double> phot::SemiAnalyticalModel::fborder_corr_angulo_flat_lateral
private

Definition at line 122 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VUVVisibility().

std::vector<std::vector<double> > phot::SemiAnalyticalModel::fborder_corr_dome
private

Definition at line 129 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VUVVisibility().

std::vector<std::vector<double> > phot::SemiAnalyticalModel::fborder_corr_flat
private
std::vector<std::vector<double> > phot::SemiAnalyticalModel::fborder_corr_flat_lateral
private

Definition at line 123 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VUVVisibility().

const TVector3 phot::SemiAnalyticalModel::fcathode_centre
private
Dims phot::SemiAnalyticalModel::fcathode_plane
private

Definition at line 109 of file SemiAnalyticalModel.h.

Referenced by detectedReflectedVisibilities(), and SemiAnalyticalModel().

double phot::SemiAnalyticalModel::fdelta_angulo_vis
private

Definition at line 139 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VISVisibility().

double phot::SemiAnalyticalModel::fdelta_angulo_vuv
private

Definition at line 113 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VUVVisibility().

const bool phot::SemiAnalyticalModel::fDoReflectedLight
private

Definition at line 136 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel().

double phot::SemiAnalyticalModel::fFieldCageTransparencyCathode
private

Definition at line 133 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel().

double phot::SemiAnalyticalModel::fFieldCageTransparencyLateral
private

Definition at line 132 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel().

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

Definition at line 98 of file SemiAnalyticalModel.h.

Referenced by opticalDetectors(), and SemiAnalyticalModel().

std::vector<std::vector<double> > phot::SemiAnalyticalModel::fGHvuvpars_dome
private

Definition at line 127 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VUVVisibility().

std::vector<std::vector<double> > phot::SemiAnalyticalModel::fGHvuvpars_flat
private
std::vector<std::vector<double> > phot::SemiAnalyticalModel::fGHvuvpars_flat_lateral
private

Definition at line 121 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VUVVisibility().

const bool phot::SemiAnalyticalModel::fIncludeAnodeReflections
private

Definition at line 137 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel().

bool phot::SemiAnalyticalModel::fIsDomePDCorr
private

Definition at line 126 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), VISVisibility(), and VUVVisibility().

bool phot::SemiAnalyticalModel::fIsFlatPDCorr
private
bool phot::SemiAnalyticalModel::fIsFlatPDCorrLat
private

Definition at line 120 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), VISVisibility(), and VUVVisibility().

const larg4::ISTPC phot::SemiAnalyticalModel::fISTPC
private

Definition at line 99 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel().

const size_t phot::SemiAnalyticalModel::fNOpDets
private
const int phot::SemiAnalyticalModel::fNTPC
private

Definition at line 100 of file SemiAnalyticalModel.h.

Referenced by isOpDetInSameTPC(), and SemiAnalyticalModel().

const std::vector<OpticalDetector> phot::SemiAnalyticalModel::fOpDetector
private
double phot::SemiAnalyticalModel::fplane_depth
private
double phot::SemiAnalyticalModel::fradius
private
const bool phot::SemiAnalyticalModel::fUseXeAbsorption
private

Definition at line 155 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VUVAbsorptionLength().

std::vector<double> phot::SemiAnalyticalModel::fvis_distances_r_dome
private

Definition at line 151 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VISVisibility().

std::vector<double> phot::SemiAnalyticalModel::fvis_distances_r_flat
private

Definition at line 143 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VISVisibility().

std::vector<double> phot::SemiAnalyticalModel::fvis_distances_r_flat_lateral
private

Definition at line 147 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VISVisibility().

std::vector<double> phot::SemiAnalyticalModel::fvis_distances_x_dome
private

Definition at line 150 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VISVisibility().

std::vector<double> phot::SemiAnalyticalModel::fvis_distances_x_flat
private

Definition at line 142 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VISVisibility().

std::vector<double> phot::SemiAnalyticalModel::fvis_distances_x_flat_lateral
private

Definition at line 146 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VISVisibility().

std::vector<std::vector<std::vector<double> > > phot::SemiAnalyticalModel::fvispars_dome
private

Definition at line 152 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VISVisibility().

std::vector<std::vector<std::vector<double> > > phot::SemiAnalyticalModel::fvispars_flat
private

Definition at line 144 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VISVisibility().

std::vector<std::vector<std::vector<double> > > phot::SemiAnalyticalModel::fvispars_flat_lateral
private

Definition at line 148 of file SemiAnalyticalModel.h.

Referenced by SemiAnalyticalModel(), and VISVisibility().

double phot::SemiAnalyticalModel::fvuv_absorption_length
private

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