JAcoustics/JEvD.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <map>
#include <memory>

#include "km3net-dataformat/definitions/module_status.hh"

#include "TROOT.h"
#include "TApplication.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TH2D.h"
#include "TText.h"
#include "TMarker.h"

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

#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JGeometry.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JEvt.hh"
#include "JAcoustics/JEvtToolkit.hh"
#include "JAcoustics/JSupport.hh"

#include "JROOT/JStyle.hh"
#include "JROOT/JCanvas.hh"

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

#include "JSystem/JKeypress.hh"

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

namespace {

  /**
   * Wild card character for file name substition.
   */
  const char WILDCARD = '%';
}


/**
 * \file
 *
 * Program to display hit probabilities.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  
  JMultipleFileScanner<JEvt> inputFile;
  JLimit_t&                numberOfEvents = inputFile.getLimit();
  string                   detectorFile;
  string                   outputFile;
  JSoundVelocity           V              = getSoundVelocity;    // default sound velocity
  JCanvas                  canvas;
  size_t                   prescale;
  double                   ratio;
  bool                     batch;
  size_t                   Nmin;
  int                      debug;


  try { 

    JParser<> zap("Program to display strings.");
    
    zap['w'] = make_field(canvas,          "size of canvas <nx>x<ny> [pixels]")  = JCanvas(1200, 600);
    zap['f'] = make_field(inputFile,       "input file (output of JKatoomba[.sh])");
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['a'] = make_field(detectorFile);
    zap['o'] = make_field(outputFile,      "graphics output file name")          = MAKE_STRING("A_" << WILDCARD << ".gif");
    zap['V'] = make_field(V,               "sound velocity")                     = JPARSER::initialised();
    zap['M'] = make_field(getMechanics,    "mechanics data")                     = JPARSER::initialised();
    zap['P'] = make_field(prescale,        "prescale")                           = 1;
    zap['R'] = make_field(ratio,           "aspect ratio")                       = 0.1;
    zap['B'] = make_field(batch,           "batch processing");
    zap['N'] = make_field(Nmin,            "minimum number of transmissions")    = 0;
    zap['d'] = make_field(debug)               = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }

  if (batch  && outputFile == "") {
    FATAL("Missing output file name " << outputFile << " in batch mode." << endl);
  }

  if (!batch && outputFile == "") {
    outputFile = MAKE_STRING(WILDCARD << ".gif");
  }

  if (outputFile.find(WILDCARD) == string::npos) {
    FATAL("Output file name " << outputFile << " has no wild card '" << WILDCARD << "'" << endl);
  }

  JDetector detector;

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

  const JLocationRouter router  (detector);
  const JGeometry       geometry(detector);

  V.set(detector.getUTMZ());

  map<int, JStatus> status;
  
  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
    status[module->getID()] = module->getStatus();
  }

  // ROOT

  gROOT->SetBatch(batch);

  TApplication* tp = new TApplication("user", NULL, NULL);
  TCanvas*      cv = new TCanvas("display", "", canvas.x, canvas.y);

  unique_ptr<TStyle> gStyle(new JStyle("gplot", cv->GetWw(), cv->GetWh()));

  gROOT->SetStyle("gplot");
  gROOT->ForceStyle();

  const size_t NY = (size_t) (sqrt(geometry.size()) * 0.75 + 0.25);
  const size_t NX = (size_t) (geometry.size() + NY - 1) / NY;

  const size_t NUMBER_OF_PADS = NX * NY;

  cv->SetFillStyle(4000);
  cv->SetFillColor(kWhite);

  TPad* p1 = new TPad("p1", NULL, 0.0, 0.00, 1.0, 0.95);
  TPad* p2 = new TPad("p2", NULL, 0.0, 0.95, 1.0, 1.00);

  p1->Divide(NX, NY, 0.01, 0.001);

  p1->Draw();
  p2->Draw();

  double zmin = 0.0;
  double zmax = 0.0;

  for (JGeometry::const_iterator i = geometry.begin(); i != geometry.end(); ++i) {

    if (i->second.getZ() < zmin) {
      zmin = i->second.getZ();
    }

    if (i->second.rbegin()->getHeight() > zmax) {
      zmax = i->second.rbegin()->getHeight();
    }
  }

  zmin -=  2.0;
  zmax += 10.0;
  
  double xmin = -ratio * (zmax - zmin);
  double xmax = +ratio * (zmax - zmin);

  TH2D   H2[NUMBER_OF_PADS];
  TText  T2[NUMBER_OF_PADS];

  for (size_t i = 0; i != geometry.size(); ++i) {

    H2[i] = TH2D(MAKE_CSTRING(FILL(4,'0') << geometry.at(i).first << FILL()), NULL, 100, xmin, xmax, 100, zmin, zmax);

    DEBUG("Pad " << setw(3) << i + 1 << ' ' << H2[i].GetName() << endl);
    
    H2[i].GetXaxis()->SetTitle("x/y [m]");
    H2[i].GetYaxis()->SetTitle("z [m]");

    H2[i].GetXaxis()->CenterTitle(true);
    H2[i].GetYaxis()->CenterTitle(true);

    H2[i].SetStats(kFALSE);

    T2[i] = TText(0.0, zmax + 0.05 * (zmax - zmin), MAKE_CSTRING(FILL(4,'0') << geometry.at(i).first << FILL()));
    
    T2[i].SetTextAlign(21);
    T2[i].SetTextFont(42);
    T2[i].SetTextSize(0.05);

    p1->cd(i+1);

    H2[i].Draw("AXIS");
    T2[i].Draw();
  }
  
  typedef JTreeScanner<JEvent, JEvent::JEvaluator>  JTreeScanner_t;

  JTreeScanner_t in(inputFile);

  JTreeScanner_t::iterator p = in.begin();

  while (inputFile.hasNext()) {

    cout << "\revent: " << setw(8) << inputFile.getCounter() << flush;

    const JEvt*  evt   = inputFile.next();
    const JModel model = getModel(*evt);

    if (batch && inputFile.getCounter()%prescale != 0) {
      continue;
    }

    if (Nmin != 0) {

      for ( ; p != in.end() && p-> begin()->getToA() <  evt->UNIXTimeStart - 0.5; ++p) {}

      JTreeScanner_t::iterator q = p;

      for ( ; q != in.end() && q->rbegin()->getToA() <= evt->UNIXTimeStop  + 0.5; ++q) {}

      if (p != q) {

        map<int, size_t> buffer;

        for (JTreeScanner_t::iterator evt = p; evt != q; ++evt) {
          for (JEvent::const_iterator i = evt->begin(); i != evt->end(); ++i) {
            buffer[i->getID()] += 1;
          }
        }

        for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {

          module->setStatus(status[module->getID()]);

          if (buffer[module->getID()] < Nmin) {
            module->set(module->getFloor() == 0 ? HYDROPHONE_DISABLE : PIEZO_DISABLE);
          }
        }
      }
    }

    // graphics

    string utc = JDateAndTime((time_t) evt->UNIXTimeStart);

    TText title(0.5, 0.5, utc.c_str());

    title.SetTextAlign(22);
    title.SetTextFont(42);
    title.SetTextSize(0.3);

    p2->cd();

    title.Draw();
    
    vector<TMarker> X[NUMBER_OF_PADS];
    vector<TMarker> Y[NUMBER_OF_PADS];

    for (JEvt::const_iterator i = evt->begin(); i != evt->end(); ++i) {

      if (geometry.has(i->id)) {

        const int index = geometry.getIndex(i->id);

        const JGEOMETRY::JString& string     = geometry    [i->id];
        const JMODEL   ::JString& parameters = model.string[i->id];

        X[index].push_back(TMarker(0.0, string.getZ(), kFullCircle));
        Y[index].push_back(TMarker(0.0, string.getZ(), kFullCircle));
          
        for (size_t floor = 1; floor < string.size(); ++floor) {

          const JLocation   location(i->id, floor);
          const JPosition3D position = string.getPosition(parameters, floor);

          if (router.hasLocation(location) && !router.getModule(location).has(PIEZO_DISABLE)) {
            X[index].push_back(TMarker(position.getX() - string.getX(), position.getZ(), kFullCircle));
            Y[index].push_back(TMarker(position.getY() - string.getY(), position.getZ(), kFullCircle));
          }
        }
      }
    }

    // draw

    for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {

      p1->cd(i+1);

      for (auto& m1 : X[i]) { static_cast<TAttMarker&>(m1) = TAttMarker(kBlack, kFullCircle, 0.7/sqrt(NY)); m1.Draw(); }
      for (auto& m1 : Y[i]) { static_cast<TAttMarker&>(m1) = TAttMarker(kRed,   kFullCircle, 0.7/sqrt(NY)); m1.Draw(); }
    }
      
    cv->Update();

    // action

    if (batch) {

      cv->SaveAs(replace(outputFile, WILDCARD, MAKE_STRING(FIXED(12,0,'0') << evt->UNIXTimeStart << FILL())).c_str());

    } else {

      static int count = 0;

      if (count++ == 0) {
        cout << endl << "Type '?' for possible options." << endl;
      }

      for (bool user = true; user; ) {
        
        cout << "\n> " << flush;

        switch (JKeypress(true).get()) {

        case '?':
          cout << endl;
          cout << "possible options: " << endl;
          cout << 'q' << " -> " << "exit application"                            << endl;
          cout << 'u' << " -> " << "update canvas"                               << endl;
          cout << 's' << " -> " << "save graphics to file"                       << endl;
          cout << 'r' << " -> " << "rewind input file"                           << endl;
          cout << 'R' << " -> " << "switch to ROOT mode (quit ROOT to continue)" << endl;
          cout << 'p' << " -> " << "print event information"                     << endl;
          cout << ' ' << " -> " << "next event (as well as any other key)"       << endl;
          break;

        case 'q':
          cout << endl;
          return 0;

        case 'u':
          cv->Update();
          break;

        case 's':
          cv->SaveAs(replace(outputFile, WILDCARD, MAKE_STRING(FIXED(12,0,'0') << evt->UNIXTimeStart << FILL())).c_str());
          break;

        case 'R':
          tp->Run(kTRUE);
          break;

        case 'p':
          cout << endl << *evt << endl;
          break;

        case 'r':
          inputFile.rewind();
          
        default:
          user  = false;
          break;
        }
      }
    }
  }
  cout << endl;
}