SDetector/Scintillator.cc

Go to the documentation of this file.

#include <iostream>
#include <sstream>
#include <vector>
#include <boost/lexical_cast.hpp>

#include <sdet/Station.h>
#include <sdet/Scintillator.h>
#include <det/Detector.h>
#include <sdet/SDetector.h>
#include <utl/ErrorLogger.h>
#include <utl/Math.h>
#include <utl/Point.h>
#include <utl/Line.h>
#include <utl/GeometryUtilities.h>
#include <utl/UTMPoint.h>
#include <utl/Reader.h>
#include <utl/CoordinateSystem.h>
#include <utl/AugerCoordinateSystem.h>
#include <fwk/LocalCoordinateSystem.h>

using namespace sdet;
using namespace std;
using namespace utl;


double
Scintillator::GetBarLength()
  const
{
  if (!fBarLength.IsValid()) {
    GetScintillatorData(fBarLength.Get(), "barLength", "scintillator", "scintillator bar length");
    fBarLength.SetValid();
  }
  return fBarLength.Get();
}


double
Scintillator::GetBarWidth()
  const
{
  if (!fBarWidth.IsValid()) {
    GetScintillatorData(fBarWidth.Get(), "barWidth", "scintillator", "scintillator bar width");
    fBarWidth.SetValid();
  }
  return fBarWidth.Get();
}


double
Scintillator::GetBarThickness()
  const
{
  if (!fBarThickness.IsValid()) {
    GetScintillatorData(fBarThickness.Get(), "barThickness", "scintillator", "scintillator bar thickness");
    fBarThickness.SetValid();
  }
  return fBarThickness.Get();
}


unsigned int
Scintillator::GetNBars()
  const
{
  if (!fNBars.IsValid()) {
    const det::ManagerRegister& m = det::Detector::GetInstance().GetSManagerRegister();
    m.GetDataOrThrow(fNBars.Get(), "scintillator", "barCount", indexMap);
    fNBars.SetValid();
  }
  return fNBars.Get();
}


double
Scintillator::GetGap()
  const
{
  if (!fGap.IsValid()) {
    GetScintillatorData(fGap.Get(), "gap", "scintillator", "gap between sides of scintillator");
    fGap.SetValid();
  }
  return fGap.Get();
}


double
Scintillator::GetFiberLength()
  const
{
  if (!fFiberLength.IsValid()) {
    GetScintillatorData(fFiberLength.Get(), "fiberLength", "scintillator",
      "length of on fiber routed through two scintillator bars");
    fFiberLength.SetValid();
  }
  return fFiberLength.Get();
}


double
Scintillator::GetFiberUTurnDiameter()
  const
{
  if (!fFiberUTurnDiameter.IsValid()) {
    GetScintillatorData(fFiberUTurnDiameter.Get(), "fiberUTurnDiameter", "scintillator",
      "diameter of the u-turn fibers make in their connections from one bar to the next");
    fFiberUTurnDiameter.SetValid();
  }
  return fFiberUTurnDiameter.Get();
}


double
Scintillator::GetHousingLength()
  const
{
  if (!fHousingLength.IsValid()) {
    GetScintillatorData(fHousingLength.Get(), "housingLength", "scintillator", "length of scintillator polystyrene housing");
    fHousingLength.SetValid();
  }
  return fHousingLength.Get();
}


double
Scintillator::GetHousingWidth()
  const
{
  if (!fHousingWidth.IsValid()) {
    GetScintillatorData(fHousingWidth.Get(), "housingWidth", "scintillator", "width of scintillator polystyrene housing");
    fHousingWidth.SetValid();
  }
  return fHousingWidth.Get();
}


double
Scintillator::GetHousingThickness()
  const
{
  if (!fHousingThickness.IsValid()) {
    GetScintillatorData(fHousingThickness.Get(), "housingThickness", "scintillator", "thickness of scintillator polystyrene housing");
    fHousingThickness.SetValid();
  }
  return fHousingThickness.Get();
}


