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 "evt/Head.hh"
#include "evt/Evt.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JDAQ/JDAQEvent.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);
  }



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

      if (angleA > angle_Deg) { 
        manager.Fill(*fitA, *fitB, angleB < angle_Deg);
      }
    }
  }
  STATUS(endl);

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