JTriggerTester.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <utility>

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

#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"

#include "JPhysics/JK40Rates.hh"

#include "km3net-dataformat/online/JDAQ.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTRouter.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JPMTParametersToolkit.hh"

#include "JTrigger/JTriggerParameters.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"

#include "JTools/JHashMap.hh"
#include "JTools/JWeight.hh"
#include "JTools/JQuantile.hh"

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

namespace {

  using namespace JPP;

  /**
   * Auxiliary data structure for monitoring number of photo-electrons.
   */
  struct tuple_type :
    public std::vector<JWeight>
  {
    /**
     * Default constructor.
     */
    tuple_type() :
      std::vector<JWeight> (3, JWeight())
    {}
  };


  /**
   * Hash for PMT identifier
   */
  struct JHashEvaluator_t {
    /**
     * Get hash value.
     *
     * \param  pmt             PMT identifier
     * \return                 hash value
     */
    inline int operator()(const JPMTIdentifier& pmt) const
    {
      return pmt.getID() * KM3NETDAQ::NUMBER_OF_PMTS  +  pmt.getPMTAddress();
    }
  };
}


/**
 * \file
 * Auxiliary program to verify processing of Monte Carlo events.
 *
 * Note that the available options are identical to those of JTriggerEfficiency.cc.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  JMultipleFileScanner<Evt> inputFile;
  JLimit_t&              numberOfEvents = inputFile.getLimit();
  string                 detectorFileA;
  string                 detectorFileB;
  JTriggerParameters     parameters;
  JPMTParametersMap      pmtParameters;
  JK40Rates              rates_Hz;
  int                    debug;

  try { 

    JParser<> zap("Auxiliary program to verify processing of Monte Carlo events.");
    
    zap['f'] = make_field(inputFile,           "input file (output of detector simulation)");
    zap['n'] = make_field(numberOfEvents)                     = JLimit::max();
    zap['a'] = make_field(detectorFileA,       "detector used for converion from Monte Carlo truth to raw data.");
    zap['b'] = make_field(detectorFileB,       "detector used for conversion of raw data to calibrated data.")     = "";
    zap['@'] = make_field(parameters,          "Trigger parameters (or corresponding file name)")                  = JPARSER::initialised();
    zap['P'] = make_field(pmtParameters,       "PMT simulation data (or corresponding file name)")                 = JPARSER::initialised();
    zap['B'] = make_field(rates_Hz,            "background rates [Hz]")                                            = JPARSER::initialised();
    zap['d'] = make_field(debug,               "debug")       = 2;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }

  cout.tie(&cerr);

  if (detectorFileB == "") {
    detectorFileB = detectorFileA;
  }


  JDetector detectorA;
  JDetector detectorB;

  try {
    load(detectorFileA, detectorA);
    load(detectorFileB, detectorB);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  NOTICE("Number of modules in detector A <" << detectorFileA << ">: " << setw(4) << detectorA.size() << endl);
  NOTICE("Number of modules in detector B <" << detectorFileB << ">: " << setw(4) << detectorB.size() << endl);

  JPMTParametersMap::Throw(true);

  if (!pmtParameters.is_valid()) {
    FATAL("Invalid PMT parameters " << pmtParameters << endl);
  }

  if (pmtParameters.getQE() != 1.0) {

    NOTICE("Correct background rates with global efficiency " << pmtParameters.getQE() << endl);

    rates_Hz.correct(pmtParameters.getQE());

    NOTICE("Back ground rates: " << rates_Hz << endl);
  }

  const JPMTRouter    pmtRouter   (detectorA);
  const JModuleRouter moduleRouter(detectorB);

  parameters.set(getMaximalDistance(detectorB));

  DEBUG("Trigger:"        << endl << parameters   << endl);
  DEBUG("PMT parameters:" << endl << pmtParameters << endl); 
  
  Head header;

  try {
    header = getHeader(inputFile);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  const JPosition3D center = get<JPosition3D>(header);

  NOTICE("Apply detector offset from Monte Carlo run header (" << center << ")" << endl);

  detectorA -= center;
  detectorB -= center;

  JHashMap<int,            JWeight>                      miss;
  JHashMap<JPMTIdentifier, JWeight,    JHashEvaluator_t> lost;
  JHashMap<JPMTIdentifier, JWeight,    JHashEvaluator_t> slip;
  JHashMap<JPMTIdentifier, tuple_type, JHashEvaluator_t> npe;

  JQuantile QT;
  
  while (inputFile.hasNext()) {

    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r');
    
    Evt* event = inputFile.next();

    if (!event->mc_hits.empty()) {
      QT.put(getTimeRange(*event).getLength());
    }

    for (vector<Hit>::const_iterator hit = event->mc_hits.begin(); hit != event->mc_hits.end(); ++hit) {

      if (!pmtRouter.hasPMT(hit->pmt_id)) {

        miss[hit->pmt_id].put(getNPE(*hit));

        continue;
      }

      const JPMTIdentifier pmt = pmtRouter.getIdentifier(hit->pmt_id);

      if (!moduleRouter.hasModule(pmt.getID())) {

        lost[pmt].put(getNPE(*hit));

        continue;
      }

      const double t0 = getTime(*hit);
      const double t1 = putTime(t0, pmtRouter   .getPMT(hit->pmt_id));
      const double t2 = getTime(t1, moduleRouter.getPMT(pmt));

      if (fabs(t2 - t0) > parameters.TMaxLocal_ns) {
        slip[pmt].put(getNPE(*hit));
      }

      const JPMTParameters data = pmtParameters.getPMTParameters(pmt);
      const double         P    = getSurvivalProbability(data);

      npe[pmt][0].put(getNPE(*hit));
      npe[pmt][1].put(getNPE(*hit) * P);
      npe[pmt][2].put(getNPE(*hit) * P * data.QE);
    }
  }  
  STATUS(endl);

  
  NOTICE("Number of PMTs absent in detector A:               " << setw(6) << miss.size() << endl);

  for (JHashMap<int, JWeight>::const_iterator i = miss.begin(); i != miss.end(); ++i) {
    DEBUG(setw(5) << i->first << ' ' << FIXED(8,0) << i->second.getTotal() << endl);
  }

  NOTICE("Number of PMTs absent in detector B:               " << setw(6) << lost.size() << endl);
  
  for (JHashMap<JPMTIdentifier, JWeight>::const_iterator i = lost.begin(); i != lost.end(); ++i) {
    DEBUG(setw(8) << i->first.getModuleID() << "[" << setw(2) << i->first.getPMTAddress() << "] " << FIXED(8,0) << i->second.getTotal() << endl);
  }

  NOTICE("Number of PMTs with t0 detector A - B > " << FIXED(4,1) << parameters.TMaxLocal_ns << " [ns]: " << setw(6) << slip.size() << endl);

  for (JHashMap<JPMTIdentifier, JWeight>::const_iterator i = slip.begin(); i != slip.end(); ++i) {
    DEBUG(setw(8) << i->first.getModuleID() << "[" << setw(2) << i->first.getPMTAddress() << "] " << FIXED(8,0) << i->second.getTotal() << endl);
  }


  NOTICE("Number of true photo-electrons, passed threshold and survived QE." << endl);
  
  tuple_type total;

  for (JHashMap<JPMTIdentifier, tuple_type>::const_iterator p = npe.begin(); p != npe.end(); ++p) {  

    DEBUG(setw(8) << p->first.getModuleID() << "[" << setw(2) << p->first.getPMTAddress() << "]");

    for (size_t i = 0; i != p->second.size(); ++i) {
      DEBUG(' ' << FIXED(8,0) << p->second[i].getTotal());
    }
    DEBUG(endl);

    for (size_t i = 0; i != p->second.size(); ++i) {
      total[i].put(p->second[i].getTotal());
    }
  }

  NOTICE(setw(12) << "total");

  for (size_t i = 0; i != total.size(); ++i) {
    NOTICE(' ' << FIXED(8,0) << total[i].getTotal());
  }

  NOTICE(endl);

  NOTICE("Time range of hits [ns]: " << QT.getXmin() << " - " << QT.getXmax() << endl);
}