All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Classes | Public Member Functions | Private Member Functions | Static Private Member Functions | Private Attributes | Static Private Attributes | List of all members
SemiAnalyticalModel Class Reference

#include <SemiAnalyticalModel.h>

Classes

struct  Dims
 
struct  OpticalDetector
 

Public Member Functions

 SemiAnalyticalModel (fhicl::ParameterSet VUVHits, fhicl::ParameterSet VISHits, bool doReflectedLight=false, bool includeAnodeReflections=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

void Initialization ()
 
int 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 double 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
 
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
 
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) const
 
 if (diff<=tolerance) return true
 
 if (diff< std::fmax(std::fabs(a), std::fabs(b))*tolerance) return true
 
 if (diff< tolerance) return true
 
 if (diff< std::fmax(std::fabs(a), std::fabs(b))*tolerance) return true
 
 if (diff > tolerance) return true
 
 if (diff > std::fmax(std::fabs(a), std::fabs(b))*tolerance) return true
 

Static Private Member Functions

template<typename TReal >
static constexpr bool
 
template<typename TReal >
static constexpr bool
 
template<typename TReal >
static constexpr bool
 

Private Attributes

return false
 
const fhicl::ParameterSet fVUVHitsParams
 
const fhicl::ParameterSet fVISHitsParams
 
const larg4::ISTPC fISTPC
 
geo::GeometryCore const & fGeom
 
const int fNTPC
 
const std::vector
< geo::BoxBoundedGeo > 		fActiveVolumes
 
const TVector3 fcathode_centre
 
const TVector3 fanode_centre
 
double fplane_depth
 
double fanode_plane_depth
 
const size_t nOpDets
 
double fradius
 
Dims fcathode_plane
 
Dims fanode_plane
 
std::vector< geo::Point_t > fOpDetCenter
 
std::vector< int > fOpDetType
 
std::vector< int > fOpDetOrientation
 
std::vector< double > fOpDetLength
 
std::vector< double > fOpDetHeight
 
const int fvuv_absorption_length
 
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
 

Static Private Attributes

template<typename TReal >
static constexpr bool
 

Detailed Description

Definition at line 31 of file SemiAnalyticalModel.h.

Constructor & Destructor Documentation

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

Definition at line 25 of file SemiAnalyticalModel.cxx.

26  :
27  fVUVHitsParams(VUVHits)
28  , fVISHitsParams(VISHits)
29  , fISTPC{*(lar::providerFrom<geo::Geometry>())}
30  , fGeom{*(lar::providerFrom<geo::Geometry>())}
31  , fNTPC(fGeom.NTPC())
32  , fActiveVolumes(fISTPC.extractActiveLArVolume(fGeom))
33  , fcathode_centre{fGeom.TPC(0, 0).GetCathodeCenter().X(),
34  fActiveVolumes[0].CenterY(),
35  fActiveVolumes[0].CenterZ()}
36  , fanode_centre{fGeom.TPC(0, 0).FirstPlane().GetCenter().X(),
37  fActiveVolumes[0].CenterY(),
38  fActiveVolumes[0].CenterZ()}
39  , nOpDets(fGeom.NOpDets())
40  , fvuv_absorption_length(VUVAbsorptionLength())
41  , fDoReflectedLight(doReflectedLight)
42  , fIncludeAnodeReflections(includeAnodeReflections)
43 {
44  // initialise parameters and geometry
45  mf::LogInfo("SemiAnalyticalModel") << "Semi-analytical model initialized." << std::endl;
46  Initialization();
47 }
larg4::ISTPC::extractActiveLArVolume
static std::vector< geo::BoxBoundedGeo > extractActiveLArVolume(geo::GeometryCore const &geom)
Definition: ISTPC.cxx:51
SemiAnalyticalModel::fcathode_centre
const TVector3 fcathode_centre
Definition: SemiAnalyticalModel.h:160
geo::TPCGeo::GetCathodeCenter
Point GetCathodeCenter() const
Definition: TPCGeo.h:298
SemiAnalyticalModel::fGeom
geo::GeometryCore const & fGeom
Definition: SemiAnalyticalModel.h:157
geo::TPCGeo::FirstPlane
geo::PlaneGeo const & FirstPlane() const
Returns the first wire plane (the closest to TPC center).
Definition: TPCGeo.h:245
SemiAnalyticalModel::fanode_centre
const TVector3 fanode_centre
Definition: SemiAnalyticalModel.h:160
SemiAnalyticalModel::fVUVHitsParams
const fhicl::ParameterSet fVUVHitsParams
Definition: SemiAnalyticalModel.h:152
SemiAnalyticalModel::fVISHitsParams
const fhicl::ParameterSet fVISHitsParams
Definition: SemiAnalyticalModel.h:153
geo::PlaneGeo::GetCenter
Point GetCenter() const
Returns the centre of the wire plane in world coordinates [cm].
Definition: PlaneGeo.h:479
SemiAnalyticalModel::fISTPC
const larg4::ISTPC fISTPC
Definition: SemiAnalyticalModel.h:156
geo::GeometryCore::NOpDets
unsigned int NOpDets() const
Number of OpDets in the whole detector.
Definition: GeometryCore.cxx:244
geo::GeometryCore::NTPC
unsigned int NTPC(unsigned int cstat=0) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:2311
geo::GeometryCore::TPC
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
Definition: GeometryCore.h:2431
SemiAnalyticalModel::fActiveVolumes
const std::vector< geo::BoxBoundedGeo > fActiveVolumes
Definition: SemiAnalyticalModel.h:159
SemiAnalyticalModel::fIncludeAnodeReflections
const bool fIncludeAnodeReflections
Definition: SemiAnalyticalModel.h:201

Member Function Documentation

template<typename TReal >
static constexpr SemiAnalyticalModel::bool ( )
inlinestaticprivate

Definition at line 124 of file SemiAnalyticalModel.h.

template<typename TReal >
static constexpr SemiAnalyticalModel::bool ( )
inlinestaticprivate

Definition at line 134 of file SemiAnalyticalModel.h.

template<typename TReal >
static constexpr SemiAnalyticalModel::bool ( )
inlinestaticprivate

Definition at line 144 of file SemiAnalyticalModel.h.

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

Definition at line 189 of file SemiAnalyticalModel.cxx.

191 {
192  DetectedVisibilities.resize(nOpDets);
193  for (size_t const OpDet : util::counter(nOpDets)) {
194  if (!isOpDetInSameTPC(ScintPoint, fOpDetCenter[OpDet])) {
195  DetectedVisibilities[OpDet] = 0.;
196  continue;
197  }
198 
199  // set detector struct for solid angle function
200  const SemiAnalyticalModel::OpticalDetector op{
201  fOpDetHeight[OpDet], fOpDetLength[OpDet],
202  fOpDetCenter[OpDet], fOpDetType[OpDet], fOpDetOrientation[OpDet]};
203 
204  DetectedVisibilities[OpDet] = VUVVisibility(ScintPoint, op);;
205  }
206 }
SemiAnalyticalModel::VUVVisibility
double VUVVisibility(geo::Point_t const &ScintPoint, OpticalDetector const &opDet) const
Definition: SemiAnalyticalModel.cxx:209
SemiAnalyticalModel::fOpDetHeight
std::vector< double > fOpDetHeight
Definition: SemiAnalyticalModel.h:172
util::counter
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
SemiAnalyticalModel::fOpDetCenter
std::vector< geo::Point_t > fOpDetCenter
Definition: SemiAnalyticalModel.h:168
SemiAnalyticalModel::fOpDetLength
std::vector< double > fOpDetLength
Definition: SemiAnalyticalModel.h:171
SemiAnalyticalModel::fOpDetType
std::vector< int > fOpDetType
Definition: SemiAnalyticalModel.h:169
SemiAnalyticalModel::isOpDetInSameTPC
bool isOpDetInSameTPC(geo::Point_t const &ScintPoint, geo::Point_t const &OpDetPoint) const
Definition: SemiAnalyticalModel.cxx:688
SemiAnalyticalModel::fOpDetOrientation
std::vector< int > fOpDetOrientation
Definition: SemiAnalyticalModel.h:170
void SemiAnalyticalModel::detectedReflectedVisibilities ( std::vector< double > &  ReflDetectedVisibilities,
geo::Point_t const &  ScintPoint,
bool  AnodeMode = false 
) const

