MDetector/PMT.cc

Go to the documentation of this file.

#include <mdet/PMT.h>

#include <mdet/Pixel.h>
#include <mdet/MHierarchyInfo.h>
//
#include <utl/ErrorLogger.h>
#include <utl/RandomEngine.h>
//
#include <fwk/RandomEngineRegistry.h>
//
#include <CLHEP/Random/RandFlat.h>
//
#include <cmath>
#include <algorithm>
#include <numeric>
#include <sstream>
//
#include <boost/algorithm/minmax_element.hpp>
#include <boost/lambda/lambda.hpp>


namespace { // unnamed namespace with local helper stuff.

  struct PixelEq {
    PixelEq(const mdet::Pixel& p) : fPixel(p) { }

    bool
    operator()(const mdet::Pixel* const o)
      const
    {
       /*
        * Same object (equal addresses).
        * Here we're relying on the fact that there
        * isn't a subverted overload for operator& 
        * among Pixels and then that we're calling the
        * built-in one which gives the object-address.
        *
        * Boost's utility AddressOf could have been used
        * but it seems that using them would have been too
        * fancy more than a real concern: don't overload
        * operator&, please (at least for a mdet::Pixel). See 
        * http://www.artima.com/cppsource/YourAddress.html
        * http://www.artima.com/weblogs/viewpost.jsp?thread=102857
        */
      if (o == &fPixel)
        return true;
      // Id has to be equal
      return o->GetId() == fPixel.GetId() && 
        // and then also within the hierarchy.
        o->GetIdsMap() == fPixel.GetIdsMap();
    }
    const mdet::Pixel& fPixel;
  };


  template<class PGroup, class Container>
  struct PixelPusher {
    PixelPusher(PGroup& g, Container& c) : fGroup(g), fContainer(c) { }

    void
    operator()(int id)
    {
      const mdet::Pixel& p = fGroup.Get(id);
      fContainer.push_back(&p);
    }
    PGroup& fGroup;
    Container& fContainer;
  };

}


namespace mdet {

  const char* const PMT::kComponentName = MHierarchyInfo::kComponentsNames[5];

  const char* const PMT::kComponentId = MHierarchyInfo::kComponentsIds[5];

  /*
   * See C++ standard ISO/IEC 14882. 
   * 9.4.2(4) Static data members
   *
   * If a static data member is of const integral or const enumeration type, its declaration in the class
   * definition can specify a constant-initializer which shall be an integral constant expression (5.19). In that
   * case, the member can appear in integral constant expressions. The member shall still be defined in a namespace
   * scope if it is used in the program and the namespace scope definition shall not contain an initializer.
   *
   * So we have to initialize this constant here.
   */ 
  const double PMT::kCentralPixelSelfTalk = 100;


  PMT::PixelGroup::SizeType 
  PMT::GetRows() 
    const 
  {
    PixelGroup::SizeType tot  = fPixels.GetNumberOfComponents();
    // Work with square PMTs.
    // To the able to work with rectangular one another datum is needed in the config,
    // for instance the number of rows.
    PixelGroup::SizeType root = static_cast<PixelGroup::SizeType>(std::sqrt(double(tot)));
    if (tot != root*root) {
      std::ostringstream ss;
      ss << "The total number of pixels is not a perfect square number in "; 
      AddIdMessage(ss);
      ss << ".";
      ERROR(ss);
    }
    return root;
  }


  PMT::PixelGroup::SizeType 
  PMT::GetCols() 
    const 
  {
    PixelGroup::SizeType tot  = fPixels.GetNumberOfComponents();
    return tot / GetRows();
  }


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


