JShowerPostfit.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 "JPhysics/JGeanz.hh"
#include "JGeometry3D/JShower3E.hh"
#include "JGeometry3D/JCylinder3D.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();
  }
}


/**
 * \file
 *
 * Example program to histogram shower fit results: it handles both neutrino and muon productions.
 *
 * \author mdejong, adomi
 */
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;
  bool                    isMuon;
  bool                    wrtNeutrino;
  double                  radius;
  JRange<double>          height;                  
  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)         = JSHOWERPREFIT, JSHOWERPOINTSIMPLEX, JSHOWERPOSITIONFIT, JSHOWERDIRECTIONPREFIT, JSHOWERCOMPLETEFIT, JSHOWER_BJORKEN_Y; 
    zap['a'] = make_field(atmosphere)          = JAtmosphericMuon(90.0, 90.0);
    zap['O'] = make_field(option)              = "E", "N", "LINE", "LOGE";
    zap['I'] = make_field(isMuon);
    zap['w'] = make_field(wrtNeutrino);
    zap['r'] = make_field(radius)              = numeric_limits<double>::max();     // no containment
    zap['z'] = make_field(height)              = JRange<double>();                  // no containment
    zap['d'] = make_field(debug)               = 2;

    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
  
  cout << "APPLICATION " << application << 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);

  JCylinder3D containment(JCircle2D(JPosition2D(), radius), height.getLowerLimit(), height.getUpperLimit());

  containment.add(origin);

  NOTICE("Offset: "      << offset      << endl);
  NOTICE("Origin: "      << origin      << endl);
  NOTICE("Cylinder:    " << cylinder    << endl);
  NOTICE("Containment: " << containment << endl);

  const double EMIN_GEV = 10e3;  // MUPAGE

  TFile out(outputFile.c_str(), "recreate");

  TH1D job("job", NULL, 100, 0.5, 100.5);

  TH1D hn("hn",  "NDF (Number of Degrees of Freedom)", 250, -0.5, 249.5);  
  TH1D hq("hq",  "Fit Quality",  300,    0.0, 600.0);    
  TH1D hi("hi",  "Fit Index vs Best Resolution",  101, -0.5, 100.5);  // to see if the fit with the best quality has also the best resolution
  TH1D hd("hd",  "Distance between Reco and True position",  2000,    0.0,  100.0);    // [m]
  TH1D hz("hz",  "Longitudinal Distance between Reco and MC lepton position",  100, -200.0, 200.0);    // [m]
  TH1D ht("ht",  "Time difference between Reco and MC lepton",  100, 0, 100.0);    // [ns]
  TH1D h4D("h4D", "4D-distance between reco and true vertex", 2000, 0.0, 100.0); //[m]
  
  TH1D ha("ha",  "Angle between reco direction and true direction",  1000, 0.0, 180.0);    // [ns]         
  
  TH1D e0("e0",  "True lepton energy",  100, volume.getXmin(), volume.getXmax());       // [log(E)]
  TH1D e2("e2",  "Reconstructed energy",  100, volume.getXmin(), volume.getXmax());       // [log(E)]
  TH1D er("er",  "Ratio of reconstructed energy and true energy",  100, -5.0, +5.0);       // [log(E/E)]
  TH1D ea("ea", "Distribution of Energy error for all the events; E_{reco}-E_{true}", 100, -30, +30);
  TH1D ed3_5GeV("ed3_5GeV", "Distribution of Energy error for events with MC energy #in [3,5] GeV; E_{reco}-E_{true}", 100, -30, +30);
  TH1D ed8_11GeV("ed8_11GeV", "Distribution of Energy error for events with MC energy #in [8,11] GeV; E_{reco}-E_{true}", 100, -30, +30);
  TH1D ed15_21GeV("ed15_21GeV", "Distribution of Energy error for events with MC energy #in [15,21] GeV; E_{reco}-E_{true}", 100, -30, +30);
  TH2D ee("ee",  "; E_{true} [GeV]; E_{reco} [GeV]",
          40, 0, 4,
          40, 0, 4);
  TH1D eeres("eeres","Median reco energy as a function of true energy",40,0,4);
  TH1D eeres16("eeres16","16%Quantile reco energy as a function of true energy",40,0,4);
  TH1D eeres84("eeres84","84%Quantile reco energy as a function of true energy",40,0,4);
  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 herrorE("herrorE", "Energy Error", X.size() - 1, X.data());
  TProfile hfracE( "hfracE", "Fractional Energy Error", X.size() - 1, X.data());
  TProfile he("he","Reconstruction Efficiency", X.size() - 1, X.data());
  TProfile htheta_nu_lep("htheta_nu_lep", "Angle between neutrino and primary lepton", X.size() - 1, X.data());
  TH2D hgetln("hgetln","Kinematic angle", 40,0,4,1000,0,180);
  TH1D hln1d("hln1d","Angle between neutrino and lepton",40,0,4);
  TH2D     hb("hb", NULL, 100, 0.0, 0.5e6, 100, -100.0, +900.0);
  TH2S     hmca("hmca", NULL, X.size() - 1, X.data(), 1000, 0, 180);
  TH2D     hae("hae",NULL,40,0,4,1000,0,180);
  TH1D hmae("hmae","Mean angle deviation as function of log MC energy",40,0,4);
  TH2D hzenith("hzenith", "Reco Zenith vs MC Energy", 100, 0, 100, ny, ymin, ymax);
  TH2D hY("hY", "Reco Bjorken Y vs MC Bjorken Y", 40, 0, 1, 40, 0, 1);
  TH2D hby3_5GeV("hby3_5GeV", "Reco Bjorken Y vs MC Bjorken Y for events with Reco E #in [3, 5] GeV", 40, 0, 1, 40, 0, 1);
  TH2D hby8_10GeV("hby8_10GeV", "Reco Bjorken Y vs MC Bjorken Y for events with Reco E #in [8, 10] GeV", 40, 0, 1, 40, 0, 1);

  JQuantile Q("Angle", true);
  JQuantile O("Omega", true);
  JQuantile QE("Energy", true);
  JQuantile Qp("Pos",true);
  JQuantile Qt("Time",true);
  JQuantile Q4d("4D",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;  // Reco event
    Evt*       event = ps;  // MC event

    const time_converter converter(*event, *tev);

    job.Fill(1.0);

    double Enu  = 0.0;
    double Elep = 0.0;
    double Emax = 0.0;

    Trk neutrino;

    if(has_neutrino(*event)){
    
      neutrino = get_neutrino(*event); 
      Enu = neutrino.E;
    
    }
    
    vector<Trk>::const_iterator lepton = event->mc_trks.end();  // can be electron or muon

    for (vector<Trk>::const_iterator mc_evt = event->mc_trks.begin(); mc_evt != event->mc_trks.end(); ++mc_evt) {
        
      if (!isMuon && is_electron(*mc_evt)) {

        Elep += mc_evt->E;
        
        if (mc_evt->E > Emax) {
          lepton = mc_evt;
          Emax = mc_evt->E;
        }


      } else if(isMuon && is_muon(*mc_evt)){
        
        Elep += mc_evt->E;
        if (mc_evt->E > Emax) {
          lepton = mc_evt;
          Emax = mc_evt->E;
        }
      }
    }

    if (lepton == event->mc_trks.end()) {
      continue;
    }

    job.Fill(3.0);

    double true_BjY = (Enu - Elep) / Enu;
    
    // abscissa

    Double_t x = 0.0;
    
    if (option.find('E') != string::npos){

      if(wrtNeutrino == true && !isMuon) x = volume.getX(neutrino.E);
      else if(!wrtNeutrino && !isMuon)   x = volume.getX(lepton->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 && numberOfPrefits < (size_t) distance(evt->begin(), evt->end())) {     
      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);
      }

      JEvt::iterator best = evt->begin();  // the best shower is the one with the best quality

      JPosition position(getPosition(lepton->pos)); // by default the lepton pos is the interaction vertex
      JPointing pointing(getDirection(*lepton));
      JShowerEnergy energy(lepton->E);

      JVector3D showerBrightPoint(getPosition(lepton->pos) + offset + getPosition(lepton->dir) * geanz.getMaximum(lepton->E));

      if(!wrtNeutrino && !isMuon){
        position = JPosition(showerBrightPoint);  // wrt the em shower brightest point
      } else if(wrtNeutrino == true && !isMuon){
        position = JPosition(getPosition(neutrino));
        pointing = JPointing(getDirection(neutrino));
        energy = JShowerEnergy(Enu);
      }

      JEvt::iterator fit_with_smallest_error = best;

      if(application == JSHOWERPREFIT || application == JSHOWERPOINTSIMPLEX || application == JSHOWERPOSITIONFIT){
        fit_with_smallest_error  = position(evt->begin(), __end);
      } else if (application == JSHOWERENERGYPREFIT){
        fit_with_smallest_error  = energy(evt->begin(), __end);
      } else {
        fit_with_smallest_error  = pointing(evt->begin(), __end);
      }
      best = fit_with_smallest_error;
      
      const Double_t beta  = pointing.getAngle(*best);
      const double   Efit  = best->getE();

      if (containment.is_inside(getPosition(*best))){
      
        // selection of fit result      
        bool ok = (Efit >= Emin_GeV);
      
        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(), fit_with_smallest_error));
          Q.put(beta);

          JShower3E ta;

          if(wrtNeutrino){
            ta = getTrack(neutrino);
          } else {
            ta = getTrack(*lepton);  // MC event
          }

          JShower3E tb = getTrack(*best);  // Reco event

          double zenith = tb.getDirection().getTheta() * 180 / PI;

          // take into account the difefrent MC/Reco geometry and time
          ta.add(offset);      
          tb.sub(converter.putTime());
          double time_true = ta.getT();
          double time_reco=tb.getT();
          double distance_m;
          
          if(!isMuon && !wrtNeutrino){
            JPosition3D posDif = tb.getPosition() - JPosition3D(showerBrightPoint);
            distance_m = posDif.getLength();
            time_true +=  geanz.getMaximum(lepton->E)*getInverseSpeedOfLight();
            hd.Fill(fabs(distance_m));
            ht.Fill(fabs(time_reco - time_true));
            Qp.put(fabs(distance_m));
              Qt.put(fabs(time_reco-time_true));
            double distance4d = sqrt(posDif.getLengthSquared() + pow(getSpeedOfLight()/getIndexOfRefraction()*(time_true-time_reco),2));
            h4D.Fill(distance4d);
            Q4d.put(distance4d);
          } else {
            JPosition3D posDif = tb.getPosition() - ta.getPosition();
            distance_m = posDif.getLength();
            hd.Fill(fabs(distance_m));
            ht.Fill(fabs(time_reco - time_true));
            Qp.put(fabs(distance_m));
              Qt.put(fabs(time_reco-time_true));
            double distance4d = sqrt(posDif.getLengthSquared() + pow(getSpeedOfLight()/getIndexOfRefraction()*(time_true-time_reco),2));
            h4D.Fill(distance4d);
            Q4d.put(distance4d);
          }
        
          if (has_neutrino(*event)) {
          
            job.Fill(6.0);

            JPosition3D mc_vx = getPosition(neutrino); // mc vertex
          
            mc_vx.add(offset);

            if (cylinder.is_inside(mc_vx)) {

              job.Fill(7.0);

              hz.Fill(tb.getIntersection(mc_vx));
            }
          }
        
          if (best->getE() >= EMIN_GEV) {
            hb.Fill(JVector2D(best->getX() - origin.getX(), best->getY() - origin.getY()).getLengthSquared(), best->getZ());
          }
        
          ha.Fill(getAngle(ta.getDirection(), tb.getDirection()));  

          htheta_nu_lep.Fill(x, getAngle(getDirection(neutrino), getDirection(*lepton)));
          hgetln.Fill(log10(Enu), getAngle(getDirection(neutrino), getDirection(*lepton)));
           for (int i = 1; i <= hgetln.GetNbinsX(); ++i) {
            std::vector<double> deltaThetaValues;
            for (int j = 1; j <= hgetln.GetNbinsY(); ++j) {
              double binContent = hgetln.GetBinContent(i, j);
              // Add Delta theta values based on the bin content (number of events)
              for (int k = 0; k < binContent; ++k) {
                double deltaTheta = hgetln.GetYaxis()->GetBinCenter(j);
                deltaThetaValues.push_back(deltaTheta);
              }
            }

            // Calculate the median using TMath::Median
            double medianDeltaTheta = TMath::Median(deltaThetaValues.size(), &deltaThetaValues[0]);

            // Set the median Delta theta as a function of energy in the graph
            hln1d.SetBinContent(i, medianDeltaTheta);
          }
          hmca.Fill(x, getAngle(ta.getDirection(), tb.getDirection()));
          hae.Fill(log10(Enu),getAngle(ta.getDirection(), tb.getDirection()));
          for (int i = 1; i <= hae.GetNbinsX(); ++i) {
            std::vector<double> deltaThetaValues;
            for (int j = 1; j <= hae.GetNbinsY(); ++j) {
              double binContent = hae.GetBinContent(i, j);
              // Add Delta theta values based on the bin content (number of events)
              for (int k = 0; k < binContent; ++k) {
                double deltaTheta = hae.GetYaxis()->GetBinCenter(j);
                deltaThetaValues.push_back(deltaTheta);
              }
            }

            // Calculate the median using TMath::Median
            double medianDeltaTheta = TMath::Median(deltaThetaValues.size(), &deltaThetaValues[0]);

            // Set the median Delta theta as a function of energy in the graph
            hmae.SetBinContent(i, medianDeltaTheta);
          }
          if(application == JSHOWER_BJORKEN_Y) {

            hY.Fill(true_BjY, best->getW(5));

            if(Efit > 2.5 && Efit < 6) hby3_5GeV.Fill(true_BjY, best->getW(5));
            else if(Efit > 7.5 && Efit < 11) hby8_10GeV.Fill(true_BjY, best->getW(5));

          }
          if(wrtNeutrino){

            e0.Fill(volume.getX(Enu,  true));
            er.Fill(volume.getX(Efit) - volume.getX(Enu));
            ee.Fill(volume.getX(Enu), volume.getX(Efit));
            ea.Fill(Efit - Enu);
            hzenith.Fill(Enu, zenith);
            QE.put(log10(Efit/Enu));
            herrorE.Fill(x, volume.getX(Efit) - volume.getX(Enu));
            hfracE.Fill(x, fabs(volume.getX(Efit) - volume.getX(Enu))/volume.getX(Enu));        

            if(Enu >= 3 && Enu <= 5) ed3_5GeV.Fill(Efit - Enu);
            else if(Enu >= 8 && Enu <= 11) ed8_11GeV.Fill(Efit - Enu);
            else if(Enu >= 15 && Enu <= 21) ed15_21GeV.Fill(Efit - Enu);

          } else {

            e0.Fill(volume.getX(Elep,  true));
            er.Fill(volume.getX(Efit) - volume.getX(Elep));
            ee.Fill(volume.getX(Elep), volume.getX(Efit));
            ea.Fill(Efit - Elep);
            hzenith.Fill(Elep, zenith);
            QE.put(log10(Efit/Elep));
            herrorE.Fill(x, volume.getX(Efit) - volume.getX(Elep));
            hfracE.Fill(x, fabs(volume.getX(Efit) - volume.getX(Elep))/volume.getX(Elep));      

            if(Elep >= 3 && Elep <= 5) ed3_5GeV.Fill(Efit - Elep);
            else if(Elep >= 8 && Elep <= 11) ed8_11GeV.Fill(Efit - Elep);
            else if(Elep >= 15 && Elep <= 21) ed15_21GeV.Fill(Efit - Elep);

          }
        for (int i = 1; i <= ee.GetNbinsX(); ++i) {
            std::vector<double> Erecovalues;
             for (int j = 1; j <= ee.GetNbinsY(); ++j) {
              double binContent = ee.GetBinContent(i, j);
              // Add E values based on the bin content (number of events)
              for (int k = 0; k < binContent; ++k) {
                double Erecobin = ee.GetYaxis()->GetBinCenter(j);
                Erecovalues.push_back(Erecobin);
              }
            }
             if(Erecovalues.empty()) continue;
             eeres.SetBinContent(i, Erecovalues[int(Erecovalues.size()*0.5)]);
             eeres16.SetBinContent(i, Erecovalues[int(Erecovalues.size()*0.16)]);
             eeres84.SetBinContent(i, Erecovalues[int(Erecovalues.size()*0.84)]);}
          e2.Fill(volume.getX(Efit, true));
        
          h2.Fill(x, beta);
        }
      }
    }
    
    he.Fill(x, W);}
 
  STATUS(endl);

  NOTICE("Number of events input           " << setw(8) << right << job.GetBinContent(1) << endl);
  if(!isMuon){
    NOTICE("Number of events with electron   " << setw(8) << right << job.GetBinContent(3) << endl);
  } else{
    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, 70, 80 and 90% quantiles of space angle [deg] " << FIXED     (6,3) << Q.getQuantile(0.5) << " "<< Q.getQuantile(0.7) << " "<< Q.getQuantile(0.8)<< " "<< Q.getQuantile(0.9)<< endl);
  }
  if (QE.getCount() != 0) {
      NOTICE("Median, 70%, 80% and 90% quantiles of energy ratio log(Ereco/Etrue) " << FIXED     (6,3) << QE.getQuantile(0.5) << " "<< QE.getQuantile(0.7) << " "<< QE.getQuantile(0.8)<< " "<< QE.getQuantile(0.9) <<  endl);
  }
  if(Qp.getCount() != 0){
        NOTICE("Median distance to vertex:      " << setw(8) << Qp.getQuantile(0.5) << " "<< Qp.getQuantile(0.7)<<" "<< Qp.getQuantile(0.8)<<  " "<<Qp.getQuantile(0.9)<<endl);
        NOTICE("Median time to vertex:      " << setw(8)  << Qt.getQuantile(0.5) << " "<< Qt.getQuantile(0.7)<<" "<< Qt.getQuantile(0.8)<<  " "<<Qt.getQuantile(0.9)<<endl);
        NOTICE("Median 4D-distance to vertex:      " << setw(8)  << Q4d.getQuantile(0.5) << " "<< Q4d.getQuantile(0.7)<<" "<< Q4d.getQuantile(0.8)<<  " "<<Q4d.getQuantile(0.9)<<endl);
      }
      out.Write();
      out.Close();
}