Definition at line 318 of file SemiAnalyticalModel.cxx.

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 dimensions
327  geo::Vector_t ScintPoint_relative;
328  Dims plane_dimensions;
329  double plane_depth;
330  if (AnodeMode) {
331  plane_dimensions = fanode_plane;
332  plane_depth = fanode_plane_depth;
333  ScintPoint_relative.SetCoordinates(std::abs(ScintPoint.X() - fanode_plane_depth), std::abs(ScintPoint.Y() - fanode_centre[1]), std::abs(ScintPoint.Z() - fanode_centre[2]));
334  }
335  else {
336  plane_dimensions = fcathode_plane;
337  plane_depth = ScintPoint.X() < 0. ? -fplane_depth : fplane_depth;
338  ScintPoint_relative.SetCoordinates(std::abs(ScintPoint.X() - plane_depth), std::abs(ScintPoint.Y() - fcathode_centre[1]), std::abs(ScintPoint.Z() - fcathode_centre[2]));
339  }
340 
341  // calculate solid angle of anode/cathode from the scintillation point, orientation always = 0 (anode/cathode)
342  double solid_angle_cathode = Rectangle_SolidAngle(plane_dimensions, ScintPoint_relative, 0);
343 
344  // calculate distance and angle between ScintPoint and hotspot
345  // vast majority of hits in hotspot region directly infront of scintpoint,
346  // therefore consider attenuation for this distance and on axis GH instead of for the centre coordinate
347  double distance_cathode = std::abs(plane_depth - ScintPoint.X());
348  // calculate hits on cathode plane via geometric acceptance
349  double cathode_visibility_geo = std::exp(-1. * distance_cathode / fvuv_absorption_length) * (solid_angle_cathode / (4. * CLHEP::pi));
350 
351  // determine Gaisser-Hillas correction including border effects
352  // use flat correction
353  double r = std::hypot(ScintPoint.Y() - fcathode_centre[1], ScintPoint.Z() - fcathode_centre[2]);
354  double pars_ini[4] = {0, 0, 0, 0};
355  double s1 = 0; double s2 = 0; double s3 = 0;
356  if(fIsFlatPDCorr) {
357  pars_ini[0] = fGHvuvpars_flat[0][0];
358  pars_ini[1] = fGHvuvpars_flat[1][0];
359  pars_ini[2] = fGHvuvpars_flat[2][0];
360  pars_ini[3] = fGHvuvpars_flat[3][0];
361  s1 = interpolate( fborder_corr_angulo_flat, fborder_corr_flat[0], 0, true);
362  s2 = interpolate( fborder_corr_angulo_flat, fborder_corr_flat[1], 0, true);
363  s3 = interpolate( fborder_corr_angulo_flat, fborder_corr_flat[2], 0, true);
364  }
365  else {
366  throw cet::exception("SemiAnalyticalModel")
367  << "Error: flat optical detector VUV correction required for reflected semi-analytic hits. - configuration error in semi-analytical model." << "\n";
368  }
369 
370  // add border correction
371  pars_ini[0] = pars_ini[0] + s1 * r;
372  pars_ini[1] = pars_ini[1] + s2 * r;
373  pars_ini[2] = pars_ini[2] + s3 * r;
374  pars_ini[3] = pars_ini[3];
375 
376  // calculate corrected number of hits
377  double GH_correction = Gaisser_Hillas(distance_cathode, pars_ini);
378  const double cathode_visibility_rec = GH_correction * cathode_visibility_geo;
379 
380  // 2). detemine visibility of each PD
381  const geo::Point_t hotspot = {plane_depth, ScintPoint.Y(), ScintPoint.Z()};
382  ReflDetectedVisibilities.resize(nOpDets);
383  for (size_t const OpDet : util::counter(nOpDets)) {
384  if (!isOpDetInSameTPC(ScintPoint, fOpDetCenter[OpDet])) {
385  ReflDetectedVisibilities[OpDet] = 0.;
386  continue;
387  }
388 
389  // set detector struct for solid angle function
390  const OpticalDetector op{
391  fOpDetHeight[OpDet], fOpDetLength[OpDet],
392  fOpDetCenter[OpDet], fOpDetType[OpDet], fOpDetOrientation[OpDet]};
393 
394  ReflDetectedVisibilities[OpDet] = VISVisibility(ScintPoint, op, cathode_visibility_rec, hotspot, AnodeMode);
395  }
396 }
geo::Vector_t
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:164
SemiAnalyticalModel::Gaisser_Hillas
double Gaisser_Hillas(const double x, const double *par) const
Definition: SemiAnalyticalModel.cxx:494
SemiAnalyticalModel::fcathode_centre
const TVector3 fcathode_centre
Definition: SemiAnalyticalModel.h:160
SemiAnalyticalModel::fGHvuvpars_flat
std::vector< std::vector< double > > fGHvuvpars_flat
Definition: SemiAnalyticalModel.h:180
SemiAnalyticalModel::fOpDetHeight
std::vector< double > fOpDetHeight
Definition: SemiAnalyticalModel.h:172
SemiAnalyticalModel::fanode_centre
const TVector3 fanode_centre
Definition: SemiAnalyticalModel.h:160
SemiAnalyticalModel::Rectangle_SolidAngle
double Rectangle_SolidAngle(const double a, const double b, const double d) const
Definition: SemiAnalyticalModel.cxx:575
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
util::counter
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
pi
pdgs pi
Definition: selectors.fcl:22
SemiAnalyticalModel::interpolate
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, double x, bool extrapolate, size_t i=0) const
Definition: SemiAnalyticalModel.cxx:726
SemiAnalyticalModel::fOpDetCenter
std::vector< geo::Point_t > fOpDetCenter
Definition: SemiAnalyticalModel.h:168
SemiAnalyticalModel::fOpDetLength
std::vector< double > fOpDetLength
Definition: SemiAnalyticalModel.h:171
SemiAnalyticalModel::fborder_corr_flat
std::vector< std::vector< double > > fborder_corr_flat
Definition: SemiAnalyticalModel.h:182
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
SemiAnalyticalModel::fOpDetType
std::vector< int > fOpDetType
Definition: SemiAnalyticalModel.h:169
SemiAnalyticalModel::isOpDetInSameTPC
bool isOpDetInSameTPC(geo::Point_t const &ScintPoint, geo::Point_t const &OpDetPoint) const
Definition: SemiAnalyticalModel.cxx:688
geo::Point_t
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
r
esac echo uname r
Definition: condor_lar_pubs.sh:698
SemiAnalyticalModel::fborder_corr_angulo_flat
std::vector< double > fborder_corr_angulo_flat
Definition: SemiAnalyticalModel.h:181
SemiAnalyticalModel::fOpDetOrientation
std::vector< int > fOpDetOrientation
Definition: SemiAnalyticalModel.h:170
SemiAnalyticalModel::VISVisibility
double VISVisibility(geo::Point_t const &ScintPoint, OpticalDetector const &opDet, const double cathode_visibility, geo::Point_t const &hotspot, bool AnodeMode=false) const
Definition: SemiAnalyticalModel.cxx:399
double SemiAnalyticalModel::Disk_SolidAngle ( const double  d,
const double  h,
const double  b 
) const
private

Definition at line 508 of file SemiAnalyticalModel.cxx.

