JMuonStart.hh

Go to the documentation of this file.

#ifndef __JRECONSTRUCTION__JMUONSTART__
#define __JRECONSTRUCTION__JMUONSTART__

#include <string>
#include <istream>
#include <ostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
#include <cmath>

#include "km3net-dataformat/online/JDAQEvent.hh"
#include "km3net-dataformat/online/JDAQSummaryslice.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorSubset.hh"
#include "JDetector/JModuleRouter.hh"

#include "JTrigger/JHitL0.hh"
#include "JTrigger/JHitR1.hh"
#include "JTrigger/JBuildL0.hh"

#include "JSupport/JSummaryRouter.hh"

#include "JFit/JLine1Z.hh"
#include "JFit/JEnergyRegressor.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JModel.hh"
#include "JFit/JNPEHit.hh"
#include "JFit/JEnergy.hh"
#include "JFit/JMEstimator.hh"

#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JMuonStartParameters_t.hh"
#include "JReconstruction/JStart.hh"

#include "JPhysics/JConstants.hh"
#include "JPhysics/JK40Rates.hh"

#include "JLang/JComparator.hh"
#include "JLang/JPredicate.hh"


/**
 * \author mdejong, azegarelli, scelli
 */
namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }
  
namespace JRECONSTRUCTION {

  using JDETECTOR::JModuleRouter;
  using JDETECTOR::JModule;
  using JSUPPORT::JSummaryRouter;
  using JFIT::JRegressor;
  using JFIT::JEnergy;
  using JPHYSICS::JK40Rates;
  using JPHYSICS::JRateL1_t;
  using KM3NETDAQ::JDAQSummaryFrame;

  /**
   * Get probability of given response in optical module due to random background.
   *
   * \param  module               module
   * \param  frame                summary frame
   * \param  R_Hz                 multiples rates [Hz]
   * \param  T_ns                 time window [ns]
   * \param  top                  list of indentifiers of PMTs with hit in given module
   * \return                      probability
   */
  inline double getProbability(const JModule&            module,
                               const JDAQSummaryFrame&   frame,
                               const JRateL1_t&          R_Hz,
                               const double              T_ns,
                               const std::multiset<int>& top)
  {
    using namespace std;
    using namespace JPP;
    using namespace KM3NETDAQ;

    size_t N = 0;         // number of active PMTs
    size_t M = 0;         // multiplicity
    double R = 0.0;       // total rate

    for (size_t i = 0; i != module.size(); ++i) {

      if (getDAQStatus(frame, module, i)            &&
          getPMTStatus(frame, module, i)            &&
          frame[i].is_valid()                       &&
          !module.getPMT(i).has(PMT_DISABLE)) {

        N += 1;
        M += top.count(i);
        R += frame.getRate(i);
      }
    }

    if (N > 0)
      return JFIT::getProbability(N, M, JK40Rates(R/N, R_Hz), T_ns);
    else
      return (M == 0 ? 1.0 : 0.0);
  }


  /**
   * Auxiliary class to determine start and end position of muon trajectory.
   */
  class JMuonStart : 
    public JMuonStartParameters_t,
    public JRegressor<JEnergy>
  {
  public:
    typedef JRegressor<JEnergy>                  JRegressor_t;
    typedef JTRIGGER::JHitL0                     hit_type;
    typedef std::vector<hit_type>                buffer_type;

    using JRegressor_t::operator();
    
    /**
     * Constructor.
     *
     * \param  parameters   parameters
     * \param  router       module  router
     * \param  summary      summary router
     * \param  pdfFile      PDF file
     * \param  rates_Hz     K40 rates [Hz]
     * \param  debug        debug
     */
    JMuonStart(const JMuonStartParameters_t& parameters,
               const JModuleRouter&          router,
               const JSummaryRouter&         summary,
               const std::string&            pdfFile,
               const JK40Rates&              rates_Hz,
               const int                     debug = 0):
      JMuonStartParameters_t(parameters),
      JRegressor_t(pdfFile),
      router(router),
      summary(summary),
      rates_Hz(rates_Hz),
      debug(debug)
    {
      JRegressor_t::debug = debug;
      JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);

      if (this->getRmax() < parameters.roadWidth_m) {
        roadWidth_m = this->getRmax();
      }
    }


