JAAnet/JEvD.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>

#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/online/JDAQHeader.hh"
#include "km3net-dataformat/tools/time_converter.hh"

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

#include "JDAQ/JDAQSummarysliceIO.hh"

#include "JTrigger/JHitL0.hh"

#include "JSupport/JSupport.hh"
#include "JSupport/JSingleFileScanner.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/JMuonGandalfParameters_t.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/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


namespace {

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


  using JEEP::JFunctionAdaptor;


  /**
   * Event selector.
   *
   * The default constructor will accept all events.\n
   * A different method can dynamically be loaded from a (shared) library using class JEEP::JFunctionAdaptor.
   * For the definition of an alternative method, see e.g.\ event_selector.cc 
   */
  struct JEventSelector :
    public JFunctionAdaptor<bool, const Trk&, const Evt&>
  {
    typedef JFunctionAdaptor<bool, const Trk&, const Evt&>  function_adaptor_type;

    /**
     * Default event selection.
     *
     * \param  trk           track
     * \param  evt           event
     * \return               true
     */
    static inline bool accept(const Trk& trk, const Evt& evt)
    {
      return true;
    }

    
    /**
     * Default constructor.
     */
    JEventSelector() :
      JFunctionAdaptor<bool, const Trk&, const Evt&>(accept)
    {}


    /**
     * Check validity of function.
     *
     * \return                     true if valid function; else false
     */
    bool is_valid() const
    {
      return (libso != "" && function_adaptor_type::is_valid());
    }


    /**
     * Re-load function from shared library.
     */
    void reload()
    {
      reset();

      load(libso);
    }


    /**
     * Read event selector from input stream.
     *
     * \param  in                  input stream
     * \param  object              event selector
     * \return                     input stream
     */
    friend std::istream& operator>>(std::istream& in, JEventSelector& object)
    {
      using namespace std;

      object.reset();

      if (in >> object.libso) {

        if (object.libso.find(function_adaptor_type::SEPARATOR) == string::npos) {
          object.libso += MAKE_CSTRING(function_adaptor_type::SEPARATOR << "accept");
        }

        object.load(object.libso);
      }

      return in;
    }

  private:
    std::string libso;      //!< shared library : symbol
  };


  /**
   * 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);
    }
  }


  /**
   * Check availability of value.
   *
   * \param  trk      track
   * \param  i        index
   * \return          true if available; else false
   */
  inline bool hasW(const Trk& trk, const int i)
  {
    return (i >= 0 && i < (int) trk.fitinf.size());
  }


  /**
   * Get associated value.
   *
   * \param  trk      track
   * \param  i        index
   * \return          value
   */
  inline double getW(const Trk& trk, const int i)
  {
    return trk.fitinf.at(i);
  }


  /**
   * Get associated value.
   *
   * \param  trk      track
   * \param  i        index
   * \param  value    default value
   * \return          value
   */
  inline double getW(const Trk& trk, const int i, const double value)
  {
    if (hasW(trk,i))
      return trk.fitinf.at(i);
    else
      return value;
  }


  /**
   * Set associated value.
   *
   * \param  trk      track
   * \param  i        index
   * \param  value    value
   */
  void setW(Trk& trk, const int i, const double value)
  {
    if (i >= (int) trk.fitinf.size()) {
      trk.fitinf.resize(i + 1, 0.0);
    }

    trk.fitinf[i] = value;
  }
}


