MDetector/Scintillator.cc

Go to the documentation of this file.

#include <mdet/Scintillator.h>

#include <mdet/Module.h>
#include <mdet/Counter.h>
#include <mdet/MHierarchyInfo.h>
//
#include <sdet/Station.h>
//
#include <utl/Point.h>
#include <utl/Line.h>
#include <utl/Vector.h>
#include <utl/ErrorLogger.h>
#include <utl/RandomEngine.h>
#include <utl/TransformationMatrix.h>
#include <utl/GeometryUtilities.h>
//
#include <fwk/RandomEngineRegistry.h>
#include <fwk/LocalCoordinateSystem.h>
//
#include <CLHEP/Random/RandExponential.h>
//
#include <cmath>
//
#include <boost/function.hpp>
#include <boost/lambda/lambda.hpp>
#include <boost/logic/tribool.hpp>
#include <boost/ptr_container/ptr_vector.hpp>

#ifdef DEBUG
#  include <iostream>
#endif


namespace mdet {

  const char* const Scintillator::kComponentName = MHierarchyInfo::kComponentsNames[3];

  const char* const Scintillator::kComponentId = MHierarchyInfo::kComponentsIds[3];


  Scintillator::Scintillator(const int sId, const det::VManager::IndexMap& parentMap, const Module& parent) :
    MDetectorComponent<Scintillator>::Type(sId, parentMap),
    det::MPositionable<Scintillator>(*this),
    fModule(parent)
  {
    Register(fLength);
    Register(fWidth);
    Register(fHeight);
    Register(fLocalSoilDensity);
    Register(fDecayTime);
  }


  utl::CoordinateSystemPtr
  Scintillator::GetReferenceCoordinateSystem()
    const
  {
    return GetModule().GetLocalCoordinateSystem();
  }


  const Module&
  Scintillator::GetModule()
    const
  {
    return fModule;
  }


  double
  Scintillator::GetLength()
    const
  {
    return GetData(fLength, "length");
  }


  double
  Scintillator::GetWidth()
    const
  {
    return GetData(fWidth, "width");
  }


  double
  Scintillator::GetHeight()
    const
  {
    return GetData(fHeight, "height");
  }


  double
  Scintillator::GetLocalSoilDensity()
    const
  {
    return GetData(fLocalSoilDensity, "localSoilDensity");
  }


  double
  Scintillator::GetDecayTime()
    const
  {
    return GetData(fDecayTime, "decayTime");
  }

  double
  Scintillator::GetEpsilon()
    const
  {
    return GetData(fEpsilon, "epsilon");
  }


  double
  Scintillator::ComputeDecayDelay()
    const
  {
          // This method gives sense to the bare number decay time: it's an exponential process.
          utl::RandomEngine& rnd = fwk::RandomEngineRegistry::GetInstance().Get(fwk::RandomEngineRegistry::eDetector);
          double decayTime = CLHEP::RandExponential::shoot(& rnd.GetEngine(), GetDecayTime());

          return decayTime;
  }


  double
  Scintillator::GetDecayDelayMean()
    const
  {
    // The mean is the parameter.
    return GetDecayTime();
  }


  double
  Scintillator::GetDecayDelayStdDev()
    const
  {
    // For an exponential process the std-dev equals the mean value.
    return GetDecayDelayMean();
  }


  bool
  Scintillator::Contains(const utl::Point& p)
    const
  {
    /*
     * XXX
     * In all these methods the Fiber position may have to be
     * taken into account for exactness:
     * that is, a part of the configured volume for the scintillator
     * is, in fact, fiber, not scint.
     */
    /*
     * As of the general recommendation, this logic
     * avoid the usage of Get<Component> methods
     * from utl::Point.
     * Instead of that, it uses geometrical operations
     * such as inner product between vectors.
     */
    /*
     * Difference against the origin of the Scintillator, to
     * refer against that origin.
     */
    const utl::Vector pLocal = p - GetPosition();
    // Project each component & see if it lies within the boundary.
    return InExtent(GetWidth(),  pLocal, utl::Vector(1, 0, 0, GetLocalCoordinateSystem())) &&
           InExtent(GetLength(), pLocal, utl::Vector(0, 1, 0, GetLocalCoordinateSystem())) &&
           InExtent(GetHeight(), pLocal, utl::Vector(0, 0, 1, GetLocalCoordinateSystem()));
  }


