JMuonEnergy.hh

Go to the documentation of this file.

#ifndef __JRECONSTRUCTION__JMUONENERGY__
#define __JRECONSTRUCTION__JMUONENERGY__

#include <string>
#include <istream>
#include <ostream>
#include <iomanip>
#include <set>
#include <vector>
#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 "JReconstruction/JHitW0.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JMuonEnergyParameters_t.hh"

#include "JPhysics/JGeane.hh"


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

  using JDETECTOR::JModuleRouter;
  using JSUPPORT::JSummaryRouter;
  using JFIT::JRegressor;
  using JFIT::JEnergy;

  /**
   * Auxiliary class to to determine muon energy.
   */
  class JMuonEnergy : 
    public JMuonEnergyParameters_t,
    public JRegressor<JEnergy>
  {
  public:
    typedef JRegressor<JEnergy>                  JRegressor_t;
    typedef JTRIGGER::JHitR1                     hit_type;
    typedef std::vector<hit_type>                buffer_type;
    typedef JTOOLS::JRange<JEnergy>              JEnergyRange;

    using JRegressor_t::operator();

    /**
     * Constructor.
     *
     * \param  parameters   parameters
     * \param  router       module  router
     * \param  summary      summary router
     * \param  pdfFile      PDF file
     * \param  debug        debug
     */
    JMuonEnergy(const JMuonEnergyParameters_t& parameters,
                const JModuleRouter&           router,
                const JSummaryRouter&          summary,
                const std::string&             pdfFile,
                const int                      debug = 0):
      JMuonEnergyParameters_t(parameters),
      JRegressor_t(pdfFile),
      router(router), 
      summary(summary),
      debug(debug)
    {
      using namespace JPP;
      
      JRegressor_t::debug  = debug;
      JRegressor_t::T_ns.setRange(TMin_ns, TMax_ns);

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

      this->estimator.reset(getMEstimator(mestimator));

      if (!this->estimator.is_valid()) {
        FATAL("Invalid M-Estimator." << endl);
      }
    }
    

    /**                                                                                                 
     * Auxiliary class for energy estimation.                                                
     */
    struct JResult {
      /**                                                                                                      
       * Constructor.                                                                      
       *                                                                                     
       * \param  x               Energy [log(E/GeV)]                                                          
       * \param  chi2            chi2                                                   
       */
      JResult(const JEnergy& x    = 0.0,
              const double   chi2 = std::numeric_limits<double>::max())
      {
        this->x    = x;
        this->chi2 = chi2;
      }

      /**   
       * Type conversion.     
       * 
       * \return                 true if valid chi2; else false                                
       */
      operator bool() const
      {
        return chi2 != std::numeric_limits<double>::max();
      }

      JEnergy x;      //!< Energy                                                          
      double  chi2;   //!< chi2                                                                   
    };


    /**
     * 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;

      JEvt out;

      typedef vector<JHitL0>  JDataL0_t;
      const JBuildL0<JHitL0>  buildL0;

      const JDetector& detector = router.getReference();

      JDataL0_t dataL0;

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

      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);
        
        multiset<JDAQPMTIdentifier> top;

        vector<JNPEHit> data;

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

          JHitR1 hit(*i);
          
          hit.rotate(R);
          
          if (match(hit)) {
            top.insert(i->getPMTIdentifier());
          }
        }
        
        JRange<double> Z_m;
        
        if ( track->hasW(JSTART_LENGTH_METRES) && track->getW(JSTART_LENGTH_METRES) > 0.0 ) {
          Z_m.setLowerLimit(this->ZMin_m);
        }

        const JDetectorSubset_t subdetector(detector, getAxis(*track), roadWidth_m, Z_m);
        
        for (JDetector::const_iterator module = subdetector.begin(); module != subdetector.end(); ++module) {
          
          if (summary.hasSummaryFrame(module->getID())) {
            
            const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID());
            
            for (size_t i = 0; i != module->size(); ++i) {
              
              if (getDAQStatus(frame, *module, i)            &&
                  getPMTStatus(frame, *module, i)            &&
                  !module->getPMT(i).has(PMT_DISABLE)) {

                const JDAQPMTIdentifier id(module->getID(), i);
                
                const double rate_Hz = summary.getRate(id);
                const size_t count   = top.count(id);
                
                data.push_back(JNPEHit(this->getNPE(module->getPMT(i), rate_Hz), count));
              }
            }
          }
        }
        
        const int NDF = distance(data.begin(), data.end()) - 1;

        if( NDF >= 0 ) {
          
          // 5-point search between given limits

          const int      N  =  5;
          
          JResult result[N];
          
          for (int i = 0; i != N; ++i) {
            result[i].x = EMin_log + i * (EMax_log - EMin_log) / (N-1);
          }
          
          map<double, double> buffer;
          
          do {
            
            int j = 0;
            
            for (int i = 0; i != N; ++i) {
              
              if (!result[i]) {
                
                const JEnergy x    = result[i].x;
                const double  chi2 = (*this)(x, data.begin(), data.end());
                
                result[i].chi2 = chi2;
                buffer[chi2]   = x.getE();
              }
              
              if (result[i].chi2 < result[j].chi2) {
                j = i;
              }
            }  
            
            
            for (int i = 0; i != N; ++i) {
              DEBUG(' ' << FIXED(5,2) << result[i].x << ' ' << FIXED(9,3) << result[i].chi2);
            }
            DEBUG(endl);
            
            // squeeze range
            
            switch (j) {
              
            case 0:
              result[4] = result[1];
              result[2] = result[0];
              result[0] = JResult(2*result[2].x - result[4].x);
              break;
              
            case 1:
              result[4] = result[2];
              result[2] = result[1];
              break;
              
            case 2:
              result[0] = result[1];
              result[4] = result[3];
              break;
              
            case 3:
              result[0] = result[2];
              result[2] = result[3];
              break;
              
            case 4:
              result[0] = result[3];
              result[2] = result[4];
              result[4] = JResult(2*result[2].x - result[0].x);
              break;
            }
            
            result[1] = JResult(0.5 * (result[0].x + result[2].x));
            result[3] = JResult(0.5 * (result[2].x + result[4].x));
            
          } while (result[4].x - result[0].x > resolution);
          
          
          if (result[1].chi2 != result[3].chi2) {
            
            result[2].x    += 0.25 * (result[3].x - result[1].x) * (result[1].chi2 - result[3].chi2) / (result[1].chi2 + result[3].chi2 - 2*result[2].chi2);
            result[2].chi2  = (*this)(result[2].x, data.begin(), data.end());
            
          }
          
          const double chi2 = result[2].chi2;
          const double E    = result[2].x.getE();
          
          // calculate additional variables

          double Emin = numeric_limits<double>::max();
          double Emax = numeric_limits<double>::lowest();
          
          for (map<double, double>::const_iterator i = buffer.begin(); i != buffer.end() && i->first <= chi2 + 1.0; ++i) {
            if (i->second < Emin) { Emin = i->second; }
            if (i->second > Emax) { Emax = i->second; }
          }
          
          const double mu_range   = gWater(E);           // range of a muon with the reconstructed energy
          
          double noise_likelihood = 0.0;                 // log-likelihood of every hit being from K40   
          int    number_of_hits   = 0;                   // number of hits selected for JEnergy
          
          for (vector<JNPEHit>::const_iterator i = data.begin(); i != data.end(); ++i) {
            noise_likelihood += log10(getP(i->getY0(), i->getN()));       // probability of getting the observed multiplicity with K40 
            number_of_hits   += i->getN();                                // hit multiplicity
          }
            
          out.push_back(JFit(*track).add(JMUONENERGY));

          out.rbegin()->setE(E);

          out.rbegin()->setW(track->getW());
          out.rbegin()->setW(JENERGY_ENERGY,            E);
          out.rbegin()->setW(JENERGY_CHI2,              chi2);
          out.rbegin()->setW(JENERGY_MUON_RANGE_METRES, mu_range);
          out.rbegin()->setW(JENERGY_NOISE_LIKELIHOOD,  noise_likelihood);
          out.rbegin()->setW(JENERGY_NDF,               NDF);
          out.rbegin()->setW(JENERGY_NUMBER_OF_HITS,    number_of_hits);
          out.rbegin()->setW(JENERGY_MINIMAL_ENERGY,    Emin);
          out.rbegin()->setW(JENERGY_MAXIMAL_ENERGY,    Emax);
        }
      }

      return out;
    }

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

#endif