JMuonPDF.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 "JDetector/JPMTParametersMap.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
    {
      if (first.getPMTIdentifier() == second.getPMTIdentifier())
        return first.getT()             < second.getT();
      else
        return first.getPMTIdentifier() < second.getPMTIdentifier();
    }
  } 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;
  JPMTParametersMap        pmtParameters;
  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['F'] = make_field(pdfFile);
    zap['@'] = make_field(parameters)          = JPARSER::initialised();
    zap['P'] = make_field(pmtParameters)       = JPARSER::initialised();
    zap['d'] = make_field(debug)               = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
  //setDAQLongprint(debug >= JEEP::debug_t);
 
 
  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);
 
  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) {
 
        const JPMTIdentifier id(i->getModuleID(), i->getPMTAddress());
 
        const JPMTParameters& wip = pmtParameters.getPMTParameters(id);
 
        const int    type = wip.getType();
        const double QE   = wip.QE;
        const double R_Hz = summary.getRate(i->getPMTIdentifier(), parameters.R_Hz);
 
        JHitW0 hit(*i, type, QE, R_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>());
          
 
      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->getPhi());
            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);
}