JCompareDynamicPositions.cc

Go to the documentation of this file.

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

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

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

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"

#include "JDynamics/JDynamics.hh"

#include "JTools/JQuantile.hh"

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

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

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


/**
 * 
 * Example program to compare dynamic position calibrations.
 *
 * \author cgatius
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  JMultipleFileScanner<JACOUSTICS::JEvt> calibrationFile1;
  JMultipleFileScanner<JACOUSTICS::JEvt> calibrationFile2;
  string          detectorFile;
  string          outputFile;
  double          Tstep_s;
  int             debug;

  try {

    JParser<> zap("Example program to compare dynamic position calibrations.");

    zap['a'] = make_field(detectorFile);
    zap['f'] = make_field(calibrationFile1, "output of JKatoomba[.sh]");
    zap['F'] = make_field(calibrationFile2, "output of JKatoomba[.sh]");
    zap['o'] = make_field(outputFile) = "compare_dynamic_positions.root";
    zap['T'] = make_field(Tstep_s) = 600.0;
    zap['d'] = make_field(debug)               = 1;

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


  JDetector detector;

  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  JDynamics dynamics1(detector, 0.0);

  STATUS("loading input from file " << calibrationFile1 << "... " << flush);

  dynamics1.load(calibrationFile1);

  STATUS("OK" << endl);

  Double_t xmin1 = numeric_limits<Double_t>::max();
  Double_t xmax1 = numeric_limits<Double_t>::lowest();

  for (JDynamics::JPosition::const_iterator string = dynamics1.position.begin(); string != dynamics1.position.end(); ++string) {

    if (!string->second.empty()) {
      xmin1 = min(xmin1, string->second.getXmin());
      xmax1 = max(xmax1, string->second.getXmax());
    }
  }

  JDynamics dynamics2(detector, 0.0);

  STATUS("loading input from file " << calibrationFile1 << "... " << flush);

  dynamics2.load(calibrationFile2);

  STATUS("OK" << endl);

  Double_t xmin2 = numeric_limits<Double_t>::max();
  Double_t xmax2 = numeric_limits<Double_t>::lowest();

  for (JDynamics::JPosition::const_iterator string = dynamics2.position.begin(); string != dynamics2.position.end(); ++string) {

    if (!string->second.empty()) {
      xmin2 = min(xmin2, string->second.getXmin());
      xmax2 = max(xmax2, string->second.getXmax());
    }
  }

  Double_t xmin = numeric_limits<Double_t>::max();
  Double_t xmax = numeric_limits<Double_t>::lowest();

  //latest start time and earliest end time
  xmin = max(xmin1, xmin2);
  xmax = min(xmax1, xmax2);

  const JFormat_t format(4, 0, std::ios_base::fmtflags(), '0');

  map<int,JGraph_t> HO;
  map<int,JGraph_t> HA;
  JManager<int,TH1D> HC(new TH1D  ("M[%].comparetx", NULL, 100,  0, 2),'%', format);
  JManager<int,TH1D> HD(new TH1D  ("N[%].comparety", NULL, 100,  0, 2),'%', format);

  for (JDynamics::JPosition::const_iterator string1 = dynamics1.position.begin(); string1 != dynamics1.position.end(); ++string1) {
    for (JDynamics::JPosition::const_iterator string2 = dynamics2.position.begin(); string2 != dynamics2.position.end(); ++string2) {

      if (string1->first == string2->first) {

        TH1D* hc = HC[string1->first];
        TH1D* hd = HD[string1->first];

        for (Double_t t = xmin; t <= xmax; t = t + Tstep_s) {

          HO[string1->first].put(t, (string1->second(t).tx - string2->second(t).tx) * 1000);
          HA[string1->first].put(t, (string1->second(t).ty - string2->second(t).ty) * 1000);
          hc->Fill((string1->second(t).tx - string2->second(t).tx) * 1000);
          hd->Fill((string1->second(t).ty - string2->second(t).ty) * 1000);
        }
           
        break;
      }
    }
  }

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

  out << HC << HD;
  for (map<int, JGraph_t>::const_iterator i = HO.begin(); i != HO.end(); ++i) {
    out << JGraph(i->second, MAKE_CSTRING("H[" << i->first << "].comparetx"));
  }

  for (map<int, JGraph_t>::const_iterator i = HA.begin(); i != HA.end(); ++i) {
    out << JGraph(i->second, MAKE_CSTRING("K[" << i->first << "].comparety"));
  }
}