JORCAShowerPrefit.hh

Go to the documentation of this file.

#ifndef JORCASHOWERPREFIT_INCLUDE
#define JORCASHOWERPREFIT_INCLUDE

/*TODO-suggested
  This code erases elements from a vector, that is slow! Try to use another STL container, as list 
 */

#include <string>
#include <iostream>
#include <vector>
#include <algorithm>
#include <memory>

#include "JDAQ/JDAQTimeslice.hh"
#include "JDAQ/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 "JFit/JFitStatus.hh"
#include "JFit/JEvt.hh"
#include "JFit/JEvtToolkit.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JPoint4DEstimator.hh"
#include "JFit/JPoint3DEstimator.hh"
#include "JFit/JEstimator.hh"
#include "JFit/JPoint4D.hh"
#include "JFit/JModel.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JPoint4DRegressor.hh"
#include "JFit/JShowerParameters.hh"

#include "JGeometry3D/JVector3D.hh"
#include "JTools/JRange.hh"
#include "JLang/JSharedPointer.hh"

/**
 * \author adomi, gmaggi
 */

namespace JFIT
{
  /**
   * class to handle first step of the shower reconstruction, mainly dedicated for ORCA
   */

  class JORCAShowerPrefit 
  {
  private:
    typedef JFIT::JRegressor<JPoint4D, JSimplex> JRegressor_t;
    std::shared_ptr<JRegressor_t> fit2_;

    JLANG::JSharedPointer<const JDETECTOR::JModuleRouter> moduleRouter_;
    JFIT::JShowerPrefitParameters_t parameters_;
    
  public:

    /**
     * Default constructor
     */
    JORCAShowerPrefit()
    {}

    /**
     * Parameterized constructor
     *
     * \param moduleRouter JSharedPointer of JModuleRouter, this is built via detector file.
     * \param parameters struct that holds default-optimized parameters for the reconstruction, which can be found in $JPP_DATA.
     */

    JORCAShowerPrefit(const JLANG::JSharedPointer<const JDETECTOR::JModuleRouter> &moduleRouter,
                      const JFIT::JShowerPrefitParameters_t &parameters):
      moduleRouter_(moduleRouter),
      parameters_(parameters)
    {
      JRegressor_t::debug  = 1;
      JRegressor_t::MAXIMUM_ITERATIONS = 10000; 

      fit2_ = std::make_shared<JRegressor_t>(true);
    }

    ~JORCAShowerPrefit(){}

    /**
     * Declaration of Member function that actually performs the reconstruction
     *
     * \param timeSliceBuildL0 which is contructed via a JDAQEvent
     * \param timeSliceBuildL2 which is contructed via a JDAQEvent
     * \param OutFits non const (output) JEvt.
     */

    void
    getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL0,
            const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL2,
            JFIT::JEvt &OutFits) const;
   
    /**
     * Member function used in the predicate of std::sort, to sort hits according to ToT and W
     *
     * \param first lhs hit
     * \param second rhs hit
     */

    static bool 
    totW(const JTRIGGER::JHitR1& first, 
         const JTRIGGER::JHitR1& second)
    {
      return (first.getToT() + first.getW() > second.getToT() + second.getW());
    }
    
  };
}

#endif

/**
 * member function definition
 */

