JCompareKatoomba.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 "TProfile.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 compare 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 compare acoustic fit results.");

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

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


  if (inputFile.size() != 2u) {
    FATAL("Invalid number of input files; need 2 files for comparison." << endl);
  }

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

  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) {

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


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

  out << H2 << *H2;

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