JORCAShowerPositionFit.hh

Go to the documentation of this file.

#ifndef JORCASHOWERPOSITIONFIT_INCLUDE
#define JORCASHOWERPOSITIONFIT_INCLUDE

/*TODO-suggested
  This code erases elements from a std::vector, that is slow!, try to use another STL container
*/

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

#include "TFile.h"

#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 "JDetector/JModuleRouter.hh"
#include "JGeometry3D/JVector3D.hh"
#include "JTools/JPermutation.hh"
#include "JTools/JRange.hh"
#include "JLang/JSharedPointer.hh"

/**
 * \author adomi, gmaggi
*/

namespace JFIT
{
  /**
   * class to handle the second step of the shower reconstruction, mainly dedicated for ORCA
   */
  class JORCAShowerPositionFit
  {
  private:

    typedef JFIT::JRegressor<JFIT::JPoint4D, JFIT::JSimplex> JRegressor_t;
    std::shared_ptr<JRegressor_t> fit0_;
    std::shared_ptr<JRegressor_t> fit_;

    JLANG::JSharedPointer<const JDETECTOR::JModuleRouter> moduleRouter_;
    JShowerPositionFitParameters_t parameters_;

  public:
    /**
     * Default constructor
     */
    JORCAShowerPositionFit()
    {}

    /**
     * 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.
     * \param pdfFile       file name with PDF
     */

    JORCAShowerPositionFit(const JLANG::JSharedPointer<const JDETECTOR::JModuleRouter> &moduleRouter,
                           const JShowerPositionFitParameters_t &parameters,
                           const std::string pdfFile):
      moduleRouter_(moduleRouter),
      parameters_(parameters)
    {
      TFile* pdf_file = new TFile(pdfFile.c_str());
      TH2D* hpdf = (TH2D*)pdf_file->Get("hPDF2Dist");

      JRegressor_t::debug = 1;
      JRegressor_t::MAXIMUM_ITERATIONS = 10000;

      fit0_ = std::make_shared<JRegressor_t>(true);
      fit0_->estimator.reset(new JMEstimatorTukey(10.0));

      fit_ = std::make_shared<JRegressor_t>(hpdf, false);
    }

    ~JORCAShowerPositionFit(){}

    /**                                                                                    
     * Container of vertex hypothesis:                                                            
     * During the fit the algorithm creates a grid in position space (x, y, z, t),          
     * for each of these points a vertex hypothesis is saved: 4D-Point and LogLikelihood value.   
     */
    
    struct JVertexHypothesis
    {
      JPoint4D  Vertex;      
      double    LogLik;
    };
    

    /**                                                
     * Function to sort different vertex hypotheses by their likelihood.                  
     */    
    
    static bool sortLogLik(JVertexHypothesis A, JVertexHypothesis B)
    {
      return(A.LogLik < B.LogLik);
    }
    
    /**
     * 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 InPreFits input fits (normally made by JORCAShowerPrefit)
     * \param OutFits output fits
     */

