software/JReconstruction/JEnergyCorrection.cc

Go to the documentation of this file.

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

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TMath.h"
#include "TFitResult.h"

#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/definitions/reconstruction.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JEnergyCorrection.hh"
#include "JGeometry3D/JTrack3E.hh"
#include "JPhysics/JGeane.hh"
#include "JTools/JQuantile.hh"
#include "JROOT/JRootToolkit.hh"
#include "JGizmo/JGizmoToolkit.hh"
#include "JLang/JToken.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Utility program to determine energy correction.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  typedef JTriggeredFileScanner<JEvt>                       JTriggeredFileScanner_t;
  typedef JTriggeredFileScanner_t::multi_pointer_type       multi_pointer_type;
  typedef JToken<';'>      JToken_t;

  JTriggeredFileScanner_t  inputFile;
  JLimit_t&                numberOfEvents = inputFile.getLimit();
  string                   outputFile;  
  JQualitySorter           quality_sorter;
  string                   formula;
  vector<JToken_t>         startValues;
  vector<JToken_t>         fixedValues;
  string                   option;
  int                      debug;

  try {

    JParser<> zap("Utility program to determine energy correction.");

    zap['f'] = make_field(inputFile,   "event file(s) or histogram file");
    zap['o'] = make_field(outputFile)          = "energy_correction.root";
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['L'] = make_field(quality_sorter)      = JPARSER::initialised();
    zap['F'] = make_field(formula,     "fit formula, e.g: \"[0]+[1]*x\"")         = "";
    zap['@'] = make_field(startValues, "start values, e.g: \"p0 = GetMaximum;\"") = JPARSER::initialised();
    zap['='] = make_field(fixedValues, "fixed values, e.g: \"p0 = GetMaximum;\"") = JPARSER::initialised();
    zap['O'] = make_field(option,      "Fit option")                              = "";
    zap['d'] = make_field(debug)               =  2;

    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }


  TFile* in = NULL;
  TH2D*  h2 = NULL;

  if (inputFile.size() == 1) {

    in = TFile::Open(inputFile[0].c_str(), "READ");

    if (in != NULL && in->IsOpen()) {
      h2 = (TH2D*) in->Get("h2");
    }
  }

  if (h2 == NULL) {

    JHead head;

    try {
      head = getHeader(inputFile);
    }
    catch(const JException& error) {
      FATAL(error);
    }

    if (!head.is_valid(&JHead::cut_nu)) {
      FATAL("Missing neutrino parameters." << endl);
    }
    
    const JPosition3D center = get<JPosition3D>(head);

    const double   Emin  = head.cut_nu.Emin;
    const double   Emax  = head.cut_nu.Emax;

    const Int_t    nx    = 40;
    const Double_t xmin  = log10(Emin);
    const Double_t xmax  = log10(Emax);

    h2 = new TH2D("h2", NULL, nx, xmin, xmax, nx, xmin, xmax);

    while (inputFile.hasNext()) {

      STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);

      multi_pointer_type ps = inputFile.next();

      JEvt*  evt   = ps;
      Evt*   event = ps;

      double Emax  = 0.0;

      vector<Trk>::const_iterator muon = event->mc_trks.end();

      for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
        if (is_muon(*track)) {
          if (track->E > Emax) {
            muon = track;
          }
        }
      }

      if (muon == event->mc_trks.end()) {
        continue;
      }

      JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_event(JMUONENERGY));

      if (evt->begin() == __end) {
        continue;
      }

      JEvt::iterator p = min_element(evt->begin(), __end, quality_sorter);

      JTrack3E ta = getTrack(*muon);
      JTrack3E tb = getTrack(*p);

      ta.add(center);

      tb.setE(p->getW(JENERGY_ENERGY));
      tb.move(tb.getIntersection(ta.getPosition()), getSpeedOfLight(), gWater); 

      if (ta.getE() > 0.0 && tb.getE() > 0.0) {
        h2->Fill(log10(tb.getE()), log10(ta.getE()));
      }
    }
    STATUS(endl);
  }

  // short hands

  const Int_t    nx   = h2->GetXaxis()->GetNbins();
  const Double_t xmin = h2->GetXaxis()->GetXmin();
  const Double_t xmax = h2->GetXaxis()->GetXmax();

  const Int_t    ny   = h2->GetYaxis()->GetNbins();
  const Double_t ymin = h2->GetYaxis()->GetXmin();
  const Double_t ymax = h2->GetYaxis()->GetXmax();

  // 1D profile

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

  for (Int_t ix = 1; ix <= nx; ++ix) {

    TH1D h0("h0", NULL, ny, ymin, ymax);

    JQuantile Q;                          // start values

    for (int iy = 1; iy <= ny; ++iy) {

      const Double_t x = h2->GetYaxis()->GetBinCenter(iy);
      const Double_t y = h2->GetBinContent(ix,iy);
      const Double_t z = h2->GetBinError  (ix,iy);

      if (y > 0.0) {

        h0.SetBinContent(iy, y);
        h0.SetBinError  (iy, z);

        Q.put(x,y);
      }
    }

    if (Q.getCount() >= 3) {

      if (debug >= debug_t) {
        Q.print(cout);
      }

      TF1 f1("f1", "[0]*TMath::Gaus(x,[1],[2])");
 
      f1.SetParameter(0, Q.getWmax());
      f1.SetParameter(1, Q.getMean());
      f1.SetParameter(2, Q.getSTDev());

      TFitResultPtr result = h0.Fit(&f1, "SLLQ", "same",
                                    Q.getMean() - 3.0 * Q.getSTDev(),
                                    Q.getMean() + 3.0 * Q.getSTDev());

      if (debug >= debug_t || !result->IsValid()) {
        cout << "Bin: "
             << setw(4)    << ix                               << ' '
             << FIXED(7,3) << h2->GetXaxis()->GetBinCenter(ix) << ' '
             << FIXED(7,3) << result->Chi2()                   << '/'
             << result->Ndf()                                  << ' '
             << (result->IsValid() ? "" : "failed")            << endl;
      }

      if (result->IsValid()) {
        h1.SetBinContent(ix, f1.GetParameter(1));
        h1.SetBinError  (ix, f1.GetParError (1));
      }
    }
  }

  TF1* f1 = NULL;

  if (formula != "") {                  // fit

    f1 = new TF1("f1", formula.c_str());

    if (f1->IsZombie()) {
      FATAL("Function: " << formula << " is zombie." << endl);
    }

    try {

      for (vector<JToken_t>::const_iterator i = startValues.begin(); i != startValues.end(); ++i) {
        f1->SetParameter(getParameter(*i), getValue(*i,&h1));
      }

      for (vector<JToken_t>::const_iterator i = fixedValues.begin(); i != fixedValues.end(); ++i) {
        f1->FixParameter(getParameter(*i), getValue(*i,&h1));
      }
    }
    catch(JLANG::JParseError& error) {
      FATAL(error << endl);
    }

    if (option.find('S') == string::npos) { option += "S"; }

    TFitResultPtr result = h1.Fit(f1, option.c_str(), "same");

  } else {                              // interpolate

    Double_t x[2] = { xmin, xmax };
    Double_t y[2] = { xmin, xmax };
    Double_t z[2] = { 0.00, 0.00 };

    ostringstream os;

    os << "   "                         // extrapolate
       << "("
       << "x >= " << FIXED(5,3) << xmax
       << ")"
       << " * "
       << "(x)";
    os << "\n";

    for (Int_t ix = nx; ix >= 1; --ix) {

      x[0] = h1.GetXaxis()->GetBinCenter(ix);
      y[0] = h1.GetBinContent(ix);
      z[0] = h1.GetBinError  (ix);

      if (z[0] == 0.0 || y[1] < y[0]) {

        if (x[1] == xmax)
          continue;                     // continue
        else
          ix = 0;                       // terminate
        
        x[0] = xmin;
        y[0] = xmin;
      }

      os << " + "
         << "("
         << "x >= " << FIXED(5,3) << x[0]
         << " && "
         << "x <  " << FIXED(5,3) << x[1]
         << ")"
         << " * "
         << "(" << FIXED(5,3) << y[0] << " + " << FIXED(5,3) << (y[1] - y[0]) << " * (x - " << FIXED(5,3) << x[0] << ") / " << FIXED(5,3) << (x[1] - x[0]) << ")";
      os << "\n";

      x[1] = x[0];
      y[1] = y[0];
      z[1] = z[0];
    }

    os << " + "                         // extrapolate
       << "("
       << "x <  " << FIXED(5,3) << xmin
       << ")"
       << " * "
       << "(x)";
    os << "\n";

    string buffer = os.str();

    NOTICE(buffer << endl);

    for (string::iterator i = buffer.begin(); i != buffer.end(); ++i) {
      if (*i == '\n') {
        *i = ' ';
      }
    }

    f1 = new TF1("f1", buffer.c_str(), xmin, xmax);

    h1.GetListOfFunctions()->Add(f1);
  }


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

  out << JMeta(argc, argv);
  out << h1 << *h2;

  if (f1 != NULL) {                     // separately store energy correction
    
    f1->GetFormula()->SetName(JEnergyCorrection::getName());

    out << *(f1->GetFormula());
  }

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

  if (in != NULL) {
    in->Close();
  }
}