    /**
     * Fit function.
     *
     * \param  event         event
     * \param  in            start values
     * \return               fit results
     */
    JEvt operator()(const KM3NETDAQ::JDAQEvent& event, const JEvt& in)
    {
      using namespace std;
      using namespace JPP;

      const JBuildL0<hit_type>  buildL0;

      buffer_type dataL0;

      buildL0(event, router, true, back_inserter(dataL0));

      return (*this)(dataL0, in);
    } 


    /**
     * Fit function.
     *
     * \param  data          dataL0
     * \param  in            start values
     * \return               fit results
     */
    JEvt operator()(const buffer_type& data, const JEvt& in)
    {
      using namespace std;
      using namespace JPP;

      JEvt out;

      const JDetector& detector = router.getReference();

      const JStart start(Pmin1, Pmin2, Nmax2);  
      
      for (JEvt::const_iterator track = in.begin(); track != in.end(); ++track) {
        
        const JRotation3D           R (getDirection(*track));
        const JLine1Z               tz(getPosition (*track).rotate(R), track->getT());
        const JFIT::JModel<JLine1Z> match(tz, roadWidth_m, JRegressor_t::T_ns);
        
        map<int, multiset<int> > top;
        
        for (buffer_type::const_iterator i = data.begin(); i != data.end(); ++i) {
          
          hit_type hit(*i);
          
          hit.rotate(R);
          
          if (match(hit)) {
            top[hit.getModuleID()].insert(hit.getPMTAddress());
          }
        }
        
        struct JHit_t {
          
          double getZ() const { return z; }
          double getP() const { return p; }
          
          double z;
          double p;
        };
        
        vector<JHit_t> data;
        
        for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
          
          if (summary.hasSummaryFrame(module->getID())) {

            const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID());
            
            JPosition3D pos(module->getPosition());
            
            pos.transform(R, tz.getPosition());
            
            if (pos.getX() <= roadWidth_m) {
              
              const double z = pos.getZ()  -  pos.getX() / getTanThetaC();
              
              const double p = getProbability(*module, frame, rates_Hz.getMultiplesRates(), JRegressor_t::T_ns.getLength(), top[module->getID()]);
              
              data.push_back({ z, p });
            }
          }
        }
        
        double Zmin       = 0.0;     // relative start position [m]
        double Zmax       = 0.0;     // relative end   position [m]
        double npe_total  = 0.0;     // npe due to mip
        double npe_missed = 0.0;     // npe due to mip not detected
        
        if (!data.empty()) {
          
          sort(data.begin(), data.end(), JLANG::make_comparator(&JHit_t::getZ));
          
          vector<JHit_t>::const_iterator         track_start = start.find(data. begin(), data. end());
          vector<JHit_t>::const_reverse_iterator track_end   = start.find(data.rbegin(), data.rend());
          
          if (track_start != data. end()) { Zmin = track_start->getZ(); }
          if (track_end   != data.rend()) { Zmax = track_end  ->getZ(); }
          
          // additional features
          
          for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) {
            
            JPosition3D pos(module->getPosition());
            
            pos.transform(R, tz.getPosition());
            
            if (pos.getX() <= roadWidth_m) {
              
              const double z = pos.getZ()  -  pos.getX() / getTanThetaC();
              
              if (z >= Zmin && z <= Zmax) {
                
                for (size_t i = 0; i != module->size(); ++i) {
                  
                  const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID());
                  
                  if (getDAQStatus(frame, *module, i)      &&
                      getPMTStatus(frame, *module, i)      &&
                      frame[i].is_valid()                  &&
                      !module->getPMT(i).has(PMT_DISABLE)) {
                    
                    JPMT pmt = module->getPMT(i);
                    
                    pmt.transform(R, tz.getPosition());
                    
                    const double ya = this->getNPE(pmt, 0.0).getYA();
                    
                    npe_total += ya;

                    if (top[module->getID()].count(i) == 0) {
                      npe_missed += ya;
                    }
                  }
                }
              }
            }
          }
        }
 
        JFit fit = *track;

        // move track
          
        fit.move(Zmin, getSpeedOfLight());
          
        out.push_back(fit.add(JMUONSTART));
          
        out.rbegin()->setW(track->getW());
        out.rbegin()->setW(JSTART_NPE_MIP_TOTAL,  npe_total);
        out.rbegin()->setW(JSTART_NPE_MIP_MISSED, npe_missed);
        out.rbegin()->setW(JSTART_LENGTH_METRES,  Zmax - Zmin);
      }
      
      return out;
    }

  private:
    const JModuleRouter&       router;
    const JSummaryRouter&      summary;
    const JK40Rates            rates_Hz;
    int                        debug;  
  };
}

#endif