double
Scintillator::GetCasingThickness()
  const
{
  if (!fCasingThickness.IsValid()) {
    GetScintillatorData(fCasingThickness.Get(), "casingThickness", "scintillator", "thickness of scintillator aluminum casing");
    fCasingThickness.SetValid();
  }
  return fCasingThickness.Get();
}


double
Scintillator::GetSandwichPanelThickness()
  const
{
  if (!fSandwichPanelThickness.IsValid()) {
    GetScintillatorData(fSandwichPanelThickness.Get(), "sandwichPanelThickness", "scintillator", "thickness of the foam portion of the sandwich board of the scintillator");
    fSandwichPanelThickness.SetValid();
  }
  return fSandwichPanelThickness.Get();
}


double
Scintillator::GetRoofOffset()
  const
{
  if (!fRoofOffset.IsValid()) {
    GetScintillatorData(fRoofOffset.Get(), "roofOffset", "scintillator", "offset of the bottom of the roof from the top of the top aluminum sheet of the scintillator");
    fRoofOffset.SetValid();
  }
  return fRoofOffset.Get();
}


double
Scintillator::GetRoofThickness()
  const
{
  if (!fRoofThickness.IsValid()) {
    GetScintillatorData(fRoofThickness.Get(), "roofThickness", "scintillator", "thickness of the scintillator roof");
    fRoofThickness.SetValid();
  }
  return fRoofThickness.Get();
}


utl::Point
Scintillator::GetPosition()
  const
{
  if (!fPosition.IsValid()) {

    vector<double> coords;
    GetScintillatorData(coords, "position", "scintillator", "scintillator position");

    // Return the scintillator position in a coordinate system whose origin is
    // at the intersection of the ground and the central axis of the tank

    CoordinateSystemPtr tCS =
      det::Detector::GetInstance().GetSDetector().GetStation(fStationId).GetLocalCoordinateSystem();

    fPosition = utl::Point(coords[0], coords[1], coords[2], tCS);

  }
  return fPosition.Get();
}


double
Scintillator::GetFiberExpDecayConstant()
  const
{
  if (!fFiberExpDecayConstant.IsValid()) {
    GetScintillatorData(fFiberExpDecayConstant.Get(), "fiberExpDecayConst", "scintillator", "fiber exponential decay constant");
    fFiberExpDecayConstant.SetValid();
  }
  return fFiberExpDecayConstant.Get();
}


double
Scintillator::GetBarExpDecayConstant()
  const
{
  if (!fBarExpDecayConstant.IsValid()) {
    GetScintillatorData(fBarExpDecayConstant.Get(), "barExpDecayConst", "scintillator", "scintillator bar exponential decay constant");
    fBarExpDecayConstant.SetValid();
  }
  return fBarExpDecayConstant.Get();
}


double
Scintillator::GetAttenuationLength()
  const
{
  if (!fAttenuationLength.IsValid()) {
    GetScintillatorData(fAttenuationLength.Get(), "attenuationLength", "scintillator", "attenuation length");
    fAttenuationLength.SetValid();
  }
  return fAttenuationLength.Get();
}


double
Scintillator::GetBarLeakageExpCoefficient()
  const
{
  if (!fBarLeakageExpCoefficient.IsValid()) {
    GetScintillatorData(fBarLeakageExpCoefficient.Get(), "barLeakageExpCoefficient",
      "scintillator", "exponential coefficient describing photon leakage at ends of bar");
    fBarLeakageExpCoefficient.SetValid();
  }
  return fBarLeakageExpCoefficient.Get();
}


double
Scintillator::GetBarLeakageMaxRatio()
  const
{
  if (!fBarLeakageMaxRatio.IsValid()) {
    GetScintillatorData(fBarLeakageMaxRatio.Get(), "barLeakageMaxRatio", "scintillator",
      "maximum ratio of photon leakage at ends of bar");
    fBarLeakageMaxRatio.SetValid();
  }
  return fBarLeakageMaxRatio.Get();
}


