JShowerDirectionPrefit.hh

Go to the documentation of this file.

#ifndef JSHOWERDIRECTIONPREFIT_INCLUDE
#define JSHOWERDIRECTIONPREFIT_INCLUDE
 
#include <string>
#include <iostream>
#include <set>
#include <vector>
#include <algorithm>
#include <memory>
 
#include "km3net-dataformat/online/JDAQTimeslice.hh"
#include "km3net-dataformat/online/JDAQSummaryslice.hh"
 
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTParametersMap.hh"
 
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JBuildL0.hh"
 
#include "JSupport/JSummaryRouter.hh"
 
#include "JFit/JPMTW0.hh"
#include "JFit/JShower3EZRegressor.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JPoint4D.hh"
#include "JFit/JModel.hh"
 
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JShowerDirectionPrefitParameters_t.hh"
#include "JReconstruction/JModuleL0.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorSubset.hh"
 
#include "JDynamics/JCoverage.hh"
 
#include "JMath/JConstants.hh"
 
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JVersor3Z.hh"
#include "JGeometry3D/JDirection3D.hh"
#include "JGeometry3D/JOmega3D.hh"
#include "JGeometry3D/JRotation3D.hh"
 
/**
 * \author adomi, vcarretero
 */
namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }
  
namespace JRECONSTRUCTION {
 
