Geant4ReweightSysts.cc File Reference

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include "cetlib/filepath_maker.h"
#include "fhiclcpp/ParameterSet.h"
#include "ProtoDUNESelectionUtils.h"
#include "TFile.h"

Go to the source code of this file.

Classes
struct	systConfig

Typedefs
typedef std::vector< std::string >	string_v
typedef std::vector< int >	int_v
typedef std::vector< double >	double_v

Functions
bool	parseArgs (int argc, char **argv)
fhicl::ParameterSet	makeParameterSet (std::string file)
void	saveHists (std::vector< TH1 * > &vec)
void	makeFracErrorHist (TH1 *variation, TH1 *nominal)
int	main (int argc, char **argv)

Variables
std::string	output_file_name
std::string	fcl_file_name

Typedef Documentation

typedef std::vector<double> double_v

Definition at line 28 of file Geant4ReweightSysts.cc.

typedef std::vector<int> int_v

Definition at line 27 of file Geant4ReweightSysts.cc.

typedef std::vector<std::string> string_v

Definition at line 26 of file Geant4ReweightSysts.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 74 of file Geant4ReweightSysts.cc.

                                 {

  if (!parseArgs(argc, argv)) return 0; 

  fhicl::ParameterSet pset = makeParameterSet(fcl_file_name);

  systConfig cfg(pset);
  double tmin = cfg.TruthBinning[0];
  double tmax = cfg.TruthBinning.back();

  std::vector<TH1 *> bg_syst_hists;
  std::vector<TH1 *> bg_nom_hists;
  std::vector<TH1 *> sig_syst_hists;
  std::vector<TH1 *> sig_nom_hists;
  std::vector<TH1 *> inc_nom_hists;
  std::vector<TH1 *> inc_syst_hists;

  for (size_t i = 0; i < cfg.ChannelNames.size(); ++i) {
  
    //Backgrounds
    for (size_t j = 0; j < cfg.BackgroundTopology.size(); ++j) {
      int topo = cfg.BackgroundTopology[j];

      //Get nominal
      bg_nom_hists.push_back(
          protoana::ProtoDUNESelectionUtils::FillMCBackgroundHistogram_Pions(
              cfg.MCFileNames[i], cfg.RecoTreeName, cfg.RecoBinning,
              cfg.ChannelNames[i], cfg.BackgroundTopologyName[j], topo,
              cfg.EndZCut, tmin, tmax, false/*cfg.DoNegativeReco*/, 0));

      //Get Variations
      for (size_t k = 0; k < cfg.SystToConsider.size(); ++k) {
        std::cout << "Considering " << cfg.SystToConsider[k] << std::endl;
        for (auto const m : {-1, 1}) {
          bg_syst_hists.push_back(
              protoana::ProtoDUNESelectionUtils::FillMCBackgroundHistogram_Pions(
                  cfg.MCFileNames[i], cfg.RecoTreeName, cfg.RecoBinning,
                  cfg.ChannelNames[i], cfg.BackgroundTopologyName[j], topo,
                  cfg.EndZCut, tmin, tmax, false/*cfg.DoNegativeReco*/,
                  /*doSyst=*/m, /*systName=*/cfg.SystToConsider[k]));

          //Turn into frac errors
          makeFracErrorHist(bg_syst_hists.back(), bg_nom_hists.back());
        }
      }
    }

    //Signals
    for (size_t j = 0; j < cfg.SignalTopology.size(); ++j) {
      int topo = cfg.SignalTopology[j];
      for (size_t k = 1; k < cfg.TruthBinning.size(); ++k) {
        
        //Get Nominal
        sig_nom_hists.push_back(
            protoana::ProtoDUNESelectionUtils::FillMCSignalHistogram_Pions(
                cfg.MCFileNames[i], cfg.RecoTreeName, cfg.RecoBinning,
                cfg.ChannelNames[i], cfg.SignalTopologyName[j], topo,
                cfg.TruthBinning[k-1], cfg.TruthBinning[k], cfg.EndZCut,
                false/*cfg.DoNegativeReco*/));

        //Get Variations
        for (size_t m = 0; m < cfg.SystToConsider.size(); ++m) {
          for (auto const n : {-1, 1}) {
            std::cout << "Considering " << cfg.SystToConsider[m] << std::endl;
            sig_syst_hists.push_back(
                protoana::ProtoDUNESelectionUtils::FillMCSignalHistogram_Pions(
                    cfg.MCFileNames[i], cfg.RecoTreeName, cfg.RecoBinning,
                    cfg.ChannelNames[i], cfg.SignalTopologyName[j], topo,
                    cfg.TruthBinning[k-1], cfg.TruthBinning[k], cfg.EndZCut,
                    false/*cfg.DoNegativeReco*/,
                    /*doSyst=*/n, /*systName=*/cfg.SystToConsider[m]));

            //Turn into frac errors
            makeFracErrorHist(sig_syst_hists.back(), sig_nom_hists.back());
          }
        }
      }
    }

  }
  
  //Incident
  for (size_t j = 0; j < cfg.IncidentTopology.size(); ++j) {
    //Form nominal
    TH1 * temp_hist =
        protoana::ProtoDUNESelectionUtils::FillMCIncidentHistogram_Pions(
            cfg.IncidentMCFileNames[0], cfg.RecoTreeName, cfg.RecoBinning,
            /*cfg.ChannelNames[0],*/ cfg.IncidentTopologyName[j], cfg.IncidentTopology[j],
            0., 10000., cfg.EndZCut, false);
    for (size_t k = 1; k < cfg.IncidentMCFileNames.size(); ++k) {
      std::cout << "Trying " << k << std::endl;
      temp_hist->Add(
          protoana::ProtoDUNESelectionUtils::FillMCIncidentHistogram_Pions(
            cfg.IncidentMCFileNames[k], cfg.RecoTreeName, cfg.RecoBinning,
            /*cfg.ChannelNames[k],*/ cfg.IncidentTopologyName[j], cfg.IncidentTopology[j],
            0., 10000., cfg.EndZCut, false));
    }
    inc_nom_hists.push_back(temp_hist);

    //GetVariations
    for (size_t m = 0; m < cfg.SystToConsider.size(); ++m) {
      for (auto const n : {-1, 1}) {
        std::cout << "Considering " << cfg.SystToConsider[m] << std::endl;

        TH1 * temp_hist = 
            protoana::ProtoDUNESelectionUtils::FillMCIncidentHistogram_Pions(
                cfg.IncidentMCFileNames[0], cfg.RecoTreeName, cfg.RecoBinning,
                /*cfg.ChannelNames[0],*/ cfg.IncidentTopologyName[j], cfg.IncidentTopology[j],
                0., 10000., cfg.EndZCut, false, n, cfg.SystToConsider[m]);
        for (size_t k = 1; k < cfg.IncidentMCFileNames.size(); ++k) {
          temp_hist->Add(
              protoana::ProtoDUNESelectionUtils::FillMCIncidentHistogram_Pions(
                  cfg.IncidentMCFileNames[k], cfg.RecoTreeName, cfg.RecoBinning,
                  /*cfg.ChannelNames[k],*/ cfg.IncidentTopologyName[j], cfg.IncidentTopology[j],
                  0., 10000., cfg.EndZCut, false, n, cfg.SystToConsider[m]));
        }
        inc_syst_hists.push_back(temp_hist);

        //Turn into frac errors
        makeFracErrorHist(inc_syst_hists.back(), inc_nom_hists.back());
      }
    }
  }
    




  TFile output(output_file_name.c_str(), "RECREATE");
  output.cd();
  saveHists(bg_syst_hists/*, cfg.SystToConsider*/);
  saveHists(sig_syst_hists/*, cfg.SystToConsider*/);
  saveHists(inc_syst_hists);
  output.Close();
  return 0;
}

