JMuonPDF.cc

Go to the documentation of this file.

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

#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"

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

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

#include "JLang/JPredicate.hh"

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

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

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

#include "JFit/JLine1Z.hh"
#include "JFit/JModel.hh"

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

#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.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
 */
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;
  JMuonGandalfParameters_t parameters;
  int                      debug;

  try { 

    parameters.numberOfPrefits = 1;

    JParser<> zap("Program to perform detector calibration using reconstructed muon trajectories.");
    
    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 JMuonPDF_t    pdf(pdfFile, parameters.TTS_ns);
  const JMuonNPE_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 track = in->begin(); track != in->end(); ++track) {

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

      JRotation3D    R (getDirection(*track));
      JLine1Z        tz(getPosition (*track).rotate(R), track->getT());
      JRange<double> Z_m;

      if (track->hasW(JSTART_LENGTH_METRES)) {
        Z_m = JRange<double>(0.0, track->getW(JSTART_LENGTH_METRES))  +  JRange<double>(parameters.ZMin_m, parameters.ZMax_m);
      }

      const JModel<JLine1Z> match(tz, parameters.roadWidth_m, T_ns, Z_m);

      // 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
          
      sort(data.begin(), data.end(), compare);

      JDataW0_t::iterator __end = unique(data.begin(), data.end(), equal_to<JDAQPMTIdentifier>());
          

      double E_GeV = parameters.E_GeV;
      /*
      if (track->getE() > 0.1) {
        E_GeV = track->getE();
      }
      */

      DEBUG("line1z: "
            << FIXED(12,3) << tz.getX()     << ' '
            << FIXED(12,3) << tz.getY()     << ' '
            << FIXED(12,3) << tz.getZ()     << ' '
            << FIXED(12,3) << tz.getT()     << ' ' 
            << FIXED(12,1) << E_GeV         << ' '
            << FIXED( 8,3) << track->getQ() << endl);

      // move track to bary-center of hits
      /*
      int    N = 0;
      double z = 0.0;

      for (JDataW0_t::const_iterator hit = data.begin(); hit != __end; ++hit) {
        N += 1;
        z += hit->getZ();
      }

      if (N != 0) {
        tz.setZ(z/N, getSpeedOfLight());
      }
      */

      double Q = 0.0;

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

        const double x  = hit->getX() - tz.getX();
        const double y  = hit->getY() - tz.getY();
        const double z  = hit->getZ() - tz.getZ();
        const double R  = sqrt(x*x + y*y);

        const double t1 = tz.getT() + (z + R * getTanThetaC()) * 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  = gWater.getE(E_GeV, z);                  // correct for energy loss
        const double E  = E_GeV;
        const double dt = T_ns.constrain(hit->getT()  -  t1);

        JMuonPDF_t::result_type H1 = pdf.calculate(E, R, theta, phi, dt);
        JMuonPDF_t::result_type H0(hit->getR() * 1e-9, 0.0, T_ns);

        H1 += H0;                                                 // signal + background

        const double chi2 = H1.getChi2() - H0.getChi2();          // likelihood ratio

        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);


      double Y = 0.0;

      for (JDetector::const_iterator i = detector.begin(); i != detector.end(); ++i) {

        JPosition3D pos(i->getPosition());

        pos.transform(R, tz.getPosition());

        if (pos.getX() <= parameters.roadWidth_m) {

          JModule module = *i;

          module.transform(R, tz.getPosition());

          for (JModule::iterator pmt = module.begin(); pmt != module.end(); ++pmt) {

            const double R     = sqrt(pmt->getX()*pmt->getX() + pmt->getY()*pmt->getY());
            const double theta = pmt->getTheta();
            const double phi   = fabs(pmt->getTheta());
            const double y     = npe.calculate(1.0, R, theta, phi);   // MIP

            DEBUG("PMT: " 
                  << setw(10) << module.getID() << ':' << setw( 2) << setfill('0') << distance(module.begin(),pmt) << setfill(' ') << ' '
                  << FIXED(9,1)      << R      << ' '
                  << FIXED(6,4)      << theta  << ' '
                  << FIXED(6,4)      << phi    << ' '
                  << SCIENTIFIC(7,2) << y      << endl);

            Y += y;
          }
        }
      }

      DEBUG("npe: " << FIXED(7,1) << Y << endl);
    }
  }
  STATUS(endl);
}