JMuonPrefit.hh

Go to the documentation of this file.

#ifndef __JRECONSTRUCTION__JMUONPREFIT__
#define __JRECONSTRUCTION__JMUONPREFIT__
 
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
 
#include "km3net-dataformat/online/JDAQEvent.hh"
#include "km3net-dataformat/online/JDAQTimeslice.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"
 
#include "JTrigger/JHitR1.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JAlgorithm.hh"
#include "JTrigger/JMatch1D.hh"
#include "JTrigger/JMatch3B.hh"
 
#include "JFit/JLine1Z.hh"
#include "JFit/JLine1ZEstimator.hh"
#include "JFit/JMatrixNZ.hh"
#include "JFit/JVectorNZ.hh"
#include "JFit/JFitToolkit.hh"
 
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JMuonPrefitParameters_t.hh"
 
#include "JMath/JConstants.hh"
#include "JTools/JPermutation.hh"
 
#include "JLang/JPredicate.hh"
 
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JDetectorToolkit.hh"
 
#include "JDynamics/JCoverage.hh"
 
#include "JGeometry3D/JOmega3D.hh"
#include "JGeometry3D/JRotation3D.hh"
 
#include "Jeep/JMessage.hh"
 
 
/**
 * \author mdejong, gmaggi, azegarelli
 */
 
namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }
 
namespace JRECONSTRUCTION {
 
