JShowerPointSimplex.hh

Go to the documentation of this file.

#ifndef JSHOWERPOINTSIMPLEX_INCLUDE
#define JSHOWERPOINTSIMPLEX_INCLUDE
 
#include <string>
#include <iostream>
#include <set>
#include <vector>
#include <algorithm>
#include <memory>
#include <math.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/JHitR1.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JAlgorithm.hh"
#include "JTrigger/JMatch3G.hh"
#include "JTrigger/JBind2nd.hh"
  
#include "JFit/JFitToolkit.hh"
#include "JFit/JEstimator.hh"
#include "JFit/JPoint4D.hh"
#include "JFit/JModel.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JGandalf.hh"
#include "JFit/JPoint4DRegressor.hh"
 
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JShowerParameters.hh"
 
#include "JLang/JPredicate.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/JDetector.hh"
#include "JDetector/JDetectorSubset.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 JDYNAMICS::coverage_type;
  using JFIT::JPoint4D;
  using JFIT::JGandalf;
  using JFIT::JRegressor;
 
  /**
   * class to handle the second position fit of the shower reconstruction, mainly dedicated for ORCA
   */
  class JShowerPointSimplex :
  
    public JShowerPointSimplexParameters_t,
    public JRegressor<JPoint4D, JGandalf>
  { 
  public:
    typedef JRegressor<JPoint4D, JGandalf>       JRegressor_t;
    typedef JTRIGGER::JHitR1                     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  dataL0;
      buffer_type  dataL1;
      coverage_type coverage;
    };
  
  public:
   
    /**
     * Parameterized constructor
     *
     * \param parameters    struct that holds default-optimized parameters for the reconstruction.
     * \param debug         debug
     */ 
    JShowerPointSimplex(const JShowerPointSimplexParameters_t&  parameters,
                        const int                                debug = 0):
      JShowerPointSimplexParameters_t(parameters),
      JRegressor_t(parameters.sigma_ns)
    {
      using namespace JPP;
      
      JRegressor_t::debug  = debug;
      JRegressor_t::MAXIMUM_ITERATIONS = NMax;
      JRegressor_t::EPSILON = 1e-4;
 
      this->estimator.reset(getMEstimator(parameters.mestimator));
      this->parameters.resize(4);
      this->parameters[0] = JPoint4D::pX();
      this->parameters[1] = JPoint4D::pY();
      this->parameters[2] = JPoint4D::pZ();
      this->parameters[3] = JPoint4D::pT();
    }
    /**
     * Get input data.
     *
     * \param  router        module router
     * \param  event         event
     * \param  in            start values
     * \param  coverage      coverage
     * \return               input data
     */
    input_type getInput(const JModuleRouter& router, const KM3NETDAQ::JDAQEvent& event, const JEvt& in, const coverage_type& coverage) const
    {
      using namespace std;
      using namespace KM3NETDAQ;
      using namespace JTRIGGER; 
      const JBuildL0<hit_type> buildL0;
      const JBuildL2<hit_type> buildL2(JL2Parameters(2, TMaxLocal_ns, ctMin));
 
      input_type   input(event.getDAQEventHeader(), in, coverage);
      
      buffer_type& dataL0 = input.dataL0;
      buffer_type& dataL1 = input.dataL1;
  
      buildL0(JDAQTimeslice(event, true), router, back_inserter(dataL0));
      buildL2(JDAQTimeslice(event, false), router, back_inserter(dataL1));  // false = triggered  
 
      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(JSHOWERPOINTSIMPLEX);
 
      JEvt out;
 
      const buffer_type& dataL0 = input.dataL0;
      const buffer_type& dataL1 = input.dataL1;
 
      // 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());
 
        JFIT::JModel<JPoint4D> match(vx, DMax_m, JTimeRange(TMin_ns, TMax_ns));
 
        buffer_type data;
 
        data.reserve(dataL0.size() + dataL1.size());
 
        for (buffer_type::const_iterator i = dataL1.begin(); i != dataL1.end(); ++i) {
          
          if (match(*i)) {
            data.push_back(*i);
          }
        }
 
        buffer_type::iterator __end = data.end();     
        
        for (buffer_type::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
          if (find_if(data.begin(), __end, make_predicate(&hit_type::getModuleID, i->getModuleID())) == __end) {
            if (match(*i)) {
              data.push_back(*i);
            }
          }
        }
        
        const int NDF = getCount(data.begin(), data.end()) - this->parameters.size();
        
        if (NDF > 0) {
 
          double chi2 = (*this)(vx, data.begin(), data.end());
          
          JShower3D result(JVertex3D(this->value.getPosition(), this->value.getT()), JDirection3D());
 
          out.push_back(getFit(JHistory(shower->getHistory(), event()), result, getQuality(chi2, NDF), NDF, shower->getE()));
 
          // 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;
    }
 
    
  };
}
#endif