JShowerBjorkenY.hh

Go to the documentation of this file.

 
#ifndef JSHOWERBJORKENY_INCLUDE
#define JSHOWERBJORKENY_INCLUDE
 
#include <string>
#include <set>
#include <vector>
#include <algorithm>
#include <memory>
 
#include "km3net-dataformat/online/JDAQTimeslice.hh"
#include "km3net-dataformat/online/JDAQSummaryslice.hh"
 
#include "JTrigger/JHit.hh"
#include "JTrigger/JTimeslice.hh"
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JHitL1.hh"
#include "JTrigger/JHitR1.hh"
#include "JTrigger/JBuildL0.hh"
#include "JTrigger/JBuildL2.hh"
#include "JTrigger/JAlgorithm.hh"
#include "JTrigger/JMatch3G.hh"
#include "JTrigger/JBind2nd.hh"
 
#include "JSupport/JSummaryRouter.hh"
 
#include "JFit/JShowerBjorkenYRegressor.hh"
#include "JFit/JFitToolkit.hh"
#include "JFit/JPoint4D.hh"
#include "JFit/JModel.hh"
#include "JFit/JSimplex.hh"
#include "JFit/JShowerEH.hh"
 
#include "JAAnet/JAAnetToolkit.hh"
 
#include "JReconstruction/JHitW0.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JEnergyCorrection.hh"
#include "JReconstruction/JShowerBjorkenYParameters_t.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorSubset.hh"
 
#include "JGeometry3D/JVector3D.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"
#include "JGeometry3D/JOmega3D.hh"
#include "JGeometry3D/JTrack3D.hh"
#include "JGeometry3D/JTrack3E.hh"
#include "JGeometry3D/JShower3EY.hh"
 
 
/**
 * \author adomi
 */
namespace JRECONSTRUCTION {}
namespace JPP { using namespace JRECONSTRUCTION; }
  
namespace JRECONSTRUCTION {
 
  using JDETECTOR::JModuleRouter;
  using JSUPPORT::JSummaryRouter;
  using JFIT::JRegressor;
  using JFIT::JEnergy;
  using JFIT::JShowerEH;
  using JFIT::JSimplex;
 
