JMuonGandalf.hh

Go to the documentation of this file.

#ifndef JMUONGANDALF_INCLUDE
#define JMUONGANDALF_INCLUDE

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

#include "JDAQ/JDAQEvent.hh"
#include "JDAQ/JDAQTimeslice.hh"
#include "JDAQ/JDAQSummaryslice.hh"// <- not sure why this dependency is needed 

#include "JTrigger/JHit.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JTriggerParameters.hh"// <- not sure why this dependency is needed 

#include "JFit/JLine1Z.hh"
#include "JFit/JLine1ZEstimator.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JEvt.hh"
#include "JFit/JEvtToolkit.hh"
#include "JFit/JFitParameters.hh"
#include "JFit/JFitStatus.hh"
#include "JFit/JMuonGandalfParameters_t.hh"
#include "JFit/JLine3Z.hh"
#include "JFit/JModel.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JLine3ZRegressor.hh"
#include "JFit/JHitW0.hh"

#include "JGeometry3D/JRotation3D.hh"

/**
 * \author mdejong, gmaggi
 */

namespace JFIT
{
  /**
   * class to handle Muon Gandalf reconstruction
   * first, two other angular reconstructionS should be ran : JMuonPrefit and JMuonSimplex 
   */
  class JMuonGandalf
  {
  private:
    typedef JFIT::JRegressor<JFIT::JLine3Z, JFIT::JGandalf>  JRegressor_t;
    
    JLANG::JSharedPointer<const JDETECTOR::JModuleRouter> router_;
    JFIT::JMuonGandalfParameters_t parameters_;
    
    mutable JRegressor_t fitRegresor_;  
    
  public:
    
    /**
     * Default constructor 
     */
    JMuonGandalf()
    {}

    /**
     * Parameterized constructor 
     *
     * \param router JSharedPointer of JModuleRouter, this is built via detector file.
     * \param parameters struct that holds default-optimized parameters for the reconstruction, for ARCA and ORCA configuration. These default parameters could be found in $JPP_DATA.
     * \param pdfFile PDF file descriptor
     */
    JMuonGandalf(const JLANG::JSharedPointer<const JDETECTOR::JModuleRouter> &router,
                 const JFIT::JMuonGandalfParameters_t &parameters,
                 std::string pdfFile):
      router_(router), 
      parameters_(parameters)
    {
      JRegressor_t::debug              = 1;//see what to do with this
      JRegressor_t::T_ns.setRange(-50.0, +450.0);// ns
      JRegressor_t::Vmax_npe           =   10.0;   
      JRegressor_t::MAXIMUM_ITERATIONS = 10000;
      
      fitRegresor_ = {pdfFile, parameters.TTS_ns}; 
      
      fitRegresor_.parameters.resize(5);
      
      fitRegresor_.parameters[0] = JFIT::JLine3Z::pX();
      fitRegresor_.parameters[1] = JFIT::JLine3Z::pY();
      fitRegresor_.parameters[2] = JFIT::JLine3Z::pT();
      fitRegresor_.parameters[3] = JFIT::JLine3Z::pDX();
      fitRegresor_.parameters[4] = JFIT::JLine3Z::pDY();
    }
    
    ~JMuonGandalf(){}

    /**
     * Declaration of Member function that actually performs the reconstruction
     *
     * \param timeSlice which is contructed via a JDAQEvent
     * \param InPreFits input JEvt which is a container of prefits, normally JMuonSimplex
     * \param OutFits non const (output) JEvt.
     */
    void
    getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSlice,
            JFIT::JEvt &InPreFits,
            JFIT::JEvt &OutFits) const;

    /**
     * Compare hits by PMT identifier and time.
     *
     * \param  first         first  hit
     * \param  second        second hit
     * \return               true if first before second; else false
     */ 
    static bool 
    compare(const JTRIGGER::JHitL0& first, 
            const JTRIGGER::JHitL0& second)
    {
      if (std::equal_to<KM3NETDAQ::JDAQPMTIdentifier>()(first, second))
        return std::less<JTRIGGER::JHit>()(first, second);
      else
        return std::less<KM3NETDAQ::JDAQPMTIdentifier>()(first, second);
    }
    
    
  };
}

/**
 * Member function definition
 */
