JReconstruction/JEvD.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include <memory>
 
#include "TROOT.h"
#include "TApplication.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TH2D.h"
#include "TArrow.h"
#include "TLatex.h"
#include "TMarker.h"
 
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/online/JDAQ.hh"
#include "km3net-dataformat/tools/time_converter.hh"
 
#include "JROOT/JStyle.hh"
#include "JROOT/JCanvas.hh"
#include "JROOT/JRootToolkit.hh"
 
#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModuleRouter.hh"
 
#include "JDynamics/JDynamics.hh"
 
#include "JTrigger/JHitL0.hh"
#include "JTrigger/JBuildL0.hh"
 
#include "JSupport/JSupport.hh"
#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JSummaryFileRouter.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
 
#include "JPhysics/JPDF_t.hh"
 
#include "JFit/JLine1Z.hh"
#include "JFit/JModel.hh"
#include "JReconstruction/JHitW0.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JReconstruction/JMuonGandalfParameters_t.hh"
#include "JReconstruction/JEventSelector.hh"
 
#include "JLang/JPredicate.hh"
#include "JLang/JComparator.hh"
 
#include "JMath/JMathToolkit.hh"
#include "JSystem/JKeypress.hh"
#include "JSystem/JProcess.hh"
 
#include "Jeep/JFunctionAdaptor.hh"
#include "Jeep/JProperties.hh"
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
namespace {
 
  /**
   * Wild card character for file name substition.
   */
  const char WILDCARD = '%';
 
  /**
   * Execute command in shell.
   *
   * \param  command       command
   */
  inline void execute(const std::string& command, int debug)
  {
    using namespace std;
    using namespace JPP;
 
    JProcess process(command);
 
    istream in(process.getInputStreamBuffer());
 
    for (string buffer; getline(in, buffer); ) {
      DEBUG(buffer << endl);
    }
  }
 