509 {
510  if (b <= 0. || d < 0. || h <= 0.) return 0.;
511  const double leg2 = (b + d) * (b + d);
512  const double aa = std::sqrt(h * h / (h * h + leg2));
513  if (isApproximatelyZero(d)) { return 2. * CLHEP::pi * (1. - aa); }
514  double bb = 2. * std::sqrt(b * d / (h * h + leg2));
515  double cc = 4. * b * d / leg2;
516 
517  if (isDefinitelyGreaterThan(d, b)) {
518  try {
519  return 2. * aa *
520  (std::sqrt(1. - cc) * boost::math::ellint_3(bb, cc, noLDoublePromote()) -
521  boost::math::ellint_1(bb, noLDoublePromote()));
522  }
523  catch (std::domain_error& e) {
524  if (isApproximatelyEqual(d, b, 1e-9)) {
525  mf::LogWarning("SemiAnalyticalModel")
526  << "Elliptic Integral in Disk_SolidAngle() given parameters "
527  "outside domain."
528  << "\nbb: " << bb << "\ncc: " << cc << "\nException message: " << e.what()
529  << "\nRelax condition and carry on.";
530  return CLHEP::pi - 2. * aa * boost::math::ellint_1(bb, noLDoublePromote());
531  }
532  else {
533  mf::LogError("SemiAnalyticalModel")
534  << "Elliptic Integral inside Disk_SolidAngle() given parameters "
535  "outside domain.\n"
536  << "\nbb: " << bb << "\ncc: " << cc << "Exception message: " << e.what();
537  return 0.;
538  }
539  }
540  }
541  if (isDefinitelyLessThan(d, b)) {
542  try {
543  return 2. * CLHEP::pi -
544  2. * aa *
545  (boost::math::ellint_1(bb, noLDoublePromote()) +
546  std::sqrt(1. - cc) * boost::math::ellint_3(bb, cc, noLDoublePromote()));
547  }
548  catch (std::domain_error& e) {
549  if (isApproximatelyEqual(d, b, 1e-9)) {
550  mf::LogWarning("SemiAnalyticalModel")
551  << "Elliptic Integral in Disk_SolidAngle() given parameters "
552  "outside domain."
553  << "\nbb: " << bb << "\ncc: " << cc << "\nException message: " << e.what()
554  << "\nRelax condition and carry on.";
555  return CLHEP::pi - 2. * aa * boost::math::ellint_1(bb, noLDoublePromote());
556  }
557  else {
558  mf::LogError("SemiAnalyticalModel")
559  << "Elliptic Integral inside Disk_SolidAngle() given parameters "
560  "outside domain.\n"
561  << "\nbb: " << bb << "\ncc: " << cc << "Exception message: " << e.what();
562  return 0.;
563  }
564  }
565  }
566  if (isApproximatelyEqual(d, b)) {
567  return CLHEP::pi - 2. * aa * boost::math::ellint_1(bb, noLDoublePromote());
568  }
569  return 0.;
570 }
h
while getopts h
Definition: getStatistics.sh:13
noLDoublePromote
boost::math::policies::policy< boost::math::policies::promote_double< false > > noLDoublePromote
Definition: OpFastScintillation.cxx:160
pi
pdgs pi
Definition: selectors.fcl:22
e
do i e
Definition: updateVersionList.sh:68
double SemiAnalyticalModel::fast_acos ( double  x) const
private

Definition at line 703 of file SemiAnalyticalModel.cxx.

704 {
705  double negate = double(x < 0);
706  x = std::abs(x);
707  x -= double(x > 1.0) * (x - 1.0); // <- equivalent to min(1.0,x), but faster
708  double ret = -0.0187293;
709  ret = ret * x;
710  ret = ret + 0.0742610;
711  ret = ret * x;
712  ret = ret - 0.2121144;
713  ret = ret * x;
714  ret = ret + 1.5707288;
715  ret = ret * std::sqrt(1.0 - x);
716  ret = ret - 2. * negate * ret;
717  return negate * 3.14159265358979 + ret;
718 }
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
double SemiAnalyticalModel::Gaisser_Hillas ( const double  x,
const double *  par 
) const
private

Definition at line 494 of file SemiAnalyticalModel.cxx.

495 {
496  double X_mu_0 = par[3];
497  double Normalization = par[0];
498  double Diff = par[1] - X_mu_0;
499  double Term = std::pow((x - X_mu_0) / Diff, Diff / par[2]);
500  double Exponential = std::exp((par[1] - x) / par[2]);
501 
502  return (Normalization * Term * Exponential);
503 }
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
SemiAnalyticalModel::if ( diff<=  tolerance)
private
SemiAnalyticalModel::if ( )
private
SemiAnalyticalModel::if ( )
private
SemiAnalyticalModel::if ( )
private
SemiAnalyticalModel::if ( diff  ,
tolerance   
)
private
SemiAnalyticalModel::if ( diff  ,
std::fmax(std::fabs(a), std::fabs(b))*  tolerance 
)
private
void SemiAnalyticalModel::Initialization ( )
private

Definition at line 50 of file SemiAnalyticalModel.cxx.

