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/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 "JDetector/JModuleRouter.hh"

#include "JGeometry3D/JRotation3D.hh"

#include "Jeep/JMessage.hh"


/**
 * \author mdejong, gmaggi
 */

namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }

namespace JRECONSTRUCTION {

  using JDETECTOR::JModuleRouter;
  using JFIT::JRegressor;
  using JFIT::JLine3Z;
  using JFIT::JSimplex;


  /**
   * 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();

    /**
     * Constructor 
     * 
     * \param  parameters    parameters
     * \param  router        module router
     * \param  debug         debug
     */
    JMuonSimplex(const JMuonSimplexParameters_t& parameters,
                 const JModuleRouter&            router,
                 const int                       debug = 0) :
      JMuonSimplexParameters_t(parameters),
      JRegressor_t(parameters.sigma_ns),
      router(router)
    {
      using namespace JFIT;

      this->estimator.reset(new JMEstimatorLorentzian());

      JRegressor_t::debug              = debug;
      JRegressor_t::MAXIMUM_ITERATIONS = 10000;
    }
       

    /**
     * Fit function.
     *
     * \param  event         event
     * \param  in            start values
     * \return               fit results
     */
    JEvt operator()(const KM3NETDAQ::JDAQEvent& event, const JEvt& in)
    {
      using namespace std;
      using namespace JFIT;
      using namespace JTRIGGER;

      const JBuildL0<hit_type> buildL0;
      const JBuildL2<hit_type> buildL2(JL2Parameters(2, TMaxLocal_ns, ctMin));

      buffer_type dataL0;
      buffer_type dataL1;

      const KM3NETDAQ::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 (*this)(dataL0, dataL1, in);
    }


    /**
     * Fit function.
     *
     * \param  dataL0        L0 hit data
     * \param  dataL1        L1 hit data
     * \param  in            start values
     * \return               fit results
     */
    JEvt operator()(const buffer_type& dataL0,
                    const buffer_type& dataL1, 
                    const JEvt& in)
    {
      using namespace std;
      using namespace JFIT;
      using namespace JGEOMETRY3D;

      JEvt out;

      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, bind2nd(equal_to<JDAQModuleIdentifier>(), 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()).add(JMUONSIMPLEX), tb, getQuality(chi2, NDF), NDF)); 

          out.rbegin()->setW(track->getW());
        }
      }

      return out;
    }

        
    const JModuleRouter& router;
  };
}

#endif