JShowerFeatures.hh

Go to the documentation of this file.

#ifndef __JRECONSTRUCTION__JSHOWERFEATURES__
#define __JRECONSTRUCTION__JSHOWERFEATURES__
 
#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
 
#include <map>
 
#include "km3net-dataformat/online/JDAQEvent.hh"
#include "km3net-dataformat/online/JDAQSummaryslice.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"
 
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"
 
#include "JFit/JFitToolkit.hh"
#include "JFit/JModel.hh"
#include "JFit/JShower3EZRegressor.hh"
 
#include "JReconstruction/JHitW0.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JShowerFitParameters_t.hh"
#include "JReconstruction/JModuleL0.hh"
 
 
#include "JPhysics/JPDFTable.hh"
#include "JPhysics/JPDFToolkit.hh"
#include "JPhysics/JGeane.hh"
 
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JLocation.hh"
#include "JDetector/JPMTParametersMap.hh"
 
#include "JDynamics/JCoverage.hh"
 
#include "JSupport/JSummaryRouter.hh"
 
#include "JGeometry3D/JRotation3D.hh"
 
#include "JTools/JRange.hh"
 
#include "Jeep/JMessage.hh"
 
#include "JLang/JVectorize.hh"
 
 
/**
 * \author vcarretero
 */
 
namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }
 
namespace JRECONSTRUCTION {
 
  using KM3NETDAQ::JDAQEvent;
  using JDETECTOR::JModuleRouter;
  using JDETECTOR::JLocation;
  using JDETECTOR::JPMTIdentifier;
  using JDETECTOR::JPMTParametersMap;
  using JSUPPORT::JSummaryRouter;
  using JFIT::JRegressor;
  using JFIT::JEnergy;
  using JFIT::JShower3EZ;
  using JFIT::JAbstractMinimiser;
  using JDYNAMICS::coverage_type;
 
  /**
   * Wrapper class to add features after the final fit of shower
   *
   * Note that no values computed at previous steps are modified.
   */
  struct JShowerFeatures :
    public JShowerFitParameters_t,
    public JRegressor<JShower3EZ, JAbstractMinimiser>
  {
    typedef JRegressor<JShower3EZ, JAbstractMinimiser>        JRegressor_t;
 
    typedef JModuleL0                         module_type;
    typedef std::vector<module_type>          detector_type;
    typedef JRegressorStorage<JShower3EZ, JFIT::JAbstractMinimiser> storage_type;
    /**
     * Input data type.
     */
    struct input_type :
      public JDAQEventHeader
    {
      /**
       * Default constructor.
       */
      input_type()
      {}
 
 
      /**
       * 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  pmtParameters PMT parameters
     * \param  debug         debug
     */
    JShowerFeatures(const JShowerFitParameters_t& parameters,
                    const storage_type&           storage,
                    const JPMTParametersMap&      pmtParameters,
                    const int                     debug = 0) :
      JShowerFitParameters_t(parameters),
      JRegressor_t(storage),
      pmtParameters(pmtParameters)
    {
      using namespace JFIT;
 
      JRegressor_t::debug              = debug;
      JRegressor_t::T_ns.setRange(TMin_ns, TMax_ns);
    }
 
    /**
     * 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 JFIT;
      using namespace JGEOMETRY3D;
 
      JEvent event(JSHOWERFEATURES);
 
      JEvt out;
      
      JEvt in = input.in;
      
      in.select(numberOfPrefits, qualitySorter);
 
      if (!in.empty()) {
        in.select(JHistory::is_event(in.begin()->getHistory()));
      }
 
      // Computation of max pmt getting a hit in a DOM for a given event, snapshot                                                   
 
      struct {
        size_t nhits =  0;
        int    id    = -1;
      } snapshot;
 
      for (const auto& module : input.data) {
        if (module.size() != 0) {
          
          size_t snapshot_hits = 0;
 
          for (size_t i = 0; i != module->size(); ++i) {
            
            if (module.getStatus(i)) {
 
              const size_t nhits = std::count_if(
                                                 module.begin(),
                                                 module.end(),
                                                 [i](const JHitR0& hit) {
                                                   return hit.getPMT() == i;
                                                 }
                                                 );
 
              if (nhits > 0) snapshot_hits++;
            }
          }
          
          if (snapshot_hits > snapshot.nhits){
            snapshot.nhits  = snapshot_hits;
            snapshot.id     = module->getID();
          }
        }
      }
      
      for (JEvt::const_iterator shower = in.begin(); shower != in.end(); ++shower) {
        
        //Computation of the minimal distance from jshower vertex to a module
        struct {
          double Dmin = numeric_limits<double>::max();
          int    id   = -1;
        } vertex;
 
        for (const auto& module : input.data) {
          if (module->getFloor() != 0) {
 
            const double D = getDistance(module->getPosition(), getPosition(*shower));
    
            if (D < vertex.Dmin) {
              
              vertex.Dmin = D;
              vertex.id   = module->getID();
              
            }
          }
        }
      
        out.push_back(JFit(*shower).add(JSHOWERFEATURES));
        out.rbegin()->setW(JSHOWERFEATURES_DMIN,               vertex.Dmin);
        out.rbegin()->setW(JSHOWERFEATURES_DMIN_MODULEID,      vertex.id);
        out.rbegin()->setW(JSHOWERFEATURES_NSNAPSHOT,          snapshot.nhits);
        out.rbegin()->setW(JSHOWERFEATURES_NSNAPSHOT_MODULEID, snapshot.id);
 
      }
 
      // apply default sorter
      
      sort(out.begin(), out.end(), qualitySorter);
 
      copy(input.in.begin(), input.in.end(), back_inserter(out));
      
      return out;
    }
    
    const JPMTParametersMap& pmtParameters;
 
  };
 
}
 
#endif