JIntegrateToT.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <map>

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JDetectorCalibration.hh"

#include "JSupport/JMeta.hh"

#include "JCalibrate/JCalibrateToT.hh"
#include "JCalibrate/JParameters.hh"

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

#include "JROOT/JManager.hh"

namespace {

  /**
   * Wild card for histogram naming.\n
   * The wild card will replaced by the module identifier.
   */
  static const char* const WILDCARD = "%";
}

/**
 * \file
 *
 * Auxiliary program to calculate the noise to signal ratio in
 * the ToT distribution
 * \author acreusot
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  typedef JRange<double> JRange_t;
  string       inputFile;
  string       outputFile;
  string       detectorFile;
  JRange_t     fitRange;
  string       regexp;
  int          runNb;
  int          debug;

  try { 

    JParser<> zap("Auxiliary program to integrate time-over-threshold distributions.");
    
    zap['f'] = make_field(inputFile,         "input file (output from JCalibrateToT).");
    zap['o'] = make_field(outputFile,        "output file.")                                                   = "results.txt";
    zap['a'] = make_field(detectorFile,      "detector file.");
    zap['x'] = make_field(fitRange,          "ROOT integration range for noise.")                          = JRange_t(0.0, 8.0);    
    zap['R'] = make_field(regexp,            "regular expression for histogram name.")                         = MAKE_CSTRING(WILDCARD << _2SToT);
    zap['r'] = make_field(runNb,             "current run number");
    zap['d'] = make_field(debug,             "debug.")                                                         = 1;

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


  //----------------------------------------------------------
  // load detector file
  //----------------------------------------------------------

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

  TFile* in = TFile::Open(inputFile.c_str(), "exist");
  if (in == NULL || !in->IsOpen()) {
    FATAL("File: " << inputFile << " not opened." << endl);
  }
  ofstream ostr(outputFile.c_str(), ios::out|ios::app);

  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
    DEBUG("Module " << module->getID() << endl);
    TH2D* h2s = (TH2D*) in->Get(replace(regexp, WILDCARD, MAKE_STRING(module->getID())).c_str());
    if (h2s == NULL) {
      WARNING("No histogram " << module->getID() << "; skip calculation." << endl);
      continue;
    }

    for (int ix = 1; ix <= h2s->GetXaxis()->GetNbins(); ++ix) {
      double noise = 0;
      double signal = 0;
      const int domId = module->getID();
      const JPMTIdentifier id(domId, ix - 1);
      const string name = MAKE_STRING(id.getPMTAddress());
      for (int iy = h2s->GetNbinsY(); iy >= 1; --iy) {
        const double y = h2s->GetBinContent(ix, iy);
        signal += y;
        if ((iy >= fitRange.getLowerLimit()) && (iy <= fitRange.getUpperLimit())) noise += y;
      }
      const double ratio = (signal != 0? noise/signal : 0);
      const JModule dom = moduleRouter.getModule(domId);
      //JPMT pmt = dom.getPMT(ix - 1);
      const JParameters curPar(runNb, domId, ix - 1, -1, noise, signal, ratio > 0.5);
      ostr << curPar << endl;
    }
  }
  ostr.close();
}