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

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:
    
    /**
     * 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 JHitR1           hit_type;
      typedef vector<hit_type> buffer_type;

      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(roadWidth_m, TMaxLocal_ns);  // causality relation for showers 
      const JHitL1Comparator compare;
  
      //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));

      for (buffer_type::const_iterator i = buffer.begin(); i != buffer.end(); ++i) {

        if (find_if(dataL1.begin(), dataL1.end(), match_t(*i, TMaxLocal_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(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, 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, 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;
            }
          }     
        }

        if (NDF >= 0) {
        
          const JGEOMETRY3D::JTrack3D fitResult(vx, JGEOMETRY3D::JVersor3Z());
        
          out.push_back(getFit(JHistory(JSHOWERPREFIT), fitResult, getQuality(chi2, N), NDF));
        
          out.rbegin()->setW(13, chi2);
          out.rbegin()->setW(14, N);           
        }
      }
   
      return out;
    
    }
        
    const JModuleRouter& router;

    struct match_t {

    match_t(const JHitR1& root, const double TMax_ns) :
      root(root),
      TMax_ns(TMax_ns)
    {}

    bool operator()(const JHitR1& hit) const
    {
      return root.getModuleID() == hit.getModuleID() && fabs(root.getT() - hit.getT()) <= TMax_ns;
    }

    const JHitR1& root;
    double TMax_ns;
  };    
 
  };
}

#endif