  using KM3NETDAQ::JDAQEvent;
  using KM3NETDAQ::JDAQSummaryFrame;
  using JDETECTOR::JModuleRouter;
  using JDETECTOR::JModule;
  using JDETECTOR::JPMTParametersMap;
  using JSUPPORT::JSummaryRouter;
  using JFIT::JRegressor;
  using JFIT::JEnergy;
  using JFIT::JShower3EZ;
  using JFIT::JAbstractMinimiser;
  using JDYNAMICS::coverage_type;
 
 
  /**
   * class to handle the direction fit of the shower reconstruction, mainly dedicated for ORCA
   */
  class JShowerDirectionPrefit : 
    public JShowerDirectionPrefitParameters_t,
    public JRegressor<JShower3EZ, JAbstractMinimiser>
  {
  public:
    typedef JRegressor<JShower3EZ, JAbstractMinimiser>  JRegressor_t;
    typedef JModuleL0                                   module_type;
    typedef std::vector<module_type>                    detector_type;
 
    using JRegressor_t::operator();
    using JRegressor_t::getH1;
    using JRegressor_t::getH0;
    
    std::vector<double> Ev;
    
    /**
     * 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;
      detector_type data;
      coverage_type coverage;
    };
 
    /**
     * Compute chi2 for a set of energies
     *
     * \param  Ev                 energy vector
     * \param  data               PMTs
     * \param  R                  Rotation
     * \param  sh                 shower
     */
    std::vector<double> computeChi2(
                                    const std::vector<double>& Ev,
                                    const std::vector<JFIT::JPMTW0>& data,
                                    const JGEOMETRY3D::JRotation3D& R,
                                    const JShower3EZ& sh
                                    ) {
      using namespace JPP;
 
      std::vector<double> chi2(Ev.size(), 0.0);
 
      for (auto i = data.begin(); i != data.end(); ++i) {
 
        JPMTW0 pmt = *i;
        pmt.rotate(R);
 
        auto H1 = getH1(sh, pmt);
        auto H0 = getH0(pmt.getR());
 
        const bool hit = pmt.getN() != 0;
        const double chi2_bg = getChi2(get_value(H0), hit);
 
        for (size_t j = 0; j < Ev.size(); ++j) {
          chi2[j] += getChi2(get_value(Ev[j] * H1  +  H0), hit) - chi2_bg; ; // H1 is linear with E, -log-lik ratio
        }
      }
 
      return chi2;
    }
    
    /**
     *  Find best energy iteratively increasing the upper threshold if there is suspicion the best-fit is outside
     *
     * \param  Ev                 energy vector
     * \param  data               PMTs
     * \param  R                  Rotation
     * \param  sh                 shower
     */
    std::pair<double, double> findBestEnergy(
                                             std::vector<double> Ev,
                                             const std::vector<JFIT::JPMTW0>& data,
                                             const JGEOMETRY3D::JRotation3D& R,
                                             const JShower3EZ& sh
                                             ) {
      using namespace JPP;
 
      const int max_extensions = 5; 
      int extension_count = 0;
 
      while (true) {
 
        std::vector<double> chi2 = computeChi2(Ev, data, R, sh);
 
        auto p = std::min_element(chi2.begin(), chi2.end());
        size_t idx = p - chi2.begin();
 
        bool at_upper_edge = (idx == Ev.size() - 1);
 
        bool monotonic_tail = false;
        if (Ev.size() >= 3) {
          size_t n = Ev.size();
          monotonic_tail = (chi2[n-3] > chi2[n-2]) &&
            (chi2[n-2] > chi2[n-1]);
        }
 
        if (!(at_upper_edge && monotonic_tail)) {
          return {Ev[idx], *p};
        }
 
        if (extension_count >= max_extensions) {
          return {Ev[idx], *p};
        }
 
        extension_count++;
 
        std::vector<double> Ev_ext;
        double E_start = Ev.back();
        double E_end   = 100.0 * E_start;
 
        size_t N = Ev.size();
 
        for (size_t k = 0; k < N; ++k) {
          double t = double (k)  / (N - 1);
          Ev_ext.push_back(E_start * std::pow(E_end / E_start, t));
        }
 
        Ev = std::move(Ev_ext);  
      }
    }
    /**
     * Parameterized constructor
     *
     * \param  parameters         parameters
     * \param  storage            storage
     * \param  pmtParameters      PMT parameters
     * \param  debug              debug
     */
    JShowerDirectionPrefit(const JShowerDirectionPrefitParameters_t&   parameters,
                           const storage_type&                         storage,
                           const JPMTParametersMap&                    pmtParameters,
                           const int                                   debug = 0):
      JShowerDirectionPrefitParameters_t(parameters),
      JRegressor_t(storage),
      pmtParameters(pmtParameters)
    {
      using namespace JPP;
      
      JRegressor_t::debug  = debug;
      JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
      JRegressor_t::Vmax_npe = VMax_npe;
          
      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));
      }
    }
 
    /**
     * 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;
 
      input_type   input(event.getDAQEventHeader(), in, coverage);
 
      const JBuildL0 <JHitR0>   buildL0;
      map<int, vector<JHitR0> > data;
 
      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()));
 
          buildL0(buffer, back_inserter(data[i->getModuleID()]));
        }
      }
 
      for (const auto& module : router.getReference()) {
        if (!module.empty()) {
          input.data.push_back(module_type(module, summary.getSummaryFrame(module.getID(), R_Hz), data[module.getID()]));
        }
      }
 
      return input;
    }
 
    /**
     * Fit function.
     *
     * \param  input         input data
     * \return               fit results
     */
    JEvt operator()(const input_type& input)
    {
      using namespace std;
      using namespace JPP;
 
      JEvent event(JSHOWERDIRECTIONPREFIT);
 
      JEvt out;
 
      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) {
 
        vector<JPMTW0> data;
 
        const JPosition3D vertex(getPosition(*shower)); 
        const double      time = shower->getT(); 
        
        for (const auto& module : input.data) {
          
          JPosition3D pos(module->getPosition());
 
          pos.sub(vertex);
          
          if (pos.getLength() <= DMax_m) {
            
            const double t1 = time + pos.getLength() * getInverseSpeedOfLight() * getIndexOfRefraction();
 
            for (size_t i = 0; i != module->size(); ++i) {
 
              if (module.getStatus(i)) {
 
                struct {
 
                  bool operator()(const JHitR0& hit) const
                  {
                    return (hit.getPMT() == pmt && T_ns(hit.getT()));
                  }
                    
                  const JTimeRange T_ns;
                  const size_t     pmt;
 
                } match = { JRegressor_t::T_ns + t1, i };
 
                const JPMTIdentifier id(module->getID(), i);
 
                const JPMTParameters& wip = pmtParameters.getPMTParameters(id);
 
                const size_t ns = count_if(module.begin(), module.end(), match);
                const double QE = wip.QE;
 
                JPMT pmt = module->getPMT(i);
 
                pmt.sub(vertex);
 
                data.push_back(JPMTW0(pmt, QE, module.frame.getRate(i), ns));
                
              }
            }
          }
        }
 
        JVector3D start_dir(0,0,0);
 
        for (vector<JPMTW0>::const_iterator i = data.begin(); i != data.end(); ++i) {
          if (i->getN() != 0) start_dir.add(i->getPosition());
        }
        
        JOmega3D scan_directions(JDirection3D(start_dir),
                                 JOmega3D::range_type(0.0,PI),
                                 scanAngle_deg * PI / 180.0);
        
        const JShower3EZ sh(JVertex3D(JVector3D(0,0,0), shower->getT()), JVersor3Z(), 1);
 
        for (JOmega3D_t::const_iterator dir = scan_directions.begin(); dir != scan_directions.end(); ++dir) {
 
          const JGEOMETRY3D::JRotation3D R(*dir);
          
          std::pair<double,double> result = findBestEnergy(Ev, data, R, sh);
          
          out.push_back(getFit(JHistory(shower->getHistory(), event()),
                               JShower3D(JVertex3D(vertex,time), JDirection3D(*dir)),
                               getQuality(result.second),
                               data.size(),
                               result.first));
 
          // 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