JShowerPrefit.hh

Go to the documentation of this file.

#ifndef JSHOWERPREFIT_INCLUDE
#define JSHOWERPREFIT_INCLUDE

#include <iostream>
#include <vector>
#include <algorithm>
#include <functional>

#include "km3net-dataformat/online/JDAQTimeslice.hh"
#include "km3net-dataformat/online/JDAQEvent.hh"

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

#include "JFit/JFitToolkit.hh"
#include "JFit/JPoint4DEstimator.hh"
#include "JFit/JEstimator.hh"
#include "JFit/JPoint4D.hh"

#include "JDetector/JModuleRouter.hh"

#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JShowerParameters.hh"

#include "JLang/JComparator.hh"

#include "JGeometry3D/JTrack3D.hh"

/**
 * \author adomi, gmaggi, vcarretero
 */

namespace JRECONSTRUCTION
{

  using JLANG::JSharedPointer;
  using JDETECTOR::JModuleRouter;

  /**
   * class to handle first step of the shower reconstruction in ORCA:
   * it reconstructs the shower vertex, intended as the shower brightest point, as the barycenter of the hits
   */

  class JShowerPrefit : public JShowerPrefitParameters_t
  {
    
  public:
    
    typedef JTRIGGER::JHitR1                     hit_type;
    typedef std::vector<hit_type>                buffer_type;

    /**
     * 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
     */

    JShowerPrefit(const JShowerPrefitParameters_t&  parameters,
                  const JModuleRouter&              router,
                  const int                         debug = 0):
      JShowerPrefitParameters_t(parameters),
      router(router)
    {}
    
    /**
     * Declaration of the operator that performs the reconstruction
     *
     * \param event which is a JDAQEvent
     */