double
Scintillator::GetEffectiveRefractiveIndex()
  const
{
  if (!fEffectiveRefractiveIndex.IsValid()) {
    GetScintillatorData(fEffectiveRefractiveIndex.Get(), "effectiveRefractiveIndex", "scintillator", "effective refractive index");
    fEffectiveRefractiveIndex.SetValid();
  }
  return fEffectiveRefractiveIndex.Get();
}


double
Scintillator::GetReferencePENumber()
  const
{
  if (!fReferencePENumber.IsValid()) {
    GetScintillatorData(fReferencePENumber.Get(), "referencePENumber", "scintillator", "non attenuated average number of PEs per VMIP");
    fReferencePENumber.SetValid();
  }
  return fReferencePENumber.Get();
}


double
Scintillator::GetSigmaBarToBar()
  const
{
  if (!fSigmaBarToBar.IsValid()) {
    GetScintillatorData(fSigmaBarToBar.Get(), "sigmaBarToBar", "scintillator", "standard deviation of bar-to-bar response");
    fSigmaBarToBar.SetValid();
  }
  return fSigmaBarToBar.Get();
}


double
Scintillator::GetReferenceEnergy()
  const
{
  if (!fReferenceEnergy.IsValid()) {
    GetScintillatorData(fReferenceEnergy.Get(), "referenceEnergy", "scintillator", "average energy deposit per VMIP");
    fReferenceEnergy.SetValid();
  }
  return fReferenceEnergy.Get();
}


void
Scintillator::NotFoundAndThrow(const string& msg)
  const
{
  ostringstream err;
  err << "Did not find requested component: '" << msg
      << "' for station = " << indexMap["stationId"];
  throw NonExistentComponentException(err.str());
}


Scintillator::Scintillator(const int stationId, const unsigned int isUUB) :
  fStationId(stationId)
{
  /*
   * Adding stationId and isUUB to index map, which is passed
   * to the managers with each GetScintillatorData request (or when
   * Handle is used).
   */
  indexMap["stationId"] = boost::lexical_cast<string>(stationId);
  indexMap["isUUB"] = boost::lexical_cast<string>(isUUB);
}


std::vector<double>
Scintillator::GetDistancesToPMT(const Point& p)
  const
{
  CoordinateSystemPtr cs =
    det::Detector::GetInstance().GetSDetector().GetStation(fStationId).GetLocalCoordinateSystem();
  /*
   * Check to ensure point is inside one of the scintillator wings. If the point is within the a
   * "halo" distance (set to the width of one bar), this check will still pass (as the intention
   * of the check is just to ensure requests are not being made with absurd coordinates and has
   * no effect on the distance calculation). A more elegant management of the halo could be
   * implemented in the future.
   */
  if (!IsInsideBar(p, GetBarWidth())) {
    ostringstream err;
    err << "Point is not within the active area (bars) of the scintillator.";
    throw OutOfBoundException(err.str());
  }
  const double longLeg =
    GetBarLength() + 0.5*GetGap() -
    abs(p.GetX(cs) - GetPosition().GetX(cs)) + 0.25*kPi*GetFiberUTurnDiameter() + 0.5*GetFiberLength();
  std::vector<double> v { GetFiberLength() - longLeg, longLeg };
  return v;
}


