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 "JGeometry3D/JCylinder3D.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);
 
  const JCylinder3D cylinder(detector.begin(), detector.end());
 
  const double Z = cylinder.getZmax() * 0.5;
  
  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 + q->getY().tx2 * Z) - (p->getY().tx + p->getY().tx2 * Z)) / t1;     // [mrad/10 min]
        const double ty = 600.0e3 * ((q->getY().ty + q->getY().ty2 * Z) - (p->getY().ty + p->getY().ty2 * Z)) / 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();
}