JMuonPostfit.cc File Reference

#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 "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JTools/JConstants.hh"
#include "JTools/JQuantile.hh"
#include "JTrigger/JTimeConverter.hh"
#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"
#include "JFit/JEvt.hh"
#include "JFit/JEvtToolkit.hh"
#include "JFit/JFitParameters.hh"
#include "JGizmo/JVolume.hh"
#include "JLang/JPredicate.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Example program to histogram fit results for the muon reconstruction chain.

Author

mdejong, bofearraigh.

Definition in file JMuonPostfit.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 73 of file JMuonPostfit.cc.

{
  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;
  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['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("Reposition can [m]: " << cylinder << endl);
 
  const double EMIN_GEV = 10;  // set to low value for ORCA..
  const double E_nu_min = head.cut_nu.Emin; //neutrino minimum energy 
  const double E_nu_max = head.cut_nu.Emax; //neutrino maximum energy
 
  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 hz("hz",  NULL,  100, -200.0, 200.0);    // intersection of track with neutrino vertex [m]
  TH1D ht("ht",  NULL,  100, -100.0, 100.0);    // time between true and best track [ns]
 
  TH1D e0("e0",  NULL,  100, volume.getXmin(), volume.getXmax());     // true track energy [log(E)]
  TH1D e1("e1",  NULL,  100, volume.getXmin(), volume.getXmax());     // uncorrected energy [log(E)]
  TH1D e2("e2",  NULL,  100, volume.getXmin(), volume.getXmax());     // energy of best track[log(E)]
  TH1D er("er",  NULL,  100, -5.0, +5.0);                             // ratio of reconstructed energy and true energy [log(E/E)]
  TH2D ee("ee",  NULL,
      40, volume.getXmin(), volume.getXmax(),
      40, volume.getXmin(), volume.getXmax());  //(Etrue,Ereco)
 
  TH2D muon_angle_Emu("muon_angle_Emu", NULL, 30, E_nu_min, E_nu_max, 100, 0,90);         // energy vs angular deviation between true muon and reco muon direction
  TH2D muon_angle_Enu("muon_angle_Enu", NULL, 30, E_nu_min, E_nu_max, 100, 0,90);         // energy vs angular deviation between true muon and reco muon direction
  TH2D neutrino_angle_Enu("neutrino_angle_Enu", NULL, 30, E_nu_min, E_nu_max, 100, 0,90); // energy vs angular deviation between neutrino and reco muon direction
  TH2D MC_angle_Enu("MC_angle_Enu", NULL, 30, E_nu_min, E_nu_max, 100, 0,90);             // energy vs angular deviation between neutrino and true muon direction
 
  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);
  TProfile he("he", NULL, X.size() - 1, X.data());
 
  TH2D   ha("ha", NULL, 100, 0.0, 0.5e6, 100, -100.0, +900.0);
  TH2D   hb("hb", NULL, 100, 0.0, 0.5e6, 100, -100.0, +900.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 JTimeConverter converter(*event, *tev);
 
    job.Fill(1.0);
 
 
    double Emu  = 0.0;
    double Emax = 0.0;
 
    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>());
 
 
    // weight for efficiency determination
 
    Double_t W = 0.0;
 
    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);
      }
 
      const JPointing pointing(getDirection(*muon));
 
      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) {
 
        W = 1.0;
 
        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());
 
        if (( has_neutrino(*event) && get_neutrino(*event).dir.z >= atmosphere.dot2) || //Atmospheric muon hypothesis test. Difference in quality between best upwards and best downgoing track. 
            (!has_neutrino(*event) && best->getDZ()        >= atmosphere.dot2)) { //negative values imply a down going track, positive an up-going track.
          hu.Fill(atmosphere(evt->begin(), __end));
        }
 
        hx.Fill(max(log10(beta), hx.GetXaxis()->GetXmin()));
 
        Q.put(beta);
 
        const Trk neutrino = get_neutrino(*event);
        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);
 
          JPosition3D vertex = getPosition(neutrino);
 
          vertex.add(center);
 
          if (cylinder.is_inside(vertex)) {
 
            job.Fill(7.0);
 
            JTrack3E tc = getTrack(*best);
 
            ha.Fill(best->getX()*best->getX() + best->getY()*best->getY(), best->getZ());   //R^{2} [m^2],
            hz.Fill(tc.getIntersection(vertex));
          }
        }
 
        if (best->getE() >= EMIN_GEV) {
          hb.Fill(best->getX()*best->getX() + best->getY()*best->getY(), best->getZ());  //R^{2} [m^2], z [m]
        }
        ht.Fill(tb.getT() - ta.getT());
 
        e0.Fill(volume.getX(Emu,  true));
        e1.Fill(volume.getX(Eraw, true));
        e2.Fill(volume.getX(Efit, true));
        er.Fill(volume.getX(Efit) - volume.getX(Emu));
        ee.Fill(volume.getX(Emu), volume.getX(Efit));
        h2.Fill(x, beta);
 
        const double Enu = neutrino.E;  //true neutrino E
 
        neutrino_angle_Enu.Fill(Enu, (180/M_PI)* acos(getDirection(*best).getDot(getDirection(neutrino))));  
        muon_angle_Enu.Fill(Enu,(180/M_PI)* acos(getDirection(*best).getDot(getDirection(*muon))));
        muon_angle_Emu.Fill(Emu,(180/M_PI)* acos(getDirection(*best).getDot(getDirection(*muon))));
        MC_angle_Enu.Fill(Enu, (180/M_PI)* acos(getDirection(*muon).getDot(getDirection(neutrino))));
 
 
        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 photoelectrons 
          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()));
        }
      }
    }
 
    he.Fill(x, W);  //efficiency
  }
  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();
}