    void getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL0,
                 const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL2,
                 JFIT::JEvt &InPreFits,
                 JFIT::JEvt &OutFits) const 
    {
      
      if ( InPreFits.empty() ) return;

      const int            numberOfPrefits   = parameters_.numberOfPrefits; 
      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 double         pos_step          = parameters_.pos_step;  // [m] density of the grid  
      const JTOOLS::JRange<double> time_grid = parameters_.time_grid;
      const double         time_step         = parameters_.time_step;  // [ns] density of the grid
      
      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(3, Tmax_ns, ctMin);
      const JTRIGGER::JMatch3G<JTRIGGER::JHitL1>  match3G(Dmax_m, Tmax_ns); // CAUSALITY FOR SHOWERS         
      
            
      JEvt::iterator __end = InPreFits.end();
      if (numberOfPrefits > 0) {
        std::advance(__end = InPreFits.begin(), std::min(static_cast<std::size_t>(numberOfPrefits), InPreFits.size()));
      }
      std::partial_sort(InPreFits.begin(), __end, InPreFits.end(), qualitySorter);
      
      if(!InPreFits.empty()) {
        
        std::copy(InPreFits.begin(), __end, std::back_inserter(OutFits));
        
        JDataL0_t dataL0;
        JDataL0_t dataL0_selected;
        JDataL1_t data;
        
        buildL0(timeSliceBuildL0, *moduleRouter_, std::back_inserter(dataL0));  
        buildL2(timeSliceBuildL2,*moduleRouter_, std::back_inserter(data));

        for (JEvt::const_iterator shower = InPreFits.begin(); shower != __end; ++shower) {
          
          const JFIT::JPoint4D vertex(JFIT::getPosition(*shower), shower->getT());
          
          for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
            
            const JTRIGGER::JHitL0 hit(*i);
            const JGEOMETRY3D::JPosition3D hit_pos(hit.getPosition());
            const double D     = hit_pos.getDistance(vertex);
            const double t_res = hit.getT() - vertex.getT(hit_pos);
            const JGEOMETRY3D::JDirection3D photonDir(hit_pos - vertex.getPosition());
            const double cosT  = photonDir.getDot(hit.getDirection());
            
            if(data.size()>0){
              if(D < Dmax_m && (t_res > -40 && t_res < 40) && (cosT >= -1 && cosT <= 0.1)){
                dataL0_selected.push_back(hit);
              }
            } else {
              if(D < Dmax_m && (t_res > -50 && t_res < 50) && (cosT >= -1 && cosT <= 0.1)){
                dataL0_selected.push_back(hit);
              }
            }
            
          }
          
          if(dataL0_selected.size() != 0){
            
            double chi2 = 0.0;
            
            JFIT::JORCAShowerPositionFit::JVertexHypothesis vertex_hypothesis;
            std::vector<JFIT::JORCAShowerPositionFit::JVertexHypothesis> vertex_hypotheses;
            
            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){
                    
                    chi2 = 0.0;
                    
                    const JGEOMETRY3D::JPosition3D pos_shifted(vertex.getX() + x, vertex.getY() + y, vertex.getZ() + z);
                    const JFIT::JPoint4D point_shifted(pos_shifted, vertex.getT() + t);
                    
                    for(JDataL0_t::const_iterator i = dataL0_selected.begin(); i != dataL0_selected.end(); ++i){
                      const JHitL0 hit(*i);
                      chi2 += (*fit0_)(point_shifted, hit);
                    }
                    
                    vertex_hypothesis.LogLik = chi2;
                    vertex_hypothesis.Vertex = point_shifted;

                    vertex_hypotheses.push_back(vertex_hypothesis);

                  }
                }
              }
            }
            
            std::vector<JVertexHypothesis>::iterator __end2 = vertex_hypotheses.begin() + 35;
            std::partial_sort(vertex_hypotheses.begin(), __end2, vertex_hypotheses.end(), sortLogLik);
              
            for(std::vector<JVertexHypothesis>::iterator hp = vertex_hypotheses.begin(); hp != __end2; ++hp){   

              double simplex_step = 1;
              fit0_->step.resize(4);
              fit0_->step[0] = JFIT::JPoint4D(JVector3D(simplex_step, 0.0, 0.0), 0.0);
              fit0_->step[1] = JFIT::JPoint4D(JVector3D(0.0, simplex_step, 0.0), 0.0);
              fit0_->step[2] = JFIT::JPoint4D(JVector3D(0.0, 0.0, simplex_step), 0.0);
              fit0_->step[3] = JFIT::JPoint4D(JVector3D(0.0, 0.0, 0.0), 20);
              
              chi2 = (*fit0_)(JPoint4D(hp->Vertex), dataL0_selected.begin(), dataL0_selected.end());
              
              const JGEOMETRY3D::JVector3D fitPos0(fit0_->value);
              const JFIT::JPoint4D pt(fitPos0, fit0_->value.getT());
              
              simplex_step = 0.3;
              fit_->step.resize(4);
              fit_->step[0] = JFIT::JPoint4D(JVector3D(simplex_step, 0.0, 0.0), 0.0);
              fit_->step[1] = JFIT::JPoint4D(JVector3D(0.0, simplex_step, 0.0), 0.0);
              fit_->step[2] = JFIT::JPoint4D(JVector3D(0.0, 0.0, simplex_step), 0.0);
              fit_->step[3] = JFIT::JPoint4D(JVector3D(0.0, 0.0, 0.0), 10);
              
              chi2 = (*fit_)(pt, dataL0_selected.begin(), dataL0_selected.end());
              
              const int NDF = dataL0_selected.size() - fit_->step.size();
              
              const JGEOMETRY3D::JVector3D fitPos(fit_->value);
              const JGEOMETRY3D::JTrack3D fitResult(fitPos, JVersor3Z(0,0), fit_->value.getT());
              
              const double energy(0);
              const JFIT::JFit &outFit = JFIT::getFit(JFIT::JHistory(shower->getHistory()).add(JSHOWERPOSITIONFIT), 
                                                      fitResult, getQuality(chi2, NDF), NDF,
                                                      energy,JFIT::JFitStatus_t::OKAY); 
              
              OutFits.push_back(outFit);
              OutFits.rbegin()->setW(13, chi2);
              OutFits.rbegin()->setW(14, NDF);
            }
          } else{
            std::cout<<"Too few hits " << std::endl;
          }
        }
        std::sort(OutFits.begin(), OutFits.end(), qualitySorter);
      }

    }
    
  };
}

#endif