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;
  double                  NPE;
  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['M'] = make_field(NPE)                 = 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)              = "E", "N", "LINE", "LOGE";
    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 = getHeader(inputFile);
  }
  catch(const JException& error) {
    FATAL(error);
  }

  const JVolume     volume(head, option != "LINE");
  const JPosition3D center   = get<JPosition3D>(head);
  JCylinder3D       cylinder = get<JCylinder3D>(head);

  cylinder.add(center);

  NOTICE("Center         [m]: " << center   << endl);
  NOTICE("Reposition can [m]: " << 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 hv("hv",  NULL,  200,   -6.0,   0.0);    // log of Poisson distribution of (#hits,number photo-electrons)
  TH1D h1("h1",  NULL,  200,   -2.0,  +2.0);    // number of photo-electrons along the whole track
  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.find('E') != string::npos) {

    for (double x = volume.getXmin(); x <= volume.getXmax(); x += (volume.getXmax() - volume.getXmin()) / 20) {
      X.push_back(x);
    }

  } else {

    double x = -0.5;

    for ( ; x <=  15.5; x +=  1.0) { X.push_back(x); }
    for ( ; x <=  25.5; x +=  2.0) { X.push_back(x); }
    for ( ; x <=  50.5; x +=  5.0) { X.push_back(x); }
    for ( ; x <= 100.5; x += 10.0) { X.push_back(x); }
    for ( ; x <= 250.5; x += 20.0) { 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", true);
  JQuantile O("Omega", true);


  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;

    if (option.find('E') != string::npos)
      x = volume.getX(event->mc_trks[0].E);
    else
      x = getCount(tev->begin<JDAQTriggeredHit>(), tev->end<JDAQTriggeredHit>());

    if (!evt->empty()) {

      JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_event(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());
      const double   mip   = best->getW(JSTART_NPE_MIP, numeric_limits<double>::max());

      // selection of fit result

      bool ok = (Efit >= Emin_GeV  && mip  >= NPE);

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

          if (cylinder.is_inside(vertex)) {

            job.Fill(7.0);

            Va.Fill(best->getX()*best->getX() + best->getY()*best->getY(), 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(best->getX()*best->getX() + best->getY()*best->getY(), 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);

        if (best->hasW(JSTART_NPE_MIP)) { 
          h1.Fill(log10(best->getW(JSTART_NPE_MIP)));
        }

        if (best->hasW(JVETO_NPE) && best->hasW(JVETO_NUMBER_OF_HITS)) {

          const double npe   = best->getW(JVETO_NPE);               // number of photo-electrons 
          const int    count = best->getW(JVETO_NUMBER_OF_HITS);    // number of hits
          const double pv    = TMath::PoissonI(count, npe);         // Poisson distribution function for (#hits,npe)

          hv.Fill(max(log10(pv), hv.GetXaxis()->GetXmin()));
        }
      }
    }
  }
  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);
  }

  out.Write();
  out.Close();
}