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 "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"
#include "JMath/JMathToolkit.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

namespace {

  using namespace JFIT;

  static const int JQUALITY = -1;   //!< fit quality identifier


  /**
   * Map key.
   */
  struct JKey {
    /**
     * Constructor.
     *
     * \param  name             parameter name
     * \param  index            parameter value (= index in user data)
     */
    JKey(const char* name,
         const int   index)
    {
      this->name  = name;
      this->index = index;
    }

    
    /**
     * Less than operator.
     *
     * \param  first            first  key
     * \param  second           second key
     * \return                  true if index of first key less than that of second; else false
     */
    friend inline bool operator<(const JKey& first, const JKey& second)
    {
      return first.index < second.index;
    }
    
    std::string name;
    int         index;
  };

  
  /**
   * Make key.
   *
   * \param  PARAMETER        parameter
   * \return                  key
   */
#define make_key(PARAMETER) JKey(#PARAMETER, PARAMETER)


  /**
   * Map value.
   */
  struct JValue {
    /**
     * Default constructor.
     */
    JValue() :
      hA(NULL),
      hB(NULL),
      hC(NULL),
      hD(NULL),
      logx(false)
    {}
    
    
    /**
     * Constructor.
     *
     * \param  hA               pointer to histogram
     * \param  hB               pointer to histogram
     * \param  hC               pointer to histogram
     * \param  hD               pointer to histogram
     * \param  logx             log10(x) filling mode
     */
    JValue(TH1D* hA,
           TH1D* hB,
           TH1D* hC,
           TH1D* hD,
           bool  logx = false)
    {
      this->hA   = hA;
      this->hB   = hB;
      this->hC   = hC;
      this->hD   = hD;
      this->logx = logx;
    }

    
    /**
     * Fill histograms.
     *
     * \param  xA               first  abscissa value
     * \param  xB               second abscissa value
     * \param  option           option
     */
    void Fill(const Double_t xA, const Double_t xB, const bool option)
    {
      hA->Fill(logx ? log10(xA) : xA);
      hB->Fill(logx ? log10(xB) : xB);
      if (option)
        hC->Fill(logx ? log10(xB/xA) : xB - xA);
      else
        hD->Fill(logx ? log10(xB/xA) : xB - xA);
    }

      
    /**
     * Scale histograms.
     */
    void Scale()
    {
      if (hA->GetEntries() != 0) { hA->Scale(1.0/hA->GetEntries()); }
      if (hB->GetEntries() != 0) { hB->Scale(1.0/hB->GetEntries()); }
      if (hC->GetEntries() != 0) { hC->Scale(1.0/hC->GetEntries()); }
      if (hD->GetEntries() != 0) { hD->Scale(1.0/hD->GetEntries()); }
    }

      
    /**
     * Write histograms to file.
     *
     * \param  out              ROOT file
     */
    void Write(TFile& out)
    {
      out.WriteTObject(hA);
      out.WriteTObject(hB);
      out.WriteTObject(hC);
      out.WriteTObject(hD);
    }

  protected:
    TH1D* hA;
    TH1D* hB;
    TH1D* hC;
    TH1D* hD;
    bool  logx;
  };


  /**
   * Auxiliary class to manage set of histograms.
   */
  struct JManager :
    public std::map<JKey, JValue>
  {
    /**
     * Insert set of histograms.
     *
     * \param  key              key
     * \param  nx               number of bins
     * \param  xmin             minimal abscissa value
     * \param  xmax             maximal abscissa value
     * \param  logx             log10(x)
     */
    void insert(const JKey&    key,
                const Int_t    nx,
                const Double_t xmin,
                const Double_t xmax,
                const double   logx = false)
    {
      using namespace std;

      TH1D* hA = new TH1D(("[A]." + key.name).c_str(), NULL, nx,  xmin, xmax);
      TH1D* hB = new TH1D(("[B]." + key.name).c_str(), NULL, nx,  xmin, xmax);
      TH1D* hC = new TH1D(("[C]." + key.name).c_str(), NULL, nx, -xmax, xmax);
      TH1D* hD = new TH1D(("[D]." + key.name).c_str(), NULL, nx, -xmax, xmax);
      
      std::map<JKey, JValue>::insert(make_pair(key, JValue(hA,hB,hC,hD,logx)));
    }

    
    /**
     * Fill histograms.
     *
     * \param  fA               first  fit result
     * \param  fB               second fit result
     * \param  option           option
     */
    void Fill(const JFit& fA, const JFit& fB, const bool option)
    {
      for (iterator i = begin(); i != end(); ++i) {

        const int index = i->first.index;
        
        if (index == JQUALITY) {
          i->second.Fill(fA.getQ(),      fB.getQ(),      option);
        } else if (fA.hasW(index) && fB.hasW(index)) {
          i->second.Fill(fA.getW(index), fB.getW(index), option);
        }
      };
    }

      
    /**
     * Scale histograms.
     */
    void Scale()
    {
      for (iterator i = this->begin(); i != this->end(); ++i) {
        i->second.Scale();
      }
    }

      
    /**
     * Write histograms to file.
     *
     * \param  out              ROOT file
     */
    void Write(TFile& out)
    {
      for (iterator i = this->begin(); i != this->end(); ++i) {
        i->second.Write(out);
      }
    }

    
    /**
     * Write histograms to file.
     *
     * \param  file_name        file name
     */
    void Write(const char* file_name)
    {
      TFile out(file_name, "RECREATE");

      this->Write(out);

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

/**
 * \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)           = 0.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.E.getLowerLimit(); //neutrino minimum energy
  const double E_nu_max = head.cut_nu.E.getUpperLimit(); //neutrino maximum energy

  inputFileA.setLimit(numberOfEvents);
  inputFileB.setLimit(numberOfEvents);


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

  JManager manager;

  manager.insert(make_key(JQUALITY),                501, 0.0, 1.0e3);
  manager.insert(make_key(JGANDALF_BETA0_RAD),      501, 0.0, 2.0e-2);
  manager.insert(make_key(JGANDALF_BETA1_RAD),      501, 0.0, 2.0e-2);
  manager.insert(make_key(JGANDALF_NUMBER_OF_HITS), 500, 0.0, 500.0);
  manager.insert(make_key(JSTART_NPE_MIP),          501, 0.0, 500.0);
  manager.insert(make_key(JSTART_NPE_MIP_TOTAL),    501, 0.0, 500.0);
  manager.insert(make_key(JSTART_LENGTH_METRES),    501, 0.0, 2.0e3);
  manager.insert(make_key(JENERGY_ENERGY),          501, 1.0,   8.0, true);
  manager.insert(make_key(JENERGY_CHI2),            501, 0.0, 500.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();
}