void makeFracErrorHist	(	TH1 *	variation,
		TH1 *	nominal
	)

Definition at line 219 of file Geant4ReweightSysts.cc.

                                                       {
  // Scale by -1 to subtract
  nominal->Scale(-1.);

  // Subtract from var
  variation->Add(nominal);

  // Scale back to normal
  nominal->Scale(-1.);

  // Divide to get fractional error
  variation->Divide(nominal);  
}

fhicl::ParameterSet makeParameterSet

(

std::string

file

)

Definition at line 268 of file Geant4ReweightSysts.cc.

                                                   {
  fhicl::ParameterSet pset;

  // Configuration file lookup policy.
  char const* fhicl_env = getenv("FHICL_FILE_PATH");
  std::string search_path;

  if (fhicl_env == nullptr) {
    std::cerr << "Expected environment variable FHICL_FILE_PATH is missing" <<
                 " or empty: using \".\"\n";
    search_path = ".";
  }
  else {
    search_path = std::string{fhicl_env};
  }

  cet::filepath_first_absolute_or_lookup_with_dot lookupPolicy{search_path};
  return fhicl::ParameterSet::make(file, lookupPolicy);
}

bool parseArgs	(	int	argc,
		char **	argv
	)

Definition at line 233 of file Geant4ReweightSysts.cc.

                                       {

  bool found_fcl_file = false;
  bool found_output_file = false;

  for (int i = 1; i < argc; ++i) {
    if ((strcmp(argv[i], "--help") == 0) || (strcmp(argv[i], "-h") == 0)) {
      std::cout << "Usage: ./Geant4ReweightSysts -c <fitter_config>.fcl " <<
                   " -o <output_file>.root" << std::endl;
      std::cout << std::endl;

      return false;
    }

    else if (strcmp(argv[i], "-c") == 0) {
      fcl_file_name = argv[i+1];
      found_fcl_file = true;
    }

    else if (strcmp(argv[i], "-o") == 0) {
      output_file_name = argv[i+1];
      found_output_file = true;
    }

  }

  if (!found_fcl_file || !found_output_file) {
    std::cout << "Error: Must provide fcl file with options '-c' and '-o'" <<
                 std::endl;
    return false;
  }

  return true;
}

void saveHists

(

std::vector< TH1 * > &

vec

)

Definition at line 211 of file Geant4ReweightSysts.cc.

                                        {
  for (auto * h : vec) {
    //std::string name = h->GetName();
    //h->SetName((name + "_" + addToName).c_str());
    h->Write(/*(name + "_" + addToName).c_str()*/);
  }
}

Variable Documentation

std::string fcl_file_name

Definition at line 16 of file Geant4ReweightSysts.cc.

std::string output_file_name

Definition at line 15 of file Geant4ReweightSysts.cc.