51 {
52  // get PDS information
53  fOpDetType.reserve(nOpDets); fOpDetOrientation.reserve(nOpDets);
54  fOpDetCenter.reserve(nOpDets); fOpDetLength.reserve(nOpDets); fOpDetHeight.reserve(nOpDets);
55  for (size_t const i : util::counter(nOpDets)) {
56 
57  geo::OpDetGeo const& opDet = fGeom.OpDetGeoFromOpDet(i);
58  fOpDetCenter.push_back(opDet.GetCenter());
59 
60  if (opDet.isSphere()) { // dome PMTs
61  fOpDetType.push_back(1); // dome
62  fOpDetOrientation.push_back(0); // anode/cathode (default)
63  fOpDetLength.push_back(-1);
64  fOpDetHeight.push_back(-1);
65  }
66  else if (opDet.isBar()) {
67  fOpDetType.push_back(0); // (X)Arapucas/Bars
68  // determine orientation to get correction OpDet dimensions
69  fOpDetLength.push_back(opDet.Length());
70  if (opDet.Width() > opDet.Height()) { // laterals, Y dimension smallest
71  fOpDetOrientation.push_back(1);
72  fOpDetHeight.push_back(opDet.Width());
73  }
74  else { // anode/cathode (default), X dimension smallest
75  fOpDetOrientation.push_back(0);
76  fOpDetHeight.push_back(opDet.Height());
77  }
78  }
79  else {
80  fOpDetType.push_back(2); // disk PMTs
81  fOpDetOrientation.push_back(0); // anode/cathode (default)
82  fOpDetLength.push_back(-1);
83  fOpDetHeight.push_back(-1);
84  }
85  }
86 
87  // Load Gaisser-Hillas corrections for VUV semi-analytic hits
88  mf::LogInfo("SemiAnalyticalModel") << "Using VUV visibility parameterization";
89 
90  fIsFlatPDCorr = fVUVHitsParams.get<bool>("FlatPDCorr", false);
91  fIsFlatPDCorrLat = fVUVHitsParams.get<bool>("FlatPDCorrLat", false);
92  fIsDomePDCorr = fVUVHitsParams.get<bool>("DomePDCorr", false);
93  fdelta_angulo_vuv = fVUVHitsParams.get<double>("delta_angulo_vuv", 10);
94  fradius = fVUVHitsParams.get<double>("PMT_radius", 10.16);
95  fApplyFieldCageTransparency = fVUVHitsParams.get<bool>("ApplyFieldCageTransparency", false);
96  fFieldCageTransparencyLateral = fVUVHitsParams.get<double>("FieldCageTransparencyLateral", 1.0);
97  fFieldCageTransparencyCathode = fVUVHitsParams.get<double>("FieldCageTransparencyCathode", 1.0);
98 
99  if (!fIsFlatPDCorr && !fIsDomePDCorr && !fIsFlatPDCorrLat) {
100  throw cet::exception("SemiAnalyticalModel")
101  << "Both isFlatPDCorr/isFlatPDCorrLat and isDomePDCorr parameters are false, at least one type of parameterisation is required for the semi-analytic light simulation." << "\n";
102  }
103  if (fIsFlatPDCorr) {
104  fGHvuvpars_flat = fVUVHitsParams.get<std::vector<std::vector<double>>>("GH_PARS_flat");
105  fborder_corr_angulo_flat = fVUVHitsParams.get<std::vector<double>>("GH_border_angulo_flat");
106  fborder_corr_flat = fVUVHitsParams.get<std::vector<std::vector<double>>>("GH_border_flat");
107  }
108  if (fIsFlatPDCorrLat) {
109  fGHvuvpars_flat_lateral = fVUVHitsParams.get<std::vector<std::vector<double>>>("GH_PARS_flat_lateral");
110  fborder_corr_angulo_flat_lateral = fVUVHitsParams.get<std::vector<double>>("GH_border_angulo_flat_lateral");
111  fborder_corr_flat_lateral = fVUVHitsParams.get<std::vector<std::vector<double>>>("GH_border_flat_lateral");
112 
113  }
114  if (fIsDomePDCorr) {
115  fGHvuvpars_dome = fVUVHitsParams.get<std::vector<std::vector<double>>>("GH_PARS_dome");
116  fborder_corr_angulo_dome = fVUVHitsParams.get<std::vector<double>>("GH_border_angulo_dome");
117  fborder_corr_dome = fVUVHitsParams.get<std::vector<std::vector<double>>>("GH_border_dome");
118  }
119 
120  // Load corrections for VIS semi-analytic hits
121  if (fDoReflectedLight) {
122  mf::LogInfo("SemiAnalyticalModel") << "Using VIS (reflected) visibility parameterization";
123  fdelta_angulo_vis = fVISHitsParams.get<double>("delta_angulo_vis");
124 
125  if (fIsFlatPDCorr) {
126  fvis_distances_x_flat = fVISHitsParams.get<std::vector<double>>("VIS_distances_x_flat");
127  fvis_distances_r_flat = fVISHitsParams.get<std::vector<double>>("VIS_distances_r_flat");
128  fvispars_flat = fVISHitsParams.get<std::vector<std::vector<std::vector<double>>>>("VIS_correction_flat");
129  }
130  if (fIsDomePDCorr) {
131  fvis_distances_x_dome = fVISHitsParams.get<std::vector<double>>("VIS_distances_x_dome");
132  fvis_distances_r_dome = fVISHitsParams.get<std::vector<double>>("VIS_distances_r_dome");
133  fvispars_dome = fVISHitsParams.get<std::vector<std::vector<std::vector<double>>>>("VIS_correction_dome");
134  }
135 
136  // set cathode plane struct for solid angle function
137  fcathode_plane.h = fActiveVolumes[0].SizeY();
138  fcathode_plane.w = fActiveVolumes[0].SizeZ();
139  fplane_depth = std::abs(fcathode_centre[0]);
140  }
141 
142  // Load corrections for Anode reflections configuration
143  if (fIncludeAnodeReflections) {
144  mf::LogInfo("SemiAnalyticalModel") << "Using anode reflections parameterization";
145  fdelta_angulo_vis = fVISHitsParams.get<double>("delta_angulo_vis");
146  fAnodeReflectivity = fVISHitsParams.get<double>("AnodeReflectivity");
147 
148  if (fIsFlatPDCorr) {
149  fvis_distances_x_flat = fVISHitsParams.get<std::vector<double>>("VIS_distances_x_flat");
150  fvis_distances_r_flat = fVISHitsParams.get<std::vector<double>>("VIS_distances_r_flat");
151  fvispars_flat = fVISHitsParams.get<std::vector<std::vector<std::vector<double>>>>("VIS_correction_flat");
152  }
153 
154  if (fIsFlatPDCorrLat) {
155  fvis_distances_x_flat_lateral = fVISHitsParams.get<std::vector<double>>("VIS_distances_x_flat_lateral");
156  fvis_distances_r_flat_lateral = fVISHitsParams.get<std::vector<double>>("VIS_distances_r_flat_lateral");
157  fvispars_flat_lateral = fVISHitsParams.get<std::vector<std::vector<std::vector<double>>>>("VIS_correction_flat_lateral");
158  }
159 
160  // set anode plane struct for solid angle function
161  fanode_plane.h = fActiveVolumes[0].SizeY();
162  fanode_plane.w = fActiveVolumes[0].SizeZ();
163  fanode_plane_depth = fanode_centre[0];
164  }
165 
166 }
SemiAnalyticalModel::fvis_distances_r_dome
std::vector< double > fvis_distances_r_dome
Definition: SemiAnalyticalModel.h:215
SemiAnalyticalModel::fvis_distances_r_flat_lateral
std::vector< double > fvis_distances_r_flat_lateral
Definition: SemiAnalyticalModel.h:211
geo::GeometryCore::OpDetGeoFromOpDet
OpDetGeo const & OpDetGeoFromOpDet(unsigned int OpDet) const
Returns the geo::OpDetGeo object for the given detector number.
Definition: GeometryCore.cxx:1919
SemiAnalyticalModel::fborder_corr_dome
std::vector< std::vector< double > > fborder_corr_dome
Definition: SemiAnalyticalModel.h:193
SemiAnalyticalModel::fvispars_flat_lateral
std::vector< std::vector< std::vector< double > > > fvispars_flat_lateral
Definition: SemiAnalyticalModel.h:212
SemiAnalyticalModel::fcathode_centre
const TVector3 fcathode_centre
Definition: SemiAnalyticalModel.h:160
geo::OpDetGeo::isBar
bool isBar() const
Returns whether the detector shape is a bar (TGeoBBox).
Definition: OpDetGeo.h:195
SemiAnalyticalModel::fGHvuvpars_flat
std::vector< std::vector< double > > fGHvuvpars_flat
Definition: SemiAnalyticalModel.h:180
SemiAnalyticalModel::fGHvuvpars_dome
std::vector< std::vector< double > > fGHvuvpars_dome
Definition: SemiAnalyticalModel.h:191
SemiAnalyticalModel::fGeom
geo::GeometryCore const & fGeom
Definition: SemiAnalyticalModel.h:157
SemiAnalyticalModel::fvis_distances_r_flat
std::vector< double > fvis_distances_r_flat
Definition: SemiAnalyticalModel.h:207
SemiAnalyticalModel::fOpDetHeight
std::vector< double > fOpDetHeight
Definition: SemiAnalyticalModel.h:172
geo::OpDetGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Definition: OpDetGeo.cxx:40
SemiAnalyticalModel::fanode_centre
const TVector3 fanode_centre
Definition: SemiAnalyticalModel.h:160
SemiAnalyticalModel::fborder_corr_flat_lateral
std::vector< std::vector< double > > fborder_corr_flat_lateral
Definition: SemiAnalyticalModel.h:187
SemiAnalyticalModel::fVUVHitsParams
const fhicl::ParameterSet fVUVHitsParams
Definition: SemiAnalyticalModel.h:152
SemiAnalyticalModel::fvispars_flat
std::vector< std::vector< std::vector< double > > > fvispars_flat
Definition: SemiAnalyticalModel.h:208
geo::OpDetGeo::isSphere
bool isSphere() const
Returns whether the detector shape is a hemisphere (TGeoSphere).
Definition: OpDetGeo.h:198
SemiAnalyticalModel::fVISHitsParams
const fhicl::ParameterSet fVISHitsParams
Definition: SemiAnalyticalModel.h:153
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
util::counter
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
SemiAnalyticalModel::fvis_distances_x_flat
std::vector< double > fvis_distances_x_flat
Definition: SemiAnalyticalModel.h:206
geo::OpDetGeo::Length
double Length() const
Definition: OpDetGeo.h:81
SemiAnalyticalModel::fvis_distances_x_flat_lateral
std::vector< double > fvis_distances_x_flat_lateral
Definition: SemiAnalyticalModel.h:210
SemiAnalyticalModel::fOpDetCenter
std::vector< geo::Point_t > fOpDetCenter
Definition: SemiAnalyticalModel.h:168
SemiAnalyticalModel::fOpDetLength
std::vector< double > fOpDetLength
Definition: SemiAnalyticalModel.h:171
SemiAnalyticalModel::fborder_corr_angulo_flat_lateral
std::vector< double > fborder_corr_angulo_flat_lateral
Definition: SemiAnalyticalModel.h:186
SemiAnalyticalModel::fborder_corr_flat
std::vector< std::vector< double > > fborder_corr_flat
Definition: SemiAnalyticalModel.h:182
SemiAnalyticalModel::fGHvuvpars_flat_lateral
std::vector< std::vector< double > > fGHvuvpars_flat_lateral
Definition: SemiAnalyticalModel.h:185
SemiAnalyticalModel::fOpDetType
std::vector< int > fOpDetType
Definition: SemiAnalyticalModel.h:169
SemiAnalyticalModel::fborder_corr_angulo_dome
std::vector< double > fborder_corr_angulo_dome
Definition: SemiAnalyticalModel.h:192
SemiAnalyticalModel::fvis_distances_x_dome
std::vector< double > fvis_distances_x_dome
Definition: SemiAnalyticalModel.h:214
SemiAnalyticalModel::fActiveVolumes
const std::vector< geo::BoxBoundedGeo > fActiveVolumes
Definition: SemiAnalyticalModel.h:159
geo::OpDetGeo::Width
double Width() const
Definition: OpDetGeo.h:82
geo::OpDetGeo::Height
double Height() const
Definition: OpDetGeo.h:83
SemiAnalyticalModel::fborder_corr_angulo_flat
std::vector< double > fborder_corr_angulo_flat
Definition: SemiAnalyticalModel.h:181
SemiAnalyticalModel::fOpDetOrientation
std::vector< int > fOpDetOrientation
Definition: SemiAnalyticalModel.h:170
SemiAnalyticalModel::fvispars_dome
std::vector< std::vector< std::vector< double > > > fvispars_dome
Definition: SemiAnalyticalModel.h:216
SemiAnalyticalModel::fIncludeAnodeReflections
const bool fIncludeAnodeReflections
Definition: SemiAnalyticalModel.h:201
double SemiAnalyticalModel::interpolate ( const std::vector< double > &  xData,
const std::vector< double > &  yData,
double  x,
bool  extrapolate,
size_t  i = 0 
) const
private