pair<double, double>
Scintillator::GetDistancesToBarEnds(const Point& p)
  const
{
  CoordinateSystemPtr cs =
    det::Detector::GetInstance().GetSDetector().GetStation(fStationId).GetLocalCoordinateSystem();
  /*
   * Check to ensure point is inside one of the scintillator wings. If the point is within the a
   * "halo" distance (set to the width of one bar), this check will still pass (as the intention
   * of the check is just to ensure requests are not being made with absurd coordinates and has
   * no effect on the distance calculation). A more elegant management of the halo could be
   * implemented in the future.
   */
  if (!IsInsideBar(p, GetBarWidth())) {
    ostringstream err;
    err << "Point is not within the active area (bars) of the scintillator.";
    throw OutOfBoundException(err.str());
  }
  // distance to end of bar closer to PMT
  const double di = abs(p.GetX(cs) - GetPosition().GetX(cs)) - 0.5*GetGap();
  return make_pair(di, GetBarLength()-di);
}


bool
Scintillator::IsHit(const Point& particlePoint, const Vector& particleDirection, const double halo)
  const
{ 
  const Line particleLine(particlePoint, particleDirection);
  CoordinateSystemPtr cs =
    det::Detector::GetInstance().GetSDetector().GetStation(fStationId).GetLocalCoordinateSystem();
  const utl::Point& position = GetPosition();
  const Plane scintillatorPlane(Point(0, 0, position.GetZ(cs), cs),
                       Vector(0, 0, 1, cs));
  const Point where = Intersection(scintillatorPlane, particleLine);

  // variables needed to define the area where particle might hit the detector
  const double gapLength = GetGap();
  const double barLength = GetBarLength();
  const int barCount = GetNBars();
  const double barWidth = GetBarWidth();
  const double width = 0.5*barCount * barWidth;

  const double scintillatorPosition_X = position.GetX(cs);
  const double scintillatorPosition_Y = position.GetY(cs); 

  // Limits for the interaction area
  const double limitX1_r = scintillatorPosition_X + 0.5*gapLength - halo;
  const double limitX2_r = scintillatorPosition_X + 0.5*gapLength + barLength + halo;
  const double limitX1_l = scintillatorPosition_X - 0.5*gapLength + halo;
  const double limitX2_l = scintillatorPosition_X - 0.5*gapLength - barLength - halo;

  const double limitY1 = scintillatorPosition_Y - 0.5*width - halo;
  const double limitY2 = scintillatorPosition_Y + 0.5*width + halo;

  const double x = where.GetX(cs);
  const double y = where.GetY(cs);
  if (limitX1_r < x && x < limitX2_r && limitY1 < y && y < limitY2) // Try right panel
    return true;
  if (limitX2_l < x && x < limitX1_l && limitY1 < y && y < limitY2) // Try left panel
    return true;

  return false;
}


bool
Scintillator::IsInsideBar(const Point& p, const double halo)
  const
{
  CoordinateSystemPtr cs =
    det::Detector::GetInstance().GetSDetector().GetStation(fStationId).GetLocalCoordinateSystem();
  // check if in xy plane of either wing
  if (!IsHit(p, Vector(0, 0, 1, cs), halo))
    return false;
  // check if z-coordinate is within bar height bounds
  if (abs(p.GetZ(cs) - GetPosition().GetZ(cs)) > 0.5*GetBarThickness() + halo)
    return false;

  return true;
}


double
Scintillator::GetMaxRadius()
  const
{
  const utl::Point& position = GetPosition();
  const double x0 = position.GetX(position.GetCoordinateSystem());
  const double y0 = position.GetY(position.GetCoordinateSystem());

  const double maxX = abs(x0) + 0.5*GetHousingLength() + GetCasingThickness();
  const double maxY = abs(y0) + 0.5*GetHousingWidth() + GetCasingThickness();

  return hypot(maxX, maxY);
}


double
Scintillator::GetMaxHeight()
  const
{
  const utl::Point& position = GetPosition();
  const double z = position.GetZ(position.GetCoordinateSystem());
  return z - 0.5*GetBarThickness() - GetCasingThickness() - GetSandwichPanelThickness() + GetHousingThickness() + GetRoofOffset() + GetRoofThickness();
}


void
Scintillator::Update()
  const
{
  // In this method we clean all the variables that might change with time.
  fPosition.SetValid(false);
}