JMuonGandalf.hh

Go to the documentation of this file.

#ifndef __JRECONSTRUCTION__JMUONGANDALF__
#define __JRECONSTRUCTION__JMUONGANDALF__

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>

#include "km3net-dataformat/online/JDAQEvent.hh"
#include "km3net-dataformat/online/JDAQSummaryslice.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"

#include "JTrigger/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"

#include "JFit/JFitToolkit.hh"
#include "JFit/JLine1Z.hh"
#include "JFit/JLine3Z.hh"
#include "JFit/JModel.hh"
#include "JFit/JGandalf.hh"
#include "JFit/JLine3ZRegressor.hh"
#include "JFit/JHitW0.hh"

#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JMuonGandalfParameters_t.hh"

#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JGeane.hh"

#include "JDetector/JModuleRouter.hh"

#include "JSupport/JSummaryRouter.hh"

#include "JGeometry3D/JRotation3D.hh"

#include "JTools/JRange.hh"

#include "Jeep/JMessage.hh"


/**
 * \author mdejong, gmaggi
 */

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

namespace JRECONSTRUCTION {

  using JDETECTOR::JModuleRouter;
  using JSUPPORT::JSummaryRouter;
  using JFIT::JRegressor;
  using JFIT::JLine3Z;
  using JFIT::JGandalf;


  /**
   * Wrapper class to make final fit of muon trajectory.
   *
   * The JMuonGandalf fit uses one or more start values (usually taken from the output of JMuonSimplex).\n
   * All hits of which the PMT position lies within a set road width (JMuonGandalfParameters_t::roadWidth_m) and
   * time is within a set window (JMuonGandalfParameters_t::TMin_ns, JMuonGandalfParameters_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 PDF is accordingly evaluated, i.e.\ the normalised probability for a first hit at the given time of the hit is taken.
   * The normalisation is consistently based on the specified time window.\n
   * Note that this hit selection is unbiased with respect to the PDF of a single PMT.
   */
  struct JMuonGandalf :
    public JMuonGandalfParameters_t,
    public JRegressor<JLine3Z, JGandalf>
  {
    typedef JRegressor<JLine3Z, JGandalf>        JRegressor_t;
    typedef JTRIGGER::JHitL0                     hit_type;
    typedef std::vector<hit_type>                buffer_type;

    using JRegressor_t::operator();

    /**
     * Constructor
     *
     * \param  parameters    parameters
     * \param  router        module router
     * \param  summary       summary file router
     * \param  pdf_file      PDF file
     * \param  debug         debug
     */
    JMuonGandalf(const JMuonGandalfParameters_t& parameters,
                 const JModuleRouter&            router,
                 const JSummaryRouter&           summary,
                 const std::string&              pdf_file,
                 const int                       debug = 0) :
      JMuonGandalfParameters_t(parameters),
      JRegressor_t(pdf_file, parameters.TTS_ns),
      router (router),
      summary(summary)
    {
      using namespace JPP;

      if (this->getRmax() < roadWidth_m) {
        roadWidth_m = this->getRmax();
      }

      JRegressor_t::debug              = debug;
      JRegressor_t::T_ns.setRange(TMin_ns, TMax_ns);
      JRegressor_t::Vmax_npe           =   10.0;
      JRegressor_t::MAXIMUM_ITERATIONS = 1000;
      
      this->parameters.resize(5);
      
      this->parameters[0] = JLine3Z::pX();
      this->parameters[1] = JLine3Z::pY();
      this->parameters[2] = JLine3Z::pT();
      this->parameters[3] = JLine3Z::pDX();
      this->parameters[4] = JLine3Z::pDY();
    }


    /**
     * 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 JPP;

      const JBuildL0<hit_type>  buildL0;

      buffer_type dataL0;

      buildL0(event, router, true, back_inserter(dataL0));

      return (*this)(dataL0, in);
    }


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

      JEvt out;
      
      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());
        JRange<double>     Z_m;
      
        if (track->hasW(JSTART_LENGTH_METRES)) {
          Z_m = JZRange(ZMin_m, ZMax_m + track->getW(JSTART_LENGTH_METRES));
        }

        const JModel<JLine1Z> match(tz, roadWidth_m, JRegressor_t::T_ns, Z_m);

        // hit selection based on start value

        vector<JHitW0> buffer;
          
        for (buffer_type::const_iterator i = data.begin(); i != data.end(); ++i) {
            
          const double rate_Hz = summary.getRate(*i);

          JHitW0 hit(*i, rate_Hz);
            
          hit.rotate(R);

          if (match(hit)) {
            buffer.push_back(hit);
          }
        }
          
        // select first hit
          
        sort(buffer.begin(), buffer.end(), compare);
          
        vector<JHitW0>::iterator __end = unique(buffer.begin(), buffer.end(), equal_to<JDAQPMTIdentifier>());
          

        const int NDF = distance(buffer.begin(), __end) - this->parameters.size();

        if (NDF > 0) {

          // set fit parameters

          if (track->getE() > 0.1)
            JRegressor_t::E_GeV   = track->getE();
          else
            JRegressor_t::E_GeV   = this->JMuonGandalfParameters_t::E_GeV;

          const double chi2 = (*this)(JLine3Z(tz), buffer.begin(), __end);
            
          JTrack3D tb(this->value);
            
          tb.rotate_back(R);

          out.push_back(getFit(JHistory(track->getHistory()).add(JMUONGANDALF), tb, getQuality(chi2), NDF));

          // set additional values

          out.rbegin()->setW(track->getW());
          out.rbegin()->setW(JGANDALF_BETA0_RAD,            sqrt(this->error.getDX() * this->error.getDX()  +
                                                                 this->error.getDY() * this->error.getDY()));
          out.rbegin()->setW(JGANDALF_BETA1_RAD,            sqrt(this->error.getDX() * this->error.getDY()));
          out.rbegin()->setW(JGANDALF_CHI2,                 chi2);
          out.rbegin()->setW(JGANDALF_NUMBER_OF_HITS,       distance(buffer.begin(), __end));
          out.rbegin()->setW(JGANDALF_LAMBDA,               this->lambda);
          out.rbegin()->setW(JGANDALF_NUMBER_OF_ITERATIONS, this->numberOfIterations);
        }
      }

      return out;
    }


    /**
     * Auxiliary data structure for sorting of hits.
     */
    static const struct {
      /**
       * Compare hits by PMT identifier and time.
       *
       * \param  first         first  hit
       * \param  second        second hit
       * \return               true if first before second; else false
       */
      template<class T>
      bool operator()(const T& first, const T& second) const
      {
        using namespace std;
        using namespace JPP;
        using namespace KM3NETDAQ;

        if (equal_to<JDAQPMTIdentifier>()(first, second))
          return less<JHit>()(first, second);
        else
          return less<JDAQPMTIdentifier>()(first, second);
      }
    } compare;

    const JModuleRouter&  router;
    const JSummaryRouter& summary;
    int                   debug;
  };
}

#endif