Definition at line 726 of file SemiAnalyticalModel.cxx.

731 {
732  if (i == 0) {
733  size_t size = xData.size();
734  if (x >= xData[size - 2]) { // special case: beyond right end
735  i = size - 2;
736  }
737  else {
738  while (x > xData[i + 1])
739  i++;
740  }
741  }
742  double xL = xData[i];
743  double xR = xData[i + 1];
744  double yL = yData[i];
745  double yR = yData[i + 1]; // points on either side (unless beyond ends)
746  if (!extrapolate) { // if beyond ends of array and not extrapolating
747  if (x < xL) return yL;
748  if (x > xR) return yL;
749  }
750  const double dydx = (yR - yL) / (xR - xL); // gradient
751  return yL + dydx * (x - xL); // linear interpolation
752 }
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
double SemiAnalyticalModel::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 
) const
private

Definition at line 755 of file SemiAnalyticalModel.cxx.

761 {
762  // interpolate in x for each r bin, for angle bin k
763  const size_t nbins_r = parameters[k].size();
764  std::vector<double> interp_vals(nbins_r, 0.0);
765  {
766  size_t idx = 0;
767  size_t size = xDistances.size();
768  if (x >= xDistances[size - 2])
769  idx = size - 2;
770  else {
771  while (x > xDistances[idx + 1])
772  idx++;
773  }
774  for (size_t i = 0; i < nbins_r; ++i) {
775  interp_vals[i] = interpolate(xDistances,
776  parameters[k][i],
777  x,
778  false,
779  idx);
780  }
781  }
782  // interpolate in r
783  double border_correction = interpolate(rDistances, interp_vals, r, false);
784  return border_correction;
785 }
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
SemiAnalyticalModel::interpolate
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, double x, bool extrapolate, size_t i=0) const
Definition: SemiAnalyticalModel.cxx:726
k
pdgs k
Definition: selectors.fcl:22
r
esac echo uname r
Definition: condor_lar_pubs.sh:698
bool SemiAnalyticalModel::isOpDetInSameTPC ( geo::Point_t const &  ScintPoint,
geo::Point_t const &  OpDetPoint 
) const
private

Definition at line 688 of file SemiAnalyticalModel.cxx.

690 {
691  // method working for SBND, uBooNE, DUNE-HD 1x2x6 and DUNE-VD 1x8x6
692  // will need to be replaced to work in full DUNE geometry, ICARUS geometry
693  // check x coordinate has same sign or is close to zero, otherwise return false
694  if (((ScintPoint.X() < 0.) != (OpDetPoint.X() < 0.)) &&
695  std::abs(OpDetPoint.X()) > 10. && fNTPC == 2) { // TODO: replace with geometry service
696  return false;
697  }
698  return true;
699 }
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
double SemiAnalyticalModel::Omega_Dome_Model ( const double  distance,
const double  theta 
) const
private

Definition at line 653 of file SemiAnalyticalModel.cxx.

654 {
655  // this function calculates the solid angle of a semi-sphere of radius b,
656  // as a correction to the analytic formula of the on-axix solid angle,
657  // as we move off-axis an angle theta. We have used 9-angular bins
658  // with delta_theta width.
659 
660  // par0 = Radius correction close
661  // par1 = Radius correction far
662  // par2 = breaking distance betwween "close" and "far"
663 
664  double par0[9] = {0., 0., 0., 0., 0., 0.597542, 1.00872, 1.46993, 2.04221};
665  double par1[9] = {0, 0, 0.19569, 0.300449, 0.555598, 0.854939, 1.39166, 2.19141, 2.57732};
666  const double delta_theta = 10.; // TODO: should this be fdelta_angulo_vuv?
667  int j = int(theta/delta_theta);
668  // PMT radius
669  const double b = fradius; // cm
670  // distance form which the model parameters break (empirical value)
671  const double d_break = 5*b; //par2
672 
673  if(distance >= d_break) {
674  double R_apparent_far = b - par1[j];
675  double ratio_square = (R_apparent_far*R_apparent_far)/(distance*distance);
676  return (2*CLHEP::pi * (1 - std::sqrt(1 - ratio_square)));
677  }
678  else {
679  double R_apparent_close = b - par0[j];
680  double ratio_square = (R_apparent_close*R_apparent_close)/(distance*distance);
681  return (2*CLHEP::pi * (1 - std::sqrt(1 - ratio_square)));
682  }
683 }
icarus::distance
double distance(geo::Point_t const &point, CathodeDesc_t const &cathode)
Returns the distance of a point from the cathode.
Definition: CathodeCrossingUtils.cxx:35
pi
pdgs pi
Definition: selectors.fcl:22
double SemiAnalyticalModel::Rectangle_SolidAngle ( const double  a,
const double  b,
const double  d 
) const
private

Definition at line 575 of file SemiAnalyticalModel.cxx.

576 {
577  double aa = a / (2. * d);
578  double bb = b / (2. * d);
579  double aux = (1. + aa * aa + bb * bb) / ((1. + aa * aa) * (1. + bb * bb));
580  return 4. * fast_acos(std::sqrt(aux));
581 }
SemiAnalyticalModel::fast_acos
double fast_acos(double x) const
Definition: SemiAnalyticalModel.cxx:703
a
process_name gaushit a
Definition: decode_signalprocess_icarus.fcl:37
double SemiAnalyticalModel::Rectangle_SolidAngle ( Dims const &  o,
geo::Vector_t const &  v,
const double  OpDetOrientation 
) const
private

Definition at line 584 of file SemiAnalyticalModel.cxx.