/**
 * \file
 *
 * Program to display hit probabilities.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
  
  JSingleFileScanner<Evt>  inputFile;
  JLimit_t&                numberOfEvents = inputFile.getLimit();
  string                   pdfFile;
  string                   outputFile;
  JMuonGandalfParameters_t parameters;
  int                      application;
  JEventSelector           event_selector;
  JCanvas                  canvas;
  bool                     batch;
  double                   arrowSize;
  string                   arrowType;
  double                   arrowScale;
  int                      debug;


  try { 

    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['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)         = JMUONENERGY, JMUONSTART, JMUONGANDALF, JMUONSIMPLEX, JMUONPREFIT;
    zap['L'] = make_field(event_selector)      = JPARSER::initialised();
    zap['S'] = make_field(arrowSize)           = 0.003;
    zap['T'] = make_field(arrowType)           = "|->";
    zap['F'] = make_field(arrowScale)          = 250.0;
    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);
  }


  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;


  // ROOT

  gROOT->SetBatch(batch);

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

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

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

  const size_t NUMBER_OF_PADS = 3;

  cv->SetFillStyle(4000);
  cv->SetFillColor(kWhite);
  cv->Divide(NUMBER_OF_PADS, 1);

  const double Dmax = 1000.0;
  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 = min(-Amax, Tmin - 0.3 * Amax);
  const double ymax = max(+Amax, Tmax + 0.5 * Amax);

  const string Xlabel[NUMBER_OF_PADS] = { "R [m]", "#phi [rad]",   "z [m]"  };
  const double Xmin  [NUMBER_OF_PADS] = {  Rmin,      -PI,      -0.5 * Dmax };
  const double Xmax  [NUMBER_OF_PADS] = {  Rmax,      +PI,      +0.5 * 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, 100, Xmin[i] - Xs[i], Xmax[i] + Xs[i], 100, 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);

    cv->cd(i+1);

    H2[i].Draw("AXIS");
    G2[i].Draw("SAME");
  }


  while (inputFile.hasNext()) {

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

    Evt* evt = inputFile.next();

    if (has_reconstructed_track<JPP_RECONSTRUCTION_TYPE>(*evt, rec_stages_range(application))) {

      Trk fit = get_best_reconstructed_track<JPP_RECONSTRUCTION_TYPE>(*evt, rec_stages_range(application));

      if (!event_selector(fit, *evt)) {
        continue;
      }
      

      JDataL0_t dataL0;
      
      for (const Hit& hit : evt->hits) {
        dataL0.push_back(JHitL0(JDAQPMTIdentifier(hit.dom_id, hit.channel_id),
                                JAxis3D(getPosition(hit.pos), getDirection(hit.dir)),
                                JHit(hit.t, hit.tot)));
      }

      summary.update(JDAQChronometer(evt->det_id,
                                     evt->run_id,
                                     evt->frame_index,
                                     JDAQUTCExtended(evt->t.GetSec(), evt->t.GetNanoSec() / 16))); 

      const time_converter converter = time_converter(*evt);

      Trk muon;                                                            // Monte Carlo true muon

      if (has_muon(*evt)) {

        for (const auto& trk : evt->mc_trks) {
          if (is_muon(trk)) {
            if (trk.E > muon.E) {

              muon    = trk;
              muon.t += converter.putTime();

              setW(muon, JSTART_LENGTH_METRES, fabs(muon.len));
            }
          }
        }
      }


      bool   next_event  = false;                                          // goto next event
      bool   monte_carlo = false;                                          // show Monte Carlo true muon

      while (!next_event) {

        Trk trk;

        if (!monte_carlo) 
          trk = fit;
        else
          trk = muon;

        JRotation3D    R (getDirection(trk));
        JLine1Z        tz(getPosition (trk).rotate(R), trk.t);
        JRange<double> Z_m;
        /*
        if (hasW(trk, JSTART_LENGTH_METRES)) {
          Z_m = JRange<double>(0.0, getW(fit,JSTART_LENGTH_METRES))  +  JRange<double>(parameters.ZMin_m, parameters.ZMax_m);
        }
        */
        const 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 (trk.E > 0.1) {
          E_GeV = trk.E;
        }
        */

        // move fit to geometrical center of hits

        JRange<double> zs(make_array(data.begin(), __end, &JHitW0::getZ));

        const double   z0 = tz.getZ();
        const double   z1 = 0.5 * (zs.getLowerLimit() + zs.getUpperLimit());

        tz.setZ(z1, getSpeedOfLight());


        // graphics

        TLatex          latex [NUMBER_OF_PADS];
        vector<TArrow>  arrow [NUMBER_OF_PADS];
        vector<TMarker> marker[NUMBER_OF_PADS];

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

          latex[i].SetTextAlign(12);
          latex[i].SetTextFont(42);
          latex[i].SetTextSize(0.06);

          latex[i].SetX(H2[i].GetXaxis()->GetXmin()  +  0.05 * (H2[i].GetXaxis()->GetXmax() - H2[i].GetXaxis()->GetXmin()));
          latex[i].SetY(H2[i].GetYaxis()->GetXmax()  -  0.05 * (H2[i].GetYaxis()->GetXmax() - H2[i].GetYaxis()->GetXmin()));
        }

        if (hasW(trk, JSTART_LENGTH_METRES) && getW(trk, JSTART_LENGTH_METRES) > 0.0) {

          marker[2].push_back(TMarker(z0 - tz.getZ(), 0.0, kFullCircle));
          marker[2].push_back(TMarker(z0 - tz.getZ() + getW(trk, 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);
        }

        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();

          const double derivative = H1.getDerivativeOfChi2() - H0.getDerivativeOfChi2();

          double size = derivative * 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) {
            arrow[i].push_back(TArrow(X[i] + xs*Xs[i], dt, X[i] + xs*Xs[i], dt + size, arrowSize, arrowType.c_str()));
          }
        }

        latex[0].SetTitle(MAKE_CSTRING(""   << FILL(6,'0') << evt->run_id << FILL() <<
                                       ":"  << evt->frame_index                     << 
                                       "  " << evt->trigger_counter));

        if (!monte_carlo)
          latex[1].SetTitle(MAKE_CSTRING("Q = " << FIXED(4,0)      << -chi2 << "  " <<
                                         "E = " << SCIENTIFIC(7,1) << trk.E << " [GeV]"));

        else
          latex[1].SetTitle(MAKE_CSTRING("Q = " << FIXED(4,0)      << -chi2));

        if (is_muon(muon)) {

          if (!monte_carlo)
            latex[2].SetTitle(MAKE_CSTRING("#Delta#alpha = " << FIXED(6,2) << getAngle(getDirection(muon), getDirection(trk)) << " [deg]" <<
                                           " z = " << FIXED(6,3) << trk.dir.z));
          else
            latex[2].SetTitle("Monte Carlo");
        }


        // draw

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

          cv->cd(i+1);

          latex[i].Draw();

          for (auto& a1 : arrow[i]) {
            a1.Draw();
          }

          for (auto& m1 : marker[i]) {
            m1.Draw();
          }
        }
      
        cv->Update();


        // action

        if (batch) {

          cv->SaveAs(replace(outputFile, WILDCARD, MAKE_STRING(inputFile.getCounter())).c_str());

          next_event = true;
      
        } else {

          for (int c = 0; c == 0; ) {
        
            if (!batch) {

              static bool first = true;

              if (first) {
                cout << endl << "Type '?' for possible options." << endl;
              }

              first = false;
            }

            cout << "\n> " << flush;

            switch((c = 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 << '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 << 'p' << " -> " << "print event information"                     << endl;
              cout << ' ' << " -> " << "next event (as well as any other key)"       << endl;
              c = 0;
              break;

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

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

            case 's':
              cv->SaveAs(replace(outputFile, WILDCARD, MAKE_STRING(inputFile.getCounter())).c_str());
              c = 0;
              break;

            case 'M':
              if (is_muon(muon))
                monte_carlo = true;
              else
                ERROR(endl << "No Monte Carlo muon available." << endl);
              break;

            case 'F':
              monte_carlo = false;
              break;

            case 'r':
              inputFile.rewind();
              next_event  = true;
              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);
              c = 0;
              break;

            case 'p':
              cout << endl;
              evt->print(cout);
              if (debug >= debug_t) {
                for (const auto& trk : evt->mc_trks) {
                  cout << "MC  "; trk.print(cout); cout << endl;
                }
                for (const auto& trk : evt->trks) {
                  cout << "fit "; trk.print(cout); cout << endl;
                }
                for (const auto& hit : evt->hits) {
                  cout << "hit "; hit.print(cout); cout << endl;
                }
              }
              c = 0;
              break;
          
            default:
              next_event  = true;
              break;
            }
          }
        }
      }
    }
  }
  cout << endl;
}