JORCAShowerFit.hh

Go to the documentation of this file.

#ifndef JORCASHOWERFIT_INCLUDE
#define JORCASHOWERFIT_INCLUDE
 
/*TODO-NOTES :
 - This algorithm erases element from STL vectors, that is slow!, try to use another STL container as a list. 
 - Prefer to accesses to container elements via at() instead of [], the latter doesn't check the range.  
 - This algorithm uses several magic numbers, that should be avoided and those should clearly be defined.
 - An exact value of pi should be used. 
 - It would be healthy to check whether some numbers are nan, e.g. , when using getDot.
*/
 
#include <string>
#include <iostream>
#include <set>
#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/JShower3EZRegressor.hh"
#include "JFit/JHitW0.hh"
#include "JFit/JEvt.hh"
#include "JFit/JEvtToolkit.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JPoint4D.hh"
#include "JFit/JModel.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JPoint4DRegressor.hh"
#include "JFit/JShowerFitParameters_t.hh"
 
#include "JDAQ/JDAQTimeslice.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorSubset.hh"
 
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"
 
#include "JLang/JSharedPointer.hh"
 
 
namespace JFIT
{
  /**
   * class to handle the third step of the shower reconstruction, mainly dedicated for ORCA
   */
 
  using namespace JPP;
 
  class JORCAShowerFit
  {
  private:
 
    typedef JFIT::JRegressor<JFIT::JShower3EZ, JFIT::JSimplex>  JRegressor_t;
    mutable JRegressor_t fit_;
 
    JLANG::JSharedPointer<const JDETECTOR::JModuleRouter> moduleRouter_;
    JFIT::JShowerFitParameters_t parameters_; 
 
  public:
 
    /**
     * Default constructor
     */
    JORCAShowerFit()
    {}
 
    /**
     * 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 PDF file
     */
 
    JORCAShowerFit(const JLANG::JSharedPointer<const JDETECTOR::JModuleRouter> &moduleRouter,
                   const JFIT::JShowerFitParameters_t &parameters,
                   const std::string pdfFile):
      moduleRouter_(moduleRouter),
      parameters_(parameters)
    {
      JRegressor_t::debug  =1;
      JRegressor_t::T_ns.setRange(-parameters.Tmax_ns,parameters.Tmax_ns);
      JRegressor_t::Vmax_npe = 20.0;
      JRegressor_t::MAXIMUM_ITERATIONS = 1000;
      
      fit_ = {pdfFile};
      fit_.estimator.reset(new JMEstimatorNull());
    }
 
    ~JORCAShowerFit(){}
    
 
    /**
     * Container of shower hypothesis: 
     * During the fit the algorithm creates a grid in direction and energy space,
     * for each of these points a shower hypothesis is saved: direction, energy and LogLikelihood value.
     */
    
    struct JShowerHypothesis
    {
      JVector3D  Direction;
      double     Energy;
      double     LogLik;
    };
 
    /**
     * Function to estimate the initial em shower energy 
     */
 
    static double estimateInitialEnergy(int nHits){
 
      double Energy = exp((nHits + 26)/26);
      if(Energy > 100) Energy = 100;
 
      return Energy;
    }
 
 
    /**
     * Function to sort different shower hypotheses by their likelihood.
     */
 
    static bool sortLogLik(JShowerHypothesis A, JShowerHypothesis B)
    {  
      return(A.LogLik < B.LogLik);
    };    
 
    /**
     * Declaration of the member function that actually performs the reconstruction
     *
     * \param timeSliceBuildL0 which is contructed via a JDAQEvent
     * \param InPreFits input fits, normally produced by JORCAShowerPositionFit
     * \param OutFits output fits
     */ 
 
    void getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL0,
                 JFIT::JEvt &InPreFits,
                 JFIT::JEvt &OutFits) const;
       
  };
  
}
 
 
/**
 * member function definition
 */
 