  const char* const histogram_t   =   "histogram";     //!< draw histogram
  const char* const arrow_t       =   "arrow";         //!< draw arrow
}
 
 
/**
 * \file
 *
 * Program to display hit probabilities.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
  
  typedef JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> > JParallelFileScanner_t;
  typedef JParallelFileScanner_t::multi_pointer_type               multi_pointer_type;
  typedef JTYPELIST<JCOMPASS::JOrientation,
                    JACOUSTICS::JEvt>::typelist                    calibration_types;
  typedef JMultipleFileScanner<calibration_types>                  JCalibration_t;
 
  JParallelFileScanner_t   inputFile;
  JLimit_t&                numberOfEvents = inputFile.getLimit();
  string                   detectorFile;
  JCalibration_t           calibrationFile;
  double                   Tmax_s;
  string                   pdfFile;
  string                   outputFile;
  JMuonGandalfParameters_t parameters;
  int                      application;
  JEventSelector           event_selector;
  JCanvas                  canvas;
  bool                     batch;
  struct : JStyle::JParameters {
    double                 arrowSize   =   0.003;
    string                 arrowType   =  "|->";
    double                 arrowScale  = 250.0;
    Width_t                lineWidth   =   2;
    Style_t                lineStyle   =   1;
    int                    nbinsX      =  50;
    int                    nbinsY      = 250;
    double                 T_ns        =   0.0;
  } graphics;
  string                   option;
  int                      debug;
 
 
  try { 
 
    JProperties properties = graphics.getProperties();
 
    properties.insert(gmake_property(graphics.arrowSize));
    properties.insert(gmake_property(graphics.arrowType));
    properties.insert(gmake_property(graphics.arrowScale));
    properties.insert(gmake_property(graphics.lineWidth));
    properties.insert(gmake_property(graphics.lineStyle));
    properties.insert(gmake_property(graphics.T_ns));
 
    parameters.numberOfPrefits = 1;
 
    JParser<> zap("Program to display hit probabilities.");
    
    zap['w'] = make_field(canvas,       "size of canvas <nx>x<ny> [pixels]")  = JCanvas(1200, 600);
    zap['f'] = make_field(inputFile,    "input file (output of JXXXMuonReconstruction.sh)");
    zap['a'] = make_field(detectorFile);
    zap['+'] = make_field(calibrationFile)     = JPARSER::initialised();
    zap['T'] = make_field(Tmax_s)              = 100.0;
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['P'] = make_field(pdfFile);
    zap['o'] = make_field(outputFile,   "graphics output file name")          = MAKE_STRING("display_" << WILDCARD << ".gif");
    zap['@'] = make_field(parameters)          = JPARSER::initialised();
    zap['A'] = make_field(application)         = JMUONGANDALF, JMUONENERGY, JMUONSTART;
    zap['L'] = make_field(event_selector)      = JPARSER::initialised();
    zap['%'] = make_field(properties)          = JPARSER::initialised();
    zap['O'] = make_field(option,       "draw option")                         = arrow_t, histogram_t;
    zap['B'] = make_field(batch,        "batch processing");
    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);
  }
  unique_ptr<JDynamics> dynamics;
 
  try {
 
    dynamics.reset(new JDynamics(detector, Tmax_s));
 
    dynamics->load(calibrationFile);
  }
  catch(const exception& error) {
    if (!calibrationFile.empty()) {
      FATAL(error.what());
    }
  }
 
  const JModuleRouter router(dynamics ? dynamics->getDetector() : detector);
 
  JSummaryFileRouter  summary(inputFile, parameters.R_Hz);
 
  const JMuonPDF_t    pdf(pdfFile, parameters.TTS_ns);
 
  const JTimeRange    T_ns(parameters.TMin_ns, parameters.TMax_ns);
 
  typedef vector<JHitL0>  JDataL0_t;
  typedef vector<JHitW0>  JDataW0_t;
 
  const JBuildL0<JHitL0>  buildL0;
 
 
  Vec offset(0.0, 0.0, 0.0);
 
  try {
    offset = getOffset(getHeader(inputFile));
  } catch(const exception& error) {}
 
 
  // 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(), graphics));
 
  gROOT->SetStyle("gplot");
  gROOT->ForceStyle();
 
  const size_t NUMBER_OF_PADS = 3;
 
  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(NUMBER_OF_PADS, 1);
 
  p1->Draw();
  p2->Draw();
 
  const double Dmax = getMaximalDistance(detector);
  const double Rmin = 0.0;
  const double Rmax = min(parameters.roadWidth_m, 0.4 * Dmax);
  const double Tmin = min(parameters.TMin_ns,  -10.0);
  const double Tmax = max(parameters.TMax_ns, +100.0);
  const double Amin = 0.002 * (Tmax - Tmin);                               // minimal arrow length [ns]
  const double Amax = 0.8   * (Tmax - Tmin);                               // maximal arrow length [ns]
  const double ymin = Tmin - (option == arrow_t ? 0.2 * Amax : 0.0);
  const double ymax = Tmax + (option == arrow_t ? 0.5 * Amax : 0.0);
 
  const string Xlabel[NUMBER_OF_PADS] = { "R [m]", "#phi [rad]",   "z [m]"  };
  const double Xmin  [NUMBER_OF_PADS] = {  Rmin,      -PI,      -0.3 * Dmax };
  const double Xmax  [NUMBER_OF_PADS] = {  Rmax,      +PI,      +0.3 * Dmax };
 
  double Xs[NUMBER_OF_PADS];
 
  for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
    Xs[i] = 0.003 * (Xmax[i] - Xmin[i]) * (0.5 * NUMBER_OF_PMTS);          // x-offset arrow as function of PMT number
  }
 
  TH2D   H2[NUMBER_OF_PADS];
  TGraph G2[NUMBER_OF_PADS];
 
  for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
 
    H2[i] = TH2D(MAKE_CSTRING("h" << i), NULL, graphics.nbinsX, Xmin[i] - Xs[i], Xmax[i] + Xs[i], graphics.nbinsY, ymin, ymax);
 
    H2[i].GetXaxis()->SetTitle(Xlabel[i].c_str());
    H2[i].GetYaxis()->SetTitle("#Deltat [ns]");
 
    H2[i].GetXaxis()->CenterTitle(true);
    H2[i].GetYaxis()->CenterTitle(true);
 
    H2[i].SetStats(kFALSE);
    
    G2[i].Set(2);
 
    G2[i].SetPoint(0, H2[i].GetXaxis()->GetXmin(), 0.0);
    G2[i].SetPoint(1, H2[i].GetXaxis()->GetXmax(), 0.0);
 
    p1->cd(i+1);
 
    H2[i].Draw("AXIS");
    G2[i].Draw("SAME");
  }
 
 
  for (JTreeScanner<Evt> mc(inputFile); inputFile.hasNext(); ) {
 
    cout << "event: " << setw(8) << inputFile.getCounter() << endl;
 
    multi_pointer_type ps = inputFile.next();
 
    JDAQEvent*  tev   = ps;
    JFIT::JEvt* in    = ps;
    Evt*        event = NULL;
    
    if (dynamics) {
      dynamics->update(*tev);
    }
 
    if (mc.getEntries() != 0) {
      event = mc.getEntry(tev->getCounter());                              // Monte Carlo true information
    }
 
    in->select(JHistory::is_event(application));
 
    if (!in->empty()) {
 
      sort(in->begin(), in->end(), qualitySorter);
 
      if (!event_selector(*tev, *in, event)) {
        continue;
      }
      
 
      JDataL0_t dataL0;
      
      buildL0(*tev, router, true, back_inserter(dataL0));
 
      summary.update(*tev);
 
      JFIT::JFit muon;                                                           // Monte Carlo true muon
 
      if (event != NULL) {
 
        const time_converter converter = time_converter(*event, *tev);
 
        for (const auto& t1 : event->mc_trks) {
          if (is_muon(t1)) {
            if (t1.E > muon.getE()) {
 
              JTrack3E ta = getTrack(t1);
 
              ta.add(getPosition(offset));
              ta.add(converter.putTime());
 
              muon = getFit(0, ta, 0.0, 0, t1.E, 1);
 
              muon.setW(JSTART_LENGTH_METRES, fabs(t1.len));
            }
          }
        }
      }
 
      bool   monte_carlo = false;                                          // show Monte Carlo true muon
      size_t index       = 0;                                              // index of fit
 
      for (bool next = false; !next; ) {
 
        for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
          H2[i].Reset();
        }
 
        JFIT::JFit fit;
 
        if (!monte_carlo) 
          fit = (*in)[index];
        else
          fit = muon;
 
        JRotation3D    R (getDirection(fit));
        JLine1Z        tz(getPosition (fit).rotate(R), fit.getT());
        JRange<double> Z_m;
        /*
        if (fit.hasW(JSTART_LENGTH_METRES)) {
          Z_m = JRange<double>(0.0, fit.getW(JSTART_LENGTH_METRES))  +  JRange<double>(parameters.ZMin_m, parameters.ZMax_m);
        }
        */
        const JFIT::JModel<JLine1Z> match(tz, parameters.roadWidth_m, T_ns, Z_m);
 
        // hit selection based on fit result
 
        JDataW0_t data;
          
        for (JDataL0_t::const_iterator i = dataL0.begin(); i != dataL0.end(); ++i) {
 
          JHitW0 hit(*i, summary.getRate(i->getPMTIdentifier()));
            
          hit.rotate(R);
 
          if (match(hit)) {
            data.push_back(hit);
          }
        }
          
        // select first hit in PMT
          
        sort(data.begin(), data.end(), JHitW0::compare);
 
        JDataW0_t::iterator __end = unique(data.begin(), data.end(), equal_to<JDAQPMTIdentifier>());      
 
        double E_GeV = parameters.E_GeV;
        /*
        if (fit.getE() > 0.1) {
          E_GeV = fit.getE();
        }
        */
 
        // move fit to geometrical center of hits
 
        JRange<double> zs = JRange<double>::DEFAULT_RANGE();
 
        for (JDataW0_t::iterator hit = data.begin(); hit != __end; ++hit) {
 
          const double x  = hit->getX() - tz.getX();
          const double y  = hit->getY() - tz.getY();
          const double z  = hit->getZ();
          const double R  = sqrt(x*x + y*y);
 
          zs.include(z - R/getTanThetaC());
        }
 
        const double   z0 = tz.getZ();
        const double   z1 = 0.5 * (zs.getLowerLimit() + zs.getUpperLimit());
 
        tz.setZ(z1, getSpeedOfLight());
 
        // graphics
 
        ostringstream   os;
        vector<TArrow>  arrow [NUMBER_OF_PADS];
        vector<TMarker> marker[NUMBER_OF_PADS];
 
        if (fit.hasW(JSTART_LENGTH_METRES) && fit.getW(JSTART_LENGTH_METRES) > 0.0) {
 
          marker[2].push_back(TMarker(z0 - tz.getZ(), 0.0, kFullCircle));
          marker[2].push_back(TMarker(z0 - tz.getZ() + fit.getW(JSTART_LENGTH_METRES), 0.0, kFullCircle));
 
          static_cast<TAttMarker&>(marker[2][0]) = TAttMarker(kRed, kFullCircle, 0.7);
          static_cast<TAttMarker&>(marker[2][1]) = TAttMarker(kRed, kFullCircle, 0.7);
        }
 
        DEBUG("trk: " 
              << FIXED(7,2)  << tz.getX() << ' '
              << FIXED(7,2)  << tz.getY() << ' '
              << FIXED(7,2)  << tz.getZ() << ' '
              << FIXED(12,2) << tz.getT() << endl);
 
        double chi2 = 0;
 
        for (JDataW0_t::const_iterator hit = data.begin(); hit != __end; ++hit) {
 
          const double x  = hit->getX() - tz.getX();
          const double y  = hit->getY() - tz.getY();
          const double z  = hit->getZ() - tz.getZ();
          const double R  = sqrt(x*x + y*y);
 
          const double t1 = tz.getT() + (z + R * getTanThetaC()) * getInverseSpeedOfLight();
 
          JDirection3D dir(hit->getDX(), hit->getDY(), hit->getDZ()); // PMT orientation
 
          dir.rotate(JRotation3Z(-atan2(y,x)));                       // rotate PMT axis to x-z plane
 
          const double theta = dir.getTheta();
          const double phi   = fabs(dir.getPhi());                    // rotational symmetry of Cherenkov cone
 
          //const double E  = gWater.getE(E_GeV, z);                  // correct for energy loss
          const double E  = E_GeV;
          const double dt = T_ns.constrain(hit->getT()  -  t1);
 
          JMuonPDF_t::result_type H1 = pdf.calculate(E, R, theta, phi, dt);
          JMuonPDF_t::result_type H0(hit->getR() * 1e-9, 0.0, T_ns);
 
          if (H1.V >= parameters.VMax_npe) {
            H1 *= parameters.VMax_npe / H1.V;
          }
 
          H1   += H0;                                                 // signal + background
 
          chi2 += H1.getChi2() - H0.getChi2();
 
          DEBUG("hit: " 
                << setw(8) << hit->getModuleID() << '.' << FILL(2,'0') << (int) hit->getPMTAddress() << FILL() << ' '
                << SCIENTIFIC(8,2) << E       << ' '
                << FIXED(7,2) << R            << ' '
                << FIXED(7,4) << theta        << ' '
                << FIXED(7,4) << phi          << ' '
                << FIXED(7,3) << dt           << ' '
                << FIXED(7,3) << H1.getChi2() << ' '
                << FIXED(7,3) << H0.getChi2() << endl);
 
          const double derivative = H1.getDerivativeOfChi2() - H0.getDerivativeOfChi2();
 
          double size = derivative * graphics.arrowScale;             // size of arrow
 
          if        (fabs(size) < Amin) { 
            size = (size > 0.0 ? +Amin : -Amin);
          } else if (fabs(size) > Amax) { 
            size = (size > 0.0 ? +Amax : -Amax);
          }
 
          const double X[NUMBER_OF_PADS] = { R, atan2(y,x), z - R/getTanThetaC() };
            
          const double xs = (double) (NUMBER_OF_PMTS - 2 * hit->getPMTAddress()) / (double) NUMBER_OF_PMTS;
 
          for (size_t i = 0; i != NUMBER_OF_PADS; ++i) {
 
            TArrow a1(X[i] + xs*Xs[i], dt + graphics.T_ns, X[i] + xs*Xs[i], dt + graphics.T_ns  + size, graphics.arrowSize, graphics.arrowType.c_str());
 
            a1.SetLineWidth(graphics.lineWidth);
            a1.SetLineStyle(graphics.lineStyle);
 
            arrow[i].push_back(a1);
 
            H2[i].Fill(X[i], dt + graphics.T_ns);
          }
        }
 
        os << FILL(6,'0') << tev->getRunNumber() << ":"  << tev->getFrameIndex() << "/" << tev->getCounter() << FILL();
        os << "  Q = "           << FIXED(4,0)      << fit.getQ();
        os << "  E = "           << SCIENTIFIC(7,1) << fit.getE() << " [GeV]";
        os << "  cos(#theta) = " << FIXED(6,3)      << fit.getDZ();
 
        if      (monte_carlo)
          os << "  Monte Carlo";
        else if (muon.getStatus() >= 0)
          os << "  #Delta#alpha = " << FIXED(6,2) << getAngle(getDirection(muon), getDirection(fit)) << " [deg]";
 
        
        // draw
 
        TLatex title(0.05, 0.5, os.str().c_str());
 
        title.SetTextAlign(12);
        title.SetTextFont(42);
        title.SetTextSize(0.6);
 
        p2->cd();
 
        title.Draw();
 
        for (int i = 0; i != NUMBER_OF_PADS; ++i) {
 
          p1->cd(i+1);
 
          if (option == arrow_t) {
 
            for (auto& a1 : arrow[i]) {
              a1.Draw();
            }
 
            for (auto& m1 : marker[i]) {
              m1.Draw();
            }
          }
 
          if (option == histogram_t) {
            H2[i].Draw("SAME");
          }
        }
      
        cv->Update();
 
 
        // action
 
        if (batch) {
 
          cv->SaveAs(replace(outputFile, WILDCARD, MAKE_STRING(inputFile.getCounter())).c_str());
 
          next = true;
      
        } 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 << '+' << " -> " << "next fit"                                    << endl;
              cout << '-' << " -> " << "previous fit"                                << endl;
              cout << 'M' << " -> " << "Monte Carlo true muon information"           << endl;
              cout << 'F' << " -> " << "fit information"                             << endl;
              if (event_selector.is_valid()) {
                cout << 'L' << " -> " << "reload event selector"                       << endl;
              }
              cout << 'r' << " -> " << "rewind input file"                           << endl;
              cout << 'R' << " -> " << "switch to ROOT mode (quit ROOT to continue)" << 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(inputFile.getCounter())).c_str());
              break;
 
            case '+':
              monte_carlo = false;
              index = (index != in->size() - 1 ? index + 1 : 0);
              user = false;
              break;
 
            case '-':
              monte_carlo = false;
              index = (index != 0 ? index - 1 : in->size() - 1);
              user = false;
              break;
 
            case 'M':
              if (muon.getStatus() >= 0)
                monte_carlo = true;
              else
                ERROR(endl << "No Monte Carlo muon available." << endl);
              user = false;
              break;
 
            case 'F':
              monte_carlo = false;
              user = false;
              break;
 
            case 'L':
              if (event_selector.is_valid()) {
                execute(MAKE_STRING("make -f " << getPath(argv[0]) << "/JMakeEventSelector libs"), 3);
                event_selector.reload();
              }
              break;
 
            case 'R':
              tp->Run(kTRUE);
              break;
 
            case 'r':
              inputFile.rewind();
          
            default:
              next = true;
              user = false;
              break;
            }
          }
        }
      }
    }
  }
  cout << endl;
}