  /**
   * class to handle the direction fit of the shower reconstruction, mainly dedicated for ORCA
   */
  class JShowerBjorkenY : 
    public JShowerBjorkenYParameters_t,
    public JRegressor<JShowerEH, JSimplex>
  {
 
    typedef JRegressor<JShowerEH, JSimplex>  JRegressor_t;
    using JRegressor_t::operator();
    
  public:
 
    /**
     * Parameterized constructor
     *
     * \param parameters   struct that holds default-optimized parameters for the reconstruction, available in $JPP_DATA. 
     * \param router       module router, this is built via detector file. 
     * \param summary      summary router 
     * \param pdfFile      PDF file
     * \param correct      energy correction
     * \param debug        debug
     */
    JShowerBjorkenY(const JShowerBjorkenYParameters_t&   parameters,
                    const JModuleRouter&                 router,
                    const JSummaryRouter&                summary,
                    const std::string                    pdfFile,
                    const JEnergyCorrection&             correct,
                    const int                            debug = 0):
      JShowerBjorkenYParameters_t(parameters),
      JRegressor_t(pdfFile),
      router(router),
      summary(summary),
      correct(correct)
    {
      using namespace JPP;
      
      JRegressor_t::debug  = debug;
      JRegressor_t::T_ns.setRange(parameters.TMin_ns, parameters.TMax_ns);
      JRegressor_t::Vmax_npe = VMax_npe;
      JRegressor_t::MAXIMUM_ITERATIONS = NMax;
      
      this->estimator.reset(getMEstimator(parameters.mestimator));
    }
 
    /**
     * Declaration of the member function that actually performs the reconstruction
     *
     * \param event = JDAQEvent
     * \param in = input fits
     */ 
    JEvt operator()(const KM3NETDAQ::JDAQEvent& event, const JFIT::JEvt &in) 
    {
      using namespace std;
      using namespace JPP;
 
      typedef vector<JHitL0>  JDataL0_t;
      JBuildL0<JHitL0>        buildL0;
 
      JEvt out;
 
      const JDetector &detector = router.getReference();
 
      JDataL0_t dataL0;
 
      buildL0(JDAQTimeslice(event, true), router, back_inserter(dataL0));
  
      for (JEvt::const_iterator shower = in.begin(); shower != in.end(); ++shower) {
        
        JShower3EY sh(getShower(*shower), 0.0);
 
        multiset<JDAQPMTIdentifier> top;
 
        const JModel<JPoint4D> match(JPoint4D(sh.getPosition(), sh.getT()), roadWidth_m, JRegressor_t::T_ns);
 
        for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
          
          if (match(*i)) {
            top.insert(i->getPMTIdentifier());
          }
        }
        
        JDetectorSubset_t subdetector(detector, sh.getPosition(), roadWidth_m);
 
        const JDirection3D conversion(sh.getDirection());           
        const JRotation3D R(conversion);
        
        vector<JPMTW0> buffer;
          
        for (JDetectorSubset_t::iterator module = subdetector.begin();
             module != subdetector.end(); ++module) {
            
          const JDAQSummaryFrame& frame = summary.getSummaryFrame(module->getID(), R_Hz);
 
          JModule dom(*module); 
 
          dom.rotate(R);  
 
          for (unsigned int i = 0; i != dom.size(); ++i) {
 
            if (getDAQStatus(frame, *module, i)            &&
                getPMTStatus(frame, *module, i)            &&
                frame[i].is_valid()                        &&
                !module->getPMT(i).has(PMT_DISABLE)) {
 
              const JDAQPMTIdentifier id(module->getID(), i);
                
              const double QE      = 1.0;
              const double rate_Hz = summary.getRate(id, this->R_Hz);
              const size_t count   = top.count(id);
 
              buffer.push_back(JPMTW0(dom.getPMT(i), QE, rate_Hz, count));
            }
          }   
        }
 
        this->step.resize(2);
        this->step[0] = JShowerEH(JPoint4D(JVector3D(), 0.0), JVersor3Z(), fit_step, 0.0, 0.0);
        this->step[1] = JShowerEH(JPoint4D(JVector3D(), 0.0), JVersor3Z(), 0.0, fit_step, 0.0);
 
        double f_h = 1 - 0.681 * (std::pow(shower->getE()/0.863, -0.207));
        
        double chi2 = (*this)(JShowerEH(JVertex3D(JVector3D(0,0,0), sh.getT()),
                                        JVersor3Z(), 
                                        log10(sh.getE()),
                                        log10(f_h*sh.getE()),
                                        sh.getBjY()),
                              buffer.begin(), buffer.end());
 
        double NDF = getCount(buffer.begin(), buffer.end()) - this->step.size(); 
        
        JShower3EY result(JShower3E(JShower3D(this->value.getVertex(),
                                              this->value.getDirection()),
                                    correct(this->value.getEem() + this->value.getEh())),
                          this->value.getBy());
 
        double y = getFinalBjY(this->value.getEem(), this->value.getEh());
        
        result.rotate_back(R);
        
        result.add(sh.getPosition());
          
        out.push_back(getFit(JHistory(shower->getHistory()).add(JSHOWER_BJORKEN_Y), result, getQuality(chi2), NDF, result.getE()));
 
        out.rbegin()->setW(5, y);
        out.rbegin()->setW(6, this->value.getEem());
        out.rbegin()->setW(7, this->value.getEh());
 
      }
 
      return out;
    }
  
    const JModuleRouter&            router;
    const JSummaryRouter&           summary;
    const JEnergyCorrection&        correct;
 
    /*
     * Function to get the reconstructed Bjorken Y from reconstructed energies.
     *
     * \param E_em   reconstructed EM shower energy
     * \param E_h    reconstructed Hadronic shower energy
     */
    static inline double getFinalBjY(double E_em, double E_h){
      return E_h / (E_em + E_h);
    }
 
  };
}
 
#endif