JParramatta.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>

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

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JTreeScanner.hh"

#include "JROOT/JRootToolkit.hh"
#include "JROOT/JGraph.hh"
#include "JROOT/JManager.hh"

#include "JAcoustics/JEvt.hh"
#include "JAcoustics/JSupport.hh"

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


/**
 * \file
 *
 * Example program to plot acoustic fit results.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  JMultipleFileScanner<JEvt> inputFile;
  JLimit_t&           numberOfEvents = inputFile.getLimit();
  string              outputFile;
  int                 debug;

  try {

    JParser<> zap("Example program to plot acoustic fit results.");

    zap['f'] = make_field(inputFile,       "input file (output of JKatoomba)");
    zap['n'] = make_field(numberOfEvents)   = JLimit::max();
    zap['o'] = make_field(outputFile)       = "parramatta.root";
    zap['d'] = make_field(debug)            = 2;

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


  if (inputFile.size() < 1u || inputFile.size() > 2u) {
    FATAL("Invalid number of input files; 1 file for histograms and 2 files for comparison." << endl);
  }


  TH1D h1("chi2/NDF", NULL, 500, 0.0, 10.0);

  JGraph_t g0;
  JGraph_t g1;
  
  JManager<int, TH2D> H2(new TH2D("string[%]", NULL, 500, -20.0, +20.0, 500, -20.0, +20.0));


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

    JTreeScanner<JEvt, JEvt::JEvaluator> inA(inputFile[0]);
    JTreeScanner<JEvt, JEvt::JEvaluator> inB(inputFile[1]);

    while (inA.hasNext() && inB.hasNext()) {

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

      JEvt* pA = inA.next();
      JEvt* pB = inB.next();

      // find the same event based on the start and stop times of the events

      while (pA->UNIXTimeStop < pB->UNIXTimeStart && inA.hasNext()) { pA = inA.next(); }
      while (pB->UNIXTimeStop < pA->UNIXTimeStart && inB.hasNext()) { pB = inB.next(); }

      if (pA->UNIXTimeStart < pB->UNIXTimeStop &&
          pB->UNIXTimeStart < pA->UNIXTimeStop) {

        for (JEvt::const_iterator iA = pA->begin(); iA != pA->end(); ++iA) {
          for (JEvt::const_iterator iB = pB->begin(); iB != pB->end(); ++iB) {

            if (iA->id == iB->id) {
              H2        ->Fill((iA->tx - iB->tx)*1.0e3, (iA->ty - iB->ty)*1.0e3);  // [mrad]
              H2[iA->id]->Fill((iA->tx - iB->tx)*1.0e3, (iA->ty - iB->ty)*1.0e3);  // [mrad]
            }
          }
        }
      }
    }
  }

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

    while (inputFile.hasNext()) {

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

      const JEvt* evt = inputFile.next();

      h1.Fill(evt->chi2 / evt->ndf);

      g0.put(0.5 * (evt->UNIXTimeStart + evt->UNIXTimeStop), evt->ndf);
      g1.put(0.5 * (evt->UNIXTimeStart + evt->UNIXTimeStop), evt->chi2 / evt->ndf);

      for (JEvt::const_iterator i = evt->begin(); i != evt->end(); ++i) {
        H2       ->Fill(i->tx*1.0e3, i->ty*1.0e3);  // [mrad]
        H2[i->id]->Fill(i->tx*1.0e3, i->ty*1.0e3);  // [mrad]
      }
    }
    STATUS(endl);
  }


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

  out <<  H2
      << *H2;

  if (inputFile.size() == 1) {
    out << h1
        << JGraph(g0, "g0")
        << JGraph(g1, "g1");
  }

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