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