    JEvt operator()(const KM3NETDAQ::JDAQEvent &event) const
    {

      using namespace std;
      using namespace JPP;

      const double STANDARD_DEVIATIONS       = 3.0;

      typedef JEstimator<JPoint4D> JEstimator_t;

      const JBuildL0<hit_type> buildL0;
      const JBuildL2<hit_type> buildL2(JL2Parameters(2, TMaxLocal_ns, ctMin));
      const JMatch3G<hit_type> match3G(DMax_m, TMaxExtra_ns);  // causality relation for showers 
  
      //define empty vectors to be filled with JHitL0 and JHitL1
      buffer_type dataL0;
      buffer_type dataL1;
  
      JEvt out;   

      buffer_type buffer;

      
      buildL0(JDAQTimeslice(event, true),  router, back_inserter(buffer));  // true  => snapshot
      buildL2(JDAQTimeslice(event, false), router, back_inserter(dataL1));  // false => triggered
      
      copy(dataL1.begin(), dataL1.end(), back_inserter(dataL0));

      if(dataL1.size() <= numberOfL1){
        for (buffer_type::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
          if (find_if(dataL1.begin(), dataL1.end(), match_t(*i, TMaxExtra_ns)) == dataL1.end()) {
            dataL0.push_back(*i); 
          } 
        }
      }
      for (buffer_type::const_iterator root = dataL1.begin(); root != dataL1.end(); ++root) {

        buffer_type data(1, *root);
      
        JBinder2nd<hit_type> matching = JBind2nd(match3G, *root);
      
        for (buffer_type::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
        
          if(( root->getModuleIdentifier() != i->getModuleIdentifier() ) && matching(*i)){
            data.push_back(*i);     
          }
        }
        
        buffer_type::iterator __end1 = clusterizeWeight(data.begin() + 1, data.end(), match3G);
      
        // 4D fit
        JEstimator_t fit;
        JPoint4D  vx;
        double    chi2 = numeric_limits<double>::max();
        int       NDF  = distance(data.begin(), __end1) - JEstimator_t::NUMBER_OF_PARAMETERS;
        int       N    = getCount(data.begin(), __end1);

        if(NDF > 0){
          if(distance(data.begin(), __end1) <= factoryLimit){
          
            double ymin = numeric_limits<double>::max();
      
            buffer_type::iterator __end2 = __end1;

            for (int n = 0; n <= numberOfOutliers && distance(data.begin(), __end2) > 
                   JEstimator_t::NUMBER_OF_PARAMETERS; ++n, --__end2) {
        
              sort(data.begin() + 1, __end1, hit_type::compare);

              do {
                try {
            
                  fit(data.begin(), __end2);
                  
                  double y = getChi2(fit, data.begin(), __end2, sigma_ns);

                  if (y < ymin) {
                    ymin = y;
                    vx   = fit;
                    chi2 = ymin;
                    NDF  = distance(data.begin(), __end2) - JEstimator_t::NUMBER_OF_PARAMETERS;
                    N    = getCount(data.begin(), __end2);
                  }
                }
                catch(JException& error) { }
          
              } while (next_permutation(data.begin() + 1, __end2, __end1, hit_type::compare));
        
              ymin -= STANDARD_DEVIATIONS * STANDARD_DEVIATIONS;
            }
            
          } else {

            const int number_of_outliers = distance(data.begin(), __end1) - JEstimator_t::NUMBER_OF_PARAMETERS - 1;

            buffer_type::iterator __end2 = __end1;

            for (int n = 0; n <= number_of_outliers; ++n) {       

              try{      

                fit(data.begin(), __end2);
                vx   = fit;
                chi2 = getChi2(fit, data.begin(), __end2, sigma_ns);
                NDF  = distance(data.begin(), __end2) - JEstimator_t::NUMBER_OF_PARAMETERS;
                N    = getCount(data.begin(), __end2);
                
              }
              catch(const JException& error){ }

              double                ymax = 0;
              buffer_type::iterator imax = __end2;

              for (buffer_type::iterator i = data.begin() + 1; i != __end2; ++i) {

                double y = getChi2(fit, *i, sigma_ns);

                if (y > ymax) {
                  ymax = y;
                  imax = i;
                }
              }

              if (ymax > STANDARD_DEVIATIONS * STANDARD_DEVIATIONS) {   
                --__end2;
                swap(*imax, *__end2);     
              } else {
                break;
              }
            }   
          }
  
          out.push_back(getFit(JHistory(JSHOWERPREFIT), JTrack3D(vx, JGEOMETRY3D::JVersor3Z()), getQuality(chi2, N), NDF));
        
        }
      }

      out.select(numberOfGrids, qualitySorter);

      size_t solutions = out.size();

      for(size_t i=0; i < solutions; i++){
        for(int x = -pos_grid_m; x < pos_grid_m + pos_step_m/2.; x += pos_step_m){
          for(int y = -pos_grid_m; y < pos_grid_m + pos_step_m/2.; y += pos_step_m){
            for(int z = -pos_grid_m; z < pos_grid_m + pos_step_m/2.; z += pos_step_m){
              for(int t = -time_grid_ns; t < time_grid_ns + time_step_ns/2.; t += time_step_ns){
                if (x != 0 || y != 0 || z != 0 || t != 0) {
                  
                  out.push_back(getFit(JHistory(JSHOWERPREFIT),
                                       JTrack3D(JPoint4D(getPosition(out[i]) + JPosition3D(x,y,z), out[i].getT()+t), JVersor3Z()),
                                       0,
                                       0));
                   
                }
              }
            }
          }
        }
      }
      return out;
    
    }
        
    const JModuleRouter& router;

    /**
     * Auxiliary class to match hit to root hit.
     */
    struct match_t {
      /**
       * Constructor.
       *
       * \param  root        root hit
       * \param  TMax_ns     maximal time [ns]
       */
      match_t(const hit_type& root, const double TMax_ns) :
        root(root),
        TMax_ns(TMax_ns)
      {}

      /**
       * Test match.
       *
       * \param  hit         hit
       * \return             true if mach; else false
       */
      bool operator()(const hit_type& hit) const
      {
        return root.getModuleID() == hit.getModuleID() && fabs(root.getT() - hit.getT()) <= TMax_ns;
      }

      const hit_type& root;
      double TMax_ns;
    };
  };
}

#endif