JCombinatoricsK40.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 "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTIdentifier.hh"

#include "JCalibrate/JCalibrateK40.hh"
#include "JROOT/JRootToolkit.hh"

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

/**
 * \file
 *
 * Auxiliary program to project single PMT data from 2D histogram.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  string         detectorFile;
  string         inputFile;
  string         outputFile;
  JPMTIdentifier pmt;
  string         extension;
  int            debug;
    
  try { 

    JParser<> zap("Auxiliary program to project single PMT data from 2D histogram.");
    
    zap['a'] = make_field(detectorFile,   "detector file.");
    zap['f'] = make_field(inputFile,      "input file.");
    zap['o'] = make_field(outputFile,     "output file.")             = "k40.root";
    zap['P'] = make_field(pmt,            "PMT identifier");
    zap['e'] = make_field(extension,      "histogram name extension") = _2R, _2S, _2F;
    zap['d'] = make_field(debug,          "debug flag.")              = 1;
     
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  gErrorIgnoreLevel = kError;


  JDetector detector;

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

  const JModuleRouter router(detector);

  const JModule&      module = router.getModule(pmt.getModuleID());
  
  JCombinatorics combinatorics;

  combinatorics.configure(module.size());

  combinatorics.sort(JPairwiseComparator(module));

  if (debug >= debug_t) {

    for (size_t i = 0; i != combinatorics.getNumberOfPairs(); ++i) {

      const JCombinatorics::pair_type pair = combinatorics.getPair(i);

      cout << setw(3) << i << " -> (" << FILL(2,'0') << pair.first << "," << FILL(2,'0') << pair.second << FILL() << ")" << endl;
    }
  }

  TFile in(inputFile.c_str(), "read");

  TH2D* h2 = (TH2D*) in.Get(MAKE_CSTRING(pmt.getModuleID() << extension));

  if (h2 == NULL) {
    FATAL("Missing histogram for module " << pmt.getModuleID() << endl);
  }
  
  TH1D  ha(MAKE_CSTRING(pmt.getModuleID() << '.' << FILL(2,'0') << pmt.getTDC() << ".1D"), NULL,
           h2->GetYaxis()->GetNbins(), h2->GetYaxis()->GetXmin(), h2->GetYaxis()->GetXmax());

  TH2D  hb(MAKE_CSTRING(pmt.getModuleID() << '.' << FILL(2,'0') << pmt.getTDC() << ".2D"), NULL,
           module.size(), -0.5, module.size() - 0.5,
           h2->GetYaxis()->GetNbins(), h2->GetYaxis()->GetXmin(), h2->GetYaxis()->GetXmax());
  
  for (int i = 0; i != (int) module.size(); ++i) {

    if (i != pmt.getTDC()) {

      const Int_t ix = combinatorics.getIndex(pmt.getTDC(), i) + 1;

      for (Int_t iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy) {

        ha.SetBinContent(iy,       ha.GetBinContent(iy) + h2->GetBinContent(ix,iy));
        ha.SetBinError  (iy, hypot(ha.GetBinError  (iy),  h2->GetBinError  (ix,iy)));

        hb.SetBinContent(i + 1, iy, h2->GetBinContent(ix,iy));
        hb.SetBinError  (i + 1, iy, h2->GetBinError  (ix,iy));
      }
    }
  }

  TFile out(outputFile.c_str(), "recreate");

  out << ha << hb;
  
  out.Write();
  out.Close();  
}