JMuonPostfit.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TProfile.h"
#include "TMath.h"
 
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/online/JDAQModuleIdentifier.hh"
#include "km3net-dataformat/definitions/reconstruction.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"
#include "km3net-dataformat/tools/time_converter.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JVolume.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JPhysics/JConstants.hh"
#include "JTools/JQuantile.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JLang/JPredicate.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
namespace {
 
  using KM3NETDAQ::JDAQEvent;
 
 
  /**
   * Count number of unique identifiers in event.
   *
   * \param  __begin       begin of data
   * \param  __end         end   of data
   * \return               number of unique identifiers
   */
  template<class JHit_t>
  inline int getCount(JDAQEvent::const_iterator<JHit_t> __begin, 
                      JDAQEvent::const_iterator<JHit_t> __end)
  {
    using namespace std;
    using namespace KM3NETDAQ;
 
    typedef JDAQModuleIdentifier  JType_t;
 
    set<JType_t> buffer;
 
    for (JDAQEvent::const_iterator<JHit_t> i = __begin; i != __end; ++i) {
      buffer.insert(JType_t(*i));
    }
 
    return buffer.size();
  }
 
 
  const char* const muon_t      =  "muon";
  const char* const neutrino_t  =  "neutrino";
}
 
/**
 * \file
 *
 * Example program to histogram fit results for the muon reconstruction chain.
 * \author mdejong, bofearraigh.
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  typedef JTriggeredFileScanner<JEvt>                       JTriggeredFileScanner_t;
  typedef JTriggeredFileScanner_t::multi_pointer_type       multi_pointer_type;
 
  JTriggeredFileScanner_t inputFile;
  JLimit_t&               numberOfEvents = inputFile.getLimit();
  string                  outputFile;  
  size_t                  numberOfPrefits;
  JQualitySorter          quality_sorter;
  JAtmosphericMuon        atmosphere;
  double                  Emin_GeV;
  int                     application;
  string                  option;
  string                  primary;
  int                     debug;
 
  try {
 
    JParser<> zap("Example program to histogram fit results.");
 
    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "postfit.root";
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['N'] = make_field(numberOfPrefits)     = 1;
    zap['L'] = make_field(quality_sorter)      = JPARSER::initialised();
    zap['E'] = make_field(Emin_GeV)            = 0.0;
    zap['A'] = make_field(application)         = JMUONPREFIT, JMUONSIMPLEX, JMUONGANDALF, JMUONENERGY, JMUONSTART, JMUONPATH; 
    zap['a'] = make_field(atmosphere)          = JAtmosphericMuon(90.0, 90.0);
    zap['O'] = make_field(option)              = "LINE", "LOGE", "LINN", "LOGN";
    zap['p'] = make_field(primary)             =  muon_t, neutrino_t;
    zap['d'] = make_field(debug)               =  2;
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
  JHead head;
 
  try {
    head = getCommonHeader(inputFile);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
  JVolume     volume(head, option != "LINE");
  JPosition3D offset   = getPosition(getOffset(head));
  JPosition3D origin   = getPosition(getOrigin(head));
  JCylinder3D cylinder = getCylinder(head);
 
  cylinder.add(offset);
 
  NOTICE("Offset: "   << offset   << endl);
  NOTICE("Cylinder: " << cylinder << endl);
 
  TFile out(outputFile.c_str(), "recreate");
 
  TH1D job("job", NULL, 100, 0.5, 100.5);
 
  TH1D hn("hn",  NULL,  250,   -0.5, 249.5);    // NDF
  TH1D hq("hq",  NULL,  300,    0.0, 600.0);    // quality parameter
  TH1D hi("hi",  NULL,  101,   -0.5, 100.5);    // index of the best event in order of descending angular resolution
  TH1D hc("hc",  NULL,  200,   -1.0,  +1.0);    // dz (z-slope)
  TH1D hu("hu",  NULL,  400, -1.0e3, 1.0e3);    // quality difference between best up-going and down-going track (check of atmospheric muon)
 
  TH1D hx("hx",  NULL,  70,   -3.0,  +2.3);     // angular deviation [log(deg)]
  TH1D hd("hd",  NULL,  100,    0.0,  10.0);    // distance between true and best event at intersection point [m]
  TH1D ht("ht",  NULL,  100, -100.0, 100.0);    // time between true and best track [ns]
 
  TH1D hz0("hz0[start]",  NULL,  400, -200.0, 200.0);    // start point [m]
  TH1D hz1("hz1[end]",    NULL,  400, -200.0, 200.0);    // end   point [m]
 
  TH1D E_0("E_0",  NULL,  100, volume.getXmin(), volume.getXmax());     // true muon energy
  TH1D E_1("E_1",  NULL,  100, volume.getXmin(), volume.getXmax());     // uncorrected energy of best track
  TH1D E_2("E_2",  NULL,  100, volume.getXmin(), volume.getXmax());     // corrected   energy of best track
  TH1D E_E("E_E",  NULL,  100, -5.0, +5.0);                             // ratio of reconstructed energy and true energy [log(E/E)]
  TH2D ExE("ExE",  NULL,
           40, volume.getXmin(), volume.getXmax(),
           40, volume.getXmin(), volume.getXmax());  // (Etrue,Ereco)
 
  const Int_t    ny   = 28800;
  const Double_t ymin =   0.0;         // [deg]
  const Double_t ymax =   180.0;       // [deg]
 
  vector<double> X;
 
  if        (option == "LINN") {
 
    const double xmin = log((double)   3);
    const double xmax = log((double) 300);
 
    for (double x = xmin, dx = (xmax - xmin) / 20; x <= xmax; x += dx) {
      X.push_back((int) exp(x));
    }
 
  } else if (option == "LOGN") {
 
    const double xmin = log10((double)   3);
    const double xmax = log10((double) 400);
 
    for (double x = xmin, dx = (xmax - xmin) / 20; x <= xmax; x += dx) {
      X.push_back(x);
    }
 
  } else {
 
    for (double x = volume.getXmin(), dx = (volume.getXmax() - volume.getXmin()) / 20.0; x < volume.getXmax() + 0.5*dx; x += dx) {
      X.push_back(x);
    }
  }
 
  TH2D h2("h2", NULL, X.size() - 1, X.data(), ny, ymin, ymax);
 
  TH2D Va("Va", NULL, 
          100, 0.0, cylinder.getRadius()*cylinder.getRadius(), 
          100, cylinder.getZmin() - 50.0, cylinder.getZmax() + 50.0);
  TH2D Vb("Vb", NULL, 
          100, 0.0, cylinder.getRadius()*cylinder.getRadius(),
          100, cylinder.getZmin() - 50.0, cylinder.getZmax() + 50.0);
 
  JQuantile Q("Angle");
  JQuantile O("Omega");
 
 
  while (inputFile.hasNext()) {
 
    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);
 
    multi_pointer_type ps = inputFile.next();
 
    JDAQEvent* tev   = ps;
    JEvt*      evt   = ps;
    Evt*       event = ps;
 
    const time_converter converter(*event, *tev);
 
    job.Fill(1.0);
 
 
    double Enu  = 0.0;
    double Emu  = 0.0;
    double Emax = 0.0;
 
    if (has_neutrino(*event)) {
      Enu = get_neutrino(*event).E;
    }
 
    vector<Trk>::const_iterator muon = event->mc_trks.end();
 
    for (vector<Trk>::const_iterator track = event->mc_trks.begin(); track != event->mc_trks.end(); ++track) {
 
      if (is_muon(*track)) {
 
        Emu += track->E;
 
        if (track->E > Emax) {
          muon = track;
          Emax = track->E;
        }
      }
    }
 
    if (muon == event->mc_trks.end()) {
      continue;
    }
 
    job.Fill(3.0);
 
 
    // abscissa
 
    Double_t x = 0.0;
 
    {
      const double E = event->mc_trks[0].E;
      const int    N = getCount(tev->begin<JDAQTriggeredHit>(), tev->end<JDAQTriggeredHit>());
    
      if      (option == "LINN")
        x = N;
      else if (option == "LOGN")
        x = log10(N);
      else
        x = volume.getX(E);
    }
 
 
    if (!evt->empty()) {
 
      JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_application(application));
 
      if (evt->begin() == __end) {
        continue;
      }
 
      job.Fill(4.0);
 
      if (numberOfPrefits > 0) {
 
        JEvt::iterator __q = __end;
 
        advance(__end = evt->begin(), min(numberOfPrefits, (size_t) distance(evt->begin(), __q)));
 
        partial_sort(evt->begin(), __end, __q, quality_sorter);
 
      } else {
 
        sort(evt->begin(), __end, quality_sorter);
      }
 
      double    E = 0.0;
      JPointing pointing;
 
      if        (primary == muon_t) {
        E        = Emu;
        pointing = JPointing(getDirection(*muon));
      } else if (primary == neutrino_t) {       
        E        = Enu;
        pointing = JPointing(getDirection(get_neutrino(*event)));
      }
 
      JEvt::iterator best  = pointing(evt->begin(), __end);
      const Double_t beta  = pointing.getAngle(*best);
      const double   Efit  = best->getE();
      const double   Eraw  = best->getW(JENERGY_ENERGY, numeric_limits<double>::min());
 
      // selection of fit result
 
      bool ok = (Efit >= Emin_GeV);
 
      if (ok) {
 
        job.Fill(5.0);
 
        hn.Fill((Double_t) best->getNDF());
        hq.Fill(best->getQ());
        hi.Fill((Double_t) distance(evt->begin(), best));
        hc.Fill(best->getDZ());
 
        // atmospheric muon hypothesis test. 
 
        if (( has_neutrino(*event) && get_neutrino(*event).dir.z >= atmosphere.dot2) || // difference in quality between best upwards and best downgoing track. 
            (!has_neutrino(*event) && best->getDZ()              >= atmosphere.dot2)) { // negative (positive) values imply a down(up) going track.
          hu.Fill(atmosphere(evt->begin(), __end));
        }
 
        hx.Fill(max(log10(beta), hx.GetXaxis()->GetXmin()));
 
        Q.put(beta);
 
        JTrack3E ta = getTrack(*muon);
        JTrack3E tb = getTrack(*best);
 
        ta.add(offset);
        tb.sub(converter.putTime());
 
        static_cast<JTrack3D&>(ta).move(ta.getIntersection(tb), getSpeedOfLight());     // move to intersection point of both tracks
        static_cast<JTrack3D&>(tb).move(tb.getIntersection(ta), getSpeedOfLight());
 
        hd.Fill((tb.getPosition() - ta.getPosition()).getLength());
 
        if (has_neutrino(*event)) {
 
          job.Fill(6.0);
 
          const Trk neutrino = get_neutrino(*event);
 
          JPosition3D vertex = getPosition(neutrino);
 
          vertex.add(offset);
 
          if (cylinder.is_inside(vertex)) {
 
            job.Fill(7.0);
 
            Va.Fill(JVector2D(best->getX() - origin.getX(), best->getY() - origin.getY()).getLengthSquared(), best->getZ());   //R^{2} [m^2], z [m]
 
            JTrack3E tc = getTrack(*best);
 
            hz0.Fill(tc.getIntersection(vertex));
            hz1.Fill(best->getW(JSTART_LENGTH_METRES, 0.0) - gWater(muon->E));
          }
        }
 
        Vb.Fill(JVector2D(best->getX() - origin.getX(), best->getY() - origin.getY()).getLengthSquared(), best->getZ());   //R^{2} [m^2], z [m]
        ht.Fill(tb.getT() - ta.getT());
 
        E_0.Fill(volume.getX(E,    true));
        E_1.Fill(volume.getX(Eraw, true));
        E_2.Fill(volume.getX(Efit, true));
        E_E.Fill(volume.getX(Efit) - volume.getX(E));
        ExE.Fill(volume.getX(E),     volume.getX(Efit));
 
        h2.Fill(x, beta);
      }
    }
  }
  STATUS(endl);
 
  NOTICE("Number of events input           " << setw(8) << right << job.GetBinContent(1) << endl);
  NOTICE("Number of events with muon       " << setw(8) << right << job.GetBinContent(3) << endl);
  NOTICE("Number of events with fit        " << setw(8) << right << job.GetBinContent(4) << endl);
  NOTICE("Number of events selected        " << setw(8) << right << job.GetBinContent(5) << endl);
  NOTICE("Number of events with neutrino   " << setw(8) << right << job.GetBinContent(6) << endl);
  NOTICE("Number of events contained       " << setw(8) << right << job.GetBinContent(7) << endl);
 
  if (Q.getCount() != 0) {
 
    NOTICE("Median space angle [deg] " << FIXED   (6,3) << Q.getQuantile(0.5) << endl);
 
    for (double q : {0.50, 0.90, 0.99}) {
      NOTICE("Space angle " << FIXED(5,1) << (100.0 * q) << "% quantile [deg] " << FIXED(6,3) << Q.getQuantile(q) << endl);
    }
  }
 
  out.Write();
  out.Close();
}