JMuonSimplex.hh

Go to the documentation of this file.

#ifndef  __JRECONSTRUCTION__JMUONSIMPLEX__
#define  __JRECONSTRUCTION__JMUONSIMPLEX__
 
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
 
#include "km3net-dataformat/online/JDAQEvent.hh"
#include "km3net-dataformat/online/JDAQTimeslice.hh"
#include "km3net-dataformat/online/JDAQSummaryslice.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"
 
#include "JTrigger/JHitR1.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL2.hh"
 
#include "JFit/JFitToolkit.hh"
#include "JFit/JLine1Z.hh"
#include "JFit/JLine3Z.hh"
#include "JFit/JModel.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JLine3ZRegressor.hh"
#include "JFit/JMEstimator.hh"
 
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JMuonSimplexParameters_t.hh"
 
#include "JLang/JPredicate.hh"
 
#include "JDetector/JModuleRouter.hh"
 
#include "JDynamics/JCoverage.hh"
 
#include "JGeometry3D/JRotation3D.hh"
 
#include "Jeep/JMessage.hh"
 
 
/**
 * \author mdejong, gmaggi
 */
 
namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }
 
namespace JRECONSTRUCTION {
 
  using KM3NETDAQ::JDAQEvent;
  using JDETECTOR::JModuleRouter;
  using JFIT::JRegressor;
  using JFIT::JLine3Z;
  using JFIT::JSimplex;
  using JDYNAMICS::coverage_type;
 
 
  /**
   * Wrapper class to make intermediate fit of muon trajectory.
   *
   * The JMuonSimplex fit uses one or more start values (usually taken from the output of JMuonPrefit) and
   * produces new start values for subsequent fits (usually JMuonGandalf).\n
   * All hits of which the PMT position lies within a set road width (JMuonSimplexParameters_t::roadWidth_m) and
   * time is within a set window (JMuonSimplexParameters_t::TMin_ns, JMuonSimplexParameters_t::TMax_ns) around the Cherenkov hypothesis are taken.\n
   * In case there are multiple hits from the same PMT is the specified window,
   * the first hit is taken and the other hits are discarded.\n
   * The chi-squared is based on an M-estimator of the time residuals.\n
   */
  struct JMuonSimplex :
    public JMuonSimplexParameters_t,
    public JRegressor<JLine3Z, JSimplex>
  {
    typedef JRegressor<JLine3Z, JSimplex>        JRegressor_t;
    typedef JTRIGGER::JHitR1                     hit_type;
    typedef std::vector<hit_type>                buffer_type;      
 
    using JRegressor_t::operator();
 
    /**
     * Input data type.
     */
    struct input_type :
      public JDAQEventHeader
    {
      /**
       * Default constructor.
       */
      input_type()
      {}
 
 
      /**
       * Constructor.
       *
       * \param  header         header
       * \param  in             start values
       * \param  coverage       coverage
       */
      input_type(const JDAQEventHeader& header,
                 const JEvt&            in,
                 const coverage_type&   coverage) :
         JDAQEventHeader(header),
         in(in),
         coverage(coverage)
      {}
 
      JEvt          in;
      buffer_type   dataL0;
      buffer_type   dataL1;
      coverage_type coverage;
     };
 
 
    /**
     * Constructor 
     * 
     * \param  parameters    parameters
     * \param  debug         debug
     */
    JMuonSimplex(const JMuonSimplexParameters_t& parameters,
                 const int                       debug = 0) :
      JMuonSimplexParameters_t(parameters),
      JRegressor_t(parameters.sigma_ns)
    {
      using namespace JFIT;
 
      this->estimator.reset(new JMEstimatorLorentzian());
 
      JRegressor_t::debug              = debug;
      JRegressor_t::MAXIMUM_ITERATIONS = NMax;
    }
       
 
    /**
     * Get input data.
     *
     * \param  router        module router
     * \param  event         event
     * \param  in            start values
     * \param  coverage      coverage
     * \return               input data
     */
    input_type getInput(const JModuleRouter& router,
                        const JDAQEvent&     event,
                        const JEvt&          in,
                        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(), in, coverage);
 
      buffer_type& dataL0 = input.dataL0;
      buffer_type& dataL1 = input.dataL1;
 
      const JDAQTimeslice timeslice(event, true);
      JSuperFrame2D<JHit> buffer;
 
      for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {
 
        if (router.hasModule(i->getModuleID())) {
 
          buffer(*i, router.getModule(i->getModuleID()));
 
          if (useL0) {
            buildL0(buffer, back_inserter(dataL0));
          }
          {
            buildL2(buffer, back_inserter(dataL1));
          }
        }
      }
 
      return input;
    }
 
 
    /**
     * Fit function.
     *
     * \param  input         input data
     * \return               fit results
     */
    JEvt operator()(const input_type& input)
    {
      using namespace std;
      using namespace JFIT;
      using namespace JGEOMETRY3D;
      using namespace JLANG;
 
      JEvent event(JMUONSIMPLEX);
 
      JEvt out;
 
      const buffer_type& dataL0 = input.dataL0;
      const buffer_type& dataL1 = input.dataL1;
 
      // select start values
 
      JEvt in = input.in;
 
      in.select(numberOfPrefits, qualitySorter);
 
      if (!in.empty()) {
        in.select(JHistory::is_event(in.begin()->getHistory()));
      }
 
      buffer_type data;
 
      data.reserve(dataL0.size() + 
                   dataL1.size());
      
      for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
 
        const JRotation3D     R (getDirection(*track));
        const JLine1Z         tz(getPosition (*track).rotate(R), track->getT());
        const JModel<JLine1Z> match(tz, roadWidth_m, JTimeRange(TMin_ns, TMax_ns));
            
        data.clear();
 
        for (buffer_type::const_iterator i = dataL1.begin(); i != dataL1.end(); ++i) {
 
          hit_type hit(*i);
              
          hit.rotate(R);
 
          if (match(hit)) {
            data.push_back(hit);
          }
        }
 
        if (useL0) {
 
          buffer_type::iterator __end = data.end();     
 
          for (buffer_type::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
              
            if (find_if(data.begin(), __end, make_predicate(&hit_type::getModuleID, i->getModuleID())) == __end) {
 
              hit_type hit(*i);
              
              hit.rotate(R);
              
              if (match(hit)) {
                data.push_back(hit);
              }
            }
          }
        }
 
 
        this->step.resize(5);
        
        this->step[0] = JLine3Z(JLine1Z(JVector3D(0.5, 0.0, 0.0), 0.0));
        this->step[1] = JLine3Z(JLine1Z(JVector3D(0.0, 0.5, 0.0), 0.0));
        this->step[2] = JLine3Z(JLine1Z(JVector3D(0.0, 0.0, 0.0), 1.0));
        this->step[3] = JLine3Z(JLine1Z(JVector3D(), 0.0), JVersor3Z(0.005, 0.0));
        this->step[4] = JLine3Z(JLine1Z(JVector3D(), 0.0), JVersor3Z(0.0, 0.005));
        
        const int NDF = getCount(data.begin(), data.end()) - this->step.size();
 
        if (NDF > 0) {
 
          const double chi2 = (*this)(JLine3Z(tz), data.begin(), data.end());
              
          JTrack3D tb(this->value);
              
          tb.rotate_back(R);
 
          out.push_back(getFit(JHistory(track->getHistory(), event()), tb, getQuality(chi2, NDF), NDF)); 
 
          // set additional values
 
          out.rbegin()->setW(track->getW());
          out.rbegin()->setW(JPP_COVERAGE_ORIENTATION, input.coverage.orientation);
          out.rbegin()->setW(JPP_COVERAGE_POSITION,    input.coverage.position);
        }
      }
 
      // apply default sorter
 
      sort(out.begin(), out.end(), qualitySorter);
 
      copy(input.in.begin(), input.in.end(), back_inserter(out));
 
      return out;
    }
  };
}
 
#endif