  bool
  Scintillator::IntersectedBy(const utl::Line& l)
    const
  {
    // About top & bot:
    // At first the line should touch these both planes in one point for each
    // , at least somewhere far from the Scint.
    // The line won't touch in a single point them if is parallel to then, which
    // won't be a normal situation (that's an horizontal particle!?!?!)
    // Given the current dimensions, if the line intersects the Scintillator,
    // then it's more likely that crosses top & bot: because of that they
    // were put first in this if-chain.
    //--
    // Theoretically if the line intersects the box in one point, then
    // that should be enough to, then, intersect it in another one.
    // Practically, it was found that, due to rounding errors, when a point
    // should intersect inside, the result lays outside. So, let's check
    // for two-fold intersection.
    //--
    // See ComputeIntersectionWith for more details.
    // Note that we may short-cut before the two final planes of intersection
    // are used: I mean the planes of the two final points used while
    // computing the intersection.
    const int limit = 2;
    // For the first two, no need to use "if".
    int count = OnSide(ComputeTopPlane(), l) ;
    count    += OnSide(ComputeBottomPlane(), l);
    if (count < limit)
      count += OnSide(ComputeLeftPlane(), l);
    if (count < limit)
      count += OnSide(ComputeRightPlane(), l);
    if (count < limit)
      count += OnSide(ComputeFrontPlane(), l);
    if (count < limit)
      count += OnSide(ComputeBackPlane(), l);
    return count >= limit;
  }


  bool
  Scintillator::OnSide(const utl::Plane& sidePlane, const utl::Line& l)
    const
  {
    /*
     * Despite being const, certain operations
     * transform one geometry object to the 
     * coordinate system of the other: so
     * the scint may
     * be intertectedBy but then the computation
     * of the intersection dies.
     * It had to be put deep in the function hierarchy,
     * because repeated calls to this function from 
     * another function in the class already triggered
     * the problem. So, to avoid that, copies were made of 
     * the line. But at last, the problem was not solved by
     * these copies, but by an increase in the epsilon.
     */ 
    // First the line must cross the given plane, and then the crossing
    // point has to belong to the Scintillator.
    bool b1 = Intersects(sidePlane, l);
    // Use b1 so as to avoid further calculations (&&-shorcut).
    bool b2 = b1 && Contains(utl::Intersection(sidePlane, l));
    return b2; // b2 is enough given the && above.
  }


  utl::Segment
  Scintillator::ComputeIntersectionWith(const utl::Line& l)
    const
  {
    boost::ptr_vector<utl::Point> pCross;
    // Collect the intersecting points.
    IntersectionHelper(pCross, l, &Scintillator::ComputeTopPlane);
    IntersectionHelper(pCross, l, &Scintillator::ComputeBottomPlane);
    IntersectionHelper(pCross, l, &Scintillator::ComputeLeftPlane);
    IntersectionHelper(pCross, l, &Scintillator::ComputeRightPlane);
    IntersectionHelper(pCross, l, &Scintillator::ComputeFrontPlane);
    IntersectionHelper(pCross, l, &Scintillator::ComputeBackPlane);
    unsigned int index2ndPoint = 1; 
    // If there are only two, the 2nd point has to be the 2nd in the vector.
    if (pCross.size() > 2) {
      // In, no pun, corner cases there maybe more than two planes; for instance
      // when the particle enters through a corner. So, the criteria is to take
      // the first, and then the one further away, in a way to avoid taking
      // the two of the same corner.
      double dist = 0; // Start with 0: we're looking for a max of positives.
      for (unsigned int i = 1; i < pCross.size(); ++i) {
        double d = (pCross[i] - pCross[0]).GetMag();
        if (d > dist) {
          dist          = d;
          index2ndPoint = i;
        }
      }
    }
    /*
     *  The tail of the arrow goes closer to the ground & the end arrow closer to the scint:
     * for that we need the ground plane and a point within it.
     */
    const utl::Plane ground = ComputeGroundPlane();
    /*
     * Retrieve the two points (here it's supposed they exist, if not the
     * precondition of the method -over l & this - is violated by the call) and
     * then see the distance to the ground.
     */
    // Note that we should have two points as of the precondition that the
    // particle has to intersect the scint (the method IntersectedBy was
    // modified to look for two planes, not 1 any more).
    const double d0 = ( pCross[0] - ground.GetAnchor() ).GetMag();
    const double d1 = ( pCross[index2ndPoint] - ground.GetAnchor() ).GetMag();
    const bool head = d1 > d0; // true -> 1 for the head.
    /*
     * "If the source type is bool, the value false is converted to zero
     * and the value true is converted to one." (4.7 (4) from C++ Standard ANSI 14882).
     * Construct the pair with both points.
     *
     * The first is the head and the second is the tail.
     */
    return utl::Segment(pCross[ !head ], pCross[ head ]);
  }


