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/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/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 "JTools/JPermutation.hh"
#include "JTools/JRange.hh"

/**
 * \author adomi
 */
namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }
  
namespace JRECONSTRUCTION {
 
  using JDETECTOR::JModuleRouter;
  using JFIT::JPoint4D;
  using JFIT::JSimplex;
  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, JSimplex>

  {
 
    typedef JRegressor<JPoint4D, JSimplex> JRegressor_t;
    using JRegressor_t::operator();
    typedef JTRIGGER::JHitR1                     hit_type;
    typedef std::vector<hit_type>                buffer_type;  
    
  public:
   
    /**
     * Parameterized constructor
     *
     * \param parameters    struct that holds default-optimized parameters for the reconstruction.
     * \param router        module router, this is built via detector file.
     * \param debug         debug
     */ 
    JShowerPointSimplex(const JShowerPointSimplexParameters_t&  parameters,
                       const JModuleRouter&                     router,
                       const int                                debug = 0):
      JShowerPointSimplexParameters_t(parameters),
      JRegressor_t(parameters.sigma_ns),
      router(router)
    {
      using namespace JPP;
      
      JRegressor_t::debug  = 1;
      JRegressor_t::MAXIMUM_ITERATIONS = NMax;

      this->estimator.reset(new JMEstimatorTukey(tukey_std_dev));
    }
    
    /**
     * Declaration of the member function that actually performs the reconstruction
     *
     * \param event which is a JDAQEvent
     * \param in input fits
     */ 
    JEvt operator()(const KM3NETDAQ::JDAQEvent& event, const JFIT::JEvt &in)
    {
      using namespace std;
      using namespace JPP;
        
      const JBuildL0<hit_type> buildL0;
      const JBuildL2<hit_type> buildL2(JL2Parameters(2, TMaxLocal_ns, ctMin));
      
      buffer_type dataL0;
      buffer_type dataL1;
      
      JEvt out;
        
      buildL0(JDAQTimeslice(event, true), router, back_inserter(dataL0));
      buildL2(JDAQTimeslice(event, false), router, back_inserter(dataL1));  // false = triggered  

      for (JEvt::const_iterator shower = in.begin(); shower != in.end(); ++shower) {

        JPoint4D vx(getPosition(*shower), shower->getT());

        buffer_type data;
  
        data.reserve(dataL0.size() + dataL1.size());

        JFIT::JModel<JPoint4D> match(vx, DMax_m, JTimeRange(TMin_ns, TMax_ns));

        JRange<double>   posGrid_m(0, 0);
        JTimeRange       timeGrid_ns(0, 0);

        /*
         * part to help the reco of events with few hits:
         * in this case a bigger time window in the hit selection and
         * a grid in position and time are considered.
         */
        if(in.size() <= minPrefitsSize){
          
          match = JFIT::JModel<JPoint4D>(vx, DMax_m, TWindow_ns);

          posGrid_m   = pos_grid_m;
          timeGrid_ns = time_grid_ns;
        }

        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->step.size();

        if (NDF > 0) {

          for(double x = posGrid_m.getLowerLimit(); x <= posGrid_m.getUpperLimit(); x += pos_step_m){
            for(double y = posGrid_m.getLowerLimit(); y <= posGrid_m.getUpperLimit(); y += pos_step_m){
              for(double z = posGrid_m.getLowerLimit(); z <= posGrid_m.getUpperLimit(); z += pos_step_m){
                for(double t = timeGrid_ns.getLowerLimit(); t <= timeGrid_ns.getUpperLimit(); t += time_step_ns){
                  
                  JPoint4D vxs(JVector3D(vx.getX() + x, vx.getY() + y, vx.getZ() + z), vx.getT() + t);
              
                  this->step.resize(4);
                  this->step[0] = JPoint4D(JVector3D(simplex_step_m, 0.0, 0.0), 0.0);
                  this->step[1] = JPoint4D(JVector3D(0.0, simplex_step_m, 0.0), 0.0);
                  this->step[2] = JPoint4D(JVector3D(0.0, 0.0, simplex_step_m), 0.0);
                  this->step[3] = JPoint4D(JVector3D(), simplex_step_ns);

                  double chi2 = (*this)(vxs, data.begin(), data.end());
            
                  JShower3E sh_fit(this->value.getPosition(), JDirection3D(),
                                   this->value.getT(), shower->getE());

                  out.push_back(getFit(JHistory(shower->getHistory()).add(JSHOWERPOINTSIMPLEX), sh_fit, 
                                       getQuality(chi2, NDF), NDF, sh_fit.getE())); 
                  
                }
              }
            }
          }
        }
      }
      return out;
    }
    const JModuleRouter& router;
  };
}
#endif