JPlotK40.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
 
#include "JDetector/JModule.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDetectorSupportkit.hh"
 
#include "JCalibrate/JCalibrateK40.hh"
#include "JCalibrate/JFitK40.hh"
 
#include "Jeep/JProperties.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Auxiliary plot fit function of JFitK40.cc.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  JK40Parameters K40 = JK40Parameters::getInstance();
 
  string    detectorFile;
  string    outputFile;
  int       debug;
 
  try { 
 
    JProperties properties;
 
    properties.insert(gmake_property(K40.R));
    properties.insert(gmake_property(K40.p1));
    properties.insert(gmake_property(K40.p2));
    properties.insert(gmake_property(K40.p3));
    properties.insert(gmake_property(K40.p4));
 
    JParser<> zap("Auxiliary plot fit function of JFitK40.cc.");
    
    zap['@'] = make_field(properties)        = JPARSER::initialised();
    zap['a'] = make_field(detectorFile)                                 = "";
    zap['o'] = make_field(outputFile,        "output file.")            = "k40.root";
    zap['d'] = make_field(debug,             "debug.")                  = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  JModule module;
 
  if (detectorFile != "") {
 
    JDetector detector;
 
    try {
      load(detectorFile, detector);
    }
    catch(const JException& error) {
      FATAL(error);
    }
 
    module = detector[0];
 
  } else {
  
    module = getModule<JKM3NeT_t>(1001);
  }
 
  const JModel  model(module, K40);
 
  const int    nx   = model.getNumberOfPairs();
  const double xmin = -0.5;
  const double xmax = nx - 0.5;
 
  if (debug >= debug_t) {
    for (size_t i = 0; i != module.size(); ++i) {
      DEBUG("PMT " << FILL(2,'0') << i << ' ' << module[i].getDirection() << endl);
    }
  }
 
  TFile out(outputFile.c_str(), "recreate");
 
  TH1D h1("h1", NULL, nx, xmin, xmax);
 
  double R = 0.0;
 
  for (int ix = 1; ix <= h1.GetXaxis()->GetNbins(); ++ix) {
 
    const Double_t  x    = h1.GetXaxis()->GetBinCenter(ix);
 
    const pair_type pair = model.getPair((int) x);
    const double    rate = model.getValue(pair);
 
    R += rate;
 
    DEBUG("bin " << setw(3) << ix                                                                   << ' '
          << "(" << FILL(2,'0') << pair.first << "," << FILL(2,'0') << pair.second << ")" << FILL() << ' '
          << FIXED(6,3) << getDot(module[pair.first] .getDirection(),
                                  module[pair.second].getDirection())                               << ' '
          << FIXED(6,3) << rate << endl);
 
    h1.SetBinContent(ix, rate);
  }
 
  cout << "Total rate " << FIXED(7,3) << R << " [Hz]." << endl;
 
  out.Write();
  out.Close();
}