JShowerPositionFit.hh

Go to the documentation of this file.

#ifndef JSHOWERPOSITIONFIT_INCLUDE
#define JSHOWERPOSITIONFIT_INCLUDE
 
#include <string>
#include <iostream>
#include <set>
#include <vector>
#include <algorithm>
#include <memory>
#include <math.h>
 
#include "TMatrixDSym.h"
#include "TMatrixDSymEigen.h"
 
 
#include "km3net-dataformat/online/JDAQEvent.hh"
#include "km3net-dataformat/online/JDAQTimeslice.hh"
#include "km3net-dataformat/online/JDAQSummaryslice.hh"
  
#include "JTrigger/JHit.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JHitL1.hh"
#include "JTrigger/JHitR1.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JAlgorithm.hh"
#include "JTrigger/JMatch3G.hh"
#include "JTrigger/JBind2nd.hh"
 
#include "JSupport/JSummaryFileRouter.hh"
 
#include "JFit/JFitToolkit.hh"
#include "JFit/JEstimator.hh"
#include "JFit/JPoint4E.hh"
#include "JFit/JModel.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JShowerBrightPointRegressor.hh"
  
#include "JReconstruction/JHitW0.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JShowerParameters.hh"
  
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JShower3D.hh"
#include "JGeometry3D/JShower3E.hh"
#include "JGeometry3D/JTrack3D.hh"
#include "JGeometry3D/JTrack3E.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"
  
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTParametersMap.hh"
 
#include "JDynamics/JCoverage.hh"
 
#include "JTools/JPermutation.hh"
#include "JTools/JRange.hh"
 
 
/**
 * \author adomi, vcarretero
 */
namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }
  
namespace JRECONSTRUCTION {
 
  using KM3NETDAQ::JDAQEvent;
  using JDETECTOR::JModuleRouter;
  using JDETECTOR::JPMTParametersMap;
  using JSUPPORT::JSummaryRouter;
  using JFIT::JPoint4D;
  using JFIT::JPoint4E;
  using JFIT::JGandalf;
  using JFIT::JRegressor;
  using JDYNAMICS::coverage_type;
 
  /**
   * class to handle the second position fit of the shower reconstruction, mainly dedicated for ORCA
   */
  class JShowerPositionFit :
  
    public JShowerPositionFitParameters_t,
    public JRegressor<JPoint4E, JGandalf>
 
  {
    typedef JHitW0                               hit_type;
    typedef std::vector<hit_type>                buffer_type;
    typedef JRegressor<JPoint4E, JGandalf> JRegressor_t;
    using JRegressor_t::operator();
 
    std::vector<double> Ev;
 
  public:
    /**
     * 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;
    };
 
    /**
     * Parameterized constructor
     *
     * \param  parameters    struct that holds default-optimized parameters for the reconstruction
     * \param  storage       storage
     * \param  pmtParameters PMT parameters
     * \param  debug         debug
     */ 
    JShowerPositionFit(const JShowerPositionFitParameters_t&  parameters,
                       const storage_type&                    storage,
                       const JPMTParametersMap&               pmtParameters,
                       const int                              debug = 0):
      JShowerPositionFitParameters_t(parameters),
      JRegressor_t(storage),
      pmtParameters(pmtParameters)
    {
      using namespace JPP;
      
      JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
      JRegressor_t::Vmax_npe           = VMax_npe;
      JRegressor_t::MAXIMUM_ITERATIONS = NMax;
      JRegressor_t::EPSILON            = 1e-3;
      JRegressor_t::debug              = debug;
      
      if (Emin_GeV > Emax_GeV || En <= 1) {
        THROW(JException, "Invalid energy input " << Emin_GeV << ' ' << Emax_GeV << ' ' << En);
      }
 
      const double base = std::pow((Emax_GeV / Emin_GeV), 1.0 / (En - 1));
 
      for (int i = 0; i != En; ++i) {
        Ev.push_back(Emin_GeV * std::pow(base, i));
      }
      
      this->parameters.resize(5);
       
      this->parameters[0] = JPoint4E::pX();
      this->parameters[1] = JPoint4E::pY();
      this->parameters[2] = JPoint4E::pZ();
      this->parameters[3] = JPoint4E::pT();
      this->parameters[4] = JPoint4E::pE();
    }
 
    /**
     * 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 JPP;
      using namespace JTRIGGER;
      
      const JBuildL0<JHitL0> buildL0;
 
      input_type     input(event.getDAQEventHeader(), in, coverage);
 
      vector<JHitL0> data;
 
      buildL0(JDAQTimeslice(event, true), router, back_inserter(data));
 
      for (const auto& hit : data) {
 
        const JPMTIdentifier id(hit.getModuleID(), hit.getPMTAddress());
 
        const JPMTParameters& wip = pmtParameters.getPMTParameters(id);
 
        const int    type = wip.getType();
        const double QE   = wip.QE;
        const double R_Hz = summary.getRate(hit.getPMTIdentifier(), this->R_Hz);
 
        input.data.push_back(hit_type(hit, type, QE, 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(JSHOWERPOSITIONFIT);
      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 shower = in.begin(); shower != in.end(); ++shower) {
 
        JPoint4D vx(getPosition(*shower), shower->getT());
 
        const JFIT::JModel<JPoint4D> match(vx, DMax_m, JRegressor_t::T_ns);
 
        buffer_type buffer;
        
        for (buffer_type::const_iterator i = data.begin(); i != data.end(); ++i) {
  
          if (match(*i)) {
            buffer.push_back(*i);
          }
        }
 
        // select first hit
        
        sort(buffer.begin(), buffer.end(), JHitL0::compare);
        
        vector<hit_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
          for (vector<double>::iterator e = Ev.begin(); e != Ev.end(); ++e) {
 
            JPoint4E sh(vx, *e);
 
            double chi2 = (*this)(sh, buffer.begin(), __end);
          
            JShower3D result(JVertex3D(this->value.getPosition(), this->value.getT()), JDirection3D());
  
         
            out.push_back(getFit(JHistory(shower->getHistory(), event()), result, getQuality(chi2), NDF, this->value.getE()));
 
            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(JSHOWERPOSITIONFIT_ISOTROPIC_ENERGY, this->value.getE());  
            out.rbegin()->setW(JSHOWERPOSITIONFIT_LAMBDA,           Y.GetNrows() != 0 ? Y[0] : 0.0);
 
            // set additional values
            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;
    }
 
    const JPMTParametersMap& pmtParameters;
  };
}
 
#endif