JDrawK40.cc

Go to the documentation of this file.

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

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"

#include "JDetector/JModule.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JDetectorSupportkit.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 program to draw K40 fit function.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  JK40Parameters K40 = JK40Parameters::getInstance();

  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 program to draw K40 fit function.");
    
    zap['@'] = make_field(properties)        = JPARSER::initialised();
    zap['o'] = make_field(outputFile,        "output file.");
    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);

  double ct = -1.0;

  for (JModule::const_iterator i = module.begin(); i != module.end(); ++i) {
    for (JModule::const_iterator j = module.begin(); j != i; ++j) {

      if (getDot(i->getDirection(), j->getDirection()) > ct) {
        ct = getDot(i->getDirection(), j->getDirection());
      }
    }
  }

  const JCombinatorics_t combinatorics(module);

  cout << "Cosine minimal angle between PMT axes " << FIXED(5,3) << ct << endl;
  cout << "Default  rate [Hz] " << FIXED(6,3) << JK40Parameters::getInstance().getValue(ct) << '/' << FIXED(6,3) << JK40Parameters::getInstance().getValue(+1.0) << endl;
  cout << "Modified rate [Hz] " << FIXED(6,3) << K40                          .getValue(ct) << '/' << FIXED(6,3) << K40                          .getValue(+1.0) << endl;

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

  const Int_t    nx   = combinatorics.getNumberOfPairs();
  const Double_t xmin = -0.5;
  const Double_t xmax = nx - 0.5;

  TH1D h0("h0", NULL, 5000, -1.0, +1.0);
  TH1D h1("h1", NULL, nx, xmin, xmax);

  for (Int_t i = 1; i <= h0.GetXaxis()->GetNbins(); ++i) {

    const Double_t x = h0.GetBinCenter(i);

    h0.SetBinContent(i, K40.getValue(x));
  }

  for (Int_t i = 1; i <= h1.GetXaxis()->GetNbins(); ++i) {

    const int index = h1.GetBinCenter(i);
    
    const JCombinatorics_t::pair_type pair = combinatorics.getPair(index);

    const double x = getDot(module[pair.first].getDirection(), module[pair.second].getDirection());

    h1.SetBinContent(i, K40.getValue(x));
  }

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