void 
JFIT::JORCAShowerFit::getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSliceBuildL0,
                              JFIT::JEvt &InPreFits,
                              JFIT::JEvt &OutFits) const
{
 
  const double Dmax_m               = parameters_.Dmax_m;
  const double R_Hz                 = parameters_.R_Hz;
  const size_t numberOfPrefits = parameters_.numberOfPrefits;
 
  if ( InPreFits.empty() ) return;
 
  typedef vector<JTRIGGER::JHitL0>  JDataL0_t;
  typedef vector<JFIT::JPMTW0> JPMTW0_t;
  JBuildL0<JTRIGGER::JHitL0>  buildL0;
 
  const JDETECTOR::JDetector &detector = moduleRouter_->getReference();
 
  JEvt::iterator __end = InPreFits.end();
  std::copy(InPreFits.begin(), __end, back_inserter(OutFits));
  
  if (numberOfPrefits > 0) {
    advance(__end = InPreFits.begin(), min(numberOfPrefits, InPreFits.size()));
  }
  
  partial_sort(InPreFits.begin(), __end, InPreFits.end(), qualitySorter);
  
  JDataL0_t dataL0;
  buildL0(timeSliceBuildL0, *moduleRouter_, back_inserter(dataL0));
  
  for (JEvt::const_iterator shower = InPreFits.begin(); shower != __end; ++shower) {
    const JVector3D pos(getPosition(*shower));   
        
    JDetectorSubset_t subdetector(detector, pos, Dmax_m);
 
    JPoint4D pt(pos, shower->getT());
    JVector3D start_dir(0,0,0);
 
    // get the detector center and find the distance between the reco vertex and it.
    const JCenter3D center(detector.begin(), detector.end());
    const double radius = pos.getDistance(center);
 
    multiset<KM3NETDAQ::JDAQPMTIdentifier> top;
    
    double Nhits = 0;
    
    for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
      
      const JHitW0 hit(*i, R_Hz);
      const JPosition3D hit_pos(hit.getPosition());
 
      const double D = hit_pos.getDistance(pt.getPosition());
      const double t_res = hit.getT() - pt.getT(hit_pos);
 
      const JDirection3D photonDir(hit_pos - pt.getPosition());
      const double cosT = photonDir.getDot(hit.getDirection());
      
      if(D < Dmax_m && (t_res >= -25 && t_res <= 25) && (cosT >= -1 && cosT <= 0.1)){
        top.insert(hit.getPMTIdentifier());
        const JVector3D d(photonDir.getDX(), photonDir.getDY(), photonDir.getDZ());
        start_dir.add(d);
      }
      if(D < Dmax_m && (t_res > -30 && t_res < 20)){
        Nhits++;
      }
    }
    
    double start_E = estimateInitialEnergy(Nhits);
 
    // rotate subdetector wrt start_dir, subdetector is already subtracted by pos
 
    const JDirection3D conversion(start_dir.getX(), start_dir.getY(), start_dir.getZ());
    const JRotation3D R(conversion);
 
    vector<JPMTW0> buffer;
 
    JFIT::JORCAShowerFit::JShowerHypothesis shower_hypothesis;
    vector<JFIT::JORCAShowerFit::JShowerHypothesis> shower_hypotheses;
    
    double chi2 = 0.0;
    double max_theta = 0, max_phi = 0, scan_step = 0;
    double Nmin = 0;
    
    if(start_E > 10 && radius < 100){
      max_theta = 30;
      max_phi = 30;
      scan_step = 5;
      Nmin = 10;
    } else {
      max_theta = 35;
      max_phi = 35;
      scan_step = 5;
      Nmin = 20;
    }
 
    const JAngle3D start_dir_angles(start_dir.getX(), start_dir.getY(), start_dir.getZ());
    
    for(double E = -15; E <= 5; E += 7.5){  // scan in energy
      for(double th = -max_theta; th <= max_theta; th += scan_step){
        for(double ph = -max_phi; ph <= max_phi; ph += scan_step){
          
          buffer.clear();
          chi2 = 0.0;
          
          //todo: use precise value of pi
          const double theta  = th * PI / 180;
          const double phi    = ph * PI / 180;
          const double scan_E = start_E + E; // [GeV] 
 
          const JAngle3D dir_shifted_angles(start_dir_angles.getTheta() + theta, 
                                                         start_dir_angles.getPhi() + phi);
          const JVector3D dir_shifted(dir_shifted_angles.getDX(), 
                                      dir_shifted_angles.getDY(), 
                                      dir_shifted_angles.getDZ());
          const JDirection3D conversion(dir_shifted.getX(), 
                                                     dir_shifted.getY(), 
                                                     dir_shifted.getZ());
          const JRotation3D R(conversion);
          pt.rotate(R);
          
          for (JDetectorSubset_t::iterator module = subdetector.begin();
               module != subdetector.end(); ++module) {
      
            module->rotate(R);
            JPosition3D module_pos(module->getPosition());
      
            for (unsigned int i = 0; i != module->size(); ++i) {
              const JDAQPMTIdentifier id(module->getID(), i);
              JPMT pmt(module->getPMT(i));
              buffer.push_back(JPMTW0(pmt, R_Hz, top.count(id)));
            }
          }
          
          for(JPMTW0_t::iterator pmt = buffer.begin(); pmt != buffer.end(); ++pmt ){
            chi2 += fit_(JShower3EZ(pt, JVersor3Z(), JEnergy(log10(scan_E))), *pmt);
          }
          
          shower_hypothesis.LogLik = chi2;
          shower_hypothesis.Energy = scan_E;
          shower_hypothesis.Direction = dir_shifted;
          
          shower_hypotheses.push_back(shower_hypothesis);
 
          pt.rotate_back(R);
          
        }
      }
    }
 
    vector<JShowerHypothesis>::iterator __end2 = shower_hypotheses.begin() + Nmin;
    partial_sort(shower_hypotheses.begin(), __end2, shower_hypotheses.end(), sortLogLik);
 
    for(vector<JShowerHypothesis>::iterator hp = shower_hypotheses.begin(); hp != __end2; ++hp){
 
      fit_.step.resize(3);
      const double dir_step = 0.05;
      fit_.step[0] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(dir_step, 0.0), JEnergy());
      fit_.step[1] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(0.0, dir_step), JEnergy());
      fit_.step[2] = JShower3EZ(JPoint4D(JVector3D(), 0.0), JVersor3Z(), 0.05);
      
      const double scan_E = hp->Energy;
 
      buffer.clear();
      
      const JDirection3D startVersor(hp->Direction);
      
      const JDirection3D conversion(startVersor.getDX(), startVersor.getDY(), startVersor.getDZ());
      JRotation3D R(conversion);
      pt.rotate(R);
 
      for(JDetector::const_iterator module = detector.begin(); 
          module != detector.end(); ++module) {
        JPosition3D pos(module->getPosition());
        pos.rotate(R);
 
        if (pt.getDistance(pos) < Dmax_m) {
          for (unsigned int i = 0; i != module->size(); ++i) {
            const JDAQPMTIdentifier id(module->getID(), i);
            JPMT pmt(module->getPMT(i));
            pmt.rotate(R);
            buffer.push_back(JPMTW0(pmt, R_Hz, top.count(id)));
          }
        }
      }
      
      const int NDF = buffer.size() - fit_.step.size();
      
      chi2 = fit_(JShower3EZ(pt, JVersor3Z(), JEnergy(log10(scan_E))), buffer.begin(), buffer.end());
      
      pt.rotate_back(R);
      
      const JVector3D resPos(fit_.value);
      const JVersor3D resDir(fit_.value.getDX(), fit_.value.getDY(), fit_.value.getDZ());
      
      const JTrack3D td(resPos, resDir, shower->getT());
      double E_reco_corrected = abs(fit_.value.getE() - 4);
 
      // Now that JShowerFit uses a Param PDF, the energy has to be corrected with a linear function obtained
      // from a fit to the reconstruction error: this can be improved.
 
      const double p0 = -2.028;
      const double p1 = 0.7683;
      if(E_reco_corrected > 6){
        E_reco_corrected = (E_reco_corrected - p0) / (1 + p1);
      }
 
      JTrack3E tb(td, E_reco_corrected);
      tb.rotate_back(R);
      
      const JFit &outFit = getFit(JHistory(shower->getHistory()).add(JSHOWERCOMPLETEFIT), 
                                  tb, -chi2, NDF, E_reco_corrected,
                                  JFitStatus_t::OKAY); 
 
      OutFits.push_back(outFit);
    }
  }
  sort(OutFits.begin(), OutFits.end(), qualitySorter);
}
 
#endif