  using KM3NETDAQ::JDAQEvent;
  using JDETECTOR::JModuleRouter;
  using JFIT::JLine1Z;
  using JFIT::JEstimator;
  using JFIT::JMatrixNZ;
  using JFIT::JVectorNZ;
  using JTRIGGER::JMatch3B;
  using JTRIGGER::JMatch1D;
  using JGEOMETRY3D::JOmega3D;
  using JDYNAMICS::coverage_type;
 
 
  /**
   * Wrapper class to make pre-fit of muon trajectory.
   *
   * The JMuonPrefit fit is used to generate start values for subsequent fits (usually JMuonSimplex and JMuonGandalf).\n
   * To this end, a scan of directions is made and the time and transverse positions of the track are fitted for each direction (JFIT::JEstimator<JLine1Z>).\n
   * The directions are spaced by the parameters JMuonPrefitParameters_t::gridAngle_deg.\n
   * This angle corresponds to the envisaged angular accuracy of the result.\n
   * The probability that one of the results is less than this angle away from the correct value,
   * multiple start values should be considered (JMuonPrefitParameters_t::numberOfPrefits).\n
   * Note that the CPU time scales with the inverse of the square of this angle.\n
   * The chi-squared is based on the time residuals.\n
   */
  struct JMuonPrefit :
    public JMuonPrefitParameters_t,
    public JEstimator<JLine1Z>
  {
    typedef JEstimator<JLine1Z>                  JEstimator_t;
    typedef JTRIGGER::JHitR1                     hit_type;
    typedef std::vector<hit_type>                buffer_type;
 
    using JEstimator_t::operator();
 
    /**
     * Input data type.
     */
    struct input_type :
      public JDAQEventHeader
    {
      /**
       * Default constructor.
       */
      input_type()
      {}
 
 
      /**
       * Constructor.
       *
       * \param  header         header
       * \param  coverage       coverage
       */
      input_type(const JDAQEventHeader& header,
                 const coverage_type&   coverage) :
        JDAQEventHeader(header),
        coverage       (coverage)
      {}
 
      buffer_type   dataL0;
      buffer_type   dataL1;
      coverage_type coverage;
    };
 
 
    /**
     * Constructor 
     * 
     * \param  parameters    parameters
     * \param  debug         debug
     */
    JMuonPrefit(const JMuonPrefitParameters_t& parameters,
                const int                      debug = 0) :
      JMuonPrefitParameters_t(parameters),
      omega (parameters.gridAngle_deg * JMATH::PI/180.0),
      debug (debug)
    {}
 
 
    /**
     * Constructor 
     * 
     * \param  parameters    parameters
     * \param  omega         directions
     * \param  debug         debug
     */
    JMuonPrefit(const JMuonPrefitParameters_t& parameters,
                const JOmega3D&                omega,
                const int                      debug = 0) :
      JMuonPrefitParameters_t(parameters),
      omega (omega),
      debug (debug)
    {}
 
 
    /**
     * Get input data.
     *
     * \param  router        module router
     * \param  event         event
     * \param  coverage      coverage
     * \return               input data
     */
    input_type getInput(const JModuleRouter& router,
                        const JDAQEvent&     event,
                        const coverage_type& coverage) const
 
    {
      using namespace std;
      using namespace JTRIGGER;
      using namespace KM3NETDAQ;
 
      const JBuildL0<hit_type> buildL0;
      const JBuildL2<hit_type> buildL2(JL2Parameters(2, TMaxLocal_ns, ctMin));
 
      input_type   input(event.getDAQEventHeader(), coverage);
    
      buffer_type& dataL0 = input.dataL0;
      buffer_type& dataL1 = input.dataL1;
 
      buildL2(event, router, !useL0, back_inserter(dataL1));
 
      // 3D cluster of unique optical modules
 
      const JMatch3B<hit_type> match3B(roadWidth_m, TMaxLocal_ns);
 
      sort(dataL1.begin(), dataL1.end(), hit_type::compare);
 
      buffer_type::iterator __end = dataL1.end();
 
      __end = unique(dataL1.begin(), __end, equal_to<JDAQModuleIdentifier>());
    
      __end = clusterizeWeight(dataL1.begin(), __end, match3B);
 
      dataL1.erase(__end, dataL1.end());
 
 
      if (useL0) {
 
        buildL0(event, router, true, back_inserter(dataL0));
 
        __end = dataL0.end();
 
        for (buffer_type::iterator i = dataL0.begin(); i != __end; ) {
 
          if (match3B.count(*i, dataL1.begin(), dataL1.end()) != 0)
            ++i;
          else
            swap(*i, *--__end);
        }
 
        dataL0.erase(__end, dataL0.end());
      }
 
      return input;
    }
 
 
    /**
     * Fit function.
     *
     * \param  input         input data
     * \return               fit results
     */
    JEvt operator()(const input_type& input)
    {
      using namespace std;
      using namespace JPP;
 
      const double HIT_OFF  =  1.0e3 * sigma_ns * sigma_ns;  // [ns^2]
 
      const JMatch1D<hit_type> match1D(roadWidth_m, TMaxLocal_ns);
 
      const buffer_type& dataL0 = input.dataL0;
      const buffer_type& dataL1 = input.dataL1;
 
      data.reserve(dataL0.size() + 
                   dataL1.size());
 
      JEvent event(JMUONPREFIT);
 
      JEvt out;
 
      for (JOmega3D_t::const_iterator dir = omega.begin(); dir != omega.end(); ++dir) {
 
        const JRotation3D R(*dir);
 
 
        buffer_type::iterator __end = copy(dataL1.begin(), dataL1.end(), data.begin());
 
        for (buffer_type::iterator i = data.begin(); i != __end; ++i) {
          i->rotate(R);
        }
 
 
        // reduce data
 
        if (distance(data.begin(), __end) > NMaxHits) {
 
          buffer_type::iterator __p = data.begin();
 
          advance(__p, NMaxHits);
 
          partial_sort(data.begin(), __p, __end, cmz);
 
          __end = __p;
        }
 
 
        // 1D cluster
 
        __end = clusterizeWeight(data.begin(), __end, match1D);
 
        if (useL0) {
 
          buffer_type::iterator p = __end;                                     // begin L0 data
          buffer_type::iterator q = copy(dataL0.begin(), dataL0.end(), p);     // end   L0 data
 
          for (buffer_type::iterator i = p; i != q; ++i) {
 
            if (find_if(data.begin(), __end, make_predicate(&hit_type::getModuleID, i->getModuleID())) == __end) {
 
              i->rotate(R);
 
              if (match1D.count(*i, data.begin(), __end) != 0) {
                *p = *i;
                ++p;
              }
            }
          }
 
          __end = clusterize(__end, p, match1D);
        }
 
          
        if (distance(data.begin(), __end) <= NUMBER_OF_PARAMETERS) {
          continue;
        }
 
 
        // 1D fit
 
        JLine1Z  tz;
        double   chi2 = numeric_limits<double>::max();
        int      NDF  = distance(data.begin(), __end) - NUMBER_OF_PARAMETERS;
        int      N    = getCount(data.begin(), __end);
 
        try {
 
          sort(data.begin(), __end, hit_type::compare);
 
          (*this)(data.begin(), __end);
 
          V.set(*this, data.begin(), __end, gridAngle_deg, sigma_ns);
          Y.set(*this, data.begin(), __end);
 
          V.invert();
 
          tz   = *this;
          chi2 = getChi2(Y, V);
        }
        catch(const exception& error) {
          continue;
        }
 
 
        // outlier removal
 
        if (distance(data.begin(), __end) <= factoryLimit) {
 
          int number_of_outliers = numberOfOutliers;
 
          if (number_of_outliers > NDF - 1) {
            number_of_outliers = NDF - 1;
          }
 
          double ymin = chi2 - STANDARD_DEVIATIONS * STANDARD_DEVIATIONS;
            
          for (int n = 1; n <= number_of_outliers; ++n, ymin -= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {
              
            sort(data.begin(), __end, hit_type::compare);
              
            buffer_type::iterator __end1 = prev(__end, n);
 
            do {
 
              try {
                
                (*this)(data.begin(), __end1);
                
                V.set(*this, data.begin(), __end1, gridAngle_deg, sigma_ns);
                Y.set(*this, data.begin(), __end1);
                
                V.invert();
                
                double y = getChi2(Y, V);
                
                if (y <= -(STANDARD_DEVIATIONS * STANDARD_DEVIATIONS)) {
 
                  WARNING(endl << "JMuonPrefit: invalid chi2 " << y << endl);
 
                } else {
 
                  if (y < 0.0) {
                    y = 0.0;
                  }
 
                  if (y < ymin) {
                    ymin = y;
                    tz   = *this;
                    chi2 = ymin;
                    NDF  = distance(data.begin(), __end1) - NUMBER_OF_PARAMETERS;
                    N    = getCount(data.begin(), __end1);
                  }
                }
              }
              catch(const exception& error) {}
 
            } while (next_permutation(data.begin(), __end1, __end, hit_type::compare));
          }
 
        } else {
 
          const int number_of_outliers = NDF - 1;
 
          for (int n = 1; n <= number_of_outliers; ++n) {
              
            double ymax =  0.0;
            int    k    = -1;
                
            for (size_t i = 0; i != Y.size(); ++i) {
                  
              double y = getChi2(Y, V, i);
                  
              if (y > ymax) {
                ymax = y;
                k    = i;
              }
            }
                
            if (ymax < STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {
              break;
            }
 
            try {
 
              V.update(k, +HIT_OFF);
 
              this->update(data.begin(), __end, V);
 
              Y.set(*this, data.begin(), __end);
 
              tz   = *this;
              chi2 = getChi2(Y, V);
              NDF -= 1;
              N   -= getCount(data[k]);
            }
            catch(const exception& error) {
              break;
            }
          }
        }
 
        tz.rotate_back(R);
 
        out.push_back(getFit(event(), tz, *dir, getQuality(chi2, N, NDF), NDF));
 
        // set additional values
 
        out.rbegin()->setW(JPP_COVERAGE_ORIENTATION, input.coverage.orientation);
        out.rbegin()->setW(JPP_COVERAGE_POSITION,    input.coverage.position);
      }
 
 
      if (numberOfPrefits > 0) {
 
        JEvt::iterator __end = out.end();
 
        if (Qwatershed > 0) {                         // sort distinct maxima
 
          __end = out.begin();                        // output
 
          for (JEvt::iterator p1 = out.begin(); p1 != out.end() && (size_t) distance(out.begin(), __end) < numberOfPrefits; ) {
 
            JEvt::iterator p2 = p1;                   // p2 -> best-of-the-rest
 
            for (JEvt::iterator i = p1; i != out.end(); ++i) {
              if (qualitySorter(*i, *p2)) {
                p2 = i;
              }
            }
 
            swap(*p2, *p1);                           // put best-of-the-rest at head
 
            p2 = p1++;                                // p2 = head; p1 = next to head
 
            const JPointing pointing(*p2);            // reference direction
 
            sort(p1, out.end(), pointing);            // sort according space angle with respect to reference direction
 
            for (double Q = p2->getQ();
                 p1 != out.end() && (p1->getQ() >= p2->getQ() - Qwatershed ||   // bottom
                                     p1->getQ() <= Q);                          // descend
                 ++p1) {
              Q = p1->getQ();
            }
 
            swap(*(__end++), *p2);                    // output head
          }
 
        } else {
 
          advance(__end = out.begin(), min(numberOfPrefits, out.size()));
 
          partial_sort(out.begin(), __end, out.end(), qualitySorter);
        }
 
        if (numberOfPostfits > 0) {
 
          __end = gridify(__end, out.end(), numberOfPostfits);
        }
 
        out.erase(__end, out.end());
 
      } else {
 
        sort(out.begin(), out.end(), qualitySorter);
      }
 
      return out;
    }
 
 
    /**
     * Auxiliary data structure for sorting of hits.
     */
    static const struct cmz {
      /**
       * Sort hits according times corrected for position along z-axis.
       *
       * \param  first             first  hit
       * \param  second            second hit
       * \return                   true if first hit earlier than second hit; else false
       */
      template<class T>
      inline bool operator()(const T& first, const T& second) const
      {
        using namespace JPP;
 
        return (first .getT() * getSpeedOfLight()  -  first .getZ()  <
                second.getT() * getSpeedOfLight()  -  second.getZ());
      }
    } cmz;
 
 
    JOmega3D    omega;
    int         debug;
 
    static constexpr double STANDARD_DEVIATIONS  =  3.0;  //!< number of standard deviations for outlier removal [unit]
 
  private:
    buffer_type data;
    JMatrixNZ   V;
    JVectorNZ   Y;
  };
}
 
#endif