  PMT::PMT(int pId, const det::VManager::IndexMap& parentMap, const Module& parent) :
    MDetectorComponent<PMT>::Type(pId, parentMap),
    fPixels(*this),
    fModule(parent)
  {
    fPixels.Update(GetIdsMap());
    InitPixels();
  }

  
  void
  PMT::InitPixels() 
  {
    if (fPixels.GetNumberOfComponents()) {
      // Take the minimun element to make everything relative to this id.
      /*
       * Something like this can't be done cause the iterators doesn't model
       * asignable concept.
       *
       *  std::min_element(PixelsBegin(), PixelsEnd(),
       *  boost::bind(&Pixel::GetId, _1) < boost::bind(&Pixel::GetId, _2))->GetId()
       *./det/ComponentGroup.h:108: error: non-static reference member 
       * 'const det::ComponentGroup<mdet::PMT, mdet::Pixel, det::ParentCreator, 
       * mdet::MManagerProvider>& 
       * det::ComponentGroup<mdet::PMT, mdet::Pixel, det::ParentCreator, mdet::MManagerProvider>
       *  ::InternalConstFunctor::fContainer', can't use default assignment operator
       */
      { // scope this!
        PixelConstIterator i = PixelsBegin();
        PixelConstIterator e = PixelsEnd();
        // It's know that there's one at least.
        fFirstIdPMT = i->GetId(); // Note that we're initializing the member variable. 
        for (++i; i != e; ++i) 
          fFirstIdPMT = std::min(fFirstIdPMT, i->GetId());
      }
      // The code inside asumes at least one: it initializes neighbors;
      // which don't exist if there's only one.
      if (fPixels.GetNumberOfComponents() > 1) {
        // At last fill the neighbor information.
        const PixelGroup::SizeType cols = GetCols();
        const PixelGroup::SizeType colsLessOne = GetCols() - 1;
        const PixelGroup::SizeType rowsLessOne = GetRows() - 1;
        // To accumulate the ids to check'em.
        std::vector<int> check;
        // Preallocation.
        check.reserve(fPixels.GetNumberOfComponents());
        // Now loop over the pixs and calculate the neighbors.
        for (PixelIterator i = fPixels.Begin(); i != fPixels.End(); ++i) {
          const int id = i->GetId();
          check.push_back(id);  
          /*
           * Relativize: note that it's positive since 
           * there's no id that's less that firstId.
           */ 
          const unsigned int idOffset = id - fFirstIdPMT;
          /*
           * Here we are assuming consecutiveness on the ids.
           */
          /*
           * From the ISO/IEC 14882:2003(E) C++ Standard about / and %:
           * 5.6(4) "If both operands are nonnegative then the remainder is 
           * nonnegative; if not, the sign of the remainder is implementation-defined"
           */ 
          const unsigned int quot = idOffset / cols;
          const unsigned int mod = idOffset % cols;
          const bool firstCol = !mod;
          const bool firstRow = !quot;
          const bool lastCol = (mod == colsLessOne);
          const bool lastRow = (quot == rowsLessOne);
          // Count the booleans. Note that the count is at most 2 (there's more 
          // than 1 pixel);  and rely on the (guaranteed) true -> 1 conversion.
          int count = firstCol + firstRow + lastCol + lastRow;
          // Preallocate the neighbors container.
          if (count == 2) 
            i->fNeighbors.reserve(kNumNeighborsCorner);
          else if (count == 1)
            i->fNeighbors.reserve(kNumNeighborsSide);
          else 
            i->fNeighbors.reserve(kNumNeighborsInside);
          // Build a helper functor to do the actual filling.
          PixelPusher<PixelGroup, NeighborsContainer> pusher(fPixels, i->fNeighbors);
          /*
           * Push the neighbors in increasing id order.
           * (cols is at least one, of course), and so
           * the order in which the following checks are
           * made isn't arbitrary.
           *
           * As indicated by the constant we work with 8 
           * neighbors.
           * 
           * This is related to the way each Pixel stores
           * it's crosstalk data.
           */ 
          if (!firstRow) {
            if (!firstCol)
              pusher(id - cols -1);  // Prev col, prev row.
            pusher(id - cols);       // Same col, prev row.
            if (!lastCol)
              pusher(id - cols + 1); // Next col, prev row.
          }
          if (!firstCol)
            pusher(id - 1);          // Prev col, same row.
          if (!lastCol) 
            pusher(id + 1);          // Next col, same row.
          if (!lastRow) {
            if (!firstCol) 
              pusher(id + cols - 1); // Prev col, next row.
            pusher(id + cols);       // Same col, next row.
            if (!lastCol) 
              pusher(id + cols + 1); // Next col, next row.
          }
        } // End loop over pixels.
        //
        // Check Consecutiveness.
        //
        // It's likely that it will be already sorted 
        // (DetectorComponent sorts while retreiving data) 
        // but ensure that it's needed.
        //
        // It's known that there's one element at least.
        const bool sorted =
          std::adjacent_find(check.begin(), check.end(), boost::lambda::_1 > boost::lambda::_2) == check.end();
        if (!sorted) 
          std::sort(check.begin(), check.end()); 
        // Take the difference between all elements & overwrite the vector.
        std::adjacent_difference(check.begin(), check.end(), check.begin());
        // Remember that the first element is kept the same, so ignore it: that's why the + 1.
        // Now take the min & max: should be one if there's consecutiveness.
        typedef std::vector<int>::const_iterator It;
        std::pair<It, It> mm = boost::minmax_element(check.begin() + 1, check.end());
        if (*mm.first != 1 || *mm.second != 1) {
          std::ostringstream ss;
          ss << "There is no consecutiveness in the pixels of PMT ";
          AddIdMessage(ss);
          ss << ".";
          FATAL(ss);
        }
      }
    }
  }


