JNarrabri.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/JGraph.hh"
#include "JROOT/JRootToolkit.hh"
#include "JROOT/JManager.hh"

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

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

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

#include "JDynamics/JDynamics.hh"

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


/**
 * \file
 * Example program to check contents of acoustic events.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  JMultipleFileScanner<JACOUSTICS::JEvt> inputFile;
  JLimit_t&       numberOfEvents = inputFile.getLimit();
  string          detectorFile;
  string          outputFile;
  double          Tmax_s;
  int             debug;

  try {

    JParser<> zap("Example program to check contents of acoustic events.");

    zap['f'] = make_field(inputFile,       "output of JKatoomba[.sh]");
    zap['n'] = make_field(numberOfEvents)      = JLimit_t::max();
    zap['a'] = make_field(detectorFile);
    zap['o'] = make_field(outputFile)          = "narrabri.root";
    zap['T'] = make_field(Tmax_s);
    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 dynamics(detector, Tmax_s);

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

  dynamics.load(inputFile);

  STATUS("OK" << endl);


  map<int, JGraph_t>  ZO;
  map<int, JGraph_t>  ZA;

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

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

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

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

  JManager<int, TH2D> H2(new TH2D  ("[%].tilt", NULL, 300, -3.0, +3.0, 300, -3.0, +3.0), '%', format);
  JManager<int, TH1D> HT(new TH1D  ("[%].time", NULL, 200,  0.0,  4.0),                  '%', format);

  JManager<int, TH1D> HO(new TH1D  ("H[%].orientation", NULL, 1000, xmin, xmax), '%', format);
  JManager<int, TH1D> HA(new TH1D  ("H[%].amplitude",   NULL, 1000, xmin, xmax), '%', format);

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

    TH1D* ho = HO[string->first];
    TH1D* ha = HA[string->first];

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

      const Double_t x = HO->GetXaxis()->GetBinCenter(i);

      try {

        const JMODEL::JString tilt = string->second(x);

        ho->SetBinContent(i, tilt.getAngle());
        ha->SetBinContent(i, tilt.getLength());
      }
      catch(const exception& error) {
        //ERROR("Error string " << setw(4) << string->first << " at " << FIXED(20,5) << x << ' ' << error.what() << endl);
      }
    }
  }


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

    if (string->second.size() > 1) {

      TH2D* h2 = H2[string->first];
      TH1D* ht = HT[string->first];

      for (typename JDynamics::JPosition::function_type::const_iterator q = string->second.begin(), p = q++;  q != string->second.end(); ++p, ++q) {

        const double t1 = q->getX() - p->getX();
        
        HT->Fill(log10(t1));
        ht->Fill(log10(t1));

        const double tx = 600.0e3 * (q->getY().tx - p->getY().tx) / t1;     // [mrad/10 min]
        const double ty = 600.0e3 * (q->getY().ty - p->getY().ty) / t1;     // [mrad/10 min]
        
        if (t1 > 0 && t1 < Tmax_s) {
          H2->Fill(tx, ty);
          h2->Fill(tx, ty);
        }
      }

      for (typename JDynamics::JPosition::function_type::const_iterator i = string->second.begin(); i != string->second.end(); ++i) {
        ZO[string->first].put(i->getX(), i->getY().getAngle());
        ZA[string->first].put(i->getX(), i->getY().getLength());
      }
    }
  }

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

  out << *H2 << H2 << *HT << HT << HO << HA;

  for (map<int, JGraph_t>::const_iterator i = ZO.begin(); i != ZO.end(); ++i) {
    out << JGraph(i->second, MAKE_CSTRING("G[" << FILL(4,'0') << i->first << "].orientation"));
  }

  for (map<int, JGraph_t>::const_iterator i = ZA.begin(); i != ZA.end(); ++i) {
    out << JGraph(i->second, MAKE_CSTRING("G[" << FILL(4,'0') << i->first << "].amplitude"));
  }

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