JPostfit2F.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
 
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JDAQ/JDAQEventIO.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 "JMath/JMathToolkit.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
#include "JHistogramToolkit.hh"
 
 
/**
 * \file
 *
 * Example program to compare fit results from two files.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  typedef JTriggeredFileScanner<JEvt>                   JTriggeredFileScanner_t;
  typedef JTriggeredFileScanner_t::multi_pointer_type   multi_pointer_type;
 
  JTriggeredFileScanner_t inputFileA;
  JTriggeredFileScanner_t inputFileB;
  JLimit_t                numberOfEvents;
  string                  outputFile;  
  double                  angle_Deg;
  int                     debug;
 
  try {
 
    JParser<> zap("Example program to compare fit results from two files.");
 
    zap['a'] = make_field(inputFileA);
    zap['b'] = make_field(inputFileB);
    zap['o'] = make_field(outputFile)          = "postfit2f.root";
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['A'] = make_field(angle_Deg)           = 1.0;
    zap['d'] = make_field(debug)               =  2;
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
  JHead head;
 
  try {
    head = getHeader(inputFileA);
  }
  catch(const JException& error) {
    FATAL(error);
  }
  const double E_nu_min = head.cut_nu.Emin; //neutrino minimum energy
  const double E_nu_max = head.cut_nu.Emax; //neutrino maximum energy
 
  inputFileA.setLimit(numberOfEvents);
  inputFileB.setLimit(numberOfEvents);
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  JManager manager;
 
  manager.insert(make_key(JQUALITY),               101, 0.0, 250.0);
  manager.insert(make_key(JGANDALF_BETA0_RAD),      51, 0.0, 2.0e-2);
  manager.insert(make_key(JGANDALF_BETA1_RAD),      51, 0.0, 2.0e-2);
  manager.insert(make_key(JGANDALF_NUMBER_OF_HITS), 51, 0.0, 100);
  manager.insert(make_key(JSTART_NPE_MIP),          51, 0.0,  10.0);
  manager.insert(make_key(JSTART_NPE_MIP_TOTAL),    51, 0.0,  10.0);
  manager.insert(make_key(JSTART_LENGTH_METRES),    51, 0.0, 100.0);
  manager.insert(make_key(JENERGY_ENERGY),          51, 1.0,   5.0, true);
  manager.insert(make_key(JENERGY_CHI2),            51, 0.0, 200.0);
 
  TH1D hA("h[A]", NULL,  100,   -3.0,  +2.3);    // [log(deg)]
  TH1D hB("h[B]", NULL,  100,   -3.0,  +2.3);    // [log(deg)]
  TH2D hA_angle_E("hA_angle_E", NULL, 20, E_nu_min, E_nu_max, 100, 0,90);  //energy vs angular deviation between true and reco muon direction
  TH2D hB_angle_E("hB_angle_E", NULL, 20, E_nu_min, E_nu_max, 100, 0,90);  //energy vs angular deviation between true and reco muon direction
  TH2D h2("h2", NULL,
      100, -100.0, +100.0,                   // quality
      360, -180.0, +180.0);                  // deg
 
  while (inputFileA.hasNext() && inputFileB.hasNext()) {
 
    STATUS("event: " << setw(10) << inputFileA.getCounter() << '\r'); DEBUG(endl);
 
    multi_pointer_type psA = inputFileA.next();
    multi_pointer_type psB = inputFileB.next();
 
    // find the same corresponding Monte Carlo true event based on the event identifier.
 
    while (psA.get<Evt>()->mc_id < psB.get<Evt>()->mc_id && inputFileA.hasNext()) { psA = inputFileA.next(); }
    while (psB.get<Evt>()->mc_id < psA.get<Evt>()->mc_id && inputFileB.hasNext()) { psB = inputFileB.next(); }
 
    if (!psA.is_valid()) { continue; }
    if (!psB.is_valid()) { continue; }
 
    if (psA.get<Evt>()->mc_id == psB.get<Evt>()->mc_id) {
 
      Evt*  event = psA;
      JEvt* evtA  = psA;
      JEvt* evtB  = psB;
 
      const Trk neutrino = get_neutrino(*event);
 
      vector<Trk>::const_iterator muon = find_if(event->mc_trks.begin(), event->mc_trks.end(), is_muon);
 
      if (muon == event->mc_trks.end()) { continue; }
 
      if (evtA->empty())                { continue; }
      if (evtB->empty())                { continue; }
 
      JEvt::const_iterator fitA = evtA->begin();
      JEvt::const_iterator fitB = evtB->begin();
 
      const Double_t angleA = getAngle(getDirection(*muon), getDirection(*fitA));
      const Double_t angleB = getAngle(getDirection(*muon), getDirection(*fitB));
 
      hA.Fill(log10(angleA));
      hB.Fill(log10(angleB));
 
      h2.Fill(fitB->getQ() - fitA->getQ(), angleB - angleA);
 
      const double Enu = neutrino.E;    //true neutrino E  
      hA_angle_E.Fill(Enu, angleA);
      hB_angle_E.Fill(Enu, angleB);
 
      if (angleA > angle_Deg) { 
        manager.Fill(*fitA, *fitB, angleB < angle_Deg);
      }
    }
  }
  STATUS(endl);
 
  out.Write();
  out.Close();
}