JSTDevK40.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <array>
 
#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 "JDetector/JDetectorAddressMap.hh"
#include "JDetector/JDetectorSupportkit.hh"
 
#include "JROOT/JRootToolkit.hh"
 
#include "JCalibrate/JCalibrateK40.hh"
 
#include "JGizmo/JGizmoToolkit.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  const char* const address_t  =  "address";
  const char* const index_t    =  "index";
 
  array<string, 2> inputFile;
  string           outputFile;
  string           detectorFile;
  string           option;
  int              debug;
 
  try { 
 
    JParser<> zap;
    
    zap['f'] = make_field(inputFile,         "1st output of JMergeCalibrateK40 and 2nd output of JFitK40");
    zap['o'] = make_field(outputFile,        "output file.")             = "stdevk40.root";
    zap['a'] = make_field(detectorFile,      "detector file.");
    zap['O'] = make_field(option,            "axis label")               = address_t, index_t;
    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);
  }
 
  const JDetectorAddressMap& demo = getDetectorAddressMap(detector.getID());
 
 
  TFile* in[] = { NULL, NULL };
 
  for (int i = 0; i != 2; ++i) {
 
    in[i] = TFile::Open(inputFile[i].c_str(), "exist");
 
    if (in[i] == NULL || !in[i]->IsOpen()) {
      FATAL("File: " << inputFile[i] << " not opened." << endl);
    }
  }
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
 
    if (module->getFloor() == 0) { 
      continue;
    }
 
    TH2D* h2[] = { 
      (TH2D*) in[0]->Get(MAKE_CSTRING(module->getID() << _2R)),
      (TH2D*) in[1]->Get(MAKE_CSTRING(module->getID() << _2F))
    };
 
    if (h2[0] == NULL || h2[0]->GetEntries() == 0 ||
        h2[1] == NULL || h2[1]->GetEntries() == 0) {
      continue;
    }
 
    DEBUG("Module " << setw(10) << module->getID() << ' ' << getLabel(module->getLocation()) << endl);
 
    const JCombinatorics_t combinatorics(*module);
 
    const JModuleAddressMap memo = demo.get(module->getID());
    
    TH1D h1(MAKE_CSTRING(module->getID() << ".1D"), NULL, 
            h2[0]->GetXaxis()->GetNbins(), h2[0]->GetXaxis()->GetXmin(), h2[0]->GetXaxis()->GetXmax());
 
    TH2D hx(MAKE_CSTRING(module->getID() << ".2X"), NULL, 
            NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5, 
            NUMBER_OF_PMTS, -0.5, NUMBER_OF_PMTS - 0.5);
 
    for (int ix = 1; ix <= h2[0]->GetXaxis()->GetNbins(); ++ix) {
 
      const pair_type pair = combinatorics.getPair(ix - 1);
 
      double Y = 0.0;                                        // summed data
      double F = 0.0;                                        // summed function value
      double V = 0.0;                                        // summed standard deviation
      int    N = 0;
 
      for (int iy = 1; iy <= h2[0]->GetYaxis()->GetNbins(); ++iy) {
 
        const double y1 = h2[0]->GetBinContent(ix,iy);
        const double w1 = h2[0]->GetBinError  (ix,iy);
        const double f1 = h2[1]->GetBinContent(ix,iy);
      
        if (w1 > 0.0) {
          Y +=  y1;
          F +=  f1;
          V += (y1 - f1) / w1; 
          N += 1;
        }
      }
 
      if (N != 0) {
 
        V /= N;
 
        h1.SetBinContent(ix, V);
 
        hx.Fill(pair.first,  pair.second, V);
        hx.Fill(pair.second, pair.first,  V);
      }
 
      DEBUG(setw(3) << ix << ' '
            << "(" << FILL(2,'0') << pair.first              << "," << FILL(2,'0') << pair.second              << ")" << FILL() << ' '
            << "(" << memo.getPMTPhysicalAddress(pair.first) << "," << memo.getPMTPhysicalAddress(pair.second) << ")" << FILL() << ' '
            << FIXED(9,2) << Y << ' '
            << FIXED(9,2) << F << ' '
            << FIXED(9,2) << V << ' ' << (fabs(V) > 3.0 ? "***" : "") <<endl);
    }
 
    if (option == address_t) {
      setAxisLabels(hx, "X", memo);
      setAxisLabels(hx, "Y", memo);
    }
 
    out << h1 << hx;
  }
 
  for (int i = 0; i != 2; ++i) {
    in[i]->Close();
  }
  
  out.Write();
  out.Close();
}