JCombinatoricsK40.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>

#include "TROOT.h"
#include "TFile.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));

  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  h1(MAKE_CSTRING(pmt.getModuleID() << '.' << FILL(2,'0') << pmt.getTDC()), NULL, h2->GetYaxis()->GetNbins(), h2->GetYaxis()->GetXmin(), h2->GetYaxis()->GetXmax());
  
  for (size_t i = 0; i != module.size(); ++i) {

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

    for (Int_t iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy) {
      h1.SetBinContent(iy, h2->GetBinContent(ix,iy));
      h1.SetBinError  (iy, h2->GetBinError  (ix,iy));
    }
  }

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

  out << h1;
  
  out.Write();
  out.Close();  
}