JShowerPositionFit.hh

Go to the documentation of this file.

#ifndef JSHOWERPOSITIONFIT_INCLUDE
#define JSHOWERPOSITIONFIT_INCLUDE

#include <string>
#include <iostream>
#include <set>
#include <vector>
#include <algorithm>
#include <memory>
#include <math.h>

#include "km3net-dataformat/online/JDAQTimeslice.hh"
#include "km3net-dataformat/online/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 "JSupport/JSummaryFileRouter.hh"

#include "JFit/JFitToolkit.hh"
#include "JFit/JEstimator.hh"
#include "JFit/JPoint4D.hh"
#include "JFit/JModel.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JShowerBrightPointRegressor.hh"
  
#include "JReconstruction/JHitW0.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JShowerParameters.hh"
  
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JShower3D.hh"
#include "JGeometry3D/JShower3E.hh"
#include "JGeometry3D/JTrack3D.hh"
#include "JGeometry3D/JTrack3E.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"
  
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorSubset.hh"
  
#include "JTools/JPermutation.hh"
#include "JTools/JRange.hh"
#include "JLang/JSharedPointer.hh"

/**
 * \author adomi
 */
namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }
  
namespace JRECONSTRUCTION {
 
  using JDETECTOR::JModuleRouter;
  using JSUPPORT::JSummaryRouter;
  using JFIT::JPoint4D;
  using JFIT::JGandalf;
  using JFIT::JRegressor;

  /**
   * class to handle the second position fit of the shower reconstruction, mainly dedicated for ORCA
   */
  class JShowerPositionFit :
  
   public JShowerPositionFitParameters_t,
   public JRegressor<JPoint4D, JGandalf>

  {
 
    typedef JRegressor<JPoint4D, JGandalf> JRegressor_t;
    using JRegressor_t::operator();
    
  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 summary       summary router
     * \param pdfFile       file name with PDF
     * \param debug         debug
     */ 
    JShowerPositionFit(const JShowerPositionFitParameters_t&  parameters,
                       const JModuleRouter&                   router,
                       const JSummaryRouter&                  summary,
                       const std::string                      pdfFile,
                       const int                              debug = 0):
      JShowerPositionFitParameters_t(parameters),
      JRegressor_t(pdfFile, parameters.TTS_ns),
      router(router),
      summary(summary)
    {
      using namespace JPP;
      
      JRegressor_t::debug  = 1;
      JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
      JRegressor_t::Vmax_npe = 20.0;
      JRegressor_t::MAXIMUM_ITERATIONS = 1000;

      this->parameters.resize(4);
       
      this->parameters[0] = JPoint4D::pX();
      this->parameters[1] = JPoint4D::pY();
      this->parameters[2] = JPoint4D::pZ();
      this->parameters[3] = JPoint4D::pT();
    }
    
    /**
     * Declaration of the member function that actually performs the reconstruction
     *
     * \param event which is a JDAQEvent
     * \param in input fits
     */ 
    JEvt operator()(const KM3NETDAQ::JDAQEvent& event, const JFIT::JEvt &in)
    {
      using namespace std;
      using namespace JPP;
        
      typedef vector<JHitL0> JDataL0_t;
      const JBuildL0<JHitL0> buildL0;
        
      JEvt out;
        
      JDataL0_t dataL0;
      buildL0(JDAQTimeslice(event, true), router, back_inserter(dataL0));
  
      for (JEvt::const_iterator shower = in.begin(); shower != in.end(); ++shower) {

        JPoint4D vx(getPosition(*shower), shower->getT());

        vector<JHitW0> data;

        const JModel<JPoint4D> match(vx, roadWidth_m, JRegressor_t::T_ns);

        for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {

          const double rate_Hz = summary.getRate(*i);
          
          if (match(*i)) {
            data.push_back(JHitW0(*i, rate_Hz));
          }
        }

        // select first hit
        
        sort(data.begin(), data.end(), compare);
        
        vector<JHitW0>::iterator __end = unique(data.begin(), data.end(), equal_to<JDAQPMTIdentifier>());

        const int NDF = distance(data.begin(), __end) - this->parameters.size();

        if (NDF > 0) {
        
          // set fit parameters
          
          if (shower->getE() > 0.1)
            JRegressor_t::E_GeV   = shower->getE();
          else
            JRegressor_t::E_GeV   = this->JShowerPositionFitParameters_t::E_GeV;

          double chi2 = (*this)(vx, data.begin(), __end);
          
          JShower3E sh_fit(this->value.getPosition(), JDirection3D(),
                           this->value.getT(), shower->getE());
  
          out.push_back(getFit(JHistory(shower->getHistory()).add(JSHOWERPOSITIONFIT), sh_fit, 
                               getQuality(chi2), NDF, sh_fit.getE())); 
          
        }    
      }
      
      return out;
    }
    
    /**
     * Auxiliary data structure for sorting of hits.
     */
    static const struct {
      /**
       * Compare hits by PMT identifier and time.
       *
       * \param  first         first  hit
       * \param  second        second hit
       * \return               true if first before second; else false
       */
      template<class T>
      bool operator()(const T& first, const T& second) const
      {
        using namespace std;
        using namespace JPP;
        using namespace KM3NETDAQ;

        if (equal_to<JDAQPMTIdentifier>()(first, second))
          return less<JHit>()(first, second);
        else
          return less<JDAQPMTIdentifier>()(first, second);
      }
    } compare;

    const JModuleRouter&       router;
    const JSummaryRouter&      summary;
  };
}

#endif