585 {
586  // v is the position of the track segment with respect to
587  // the center position of the arapuca window
588 
589  // solid angle calculation depends on orientation of PD, set correct distances to use
590  double d1;
591  double d2;
592  if (OpDetOrientation == 1) {
593  // lateral PD, arapuca plane fixed in y direction
594  d1 = std::abs(v.X());
595  d2 = std::abs(v.Y());
596  }
597  else {
598  // anode/cathode PD, arapuca plane fixed in x direction [default]
599  d1 = std::abs(v.Y());
600  d2 = std::abs(v.X());
601  }
602  // arapuca plane fixed in x direction
603  if (isApproximatelyZero(d1) && isApproximatelyZero(v.Z())) {
604  return Rectangle_SolidAngle(o.h, o.w, d2);
605  }
606  if (isDefinitelyGreaterThan(d1, o.h * .5) && isDefinitelyGreaterThan(std::abs(v.Z()), o.w * .5)) {
607  double A = d1 - o.h * .5;
608  double B = std::abs(v.Z()) - o.w * .5;
609  double to_return = (Rectangle_SolidAngle(2. * (A + o.h), 2. * (B + o.w), d2) -
610  Rectangle_SolidAngle(2. * A, 2. * (B + o.w), d2) -
611  Rectangle_SolidAngle(2. * (A + o.h), 2. * B, d2) +
612  Rectangle_SolidAngle(2. * A, 2. * B, d2)) *
613  .25;
614  return to_return;
615  }
616  if ((d1 <= o.h * .5) && (std::abs(v.Z()) <= o.w * .5)) {
617  double A = -d1 + o.h * .5;
618  double B = -std::abs(v.Z()) + o.w * .5;
619  double to_return = (Rectangle_SolidAngle(2. * (o.h - A), 2. * (o.w - B), d2) +
620  Rectangle_SolidAngle(2. * A, 2. * (o.w - B), d2) +
621  Rectangle_SolidAngle(2. * (o.h - A), 2. * B, d2) +
622  Rectangle_SolidAngle(2. * A, 2. * B, d2)) *
623  .25;
624  return to_return;
625  }
626  if (isDefinitelyGreaterThan(d1, o.h * .5) && (std::abs(v.Z()) <= o.w * .5)) {
627  double A = d1 - o.h * .5;
628  double B = -std::abs(v.Z()) + o.w * .5;
629  double to_return = (Rectangle_SolidAngle(2. * (A + o.h), 2. * (o.w - B), d2) -
630  Rectangle_SolidAngle(2. * A, 2. * (o.w - B), d2) +
631  Rectangle_SolidAngle(2. * (A + o.h), 2. * B, d2) -
632  Rectangle_SolidAngle(2. * A, 2. * B, d2)) *
633  .25;
634  return to_return;
635  }
636  if ((d1 <= o.h * .5) && isDefinitelyGreaterThan(std::abs(v.Z()), o.w * .5)) {
637  double A = -d1 + o.h * .5;
638  double B = std::abs(v.Z()) - o.w * .5;
639  double to_return = (Rectangle_SolidAngle(2. * (o.h - A), 2. * (B + o.w), d2) -
640  Rectangle_SolidAngle(2. * (o.h - A), 2. * B, d2) +
641  Rectangle_SolidAngle(2. * A, 2. * (B + o.w), d2) -
642  Rectangle_SolidAngle(2. * A, 2. * B, d2)) *
643  .25;
644  return to_return;
645  }
646 
647  return 0.;
648 }
SemiAnalyticalModel::Rectangle_SolidAngle
double Rectangle_SolidAngle(const double a, const double b, const double d) const
Definition: SemiAnalyticalModel.cxx:575
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
lib.dedx.A
float A
Definition: dedx.py:137
double 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 399 of file SemiAnalyticalModel.cxx.

402 {
403  // set correct plane_depth
404  double plane_depth;
405  if (AnodeMode) plane_depth = fanode_plane_depth;
406  else plane_depth = ScintPoint.X() < 0. ? -fplane_depth : fplane_depth;
407 
408  // calculate visibility of the optical
409  // detector from the hotspot using solid angle:
410 
411  geo::Vector_t const emission_relative = hotspot - opDet.OpDetPoint;
412 
413  // calculate distances and angles for application of corrections
414  // distance from hotspot to optical detector
415  const double distance_vis = emission_relative.R();
416  // angle between hotspot and optical detector
417  double cosine_vis;
418  if (opDet.orientation == 1) { // lateral
419  cosine_vis = std::abs(emission_relative.Y()) / distance_vis;
420  }
421  else { // anode/cathode (default)
422  cosine_vis = std::abs(emission_relative.X()) / distance_vis;
423  }
424  const double theta_vis = fast_acos(cosine_vis) * 180. / CLHEP::pi;
425 
426  // calculate solid angle of optical channel
427  double solid_angle_detector = 0.;
428  // ARAPUCAS/Bars (rectangle)
429  if (opDet.type == 0) {
430  // get hotspot coordinates relative to opDet
431  geo::Vector_t const abs_emission_relative{std::abs(emission_relative.X()),
432  std::abs(emission_relative.Y()),
433  std::abs(emission_relative.Z())};
434  solid_angle_detector = Rectangle_SolidAngle(Dims{opDet.h, opDet.w}, abs_emission_relative, opDet.orientation);
435  }
436  // PMTS (dome)
437  else if (opDet.type == 1) {
438  solid_angle_detector = Omega_Dome_Model(distance_vis, theta_vis);
439  }
440  // PMTs (disk)
441  else if (opDet.type == 2) {
442  const double zy_offset = std::sqrt(emission_relative.Y() * emission_relative.Y() +
443  emission_relative.Z() * emission_relative.Z());
444  const double x_distance = std::abs(emission_relative.X());
445  solid_angle_detector = Disk_SolidAngle(zy_offset, x_distance, fradius);
446  }
447  else {
448  throw cet::exception("SemiAnalyticalModel")
449  << "Error: Invalid optical detector shape requested - configuration error in semi-analytical model, only rectangular, dome or disk optical detectors are supported." << "\n";
450  }
451 
452  // calculate number of hits via geometeric acceptance
453  double visibility_geo = (solid_angle_detector / (2. * CLHEP::pi)) *
454  cathode_visibility; // 2*pi due to presence of reflective foils
455 
456  // determine correction factor, depending on PD type
457  const size_t k = (theta_vis / fdelta_angulo_vis); // off-set angle bin
458  double r = std::hypot(ScintPoint.Y() - fcathode_centre[1], ScintPoint.Z() - fcathode_centre[2]);
459  double d_c = std::abs(ScintPoint.X() - plane_depth); // distance to cathode
460  double border_correction = 0;
461  // flat PDs
462  if ((opDet.type == 0 || opDet.type == 2) && (fIsFlatPDCorr || fIsFlatPDCorrLat)) {
463  // cathode/anode case
464  if (opDet.orientation == 0 && fIsFlatPDCorr) border_correction = interpolate2(fvis_distances_x_flat, fvis_distances_r_flat, fvispars_flat, d_c, r, k);
465  // laterals case
466  else if (opDet.orientation == 1 && fIsFlatPDCorrLat) border_correction = interpolate2(fvis_distances_x_flat_lateral, fvis_distances_r_flat_lateral, fvispars_flat_lateral, d_c, r, k);
467  else {
468  throw cet::exception("SemiAnalyticalModel")
469  << "Error: Invalid optical detector shape requested or corrections are missing - configuration error in semi-analytical model." << "\n";
470  }
471  }
472  // dome PDs
473  else if (opDet.type == 1 && fIsDomePDCorr) border_correction = interpolate2(fvis_distances_x_dome, fvis_distances_r_dome, fvispars_dome, d_c, r, k);
474  else {
475  throw cet::exception("SemiAnalyticalModel")
476  << "Error: Invalid optical detector shape requested or corrections are missing - configuration error in semi-analytical model." << "\n";
477  }
478 
479  // apply anode reflectivity factor
480  if (AnodeMode) border_correction = border_correction * fAnodeReflectivity;
481 
482  // apply field cage transparency factor
483  if (fApplyFieldCageTransparency) {
484  if (opDet.orientation == 1) border_correction = border_correction * fFieldCageTransparencyLateral;
485  else if (opDet.orientation == 0) border_correction = border_correction * fFieldCageTransparencyCathode;
486  }
487 
488  return border_correction * visibility_geo / cosine_vis;
489 }
SemiAnalyticalModel::fvis_distances_r_dome
std::vector< double > fvis_distances_r_dome
Definition: SemiAnalyticalModel.h:215
geo::Vector_t
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:164
SemiAnalyticalModel::fvis_distances_r_flat_lateral
std::vector< double > fvis_distances_r_flat_lateral
Definition: SemiAnalyticalModel.h:211
SemiAnalyticalModel::interpolate2
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) const
Definition: SemiAnalyticalModel.cxx:755
SemiAnalyticalModel::fvispars_flat_lateral
std::vector< std::vector< std::vector< double > > > fvispars_flat_lateral
Definition: SemiAnalyticalModel.h:212
SemiAnalyticalModel::Omega_Dome_Model
double Omega_Dome_Model(const double distance, const double theta) const
Definition: SemiAnalyticalModel.cxx:653
SemiAnalyticalModel::fcathode_centre
const TVector3 fcathode_centre
Definition: SemiAnalyticalModel.h:160
SemiAnalyticalModel::fvis_distances_r_flat
std::vector< double > fvis_distances_r_flat
Definition: SemiAnalyticalModel.h:207
SemiAnalyticalModel::fvispars_flat
std::vector< std::vector< std::vector< double > > > fvispars_flat
Definition: SemiAnalyticalModel.h:208
SemiAnalyticalModel::fast_acos
double fast_acos(double x) const
Definition: SemiAnalyticalModel.cxx:703
SemiAnalyticalModel::Rectangle_SolidAngle
double Rectangle_SolidAngle(const double a, const double b, const double d) const
Definition: SemiAnalyticalModel.cxx:575
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
SemiAnalyticalModel::fvis_distances_x_flat
std::vector< double > fvis_distances_x_flat
Definition: SemiAnalyticalModel.h:206
pi
pdgs pi
Definition: selectors.fcl:22
SemiAnalyticalModel::fvis_distances_x_flat_lateral
std::vector< double > fvis_distances_x_flat_lateral
Definition: SemiAnalyticalModel.h:210
SemiAnalyticalModel::fvis_distances_x_dome
std::vector< double > fvis_distances_x_dome
Definition: SemiAnalyticalModel.h:214
k
pdgs k
Definition: selectors.fcl:22
SemiAnalyticalModel::Disk_SolidAngle
double Disk_SolidAngle(const double d, const double h, const double b) const
Definition: SemiAnalyticalModel.cxx:508
r
esac echo uname r
Definition: condor_lar_pubs.sh:698
SemiAnalyticalModel::fvispars_dome
std::vector< std::vector< std::vector< double > > > fvispars_dome
Definition: SemiAnalyticalModel.h:216
int SemiAnalyticalModel::VUVAbsorptionLength ( ) const
private