  template<class Container>
  void
  Scintillator::IntersectionHelper(Container& v, const utl::Line& l,
                                   utl::Plane (Scintillator::* ComputePlane)() const)
    const
  {
    // Compute the plane calling the pointed method.
    const utl::Plane p = (this->*ComputePlane)();
    if (OnSide(p, l)) {
      /*
       * The container has to (and actually does, see
       * ComputeIntersectionWith) take care of deletion.
       */
      v.push_back(new utl::Point(utl::Intersection(p, l)));
    }
  }


  utl::Segment
  Scintillator::ComputeUndergroundIntersectionOf(const utl::Line& l)
    const
  {
    // With the ground...
    const utl::Plane ground  = ComputeGroundPlane();
    const utl::Point pGround = utl::Intersection(ground, l);
    // With the scint...
    const utl::Segment track = ComputeIntersectionWith(l);
    // Now the segment intersecting the ground / soil goes from the
    // point at the ground until it reaches the scintillator.
    return utl::Segment(pGround, track.GetBegin());
  }


  double
  Scintillator::ComputeUndergroundGrammageFor(const utl::Line& l)
    const
  {
    // Simple density parametrization.
    return GetLocalSoilDensity() * ComputeUndergroundIntersectionOf(l).GetLength();
  }


  utl::Plane
  Scintillator::ComputeGroundPlane()
    const
  {
    utl::CoordinateSystemPtr cs = GetModule().GetCounter().GetAssociatedTank().GetLocalCoordinateSystem();
    const utl::Point topAnchor(0, 0, 0, cs);
    const utl::Vector normalZ(0, 0, 1, cs);
    return utl::Plane(topAnchor, normalZ);
  }


//  double
//  Scintillator::ComputeEffectiveArea(const utl::Line& /*dir*/)
/*    const
  {
    // TODO Simple calculation by now: just the 2-D area.
    // No angle involved. No height.
    return GetLength() * GetWidth();
  }
*/


  double
  Scintillator::GetArea()
    const
  {
    return GetLength() * GetWidth();
  }


  /*
   * All these methods for computing different planes define the extent of the scintillator.
   */

  // This method also defines how the extent is taken.
  bool
  Scintillator::InExtent(double ext, const utl::Vector& pLocal, const utl::Vector& normal)
    const
  {
    const double c = normal * pLocal;
    /*
     * Comparing with these bounds (having the upper coming from one of width, length and height)
     * defines that the origin is in the border).
     * See also the different plane computing methods.
     *
     * In addition to the operator comparison, we resort to a tolerance one also, against the
     * borders.
     */
    const double ext2 = ext / 2;
    // With the abs of c we compare both sides.
    return (-ext2 <= c && c <= ext2) || IsCloseTo(std::fabs(c), ext2);
  }