  NeighborConstIterator 
  PMT::NeighborsBegin(const Pixel& p) 
    const 
  {
    return boost::make_indirect_iterator(p.fNeighbors.begin());
  }

  
  NeighborConstIterator 
  PMT::NeighborsEnd(const Pixel& p) 
    const 
  {
    return boost::make_indirect_iterator(p.fNeighbors.end());
  }


  bool 
  PMT::IsNeighbor(const Pixel& src, const Pixel& dst) 
    const
  {
    PixelEq e(dst); 
    return IsNeighbor(src, dst, e) != src.fNeighbors.end();
  }


  template<class Predicate>
  NeighborsContainer::const_iterator 
  PMT::IsNeighbor(const Pixel& src, const Pixel& /*dst*/, const Predicate& p) 
    const
  {
    // Look for the destination in the neighbors of src.
    return std::find_if(src.fNeighbors.begin(), src.fNeighbors.end(), p);
  }


  const Pixel& 
  PMT::ComputePulseDestination(const Pixel& src) 
    const
  {
    utl::RandomEngine& eng = fwk::RandomEngineRegistry::GetInstance().Get(fwk::RandomEngineRegistry::eDetector);
    // We shoot an uniform random number between 0 and 1 (that's why we are working with xtalk proportion).
    const double rnd = CLHEP::RandFlat::shoot(& eng.GetEngine());
    /*
     * We iterate over the neighbors, and in this variable we represent the right side limit
     * of an internal: if the random number lays to the left: that's the pixel (note that
     * laying to the left means that lays to the right of the previous right limit).
     */  
    double rightSide = 0;
    for (NeighborConstIterator nIt = NeighborsBegin(src), end = NeighborsEnd(src); nIt != end; ++nIt) {
      // Increase the limit with the crosstalk proportion of the current neighbor:
      rightSide += GetCrossTalkProportion(src, *nIt);
      /*
       * If rand lays between "previous rightSize" (it can't happen to be to the left of the
       * "previous rightSize" because in that case we would have exited on the previous loop)
       * and rightSide, that's the detination neighbor.
       */
      if (rnd < rightSide)
        return *nIt;
    }
    // Stayed in the same pixel.
    return src;
  }


  double
  PMT::GetCrossTalkNormalizationFactor(const mdet::Pixel& pix) 
    const
  {
    // Lazy evaluation...
    if (!pix.fCrossTalkNormalizationFactor.IsValid()) {
      const CrossTalkContainer* xt = nullptr;
      switch (pix.fNeighbors.size()) {
      case kNumNeighborsInside:
        xt = &pix.GetNeighborsCrossTalkInside();
        break;
      case kNumNeighborsSide:
        xt = &pix.GetNeighborsCrossTalkSide();
        break;
      case kNumNeighborsCorner:
        xt = &pix.GetNeighborsCrossTalkCorner();
        break;
      }
      // Assume xt != 0 (see mdet::PMT::GetCrossTalk).
      pix.fCrossTalkNormalizationFactor = std::accumulate(xt->begin(), xt->end(), PMT::kCentralPixelSelfTalk);
    }
    return pix.fCrossTalkNormalizationFactor.Get();
  }


  double 
  PMT::GetCrossTalkProportion(const Pixel& src, const Pixel& dst) 
    const
  {
    return GetCrossTalk(src, dst) / GetCrossTalkNormalizationFactor(src);
  }


