JPlotPMTSystematics.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
 
#include "km3net-dataformat/online/JDAQ.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
 
#include "JTools/JRange.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 * Auxiliary application to convert PMT parameters text file to ROOT histograms.
 * \author mkarel
 */
 
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  typedef JRange<double>              JRange_t;
  typedef multimap<string, JRange_t>  map_type;
 
  string    detectorFile;
  string    inputFile;
  string    outputFile;
  map_type  parameters;
  int       debug;
 
  try { 
 
    JParser<> zap;
 
    zap['a'] = make_field(detectorFile);
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)      = "PMTsystematics.root";
    zap['p'] = make_field(parameters);
    zap['d'] = make_field(debug)           = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  using namespace KM3NETDAQ;
 
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  if (detector.empty()) {
    FATAL("Empty detector." << endl);
  }
 
  const int number_of_strings = getNumberOfStrings(detector);
  const int number_of_floors  = getNumberOfFloors (detector);
 
  const int ny = number_of_strings * number_of_floors;
 
 
  TFile* in = TFile::Open(inputFile.c_str(), "read");
 
  if (in == NULL) {
    FATAL("No data in inputfile " << inputFile << endl);
  }
 
  TFile out(outputFile.c_str(),"recreate");
 
  for (map_type::const_iterator i = parameters.begin(); i != parameters.end(); ++i) {
 
    const string   key   = i->first;
    const JRange_t range = i->second;
 
    TH1D h1(MAKE_CSTRING(key << ".1D"),
            NULL,
            100, range.first, range.second);
 
    TH2D h2(MAKE_CSTRING(key << ".2D"),
            NULL,
            NUMBER_OF_PMTS,  -0.5, NUMBER_OF_PMTS  - 0.5,
            ny,              -0.5, ny              - 0.5);
  
    // fill
 
    for (size_t dom = 0; dom != detector.size(); ++dom) {
  
      const JModule& module = detector.at(dom);
 
      TH1D* p1 = (TH1D*) in->Get(MAKE_CSTRING(module.getID() << ".1" << key));
  
      if (p1 == NULL) {
        continue;
      }
  
      const int iy = (module.getString() - 1) * number_of_floors  +  module.getFloor();
 
      for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
  
        const double value = p1->GetBinContent(pmt+1);
  
        h1.Fill(value);
 
        h2.SetBinContent(pmt+1, iy, value);
      }
    } 
 
    // axis labels
 
    h1.GetXaxis()->SetTitle(key.c_str());
    h1.GetYaxis()->SetTitle("PMTs");
 
    h2.GetXaxis()->SetTitle("PMT DAQ index");
 
    for (size_t dom = 0; dom != detector.size(); ++dom) {
 
      const JModule& module = detector.at(dom);
 
      h2.GetYaxis()->SetBinLabel(dom+1, MAKE_CSTRING("" 
                                                     << setw(3) << setfill('0') << module.getString() << ' ' 
                                                     << setw(2) << setfill('0') << module.getFloor()));
    }
 
    h2.GetZaxis()->SetTitle(key.c_str());
    
    // write
    h2.Write();
    h1.Write();
  }
 
  out.Close();
}