JReconstruction/JShowerNPE.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>

#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"

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

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

#include "JSupport/JSupport.hh"
#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JSummaryFileRouter.hh"

#include "JPhysics/JNPE_t.hh"
#include "JPhysics/JPDF_t.hh"

#include "JFit/JPoint4D.hh"
#include "JFit/JModel.hh"
#include "JFit/JFitToolkit.hh"

#include "JReconstruction/JHitW0.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JShowerFitParameters_t.hh"

#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

#include "JTools/JResult.hh"

namespace {

  /**
   * Auxiliary data structure for sorting of hits.
   */
  static const struct {
    /**
     * Compare hits by PMT identifier and time.
     *
     * \param  first         first  hit
     * \param  second        second hit
     * \return               true if first before second; else false
     */
    template<class T>
    bool operator()(const T& first, const T& second) const
    {
      using namespace std;
      using namespace JPP;
      using namespace KM3NETDAQ;

      if (equal_to<JDAQPMTIdentifier>()(first, second))
        return less<JHit>()(first, second);
      else
        return less<JDAQPMTIdentifier>()(first, second);
    }
  } compare;
}


/**
 * \file
 *
 * Program to evaluate hit probabilities.
 *
 * \author mdejong and adomi
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
  
  typedef JParallelFileScanner< JTypeList<JDAQEvent, JEvt> >       JParallelFileScanner_t;
  typedef JParallelFileScanner_t::multi_pointer_type               multi_pointer_type;

  JParallelFileScanner_t   inputFile;
  JLimit_t&                numberOfEvents = inputFile.getLimit();
  string                   detectorFile;
  string                   pdfFile;
  JShowerFitParameters_t   parameters;
  int                      debug;

  try { 

    parameters.numberOfPrefits = 1;

    JParser<> zap("Program to evaluate hit probabilities.");
    
    zap['f'] = make_field(inputFile);
    zap['a'] = make_field(detectorFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['P'] = make_field(pdfFile);
    zap['@'] = make_field(parameters)          = JPARSER::initialised();
    zap['d'] = make_field(debug)               = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }

  //setDAQLongprint(debug >= JEEP::debug_t);

  cout.tie(&cerr);

  if (parameters.numberOfPrefits != 1) {
    WARNING("Number of prefits " << parameters.numberOfPrefits << " != " << 1 << endl);
  }

  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  const JModuleRouter router(detector);

  JSummaryFileRouter  summary(inputFile, parameters.R_Hz);

  const JShowerNPE_t    npe(pdfFile);
  
  const JTimeRange    T_ns(parameters.TMin_ns, parameters.TMax_ns);

  typedef vector<JHitL0>  JDataL0_t;
  typedef vector<JHitW0>  JDataW0_t;

  const JBuildL0<JHitL0>  buildL0;

  while (inputFile.hasNext()) {

    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);

    multi_pointer_type ps = inputFile.next();

    JDAQEvent* tev = ps;
    JEvt*      in  = ps;

    DEBUG("event: " << *tev << endl);

    summary.update(*tev);

    in->select(parameters.numberOfPrefits, qualitySorter);

    JDataL0_t dataL0;
      
    buildL0(*tev, router, true, back_inserter(dataL0));

    for (JEvt::const_iterator shower = in->begin(); shower != in->end(); ++shower) {

      DEBUG("shower: " << *shower << endl);

      JRotation3D    R (getDirection(*shower));
      JPoint4D       vx(getPosition(*shower).rotate(R), shower->getT());

      const JModel<JPoint4D> match(vx, parameters.roadWidth_m, T_ns);

      // hit selection based on start value

      JDataW0_t data;
          
      for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
            
        double rate_Hz = summary.getRate(*i);

        if (rate_Hz <= 0.0) {
          rate_Hz = summary.getRate();
        }

        JHitW0 hit(*i, rate_Hz);
            
        hit.rotate(R);

        if (match(hit)) {
          data.push_back(hit);
        }
      }
          
      // select first hit in PMT
          
      sort(data.begin(), data.end(), compare);

      JDataW0_t::iterator __end = unique(data.begin(), data.end(), equal_to<JDAQPMTIdentifier>());
          
      DEBUG("point4d: "
            << FIXED(12,3) << vx.getX()     << ' '
            << FIXED(12,3) << vx.getY()     << ' '
            << FIXED(12,3) << vx.getZ()     << ' '
            << FIXED(12,3) << vx.getT()     << ' ' 
            << FIXED( 8,3) << shower->getQ() << endl);

      double Q = 0.0;

      for (JDataW0_t::const_iterator hit = data.begin(); hit != __end; ++hit) {

        JPosition3D D(hit->getPosition());
        D.sub(vx);

        const double x  = hit->getX() - vx.getX();
        const double y  = hit->getY() - vx.getY();
        const double z  = hit->getZ() - vx.getZ();
        const double cd = z/D.getLength();  // Delta = angle between shower direction and PMT position

        const double t1 = vx.getT() + (D.getLength() * getIndexOfRefraction() * getInverseSpeedOfLight());
        
        JDirection3D u(hit->getDX(), hit->getDY(), hit->getDZ());  // PMT orientation

        u.rotate(JRotation3Z(-atan2(y,x)));                        // rotate PMT axis to x-z plane

        const double theta = u.getTheta();
        const double phi   = fabs(u.getPhi());                     // rotational symmetry of Cherenkov cone
        
        const double E  = shower->getE();
        const double dt = T_ns.constrain(hit->getT()  -  t1);

        double H1 = npe.calculate(E, D.getLength(), cd, theta, phi);
        double H0 = hit->getR() * 1e-9 * T_ns.getLength();

        double  Vmax_npe = 20.0;
        if (H1 >= Vmax_npe) {
          H1 *= Vmax_npe / H1;
        }

        H1 += H0;                                                // signal + background

        const double chi2 = getChi2(H1, true);                   // -log(lik)    

        DEBUG("hit: " 
              << setw(10) << hit->getModuleID() << ':' << setw( 2) << setfill('0') << hit->getPMTAddress() << setfill(' ') << ' '
              << FIXED(12,1) << E           << ' '
              << FIXED( 9,1) << R           << ' '
              << FIXED( 6,4) << theta       << ' '
              << FIXED( 6,4) << phi         << ' '
              << FIXED( 8,3) << dt          << ' '
              << FIXED(12,3) << chi2 << endl);

        Q += getQuality(chi2);
      }

      DEBUG("quality: " << FIXED(8,3) << Q << ' ' << distance(data.begin(), __end) << endl);
    }
  }
  STATUS(endl);
}