JPMTStatus.cc

Go to the documentation of this file.

#include <string>
#include <iostream>

#include "TROOT.h"
#include "TFile.h"
#include "TH2D.h"

#include "km3net-dataformat/definitions/pmt_status.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JCalibrate/JCalibrateToT.hh"
#include "JTools/JRange.hh"
#include "JLang/JLangToolkit.hh"
#include "JSupport/JMeta.hh"

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


/**
 * \file
 * Auxiliary program to update the PMT status in the detector file based on the fraction of short time-over-threshold values.
 *
 * \author acreusot, mdejong
 */
int main(const int argc, const char * const argv[])
{
  using namespace std;
  using namespace JPP;

  typedef JRange<double> JRange_t;

  string       inputFile;
  string       detectorFile;
  JRange_t     range;
  double       fraction;
  bool         overwriteDetector;
  int          debug;

  try {

    JParser<> zap("Auxiliary program to update the PMT status in the detector file based on the fraction of short time-over-threshold values.");

    zap['f'] = make_field(inputFile,         "input file (output from JCalibrateToT).");
    zap['a'] = make_field(detectorFile,      "detector file.");
    zap['x'] = make_field(range,             "integration range for noise.")                  = JRange_t(0.0, 8.0);
    zap['t'] = make_field(fraction,          "maximal fraction of signal allowed for noise.") = 0.5;
    zap['A'] = make_field(overwriteDetector, "overwrite detector file.");
    zap['d'] = make_field(debug,             "debug.")                                        = 1;

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


  JDetector detector;

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

  detector.comment.add(JMeta(argc,argv));

  if (detector.setToLatestVersion()) {
    NOTICE("Set detector version to " << detector.getVersion() << endl);
  }


  gErrorIgnoreLevel = kError;

  TFile* in = TFile::Open(inputFile.c_str(), "exist");

  if (in == NULL || !in->IsOpen()) {
    FATAL("File: " << inputFile << " not opened." << endl);
  }


  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {

    if (!module->empty()) {

      TH2D* h2s = (TH2D*) in->Get(MAKE_CSTRING(module->getID() << _2SToT));

      if (h2s != NULL) {

        for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {

          double noise  = 0.0;
          double signal = 0.0;

          for (int iy = 1; iy <= h2s->GetYaxis()->GetNbins(); ++iy) {

            const double x = h2s->GetYaxis()->GetBinCenter(iy);
            const double y = h2s->GetBinContent(ix, iy);

            signal += y;

            if (range(x)) {
              noise += y;
            }
          }

          const bool disable = noise > fraction * signal;

          if (debug >= debug_t || disable) {
            cout << "PMT " << FILL(10,'0') << module->getID() << "." << FILL(2,'0') << (ix - 1) << FILL()
                 << " noise/signal " << FIXED(7,0) << noise << "/" << FIXED(7,0) << signal << (disable ? " *" : "") << endl;
          }

          if (overwriteDetector) {
            if (disable) {
              module->getPMT(ix - 1).getStatus().set(PMT_DISABLE);
            }
          }
        }

      } else {

        WARNING("No histogram for module " << module->getID() << "; skip." << endl);
      }
    }
  }

  in->Close();
  
  if (overwriteDetector) {
  
    NOTICE("Store PMT status on file " << detectorFile << endl);

    store(detectorFile, detector);
  }
}