void
JFIT::JMuonGandalf::getJEvt(const KM3NETDAQ::JDAQTimeslice &timeSlice,
                            JFIT::JEvt &InPreFits,
                            JFIT::JEvt &OutFits) const
{
  if ( InPreFits.empty() ) return;
  
  const bool reprocess          = parameters_.reprocess;
  double roadWidth_m            = parameters_.roadWidth_m;
  const double R_Hz             = parameters_.R_Hz; 
  const size_t numberOfPrefits  = parameters_.numberOfPrefits; 
  const double E_GeV            = parameters_.E_GeV; 

  typedef std::vector<JTRIGGER::JHitL0>  JDataL0_t;
  typedef std::vector<JHitW0>  JDataW0_t;

  const JTRIGGER::JBuildL0<JTRIGGER::JHitL0>  buildL0;
  
  if (fitRegresor_.getRmax() < roadWidth_m) {
    roadWidth_m = fitRegresor_.getRmax();
  }

  JEvt::iterator __end = InPreFits.end();

  if (reprocess) {
    __end = std::partition( InPreFits.begin(), __end, 
                            JFIT::JHistory::is_not_event(JFIT::JFitApplication_t::JMUONGANDALF)
                            );
  }

  if (InPreFits.begin() != __end) {

    std::copy(InPreFits.begin(), __end, std::back_inserter(OutFits));

    if (numberOfPrefits > 0) {
      std::advance(__end = InPreFits.begin(), std::min(numberOfPrefits, OutFits.size()));
    }

    std::partial_sort(InPreFits.begin(), __end, InPreFits.end(), qualitySorter);

    __end = std::partition(InPreFits.begin(), __end, 
                           JFIT::JHistory::is_event(InPreFits.begin()->getHistory())
                           );

    JDataL0_t dataL0;

    buildL0(timeSlice, *router_, std::back_inserter(dataL0));

    if (dataL0.size() >= fitRegresor_.parameters.size()) {

      for (JEvt::const_iterator track = InPreFits.begin(); track != __end; ++track) {

        const JGEOMETRY3D::JRotation3D     R (JFIT::getDirection(*track));
        const JFIT::JLine1Z         tz(getPosition (*track).rotate(R), track->getT());
        const JFIT::JModel<JFIT::JLine1Z> match(tz, roadWidth_m, 
                                                JFIT::JRegressor<JFIT::JLine3Z, JFIT::JGandalf>::T_ns);

        //hit selection based on start value

        JDataW0_t data;

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

          JFIT::JHitW0 hit(*i, R_Hz);

          hit.rotate(R);

          if (match(hit)) {
            data.push_back(hit);
          }
        }

        //select first hit

        std::sort(data.begin(), data.end(), compare);

        JDataW0_t::iterator __end = std::unique( data.begin(), data.end(), 
                                                 std::equal_to<KM3NETDAQ::JDAQPMTIdentifier>()
                                                 );

        const int NDF = std::distance(data.begin(), __end) - fitRegresor_.parameters.size();

        if (NDF >= 0) {
          
          //set fit parameters  

          if (track->getE() > 0.1)
            fitRegresor_.E_GeV   = track->getE();
          else
            fitRegresor_.E_GeV   = E_GeV;

            const double chi2 = fitRegresor_(JLine3Z(tz), data.begin(), __end);

          JGEOMETRY3D::JTrack3D tb(fitRegresor_.value);

          tb.rotate_back(R);

          const double energy(0);
          OutFits.push_back(JFIT::getFit(JFIT::JHistory(track->getHistory()).add(JFIT::JFitApplication_t::JMUONGANDALF), tb, 
                                     JFIT::getQuality(chi2), NDF,energy,JFIT::JFitStatus_t::OKAY));

          //set additional weights

          OutFits.rbegin()->setW(JGANDALF_BETA0_RAD,sqrt(fitRegresor_.error.getDX() * fitRegresor_.error.getDX() +
                                                     fitRegresor_.error.getDY() * fitRegresor_.error.getDY()));
          OutFits.rbegin()->setW(JGANDALF_BETA1_RAD,sqrt(fitRegresor_.error.getDX() * fitRegresor_.error.getDY()));
          OutFits.rbegin()->setW(JGANDALF_CHI2,                 chi2);
          OutFits.rbegin()->setW(JGANDALF_NUMBER_OF_HITS,       distance(data.begin(), __end));
          OutFits.rbegin()->setW(JGANDALF_LAMBDA,               fitRegresor_.lambda);
          OutFits.rbegin()->setW(JGANDALF_NUMBER_OF_ITERATIONS, fitRegresor_.numberOfIterations);
        }
      }
    }
    //apply default sorter

    std::sort(OutFits.begin(), OutFits.end(), qualitySorter);
  }
    
}

#endif