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