Definition at line 169 of file SemiAnalyticalModel.cxx.

170 {
171  // determine LAr absorption length in cm
172  std::map<double, double> abs_length_spectrum = lar::providerFrom<detinfo::LArPropertiesService>()->AbsLengthSpectrum();
173  std::vector<double> x_v, y_v;
174  for (auto elem : abs_length_spectrum) {
175  x_v.push_back(elem.first);
176  y_v.push_back(elem.second);
177  }
178  int vuv_absorption_length = std::round(interpolate(x_v, y_v, 9.7, false)); // 9.7 eV: peak of VUV emission spectrum // TO DO UNHARDCODE FOR XENON
179  if (vuv_absorption_length <= 0){
180  throw cet::exception("SemiAnalyticalModel")
181  << "Error: VUV Absorption Length is 0 or negative.\n";
182  }
183  return vuv_absorption_length;
184 }
AbsLengthSpectrum
BEGIN_PROLOG AbsLengthSpectrum
Definition: larproperties.fcl:21
SemiAnalyticalModel::interpolate
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, double x, bool extrapolate, size_t i=0) const
Definition: SemiAnalyticalModel.cxx:726
double SemiAnalyticalModel::VUVVisibility ( geo::Point_t const &  ScintPoint,
OpticalDetector const &  opDet 
) const
private

Definition at line 209 of file SemiAnalyticalModel.cxx.

210 {
211  // distance and angle between ScintPoint and OpDetPoint
212  geo::Vector_t const relative = ScintPoint - opDet.OpDetPoint;
213  const double distance = relative.R();
214  double cosine;
215  if (opDet.orientation == 1) cosine = std::abs(relative.Y()) / distance;
216  else cosine = std::abs(relative.X()) / distance;
217  const double theta = fast_acos(cosine) * 180. / CLHEP::pi;
218 
219  double solid_angle = 0.;
220  // ARAPUCAS/Bars (rectangle)
221  if (opDet.type == 0) {
222  // get scintillation point coordinates relative to arapuca window centre
223  geo::Vector_t const abs_relative{std::abs(relative.X()), std::abs(relative.Y()), std::abs(relative.Z())};
224  solid_angle = Rectangle_SolidAngle(Dims{opDet.h, opDet.w}, abs_relative, opDet.orientation);
225  }
226  // PMTs (dome)
227  else if (opDet.type == 1) {
228  solid_angle = Omega_Dome_Model(distance, theta);
229  }
230  // PMTs (disk)
231  else if (opDet.type == 2) {
232  const double zy_offset = std::sqrt(relative.Y() * relative.Y() + relative.Z() * relative.Z());
233  const double x_distance = std::abs(relative.X());
234  solid_angle = Disk_SolidAngle(zy_offset, x_distance, fradius);
235  }
236  else {
237  throw cet::exception("SemiAnalyticalModel")
238  << "Error: Invalid optical detector shape requested - configuration error in semi-analytical model, only rectangular, dome or disk optical detectors are supported." << "\n";
239  }
240 
241  // calculate visibility by geometric acceptance
242  // accounting for solid angle and LAr absorbtion length
243  double visibility_geo = std::exp(-1. * distance / fvuv_absorption_length) * (solid_angle / (4 * CLHEP::pi));
244 
245  // apply Gaisser-Hillas correction for Rayleigh scattering distance
246  // and angular dependence offset angle bin
247  const size_t j = (theta / fdelta_angulo_vuv);
248 
249  // determine GH parameters, accounting for border effects
250  // radial distance from centre of detector (Y-Z)
251  double r = std::hypot(ScintPoint.Y() - fcathode_centre[1], ScintPoint.Z() - fcathode_centre[2]);
252 
253  double pars_ini[4] = {0, 0, 0, 0};
254  double s1 = 0; double s2 = 0; double s3 = 0;
255  // flat PDs
256  if ((opDet.type == 0 || opDet.type == 2) && (fIsFlatPDCorr || fIsFlatPDCorrLat)){
257  if (opDet.orientation == 1 && 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 = interpolate( fborder_corr_angulo_flat_lateral, fborder_corr_flat_lateral[0], theta, true);
263  s2 = interpolate( fborder_corr_angulo_flat_lateral, fborder_corr_flat_lateral[1], theta, true);
264  s3 = interpolate( fborder_corr_angulo_flat_lateral, fborder_corr_flat_lateral[2], theta, true);
265 
266  }
267  else if (opDet.orientation == 0 && fIsFlatPDCorr) { // cathode/anode
268  pars_ini[0] = fGHvuvpars_flat[0][j];
269  pars_ini[1] = fGHvuvpars_flat[1][j];
270  pars_ini[2] = fGHvuvpars_flat[2][j];
271  pars_ini[3] = fGHvuvpars_flat[3][j];
272  s1 = interpolate( fborder_corr_angulo_flat, fborder_corr_flat[0], theta, true);
273  s2 = interpolate( fborder_corr_angulo_flat, fborder_corr_flat[1], theta, true);
274  s3 = interpolate( fborder_corr_angulo_flat, fborder_corr_flat[2], theta, true);
275  }
276  else {
277  throw cet::exception("SemiAnalyticalModel")
278  << "Error: flat optical detectors are found, but parameters are missing - configuration error in semi-analytical model." << "\n";
279  }
280  }
281  // dome PDs
282  else if (opDet.type == 1 && fIsDomePDCorr) {
283  pars_ini[0] = fGHvuvpars_dome[0][j];
284  pars_ini[1] = fGHvuvpars_dome[1][j];
285  pars_ini[2] = fGHvuvpars_dome[2][j];
286  pars_ini[3] = fGHvuvpars_dome[3][j];
287  s1 = interpolate( fborder_corr_angulo_dome, fborder_corr_dome[0], theta, true);
288  s2 = interpolate( fborder_corr_angulo_dome, fborder_corr_dome[1], theta, true);
289  s3 = interpolate( fborder_corr_angulo_dome, fborder_corr_dome[2], theta, true);
290  }
291  else {
292  throw cet::exception("SemiAnalyticalModel")
293  << "Error: Invalid optical detector shape requested or corrections are missing - configuration error in semi-analytical model." << "\n";
294  }
295 
296  // add border correction to parameters
297  pars_ini[0] = pars_ini[0] + s1 * r;
298  pars_ini[1] = pars_ini[1] + s2 * r;
299  pars_ini[2] = pars_ini[2] + s3 * r;
300  pars_ini[3] = pars_ini[3];
301 
302  // calculate correction
303  double GH_correction = Gaisser_Hillas(distance, pars_ini);
304 
305  // apply field cage transparency factor
306  if (fApplyFieldCageTransparency) {
307  if (opDet.orientation == 1) GH_correction = GH_correction * fFieldCageTransparencyLateral;
308  else if (opDet.orientation == 0) GH_correction = GH_correction * fFieldCageTransparencyCathode;
309  }
310 
311  // determine corrected visibility of photo-detector
312  return GH_correction * visibility_geo / cosine;
313 }
geo::Vector_t
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:164
SemiAnalyticalModel::Gaisser_Hillas
double Gaisser_Hillas(const double x, const double *par) const
Definition: SemiAnalyticalModel.cxx:494
SemiAnalyticalModel::fborder_corr_dome
std::vector< std::vector< double > > fborder_corr_dome
Definition: SemiAnalyticalModel.h:193
SemiAnalyticalModel::Omega_Dome_Model
double Omega_Dome_Model(const double distance, const double theta) const
Definition: SemiAnalyticalModel.cxx:653
SemiAnalyticalModel::fcathode_centre
const TVector3 fcathode_centre
Definition: SemiAnalyticalModel.h:160
SemiAnalyticalModel::fGHvuvpars_flat
std::vector< std::vector< double > > fGHvuvpars_flat
Definition: SemiAnalyticalModel.h:180
SemiAnalyticalModel::fGHvuvpars_dome
std::vector< std::vector< double > > fGHvuvpars_dome
Definition: SemiAnalyticalModel.h:191
SemiAnalyticalModel::fborder_corr_flat_lateral
std::vector< std::vector< double > > fborder_corr_flat_lateral
Definition: SemiAnalyticalModel.h:187
SemiAnalyticalModel::fast_acos
double fast_acos(double x) const
Definition: SemiAnalyticalModel.cxx:703
SemiAnalyticalModel::Rectangle_SolidAngle
double Rectangle_SolidAngle(const double a, const double b, const double d) const
Definition: SemiAnalyticalModel.cxx:575
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
icarus::distance
double distance(geo::Point_t const &point, CathodeDesc_t const &cathode)
Returns the distance of a point from the cathode.
Definition: CathodeCrossingUtils.cxx:35
pi
pdgs pi
Definition: selectors.fcl:22
SemiAnalyticalModel::interpolate
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, double x, bool extrapolate, size_t i=0) const
Definition: SemiAnalyticalModel.cxx:726
SemiAnalyticalModel::fborder_corr_angulo_flat_lateral
std::vector< double > fborder_corr_angulo_flat_lateral
Definition: SemiAnalyticalModel.h:186
SemiAnalyticalModel::fborder_corr_flat
std::vector< std::vector< double > > fborder_corr_flat
Definition: SemiAnalyticalModel.h:182
SemiAnalyticalModel::fGHvuvpars_flat_lateral
std::vector< std::vector< double > > fGHvuvpars_flat_lateral
Definition: SemiAnalyticalModel.h:185
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
SemiAnalyticalModel::fborder_corr_angulo_dome
std::vector< double > fborder_corr_angulo_dome
Definition: SemiAnalyticalModel.h:192
SemiAnalyticalModel::Disk_SolidAngle
double Disk_SolidAngle(const double d, const double h, const double b) const
Definition: SemiAnalyticalModel.cxx:508
r
esac echo uname r
Definition: condor_lar_pubs.sh:698
SemiAnalyticalModel::fborder_corr_angulo_flat
std::vector< double > fborder_corr_angulo_flat
Definition: SemiAnalyticalModel.h:181

