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 "TMatrixDSym.h"
#include "TMatrixDSymEigen.h"
 
#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 "JReconstruction/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 "JDynamics/JCoverage.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 KM3NETDAQ::JDAQEvent;
  using JDETECTOR::JModuleRouter;
  using JSUPPORT::JSummaryRouter;
  using JFIT::JRegressor;
  using JFIT::JLine3Z;
  using JFIT::JGandalf;
  using JFIT::JTimeRange;
  using JDYNAMICS::coverage_type;
 
 
  /**
   * 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 JHitW0                               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   data;
      coverage_type coverage;
    };
 
 
    /**
     * Constructor
     *
     * \param  parameters    parameters
     * \param  storage       storage
     * \param  debug         debug
     */
    JMuonGandalf(const JMuonGandalfParameters_t& parameters,
                 const storage_type&             storage,
                 const int                       debug = 0) :
      JMuonGandalfParameters_t(parameters),
      JRegressor_t(storage)
    {
      if (this->getRmax() < roadWidth_m) {
        roadWidth_m = this->getRmax();
      }
      
      JRegressor_t::debug              = debug;
      JRegressor_t::Vmax_npe           = VMax_npe;
      JRegressor_t::MAXIMUM_ITERATIONS = NMax;
 
      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();
    }
 
 
    /**
     * Get input data.
     *
     * \param  router        module router
     * \param  summary       summary data
     * \param  event         event
     * \param  in            start values
     * \param  coverage      coverage
     * \return               input data
     */
    input_type getInput(const JModuleRouter&  router,
                        const JSummaryRouter& summary,
                        const JDAQEvent&      event,
                        const JEvt&           in,
                        const coverage_type&  coverage) const
    {
      using namespace std;
      using namespace JTRIGGER;
 
      const JBuildL0<JHitL0>  buildL0;
 
      input_type     input(event.getDAQEventHeader(), in, coverage);
 
      vector<JHitL0> data;
 
      buildL0(event, router, true, back_inserter(data));
 
      for (const auto& hit : data) {
        input.data.push_back(hit_type(hit, summary.getRate(hit.getPMTIdentifier(), this->R_Hz)));
      }
 
      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;
 
      JEvent event(JMUONGANDALF);
 
      JEvt out;
      
      const buffer_type& data = input.data;
 
      // select start values
 
      JEvt in = input.in;
 
      in.select(numberOfPrefits, qualitySorter);
 
      if (!in.empty()) {
        in.select(JHistory::is_event(in.begin()->getHistory()));
      }
 
      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->getW(JSTART_LENGTH_METRES, 0.0) > 0.0) {
          Z_m = JZRange(track->getW(JSTART_ZMIN_M) + ZMin_m,
                        track->getW(JSTART_ZMAX_M) + ZMax_m);
        }
 
        const JModel<JLine1Z> match(tz, roadWidth_m, T_ns, Z_m);
 
        // hit selection based on start value
 
        buffer_type buffer;
          
        for (buffer_type::const_iterator i = data.begin(); i != data.end(); ++i) {
 
          hit_type hit(*i);
            
          hit.rotate(R);
 
          if (match(hit)) {
            buffer.push_back(hit);
          }
        }
          
        // select first hit
          
        sort(buffer.begin(), buffer.end(), JHitL0::compare);
          
        buffer_type::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 (this->JMuonGandalfParameters_t::E_GeV > 0.1)
            JRegressor_t::E_GeV  =  this->JMuonGandalfParameters_t::E_GeV;
          else
            JRegressor_t::E_GeV  =  track->getE();
 
          const double chi2 = (*this)(JLine3Z(tz), buffer.begin(), __end);
 
          // check error matrix
 
          bool status = true;
 
          for (size_t i = 0; i != this->V.size(); ++i) {
            if (std::isnan(this->V(i,i)) || this->V(i,i) <= 0.0) {
              status = false;
            }
          }
 
          if (status) {
 
            JTrack3D tb(this->value);
            
            tb.rotate_back(R);
 
            out.push_back(getFit(JHistory(track->getHistory(), event()), tb, getQuality(chi2), NDF));
 
            // set additional values
 
            const size_t N = this->V.size();
 
            TMatrixDSym  M(N);
 
            for (size_t row = 0; row != N; ++row) {
              for (size_t col = 0; col != N; ++col) {
                M(row,col) = this->V(row,col);
              }
            }
 
            const TMatrixDSymEigen E(M);
            const TVectorD&        Y = E.GetEigenValues();
 
            out.rbegin()->setV(this->V.size(), this->V);
 
            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_NUMBER_OF_HITS,       distance(buffer.begin(), __end));
            out.rbegin()->setW(JGANDALF_LAMBDA,               Y.GetNrows() != 0 ? Y[0] : 0.0);
            out.rbegin()->setW(JGANDALF_NUMBER_OF_ITERATIONS, this->numberOfIterations);
            out.rbegin()->setW(JPP_COVERAGE_ORIENTATION,      input.coverage.orientation);
            out.rbegin()->setW(JPP_COVERAGE_POSITION,         input.coverage.position);
          }
        }
      }
 
      {
        for (JEvt::iterator i = out.begin(); i != out.end(); ++i) {
 
          double Q = numeric_limits<double>::max();
 
          for (JEvt::iterator p = out.begin(); p != out.end(); ++p) {
            if (p != i) {
              if (getDot(getDirection(*i), getDirection(*p)) <= cosLR) {
                if (i->getQ() - p->getQ() < Q) {
                  Q = i->getQ() - p->getQ();
                }
              }
            }
          }
 
          i->setW(JGANDALF_LIKELIHOOD_RATIO, Q);
        }
      }
 
      // apply default sorter
 
      sort(out.begin(), out.end(), qualitySorter);
 
      copy(input.in.begin(), input.in.end(), back_inserter(out));
 
      return out;
    }
  };
}
 
#endif