void
JFIT::JORCAShowerPrefit::getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL0,
                                 const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL1,
                                 JFIT::JEvt &OutFits) const
{
  const double                   Tmax_ns   = parameters_.Tmax_ns;  // [ns] 
  const double                   ctMin     = parameters_.ctMin;  // optimal value
  const double                   Dmax_m    = parameters_.Dmax_m;  // [m] optimal value 
  const JTOOLS::JRange<double>   pos_grid  = parameters_.pos_grid;
  const JTOOLS::JRange<double>   time_grid = parameters_.time_grid;
  const double                   pos_step  = parameters_.pos_step;  // [m] density of the grid
  const double                   time_step = parameters_.time_step;  // [ns] density of the grid
  const double                   sigma_ns  = parameters_.sigma_ns;
  const int                      numberOfOutliers = parameters_.numberOfOutliers;
  const bool                     useL0     = parameters_.useL0;
  
  const int    FACTORY_LIMIT             = 40;
  const double STANDARD_DEVIATIONS       = 1.0;
  const int    NUMBER_OF_L1HITS          = (useL0 ? 1 : 4);

  typedef std::vector<JTRIGGER::JHitL0> JDataL0_t;
  typedef std::vector<JTRIGGER::JHitL1> JDataL1_t;

  const JTRIGGER::JBuildL0<JTRIGGER::JHitL0> buildL0;
  const JTRIGGER::JBuildL2<JTRIGGER::JHitL1> buildL2(2, Tmax_ns, ctMin);
  const JTRIGGER::JMatch3G<JTRIGGER::JHitL1> match3G(Dmax_m, Tmax_ns); // CAUSALITY FOR SHOWERS
  const JFIT::JHitL1Comparator compare;
  
  //define empty vectors to fill with JHitL0 and JHitL1
  JDataL0_t dataL0;
  JDataL0_t dataL0_selected;
  JDataL1_t data;
  
  buildL0(timeSliceBuildL0,  *moduleRouter_, back_inserter(dataL0));
  buildL2(timeSliceBuildL1, *moduleRouter_, back_inserter(data));
  
  JDataL1_t::iterator __end = clusterizeWeight(data.begin(), data.end(), match3G);
  
  if(std::distance(data.begin(), __end) >= NUMBER_OF_L1HITS){
    
    JDataL1_t dataL1;
    std::copy(data.begin(), data.end(), std::back_inserter(dataL1));
    JDataL1_t::iterator __end1 = dataL1.end();
    
    __end1 = clusterizeWeight(dataL1.begin(), __end1, match3G);
    
    if (useL0) {
      
      dataL1.erase(__end1, dataL1.end()); 
      dataL1.reserve(dataL1.size() + dataL0.size());
      __end1 = dataL1.end();
      
      for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
        
        const JTRIGGER::JHitL1 hit(*i);
        
        if (std::count_if(dataL1.begin(), __end1, JTRIGGER::JBind2nd(match3G,hit)) >= NUMBER_OF_L1HITS){
          dataL1.push_back(hit);            
        }
      }
      
      std::sort(__end1, dataL1.end(), compare);
      __end1 = clusterize(__end1, dataL1.end(), match3G);

    }
    
    // 4D fit
    JFIT::JEstimator<JFIT::JPoint4D> fit;
    JFIT::JPoint4D vx;
    double   chi2 = std::numeric_limits<double>::max();
    int      NDF  = 0;
    int      N    = 0;
    
    if(std::distance(dataL1.begin(), __end1) <= FACTORY_LIMIT){
      
      double ymin = std::numeric_limits<double>::max();
      
      JDataL1_t::iterator __end2 = __end1;
      
      for (int n = 0; n <= numberOfOutliers && std::distance(dataL1.begin(), __end2) > 
             JFIT::JEstimator<JFIT::JPoint4D>::NUMBER_OF_PARAMETERS; ++n, --__end2) {
                
        std::sort(dataL1.begin(), __end1, compare);
        
        do {
          try {
            
            fit(dataL1.begin(), __end2);
            double y = JFIT::getChi2(fit, dataL1.begin(), __end2, sigma_ns);
            
            if (y <= 0.0) {
              std::cout << "chi2(1) " << y << std::endl; // <- it was used WARNING, which gives a bug.
            } else if (y < ymin) {
              ymin = y;
              vx   = fit;
              NDF  = std::distance(dataL1.begin(), __end2) - JFIT::JEstimator<JFIT::JPoint4D>::NUMBER_OF_PARAMETERS;
              N    = JFIT::getCount(dataL1.begin(), __end2);
              chi2 = y;
            }
          }
          catch(JLANG::JException& error) {}
          
        } while (JTOOLS::next_permutation(dataL1.begin(), __end2, __end1, compare));
        
        ymin -= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS;
      }
      
    } else {
      
      try{      
        fit(dataL1.begin(), __end1);
        NDF  = std::distance(dataL1.begin(), __end1) - JEstimator<JFIT::JPoint4D>::NUMBER_OF_PARAMETERS;
        N    = JFIT::getCount(dataL1.begin(), __end1);
        vx   = fit;
        chi2 = JFIT::getChi2(fit, dataL1.begin(), __end1, sigma_ns);
      }
      catch(const JLANG::JException& error){}
    }
    
    if (NDF >= 0) {
            
      for(double x = pos_grid.getLowerLimit(); x <= pos_grid.getUpperLimit(); x += pos_step){
        for(double y = pos_grid.getLowerLimit(); y <= pos_grid.getUpperLimit(); y += pos_step){
          for(double z = pos_grid.getLowerLimit(); z <= pos_grid.getUpperLimit(); z += pos_step){
            for(double t = time_grid.getLowerLimit(); t <= time_grid.getUpperLimit(); t += time_step){

              const JGEOMETRY3D::JVector3D pos_shifted(vx.getX() + x, vx.getY() + y, vx.getZ() + z);
              const JFIT::JPoint4D point_shifted(pos_shifted, vx.getT() + t);
              
              dataL0_selected.clear();
              
              //select a subcluster of L0s
              const JFIT::JModel<JFIT::JPoint4D> match(vx, 80, JTOOLS::JTimeRange(-100, +100));
              
              for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {       
                
                const JHitL0 hit(*i);
                
                if(match(hit)){
                  dataL0_selected.push_back(hit);
                }
              }
              
              double simplex_step = 2;  // [m]
              fit2_->step.resize(4);
              fit2_->step[0] = JFIT::JPoint4D(JVector3D(simplex_step, 0.0, 0.0), 0.0);
              fit2_->step[1] = JFIT::JPoint4D(JVector3D(0.0, simplex_step, 0.0), 0.0);
              fit2_->step[2] = JFIT::JPoint4D(JVector3D(0.0, 0.0, simplex_step), 0.0);
              fit2_->step[3] = JFIT::JPoint4D(JVector3D(0.0, 0.0, 0.0), 20);
       
              chi2 = (*fit2_)(point_shifted, dataL0_selected.begin(), dataL0_selected.end());
       
              const JGEOMETRY3D::JVector3D fitPos(fit2_->value);
              const JGEOMETRY3D::JTrack3D fitResult(fitPos, JGEOMETRY3D::JVersor3Z(), fit2_->value.getT());

              const double energy(0);
              const JFIT::JFit &outFit = JFIT::getFit(JFIT::JHistory(JSHOWERPREFIT), 
                                                      fitResult, getQuality(chi2, (useL0 ? N : NDF)),
                                                      NDF, energy,JFIT::JFitStatus_t::OKAY );
                                              
              OutFits.push_back(outFit);
              OutFits.rbegin()->setW(13, chi2);
              OutFits.rbegin()->setW(14, NDF);
              
            }
          }
        }
      }
    }
    std::sort(OutFits.begin(), OutFits.end(), qualitySorter);
  }else{
    std::cout<<"Too few hits " << std::endl;// <- it was used DEBUG, which gives a bug.
  }
  
}