JTriggerTester.cc File Reference

Auxiliary program to verify processing of Monte Carlo events. More...

#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"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Auxiliary program to verify processing of Monte Carlo events.

Note that the available options are identical to those of JTriggerEfficiency.cc.

Author

mdejong

Definition in file JTriggerTester.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 84 of file JTriggerTester.cc.

{
  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.getMin() << " - " << QT.getMax() << endl);
}