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 "JDynamics/JCoverage.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 "JReconstruction/JModuleL0.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 KM3NETDAQ::JDAQEvent;
  using KM3NETDAQ::JDAQSummaryFrame;
  using JDETECTOR::JModuleRouter;
  using JDETECTOR::JModule;
  using JSUPPORT::JSummaryRouter;
  using JFIT::JRegressor;
  using JFIT::JEnergy;
  using JPHYSICS::JK40Rates;
  using JDYNAMICS::coverage_type;
 
 
  /**
   * 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 JModuleL0                            module_type;
    typedef std::vector<module_type>             detector_type;
 
    using JRegressor_t::operator();
    
    /**
     * Input data type.
     */
    struct input_type :
      public JDAQEventHeader
    {
      /**
       * Default constructor.
       */
      input_type() :
        JDAQEventHeader()
      {}
 
 
      /**
       * Constructor.
       *
       * \param  header         header
       * \param  in             start values
       * \param  coverage       coverage
       */
      input_type(const JDAQEventHeader& header,
                 const JEvt&            in,
                 const coverage_type&   coverage) :
        JDAQEventHeader(header),
        in(in),
        coverage(coverage)
      {}
 
      JEvt          in;
      detector_type data;
      coverage_type coverage;
    };
 
 
    /**
     * Constructor.
     *
     * \param  parameters    parameters
     * \param  storage       storage
     * \param  rates_Hz      K40 rates [Hz]
     * \param  debug         debug
     */
    JMuonStart(const JMuonStartParameters_t& parameters,
               const storage_type&           storage,
               const JK40Rates&              rates_Hz,
               const int                     debug = 0):
      JMuonStartParameters_t(parameters),
      JRegressor_t(storage),
      rates_Hz(rates_Hz)
    {
      if (this->getRmax() < parameters.roadWidth_m) {
        roadWidth_m = this->getRmax();
      }
 
      JRegressor_t::debug = debug;
    }
 
 
    /**
     * Get input data.
     *
     * \param  router        module router
     * \param  summary       summary data
     * \param  event         event
     * \param  in            start values
     * \param  coverage      coverage
     * \return               input data
     */
    input_type getInput(const JModuleRouter&  router,
                        const JSummaryRouter& summary,
                        const JDAQEvent&      event,
                        const JEvt&           in,
                        const coverage_type&  coverage) const
    {
      using namespace std;
      using namespace JTRIGGER;
 
      input_type   input(event.getDAQEventHeader(), in, coverage);
 
      const JBuildL0 <JHitR0>   buildL0;
      map<int, vector<JHitR0> > data;
 
      const JDAQTimeslice timeslice(event, true);
 
      JSuperFrame2D<JHit> buffer;
 
      for (JDAQTimeslice::const_iterator i = timeslice.begin(); i != timeslice.end(); ++i) {
 
        if (router.hasModule(i->getModuleID())) {
 
          buffer(*i, router.getModule(i->getModuleID()));
 
          buildL0(buffer, back_inserter(data[i->getModuleID()]));
        }
      }
 
      for (const auto& module : router.getReference()) {
        if (!module.empty()) {
          input.data.push_back(module_type(module, summary.getSummaryFrame(module.getID(), R_Hz), data[module.getID()]));
        }
      }
 
      return input;
    }
 
 
    /**
     * Fit function.
     *
     * \param  input         input data
     * \return               fit results
     */
    JEvt operator()(const input_type& input)
    {
      using namespace std;
      using namespace JPP;
 
      JEvent event(JMUONSTART);
 
      JEvt out;
 
      // select start values
 
      JEvt in  = input.in;
 
      sort(in.begin(), in.end(), qualitySorter);
 
      if (numberOfPrefits != 0 && numberOfPrefits < in.size()) {
 
        JEvt::iterator __end = next(in.begin(), numberOfPrefits);
 
        if (numberOfPostfits > 0) {
 
          __end = gridify(__end, in.end(), numberOfPostfits);
        }
 
        in.erase(__end, in.end());
      }
 
      if (!in.empty()) {
        in.select(JHistory::is_event(in.begin()->getHistory()));
      }
 
      const JStart start(Pmin1, Pmin2, Nmax2);  
      
      struct JHit_t {
          
        double getZ()  const { return z; }
        double getP()  const { return p; }
          
        double z;
        double p;
        double y1;
        double yx;
      };
        
      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());
        
        vector<JHit_t> data;
 
        for (const auto& module : input.data) {
 
          JPosition3D pos(module->getPosition());
            
          pos.transform(R, tz.getPosition());
            
          if (pos.getX() <= roadWidth_m) {
 
            const double z = pos.getZ()  -  pos.getX() / getTanThetaC();
            const double t = tz .getT()  + (pos.getZ() + pos.getX() * getKappaC()) * getInverseSpeedOfLight();
 
            const double p = module.getProbability(rates_Hz.getMultiplesRates(), T_ns + t);
 
            double y1 = 0.0;
            double yx = 0.0;
 
            for (size_t i = 0; i != module->size(); ++i) {
 
              if (module.getStatus(i)) {
 
                JPMT pmt = module->getPMT(i);
                    
                pmt.transform(R, tz.getPosition());
                    
                const double npe = this->getY1(pmt);
 
                y1 += npe;
 
                if (count_if(module.begin(), module.end(), make_predicate(&JHitR0::getPMT, (JDAQHit::JPMT_t) i)) == 0) {
                  yx += npe;
                }
              }
            }
 
            data.push_back({ z, p, y1, yx });
          }
        }
 
        double Zmin = 0.0;     // start position
        double Zmax = 0.0;     // end   position
        double y1   = 0.0;     // npe due to MIP between start and end position
        double yx   = 0.0;     // npe due to MIP between start and end position that are not detected
        double pb   = 0.0;     // background probability
        
        if (!data.empty()) {
          
          sort(data.begin(), data.end(), make_comparator(&JHit_t::getZ));
 
          vector<JHit_t>::const_iterator         q1 = start.find(data. begin(), data. end());
          vector<JHit_t>::const_reverse_iterator q2 = start.find(data.rbegin(), data.rend());
 
          if (q1 != data.end() && q2 != data.rend()) {
 
            vector<JHit_t>::const_iterator         p1 = q1;  if (p1 != data. begin()) { --p1; }
            vector<JHit_t>::const_reverse_iterator p2 = q2;  if (p2 != data.rbegin()) { --p2; }
 
            Zmin = 0.5 * (p1->getZ() + q1->getZ());
            Zmax = 0.5 * (p2->getZ() + q2->getZ());
 
            for (vector<JHit_t>::const_iterator i = q1; i != q2.base(); ++i) {
              y1 += i->y1;
              yx += i->yx;
              pb += log(i->p > Pmin ? i->p : Pmin);
            }
          }
        }
 
        JFit fit = *track;
          
        fit.push_back(event());
 
        // move track
 
        fit.move(Zmin, getSpeedOfLight());
          
        out.push_back(fit);
 
        // set additional values
 
        out.rbegin()->setW(track->getW());
        out.rbegin()->setW(JSTART_NPE_MIP_TOTAL,     y1);
        out.rbegin()->setW(JSTART_NPE_MIP_MISSED,    yx);
        out.rbegin()->setW(JSTART_LENGTH_METRES,     Zmax - Zmin);
        out.rbegin()->setW(JSTART_BACKGROUND_LOGP,   pb);
        out.rbegin()->setW(JPP_COVERAGE_ORIENTATION, input.coverage.orientation);
        out.rbegin()->setW(JPP_COVERAGE_POSITION,    input.coverage.position);
      }
 
      // apply default sorter
 
      sort(out.begin(), out.end(), qualitySorter);
 
      copy(input.in.begin(), input.in.end(), back_inserter(out));
      
      return out;
    }
 
  private:
    const JK40Rates  rates_Hz;
  };
}
 
#endif