  utl::Plane
  Scintillator::ComputeTopPlane()
    const
  {
    return ComputePlaneHelper(GetHeight() / 2, 0, 0, 1);
  }


  utl::Plane
  Scintillator::ComputeBottomPlane()
    const
  {
    return ComputePlaneHelper(- GetHeight() / 2, 0, 0, 1);
  }


  utl::Plane
  Scintillator::ComputeBackPlane()
    const
  {
    return ComputePlaneHelper(GetLength() / 2, 0, 1, 0);
  }


  utl::Plane
  Scintillator::ComputeFrontPlane()
    const
  {
    return ComputePlaneHelper(- GetLength() / 2, 0, 1, 0);
  }


  utl::Plane
  Scintillator::ComputeRightPlane()
    const
  {
    return ComputePlaneHelper(GetWidth() / 2, 1, 0, 0);
  }


  utl::Plane
  Scintillator::ComputeLeftPlane()
    const
  {
    return ComputePlaneHelper(- GetWidth() / 2, 1, 0, 0);
  }


  utl::Plane 
  Scintillator::ComputePlaneHelper(double ext, double x, double y, double z)
    const
  {
    const utl::Point a(x * ext, y * ext, z * ext, GetLocalCoordinateSystem());
    const utl::Vector n(x, y, z, GetLocalCoordinateSystem());
    return utl::Plane(a, n);
  }


  bool
  Scintillator::IsCloseTo(double a, double b)
    const
  {
    return std::fabs(a - b) <=  GetEpsilon() * std::max(std::fabs(a), std::fabs(b));
  }


  bool
  Scintillator::Intersects(const utl::Plane& p, const utl::Line& l)
    const
  {
    /*
     *  double to bool conversion, implicit comparission with zero.
     *
     *  This comparison is also used in the following code:
     *
     * ./Modules/HybridGeometryFinderOG/HybridGeometryFinder.cc:213:
     *    if(my_sdp.GetMag() == 0 ) {
     * ./Modules/FdSDPFinderOG/FdSDPFinder.cc:306:
     *    if ( crossprod.GetMag() != 0. )
     * ./Documentation/ExampleApplications/HReconstruction/UserModule.cc:245:
     *    if(sdp.GetMag() == 0 )
     * ./Documentation/ExampleApplications/HReconstruction/UserModule.cc:1231:
     *    if (core_station.GetMag() != 0.0 )
     *
     *  Explicit comparison isn't necessary: see Stroustrup's book (he didn't write
     *  a single one, but I bet you know what I'm talking about).
     *  C.6.2.5 Boolean Conversions (p. 835)
     *  Pointers, integral, and floating point values can be implicitly converted to bool (§4.2).
     *  A nonzero value converts to true; a zero value converts to false .
     *
     * (ok, let's say it: Stroustrup, Bjarne, "The C++ Programming Language", Special Edition.)
     *
     * 6/8/2010 - Change the comparisson == 0 (implicit or not), to a more proper floating point
     * comparison.
     *
     * If the direction is perpendicual to the normal (ie dot product == 0) then the line
     * is parallel to the plane.
     */
    double d = p.GetNormal() * l.GetDirection();
    return ! IsCloseTo(d, 0);
  }


  // At last, there was already a method for this in utl.
  /*
  utl::Point
  Scintillator::ComputeIntersection(const utl::Plane& p, const utl::Line& l)
  {
    // See:
    // http://softsurfer.com/Archive/algorithm_0104/algorithm_0104B.htm
    // http://en.wikipedia.org/wiki/Line-plane_intersection
    // Calculate the parameter that solves the intesection between plane & line.
    double alfa =  (p.GetNormal() * (p.GetAnchor() - l.GetAnchor())) /
    // quotient    ____________________________________________________
                        (p.GetNormal() * l.GetDirection());
    // Now use the parameter to compute the specific point & we're done!
    // XXX Will the "computational representation" of this point also
    // be in the Plane object (it should)?
    return l.Propagate(alfa);
  }
  */

}