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/JDetectorToolkit.hh"
#include "JTools/JCombinatorics.hh"
#include "JROOT/JRootToolkit.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;

  JFitK40Parameters fitk40;               // setting of internal fit parameters

  string    outputFile;
  int       debug;

  try { 

    JProperties properties;

    properties.insert(gmake_property(fitk40.Rate_Hz));
    properties.insert(gmake_property(fitk40.p1));
    properties.insert(gmake_property(fitk40.p2));
    properties.insert(gmake_property(fitk40.p3));
    properties.insert(gmake_property(fitk40.p4));
 
    JParser<> zap("Auxiliary plot fit function of JFitK40.cc.");
    
    zap['@'] = make_field(properties)        = JPARSER::initialised();
    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);
  }


  const JModule module = getModule<JKM3NeT_t>(1001);

  const JCombinatorics combinatorics(module.size());

  const int    nx   = combinatorics.getNumberOfPairs();
  const double xmin = -0.5;
  const double xmax = nx - 0.5;


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

  TH1D h1("h1", NULL, nx, xmin, xmax);

  JFitK40 fit(module,
              h1.GetXaxis()->GetXmin(),
              h1.GetXaxis()->GetXmax(),
              -20.0,
              +20.0,
              true);

  fit.setModelParameters(fitk40.getModelParameters());

  fit.sort(JPairwiseComparator(module));

  for (int ix = 1; ix <= h1.GetXaxis()->GetNbins(); ++ix) {

    const Double_t x = h1.GetXaxis()->GetBinCenter(ix);

    const JFitK40::pair_type pair = fit.getPair((int) x);
    const double             rate = fit.getValue(pair);
 
    h1.SetBinContent(ix, rate);
  }

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