Member Data Documentation

template<typename TReal >
constexpr SemiAnalyticalModel::bool
inlinestaticprivate
Initial value:
{
TReal diff = std::fabs(a - b)

Definition at line 111 of file SemiAnalyticalModel.h.

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

Definition at line 159 of file SemiAnalyticalModel.h.

return SemiAnalyticalModel::false
private

Definition at line 115 of file SemiAnalyticalModel.h.

const TVector3 SemiAnalyticalModel::fanode_centre
private

Definition at line 160 of file SemiAnalyticalModel.h.

Dims SemiAnalyticalModel::fanode_plane
private

Definition at line 167 of file SemiAnalyticalModel.h.

double SemiAnalyticalModel::fanode_plane_depth
private

Definition at line 161 of file SemiAnalyticalModel.h.

double SemiAnalyticalModel::fAnodeReflectivity
private

Definition at line 204 of file SemiAnalyticalModel.h.

bool SemiAnalyticalModel::fApplyFieldCageTransparency
private

Definition at line 195 of file SemiAnalyticalModel.h.

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

Definition at line 192 of file SemiAnalyticalModel.h.

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

Definition at line 181 of file SemiAnalyticalModel.h.

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

Definition at line 186 of file SemiAnalyticalModel.h.

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

Definition at line 193 of file SemiAnalyticalModel.h.

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

Definition at line 182 of file SemiAnalyticalModel.h.

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

Definition at line 187 of file SemiAnalyticalModel.h.

const TVector3 SemiAnalyticalModel::fcathode_centre
private

Definition at line 160 of file SemiAnalyticalModel.h.

Dims SemiAnalyticalModel::fcathode_plane
private

Definition at line 166 of file SemiAnalyticalModel.h.

double SemiAnalyticalModel::fdelta_angulo_vis
private

Definition at line 203 of file SemiAnalyticalModel.h.

double SemiAnalyticalModel::fdelta_angulo_vuv
private

Definition at line 177 of file SemiAnalyticalModel.h.

const bool SemiAnalyticalModel::fDoReflectedLight
private

Definition at line 200 of file SemiAnalyticalModel.h.

double SemiAnalyticalModel::fFieldCageTransparencyCathode
private

Definition at line 197 of file SemiAnalyticalModel.h.

double SemiAnalyticalModel::fFieldCageTransparencyLateral
private

Definition at line 196 of file SemiAnalyticalModel.h.

geo::GeometryCore const& SemiAnalyticalModel::fGeom
private

Definition at line 157 of file SemiAnalyticalModel.h.

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

Definition at line 191 of file SemiAnalyticalModel.h.

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

Definition at line 180 of file SemiAnalyticalModel.h.

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

Definition at line 185 of file SemiAnalyticalModel.h.

const bool SemiAnalyticalModel::fIncludeAnodeReflections
private

Definition at line 201 of file SemiAnalyticalModel.h.

bool SemiAnalyticalModel::fIsDomePDCorr
private

Definition at line 190 of file SemiAnalyticalModel.h.

bool SemiAnalyticalModel::fIsFlatPDCorr
private

Definition at line 179 of file SemiAnalyticalModel.h.

bool SemiAnalyticalModel::fIsFlatPDCorrLat
private

Definition at line 184 of file SemiAnalyticalModel.h.

const larg4::ISTPC SemiAnalyticalModel::fISTPC
private

Definition at line 156 of file SemiAnalyticalModel.h.

const int SemiAnalyticalModel::fNTPC
private

Definition at line 158 of file SemiAnalyticalModel.h.

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

Definition at line 168 of file SemiAnalyticalModel.h.

std::vector<double> SemiAnalyticalModel::fOpDetHeight
private

Definition at line 172 of file SemiAnalyticalModel.h.

std::vector<double> SemiAnalyticalModel::fOpDetLength
private

Definition at line 171 of file SemiAnalyticalModel.h.

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

Definition at line 170 of file SemiAnalyticalModel.h.

std::vector<int> SemiAnalyticalModel::fOpDetType
private

Definition at line 169 of file SemiAnalyticalModel.h.

double SemiAnalyticalModel::fplane_depth
private

Definition at line 161 of file SemiAnalyticalModel.h.

double SemiAnalyticalModel::fradius
private

Definition at line 165 of file SemiAnalyticalModel.h.

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

Definition at line 215 of file SemiAnalyticalModel.h.

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

Definition at line 207 of file SemiAnalyticalModel.h.

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

Definition at line 211 of file SemiAnalyticalModel.h.

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

Definition at line 214 of file SemiAnalyticalModel.h.

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

Definition at line 206 of file SemiAnalyticalModel.h.

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

Definition at line 210 of file SemiAnalyticalModel.h.

const fhicl::ParameterSet SemiAnalyticalModel::fVISHitsParams
private

Definition at line 153 of file SemiAnalyticalModel.h.

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

Definition at line 216 of file SemiAnalyticalModel.h.

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

Definition at line 208 of file SemiAnalyticalModel.h.

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

Definition at line 212 of file SemiAnalyticalModel.h.

const int SemiAnalyticalModel::fvuv_absorption_length
private

Definition at line 174 of file SemiAnalyticalModel.h.

const fhicl::ParameterSet SemiAnalyticalModel::fVUVHitsParams
private

Definition at line 152 of file SemiAnalyticalModel.h.

const size_t SemiAnalyticalModel::nOpDets
private

Definition at line 164 of file SemiAnalyticalModel.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.5