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