  double 
  PMT::GetCrossTalk(const Pixel& src, const Pixel& dst) 
    const 
  {
    /*
     * (Almost) Repeated code: see constructor.
     */
    const unsigned int idOffset = src.GetId() - fFirstIdPMT;
    const unsigned int quot = idOffset / GetCols();
    const unsigned int mod = idOffset % GetCols();
    // no need const bool firstCol = mod  == 0;
    const bool firstRow = !quot;
    const bool lastCol = (mod == GetCols() - 1);
    const bool lastRow = (quot == GetCols() - 1);
    //----
    //----
    PixelEq e(dst); // to be used twice.
    NeighborsContainer::const_iterator i = IsNeighbor(src, dst, e);
    if (i != src.fNeighbors.end()) {
      // take id offset wrt the first neighbor.
      NeighborsContainer::difference_type o = i - src.fNeighbors.begin();
      // Select the kind of cross-talk.
      // Use pointer instead of reference to avoid the need to initialize
      // it here. Note that we are taking the address of the return value,
      // which is returned by const reference (it's not a temporary).
      const CrossTalkContainer* xt = nullptr;
      // Switch over to chose the proper cross-talk.
      switch (src.fNeighbors.size()) {
      /*
       * See the declaration of these constants for 
       * the specification of how the coefficients
       * are supposed to be laid out in the CrossTalkContainer.
       *
       * Note that we use integral literals in place-of
       * NeighborsContainer::difference_type. We may
       * assume that there's an implicit conversion
       * available.
       */ 
      case kNumNeighborsInside:
        xt = & src.GetNeighborsCrossTalkInside();
        // the index remains the same.
        break;
      case kNumNeighborsSide:
        xt = & src.GetNeighborsCrossTalkSide();
        // XXX What an ugly piece of code!!!
        // Recall we're in a side.
        if (lastCol) {
          if (o == 0)
            o = 1;
          else if (o == 1)
            o = 0;
          // the rest stay the same.
        } else if (lastRow) {
          if (o != 4) {
            if (o == 3)
              o = 0;
            else 
              ++o;
          }
        } else if (firstRow) {
          if (o != 0) {
            if (o == 1)
              o = 4;
            else
              --o;
          }
        } else {
          // firstCol
          if (o == 3)
            o = 4;
          else if (o == 4)
            o = 3;
          // the rest stay the same.
        }
        // End of ugly code!
        break;
      case kNumNeighborsCorner:
        xt = & src.GetNeighborsCrossTalkCorner();
        // Recall we're in a corner.
        if (lastCol) {
          if (lastRow) {
            o = (o + 1) % kNumNeighborsCorner;
          }
          // if firstRow the indexes remain the same.
        } else {
          // firstCol
          if (firstRow ) {
            if (o != 0) {
              // if not 2, then is o == 1.
              o = o == 2 ? 1 : 2;
            }
          } else {
            // lastRow
            if (o != 1) {
              o = o == 0 ? 2 : 0;
            }
          }
        }
        break;
      default:
        std::ostringstream ss;
        ss << "The neighbors list for the Pixel hasn't got the "
              "expected size, it has " 
           << src.fNeighbors.size()
           << " in ";
        AddIdMessage(ss);
        ss << '.';
        FATAL(ss);
        return 0; // Fatal error happened: return 0.
      }
      if (xt->size() != src.fNeighbors.size()) {
        std::ostringstream ss;
        ss << "The neighbors coefficients list for the Pixel hasn't got the "
              "same length that the cross-talk one, it has " 
           << src.fNeighbors.size()
           << " and the second has "
           << xt->size()
           << " in ";
        AddIdMessage(ss);
        ss << '.';
        FATAL(ss);
        return 0; // idem.
      }
      return (*xt)[o];
    } else {
      // It can be the pixel itself (use the convention constant)...
      if (e(&src))
        return kCentralPixelSelfTalk;
      return 0; // ..or a non-first neighbor pixel: there's no cross-talk.
    }
  }


  void 
  PMT::Update(bool invalidateData, bool invalidateComponents)
  {
    MDetectorComponent<PMT>::Type::Update(invalidateData, invalidateComponents);
    fPixels.Update(GetIdsMap(), invalidateData, invalidateComponents);
    // After update pixels, we may re-init them; if no component invalidation then
    // we skip this part which initializes fFirstIdPMT (which won't change if no
    // component invalidation is performed).
    if (invalidateComponents)
      InitPixels();
  }

}