JPlotEvtWeightFactorMeffRatio.cc

Go to the documentation of this file.

#include <iostream>
#include <string>
 
#include "TH1.h"
#include "TFile.h"
#include "TGraph.h"
#include "TGraph2D.h"
#include "TMultiGraph.h"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
#include "Jeep/JProperties.hh"
 
#include "JAAnet/JEvtWeightFactorMeffRatio.hh"
 
#include "JSirene/JVisibleEnergyToolkit.hh"
 
 
/**
 * \file
 * Example program for plotting effective mass ratios.
 *
 * \author bjjung
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  typedef double (JEvtWeightFactorMeffRatio::*pFunction1D)(const double) const;
  typedef double (JEvtWeightFactorMeffRatio::*pFunction2D)(const double,
                                                           const double) const;
  
  
  string                    outputFile;
  
  JEvtWeightFactorMeffRatio factor;
  
  int                       debug;
  
  try {
        
    JParser<> zap;
    
    zap['o'] = make_field(outputFile);
    zap['@'] = make_field(factor)          = JPARSER::initialised();
    zap['d'] = make_field(debug)           = 1;
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  TFile out(outputFile.c_str(), "recreate");  
 
  const JProperties properties = factor.getProperties();
 
  const bool logE   = properties.getValue<bool>("logE");
  const int  option = properties.getValue<int> ("option");
  
  const TH1& hMeff = factor.getHistogram();
 
  const int Nbins = hMeff.GetNbinsX();
 
  if (option != (int) JEvtWeightFactorMeffRatio::EVIS_LEADING_LEPTON_CONTR &&
      option != (int) JEvtWeightFactorMeffRatio::EVIS_HADRONIC_CONTR) {
 
    pFunction1D mfp;
    
    switch (option) { 
    case (int) JEvtWeightFactorMeffRatio::ENERGY_INITIAL_STATE :
      mfp = &JEvtWeightFactorMeffRatio::getMeffRatio1;
      break;
    case (int) JEvtWeightFactorMeffRatio::ENERGY_NEUTRINO:
      mfp = &JEvtWeightFactorMeffRatio::getMeffRatio2;
      break;
    case (int) JEvtWeightFactorMeffRatio::EVIS:
      mfp = &JEvtWeightFactorMeffRatio::getMeffRatio3;
      break;
    default:
      THROW(JValueOutOfRange, "JEvtWeightFactorRatio::configure(): Invalid 1D option " << option);
    }
 
    TGraph g(Nbins);
    g.SetName(hMeff.GetName());
  
    for (Int_t i = 0; i < Nbins; ++i) {
      
      const double x = hMeff.GetBinCenter(i+1);
      const double y = (factor.*mfp)(logE ? pow(10.0, x) : x);
 
      if (isfinite(y)) {
        g.SetPoint(i, x, y);
      } else {
        WARNING("Meff ratio not finite for x = " << SCIENTIFIC(10,3) << x << endl);
      }
    }
 
    g.Write();
    
  } else {
    
    pFunction2D mfp;
 
    switch (option) { 
    case (int) JEvtWeightFactorMeffRatio::EVIS_LEADING_LEPTON_CONTR:
      mfp = &JEvtWeightFactorMeffRatio::getMeffRatio4;
      break;
    case (int) JEvtWeightFactorMeffRatio::EVIS_HADRONIC_CONTR:
      mfp = &JEvtWeightFactorMeffRatio::getMeffRatio5;
      break;
    default:
      THROW(JValueOutOfRange, "JEvtWeightFactorRatio::configure(): Invalid 2D option " << option);
    }
 
    const Int_t N = Nbins * Nbins * (Nbins + 1) / 2;
 
    TGraph2D g(N);
    g.SetName(hMeff.GetName());
 
    Int_t k = 0;
    
    for (Int_t i = 0; i < Nbins; ++i) {
 
      const double x = hMeff.GetBinCenter(i+1);
      
      for (Int_t j = 0; j < i; ++j, ++k) {
        
        const double y = hMeff.GetBinCenter(j+1);
        const double z = (logE ?
                          (factor.*mfp)(pow(10.0, x), pow(10.0, y)) :
                          (factor.*mfp)(x, y));
        
        g.SetPoint(k, x, y, z);
      }
    }
 
    g.Write();
  }
  
  out.Write();
  out.Close();
  
  return 0;
}