JPostfit2F.cc File Reference

Example program to compare fit results from two files. More...

#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"

Go to the source code of this file.

Macros
#define	make_key(PARAMETER)   JKey(#PARAMETER, PARAMETER)
	Make key. More...

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

Detailed Description

Example program to compare fit results from two files.

Author

mdejong

Definition in file JPostfit2F.cc.

Macro Definition Documentation

#define make_key

(

PARAMETER

)

JKey(#PARAMETER, PARAMETER)

Make key.

Parameters

PARAMETER

parameter

Returns

key

Definition at line 72 of file JPostfit2F.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 267 of file JPostfit2F.cc.

{
  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();
}