res_plots.C

#include "FilePropertiesManager.hh"
#include "HistUtils.hh"
#include "UniverseMaker.hh"


// Script intended to help with choosing binning for kinematic variables

// The anticipated POT to use when scaling the MC prediction in the expected
// reco events plot. This will help ensure that all choices of reco binning
// are informed by the expected statistical uncertainties when the full dataset
// is analyzed.
constexpr double EXPECTED_POT = 6.790e20; // Full dataset for Runs 1-3

// Number of true bins to use when plotting true distributions in a given
// reco bin
constexpr int DEFAULT_TRUE_BINS = 100;

// ROOT integer code for Arial font
constexpr int FONT_STYLE = 62; // Arial

// When bins have zero content, set them to this very small value so that the
// colz style will still paint them
constexpr double REALLY_SMALL = 1e-11;


void make_res_plots( const std::string& branchexpr,
  const std::string& variable_title, const std::string& selection,
  const std::set<int>& runs, std::vector<double> bin_low_edges,
  bool show_bin_plots = true,
  bool show_smear_numbers = false,
  int num_true_bins = DEFAULT_TRUE_BINS,
  const std::string& mc_branchexpr = "",
  const std::string& signal_cuts = "mc_is_signal",
  const std::string& mc_event_weight = DEFAULT_MC_EVENT_WEIGHT )
{
  // Get the outer edges of the reco-space bins. This will be used to set the
  // plot range for the true-space histograms.
  double xmin = bin_low_edges.front();
  double xmax = bin_low_edges.back();

  // If the user hasn't explicitly specified a branch expression for the
  // true quantity, assume that it's the same as the reco quantity but
  // with the prefix "mc_" added.
  std::string true_branchexpr = mc_branchexpr;
  if ( true_branchexpr.empty() ) {
    true_branchexpr = "mc_" + branchexpr;
  }

  // Get access to the singleton utility class that manages the processed
  // ntuple files
  const FilePropertiesManager& fpm = FilePropertiesManager::Instance();

  // Make a TChain to process the CV numu ntuples from the requested run(s).
  // For the resolution studies, this is all we really need. Add the
  // appropriate ntuples to the TChain. Also tally the total simulated
  // POT for later scaling purposes.
  TChain chain( "stv_tree" );
  double total_simulated_POT = 0.;

  const auto& ntuple_map = fpm.ntuple_file_map();
  for ( const auto& run : runs ) {
    const auto& ntuple_files = ntuple_map.at( run )
      .at( NtupleFileType::kNumuMC );
    for ( const auto& file_name : ntuple_files ) {
      chain.Add( file_name.c_str() );

      TFile temp_file( file_name.c_str(), "read" );
      TParameter<float>* temp_pot = nullptr;
      temp_file.GetObject( "summed_pot", temp_pot );
      double pot = temp_pot->GetVal();
      total_simulated_POT += pot;
    }
  }

  // Dummy counter used to ensure that each histogram generated by this
  // function has a unique name to use with TTree::Draw()
  static int hist_count = 0;

  if ( show_bin_plots ) {
    for ( size_t b = 1u; b < bin_low_edges.size(); ++b ) {
      ++hist_count;

      TCanvas* c = new TCanvas;
      std::string true_hist_name = "true_hist" + std::to_string( hist_count );

      TH1D* true_hist = new TH1D( true_hist_name.c_str(),
        ("true events in " + variable_title + " reco bin "
        + std::to_string(b) + "; " + variable_title + "; events").c_str(),
        num_true_bins, xmin, xmax );

      double reco_bin_min = bin_low_edges.at( b - 1 );
      double reco_bin_max = bin_low_edges.at( b );
      std::string cuts = mc_event_weight + " * (is_mc && " + signal_cuts
        + " && " + selection + " && " + branchexpr + " >= "
        + std::to_string( reco_bin_min ) + " && " + branchexpr
        + " < " + std::to_string( reco_bin_max ) + ')';

      chain.Draw( (true_branchexpr + " >> " + true_hist_name).c_str(),
        cuts.c_str(), "goff" );

      true_hist->SetStats( false );
      true_hist->SetLineWidth( 2 );
      true_hist->SetLineColor( kBlack );

      true_hist->Draw( "hist pe" );

      // Prepare vertical lines to draw on the plot. These will show the
      // boundaries of the reco bin in true space
      double max_for_lines = std::numeric_limits<double>::max();

      TLine* line_bin_min = new TLine( reco_bin_min, 0.,
        reco_bin_min, max_for_lines );

      TLine* line_bin_max = new TLine( reco_bin_max, 0.,
        reco_bin_max, max_for_lines );

      line_bin_min->SetLineColor( kRed );
      line_bin_min->SetLineWidth( 2 );
      line_bin_min->SetLineStyle( 1 );
      line_bin_min->Draw( "same" );

      line_bin_max->SetLineColor( kRed );
      line_bin_max->SetLineWidth( 2 );
      line_bin_max->SetLineStyle( 1 );
      line_bin_max->Draw( "same" );

    } // loop over reco bins

  } // show bin plots

  // Also get the total number of reco bins for the 2D smearing plot
  int num_reco_bins = bin_low_edges.size() - 1u;

  // Compute the smearing matrix for a choice of true bins that exactly
  // match the ones in reco space.
  std::string smear_hist_name = "smear_hist" + std::to_string( hist_count );
  TH2D* smear_hist = new TH2D( smear_hist_name.c_str(),
    ("smearing matrix for " + variable_title + "; true " + variable_title
    + "; reco " + variable_title).c_str(), num_reco_bins, bin_low_edges.data(),
    num_reco_bins, bin_low_edges.data() );

  std::string smear_expr = branchexpr + " : " + true_branchexpr
    + " >> " + smear_hist_name;

  std::string smear_cuts = mc_event_weight + " * (is_mc && " + signal_cuts
    + " && " + selection + ')';

  chain.Draw( smear_expr.c_str(), smear_cuts.c_str(), "goff" );

  // Before renormalizing the smearing matrix histogram, take a projection
  // along the reco (y) axis. This will show the expected number of signal
  // events in each reco bin according to our central value MC model. Reco bins
  // should be chosen to have sufficient expected statistics in addition to
  // small smearing.
  TH1D* expected_reco_hist = smear_hist->ProjectionY();

  // Scale the expected reco bin counts to the POT analyzed for the full
  // dataset. Also set the bin stat uncertainties to the square root of their
  // contents. This is not correct for getting the MC statistical uncertainties
  // (which should use the sum of the squares of the weights to get the
  // variance), but we're less interested in those. Primarily we'd like to know
  // what the anticipated statistical uncertainties on the *measurement* will
  // be. We can estimate that by choosing the bin errors in this way. This will
  // help in the effort to choose suitable bins for reporting the final result.
  expected_reco_hist->Scale( EXPECTED_POT / total_simulated_POT );
  for ( int eb = 0; eb <= num_reco_bins + 1; ++eb ) {
    double bin_events = expected_reco_hist->GetBinContent( eb );
    double bin_stat_err = std::sqrt( std::max(0., bin_events) );
    expected_reco_hist->SetBinError( eb, bin_stat_err );
  }

  expected_reco_hist->SetStats( false );
  expected_reco_hist->SetLineColor( kBlack );
  expected_reco_hist->SetLineWidth( 2 );

  std::stringstream temp_ss;
  temp_ss << "expected reco bin counts (" << EXPECTED_POT << " POT);"
    << " reco " << variable_title << "; events";

  expected_reco_hist->SetTitle( temp_ss.str().c_str() );

  TCanvas* c_expected = new TCanvas;
  expected_reco_hist->Draw( "hist e" );

  // Normalize the smearing matrix elements so that a sum over all reco bins
  // (including the under/overflow bins) yields a value of one. This means that
  // every selected signal event must end up somewhere in reco space.
  int num_bins_x = smear_hist->GetXaxis()->GetNbins();
  int num_bins_y = smear_hist->GetYaxis()->GetNbins();

  // Loop over the true (x) bins. Include the underflow (index zero) and
  // overflow (index num_bins_x + 1) bins.
  for ( int bx = 0; bx <= num_bins_x + 1; ++bx ) {

    // For the current true (x) bin, compute the sum of all reco (y) bins.
    double y_sum = 0.;
    for ( int by = 0; by <= num_bins_y + 1; ++by ) {
      y_sum += smear_hist->GetBinContent( bx, by );
    }

    // Normalize each of the reco (y) bins so that the sum over y is unity.
    for ( int by = 0; by <= num_bins_y + 1; ++by ) {

      // To avoid dividing by zero, set the bin content to zero if the sum of
      // the reco (y) bins is not positive.
      if ( y_sum <= 0. ) {
        //smear_hist->SetBinContent( bx, by, REALLY_SMALL );
        smear_hist->SetBinContent( bx, by, 0. );
      }
      else {
        // Otherwise, normalize in the usual way
        double bc = smear_hist->GetBinContent( bx, by );

        double content = std::max( bc / y_sum, REALLY_SMALL );

        smear_hist->SetBinContent( bx, by, content );
      }
    } // loop over reco (y) bins

  } // loop over true (x) bins

  // Smearing matrix histogram style options
  smear_hist->GetXaxis()->SetTitleFont( FONT_STYLE);
  smear_hist->GetYaxis()->SetTitleFont( FONT_STYLE );
  smear_hist->GetXaxis()->SetTitleSize( 0.05 );
  smear_hist->GetYaxis()->SetTitleSize( 0.05 );
  smear_hist->GetXaxis()->SetLabelFont( FONT_STYLE );
  smear_hist->GetYaxis()->SetLabelFont( FONT_STYLE );
  smear_hist->GetZaxis()->SetLabelFont( FONT_STYLE );
  smear_hist->GetZaxis()->SetLabelSize( 0.03 );
  smear_hist->GetXaxis()->CenterTitle();
  smear_hist->GetYaxis()->CenterTitle();
  smear_hist->GetXaxis()->SetTitleOffset( 1.2 );
  smear_hist->GetYaxis()->SetTitleOffset( 1.1 );
  smear_hist->SetStats( false );
  smear_hist->SetMarkerSize( 1.8 ); // text size
  smear_hist->SetMarkerColor( kWhite ); // text color

  // Draw the smearing matrix plot
  TCanvas* c_smear = new TCanvas;
  c_smear->SetBottomMargin( 0.15 );
  c_smear->SetLeftMargin( 0.13 );

  if ( show_smear_numbers ) {
    // Round all numbers to this precision when rendering them
    gStyle->SetPaintTextFormat( "4.2f" );

    smear_hist->Draw("text colz");
  }
  else {
    smear_hist->Draw( "colz" );
  }

  // For each true bin, print the fraction of events that are reconstructed
  // in the correct corresponding reco bin.
  for ( int bb = 1; bb <= num_reco_bins; ++bb ) {
    std::cout << "bin #" << bb << ": "
      << expected_reco_hist->GetBinLowEdge( bb ) << ", "
      << smear_hist->GetBinContent(bb, bb) << '\n';
  }

}

// Overloaded version that uses a fixed number of equal-width bins
void make_res_plots( const std::string& branchexpr,
  const std::string& variable_title, const std::string& selection,
  const std::set<int>& runs,
  double xmin, double xmax, int Nbins,
  bool show_bin_plots = true,
  bool show_smear_numbers = false,
  int num_true_bins = DEFAULT_TRUE_BINS,
  const std::string& mc_branchexpr = "",
  const std::string& signal_cuts = "mc_is_signal",
  const std::string& mc_event_weight = DEFAULT_MC_EVENT_WEIGHT )
{
  auto low_edges = get_bin_low_edges( xmin, xmax, Nbins );
  return make_res_plots( branchexpr, variable_title, selection, runs,
    low_edges, show_bin_plots, show_smear_numbers, num_true_bins,
    mc_branchexpr, signal_cuts, mc_event_weight );
}

void make_res_plots( const std::string& rmm_config_file_name,
  const std::set<int>& runs,
  const std::string& universe_branch_name = "TunedCentralValue_UBGenie",
  size_t universe_index = 0u,
  bool show_smear_numbers = false )
{
  const std::string variable_title = "bin";

  // Create a UniverseMaker object that will handle the actual
  // calculation of the smearing matrix
  UniverseMaker rmm( rmm_config_file_name );

  // Get access to the singleton utility class that manages the processed
  // ntuple files
  const FilePropertiesManager& fpm = FilePropertiesManager::Instance();

  // TODO: Reduce code duplication for the POT tallying

  // Add the appropriate CV numu ntuples from the requested run(s) to the
  // TChain owned by the UniverseMaker object. For the resolution
  // studies, this is all we really need. Also tally the total simulated POT
  // for later scaling purposes.
  double total_simulated_POT = 0.;

  const auto& ntuple_map = fpm.ntuple_file_map();
  for ( const auto& run : runs ) {
    const auto& ntuple_files = ntuple_map.at( run )
      .at( NtupleFileType::kNumuMC );
    for ( const auto& file_name : ntuple_files ) {
      rmm.add_input_file( file_name );

      TFile temp_file( file_name.c_str(), "read" );
      TParameter<float>* temp_pot = nullptr;
      temp_file.GetObject( "summed_pot", temp_pot );
      double pot = temp_pot->GetVal();
      total_simulated_POT += pot;
    }
  }

  // Look up the MC event weights from the input files and construct the
  // response matrices in the usual way. For speed, restrict the calculation to
  // just the universe branch requested by the user (typically the CV branch).
  rmm.build_universes( { universe_branch_name } );

  // For all but the "unweighted" universe, the key used to look up
  // the map entry is "weight_" prepended to the original ntuple branch name.
  std::string universe_key = "unweighted";
  if ( universe_key != universe_branch_name ) {
    universe_key = "weight_" + universe_branch_name;
  }

  // Get access to the Universe object that stores the histograms of summed MC
  // event weights that we need
  const auto& universe = rmm.universe_map().at( universe_key )
    .at( universe_index );

  TH2D* smear_hist = dynamic_cast< TH2D* >(
    universe.hist_2d_->Clone("smear_hist") );

  // TODO: also reduce code duplication here

  // Before renormalizing the smearing matrix histogram, take a projection
  // along the reco (y) axis. This will show the expected number of signal
  // events in each reco bin according to our central value MC model. Reco bins
  // should be chosen to have sufficient expected statistics in addition to
  // small smearing.
  TH1D* expected_reco_hist = smear_hist->ProjectionY();
  int num_reco_bins = expected_reco_hist->GetNbinsX();

  // Scale the expected reco bin counts to the POT analyzed for the full
  // dataset. Also set the bin stat uncertainties to the square root of their
  // contents. This is not correct for getting the MC statistical uncertainties
  // (which should use the sum of the squares of the weights to get the
  // variance), but we're less interested in those. Primarily we'd like to know
  // what the anticipated statistical uncertainties on the *measurement* will
  // be. We can estimate that by choosing the bin errors in this way. This will
  // help in the effort to choose suitable bins for reporting the final result.
  expected_reco_hist->Scale( EXPECTED_POT / total_simulated_POT );
  for ( int eb = 0; eb <= num_reco_bins + 1; ++eb ) {
    double bin_events = expected_reco_hist->GetBinContent( eb );
    double bin_stat_err = std::sqrt( std::max(0., bin_events) );
    expected_reco_hist->SetBinError( eb, bin_stat_err );
  }

  expected_reco_hist->SetStats( false );
  expected_reco_hist->SetLineColor( kBlack );
  expected_reco_hist->SetLineWidth( 2 );

  std::stringstream temp_ss;
  temp_ss << "expected reco bin counts (" << EXPECTED_POT << " POT);"
    << " reco " << variable_title << "; events";

  expected_reco_hist->SetTitle( temp_ss.str().c_str() );

  TCanvas* c_expected = new TCanvas;
  expected_reco_hist->Draw( "hist e" );

  // Normalize the smearing matrix elements so that a sum over all reco bins
  // (including the under/overflow bins) yields a value of one. This means that
  // every selected signal event must end up somewhere in reco space.
  int num_bins_x = smear_hist->GetXaxis()->GetNbins();
  int num_bins_y = smear_hist->GetYaxis()->GetNbins();

  // Loop over the true (x) bins. Include the underflow (index zero) and
  // overflow (index num_bins_x + 1) bins.
  for ( int bx = 0; bx <= num_bins_x + 1; ++bx ) {

    // For the current true (x) bin, compute the sum of all reco (y) bins.
    double y_sum = 0.;
    for ( int by = 0; by <= num_bins_y + 1; ++by ) {
      y_sum += smear_hist->GetBinContent( bx, by );
    }

    // Normalize each of the reco (y) bins so that the sum over y is unity.
    for ( int by = 0; by <= num_bins_y + 1; ++by ) {

      // To avoid dividing by zero, set the bin content to zero if the sum of
      // the reco (y) bins is not positive.
      if ( y_sum <= 0. ) {
        //smear_hist->SetBinContent( bx, by, REALLY_SMALL );
        smear_hist->SetBinContent( bx, by, 0. );
      }
      else {
        // Otherwise, normalize in the usual way
        double bc = smear_hist->GetBinContent( bx, by );

        double content = std::max( bc / y_sum, REALLY_SMALL );

        smear_hist->SetBinContent( bx, by, content );
      }
    } // loop over reco (y) bins

  } // loop over true (x) bins

  // Smearing matrix histogram style options
  smear_hist->GetXaxis()->SetTitleFont( FONT_STYLE);
  smear_hist->GetYaxis()->SetTitleFont( FONT_STYLE );
  smear_hist->GetXaxis()->SetTitleSize( 0.05 );
  smear_hist->GetYaxis()->SetTitleSize( 0.05 );
  smear_hist->GetXaxis()->SetLabelFont( FONT_STYLE );
  smear_hist->GetYaxis()->SetLabelFont( FONT_STYLE );
  smear_hist->GetZaxis()->SetLabelFont( FONT_STYLE );
  smear_hist->GetZaxis()->SetLabelSize( 0.03 );
  smear_hist->GetXaxis()->CenterTitle();
  smear_hist->GetYaxis()->CenterTitle();
  smear_hist->GetXaxis()->SetTitleOffset( 1.2 );
  smear_hist->GetYaxis()->SetTitleOffset( 1.1 );
  smear_hist->SetStats( false );
  smear_hist->SetMarkerSize( 1.8 ); // text size
  smear_hist->SetMarkerColor( kWhite ); // text color

  // Draw the smearing matrix plot
  TCanvas* c_smear = new TCanvas;
  c_smear->SetBottomMargin( 0.15 );
  c_smear->SetLeftMargin( 0.13 );

  if ( show_smear_numbers ) {
    // Round all numbers to this precision when rendering them
    gStyle->SetPaintTextFormat( "4.2f" );

    smear_hist->Draw("text colz");
  }
  else {
    smear_hist->Draw( "colz" );
  }

  // Draw a vertical red line to separate the signal true bins from the
  // background true bins in the smearing matrix plot. Start by finding
  // where the first background bin is. Here we assume that the full set
  // of signal bins comes before any background bins.
  const auto& true_bins = rmm.true_bins();
  size_t num_true_bins = true_bins.size();
  size_t first_bkgd_bin_idx = num_true_bins;
  for ( size_t t = 0u; t < num_true_bins; ++t ) {
    const auto& tbin = true_bins.at( t );
    if ( tbin.type_ == TrueBinType::kBackgroundTrueBin ) {
      first_bkgd_bin_idx = t;
      break;
    }
  }

  // We've found the bin index. Now draw the line to indicate the
  // signal/background boundary in true space
  TLine* bkgd_line = new TLine( first_bkgd_bin_idx, 0.,
    first_bkgd_bin_idx, num_reco_bins );

  bkgd_line->SetLineColor( kRed );
  bkgd_line->SetLineWidth( 3 );
  bkgd_line->SetLineStyle( 2 );
  bkgd_line->Draw( "same" );

  // For each true bin, print the fraction of events that are reconstructed
  // in the correct corresponding reco bin.
  for ( int bb = 1; bb <= num_reco_bins; ++bb ) {
    std::cout << "bin #" << bb << ": "
      << expected_reco_hist->GetBinLowEdge( bb ) << ", "
      << smear_hist->GetBinContent(bb, bb) << '\n';
  }

}

void res_plots() {

  //auto& fpm = FilePropertiesManager::Instance();
  //fpm.load_file_properties( "new_file_properties.txt" );

  //make_res_plots( "delta_alphaT * 180 / TMath::ACos(-1.)", "#delta#alpha_{T}", "sel_CCNp0pi", {1},
  //  { 0, 25., 60., 95., 120., 145., 165., 180. },
  //  false );

  //make_res_plots( "delta_pT", "#deltap_{T}", "sel_CCNp0pi", {1},
  //  { 0, 0.06, 0.12, 0.18, 0.24, 0.32, 0.4, 0.48, 0.55, 0.68,
  //    0.75, 0.9 },
  //  false );


  // deltaPT in deltaAlphaT slices
  //{
  //  { 0., { 0, 0.06, 0.12, 0.18, 0.24, 0.32, 0.4, 0.48, 0.9 } },
  //  { 45., { 0, 0.06, 0.12, 0.18, 0.24, 0.32, 0.4, 0.48, 0.55, 0.9 },
  //  { 90., { 0, 0.06, 0.12, 0.18, 0.24, 0.32, 0.4, 0.48, 0.55, 0.68, 0.9 },
  //  { 135., { 0, 0.06, 0.12, 0.18, 0.24, 0.32, 0.4, 0.5, 0.6, 0.72, 0.9 } },
  //  { 180., {} },
  //}

  //// deltaAlphaT in deltaPT slices
  //{
  //  { 0., { 0, 25., 60., 95., 120., 145., 165., 180. } },
  //  { 0.2, { 0, 25., 60., 95., 120., 145., 165., 180. } },
  //  { 0.3, { 0, 25., 60., 95., 120., 145., 165., 180. } },
  //  { >= 0.4, { 0, 25., 60., 95., 120., 145., 165., 180. } },
  //}

  //make_res_plots( "delta_pTy", "#deltap_{Ty}", "sel_CCNp0pi", {1},
  //  { -0.6, -0.45, -0.35, -0.25, -0.15, -0.075, 0, 0.075, 0.15, 0.25,
  //    0.35, 0.45, 0.6 },
  // false );

  //make_res_plots( "p3_lead_p.Mag()", "p_p", "sel_CCNp0pi", {1},
  //  { 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8,
  //    0.85, 0.9, 0.95, 1.0 },
  // false );

  //make_res_plots( "p3_mu.Mag()", "p_#mu", "sel_CCNp0pi", {1},
  //  { 0.1, 0.15, 0.175, 0.2, 0.225, 0.25, 0.275, 0.3, 0.325, 0.35, 0.375, 0.4,
  //    0.425, 0.45, 0.475, 0.5, 0.55, 0.6, 0.65,
  //    0.7, 0.75, 0.8, 0.85, 0.9,
  //    0.95, 1.0, 1.1, 1.2 },
  // false );

  make_res_plots( "p3_lead_p.CosTheta()", "#cos#theta_{p}", "sel_CCNp0pi", {1},
    { -1., -0.9, -0.75, -0.6, -0.45, -0.3, -0.15, 0.0,
      0.15, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.85, 0.9, 0.925, 0.95, 0.975, 1.0 }, false );

  // delta_pTx in delta_pTy slices
  //{
  //  { < -0.15, { -0.6, -0.45, -0.35, -0.25, -0.15, -0.075, 0, 0.075, 0.15, 0.25,
  //   0.35, 0.45, 0.6 } },
  //  { >= -0.15 && < 0.15, { -0.6, -0.45, -0.35, -0.25, -0.15, -0.075, 0, 0.075, 0.15, 0.25,
  //   0.35, 0.45, 0.6 } },
  //  { > 0.15, { -0.6, -0.35, -0.25, -0.15, 0, 0.15, 0.25, 0.35, 0.6 } },
  //}

  //// Muon-proton opening angle
  //make_res_plots( "TMath::ACos( (p3_mu.X()*p3_lead_p.X() + "
  //  "p3_mu.Y()*p3_lead_p.Y() + p3_mu.Z()*p3_lead_p.Z()) / p3_mu.Mag()"
  //  "/ p3_lead_p.Mag() ) * 180. / TMath::ACos(-1.)",
  //  "#theta_{#mu-p}", "sel_CCNp0pi", std::set<int>{1},
  //  { 0., 30., 40., 50., 60., 70., 80., 90., 100., 110., 120., 130.,
  //    140., 150., 180. },
  //  false, false, DEFAULT_TRUE_BINS,

  //  "TMath::ACos( (mc_p3_mu.X()*mc_p3_lead_p.X() + mc_p3_mu.Y()"
  //  "*mc_p3_lead_p.Y() + mc_p3_mu.Z()*mc_p3_lead_p.Z()) / mc_p3_mu.Mag()"
  //  " / mc_p3_lead_p.Mag() ) * 180. / TMath::ACos(-1.)"
  //);

  //make_res_plots( "pn", "p_{n}", "sel_CCNp0pi", {1},
  //  { 0., 0.07, 0.14, 0.21, 0.28, 0.35,
  //    0.45, 0.54, 0.66, 0.77, 0.9 },
  //  false );

  // Muon-proton opening angle in pn slices
  //make_res_plots( "TMath::ACos( (p3_mu.X()*p3_lead_p.X() + "
  //  "p3_mu.Y()*p3_lead_p.Y() + p3_mu.Z()*p3_lead_p.Z()) / p3_mu.Mag()"
  //  "/ p3_lead_p.Mag() ) * 180. / TMath::ACos(-1.)",
  //  "#theta_{#mu-p}", "sel_CCNp0pi && pn >= 0.45", std::set<int>{1},
  //  { 0., 30., 45., 60., 75., 90., 105., 120., 135.,
  //    150., 180. },
  //  false, false, DEFAULT_TRUE_BINS,
  //  "TMath::ACos( (mc_p3_mu.X()*mc_p3_lead_p.X() + mc_p3_mu.Y()"
  //  "*mc_p3_lead_p.Y() + mc_p3_mu.Z()*mc_p3_lead_p.Z()) / mc_p3_mu.Mag()"
  //  " / mc_p3_lead_p.Mag() ) * 180. / TMath::ACos(-1.)"
  //);

  // Muon-proton opening angle in pn slices
  //{
  //  { 0. <= pn < 0.21, { 0., 60., 70., 80., 90., 100., 110., 120., 130.,
  //    140., 150., 180. } },
  //  { 0.21 <= pn < 0.45, { 0., 45., 60., 75., 90., 100., 110., 120., 130.,
  //    140., 150., 180. },
  //  { pn > 0.45, { 0., 30., 45., 60., 75., 90., 105., 120., 135.,
  //    150., 180. } },

  //  "TMath::ACos( (mc_p3_mu.X()*mc_p3_lead_p.X() + "
  //  "mc_p3_mu.Y()*mc_p3_lead_p.Y() + mc_p3_mu.Z()*mc_p3_lead_p.Z())"
  //  " / mc_p3_mu.Mag() / mc_p3_lead_p.Mag() )"
  //);

  //make_res_plots( "delta_pTx", "#deltap_{Tx}",
  //  "sel_CCNp0pi && mc_delta_pTy > 0.15", {1},
  // { -0.6, -0.35, -0.25, -0.15, 0, 0.15, 0.25,
  //   0.35, 0.6 },
  // false );


  //make_res_plots( "delta_alphaT * 180 / TMath::ACos(-1.)", "#delta#alpha_{T}",
  //  "sel_CCNp0pi && mc_delta_pT >= 0.4", {1},
  //  { 0, 25., 60., 95., 120., 145., 165., 180. },
  //  false );

  //make_res_plots( "myconfig_mcc9_2D_proton.txt", {1} );

  // Muon momentum
  //make_res_plots( "p3_mu.Mag()", "p_{#mu}", "sel_CCNp0pi", {1},
  //  { 0.1, 0.17, 0.2, 0.23, 0.26, 0.29, 0.32, 0.35, 0.38, 0.42,
  //    0.45, 0.48, 0.51, 0.55, 0.59, 0.64, 0.69, 0.74, 0.79,
  //    0.84, 0.89, 0.94, 1.0, 1.1, 1.2 },
  //  false );

  // Leading proton momentum
  //make_res_plots( "p3_lead_p.Mag()", "p_{lead p}",
  // "sel_CCNp0pi", {1},
  //{ 0.250, 0.325, 0.4, 0.45, 0.5, 0.550, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85,
  //  0.9, 1.0 },
  //false );

//  make_res_plots( "p3_lead_p.CosTheta()", "cos#theta_{p}",
//    "sel_CCNp0pi && mc_p3_lead_p.Mag() >= 0.75 && mc_p3_lead_p.Mag() < 0.8",
//    {1},
//{ -1, 0.55, 0.7, 0.8, 0.87, 1.0 }
//    , false );

/////// END NEW STUFF

  //make_res_plots( "p3_lead_p.CosTheta()", "cos#theta_{p}",
  // "sel_CCNp0pi", {1},
  // { -1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.5, 0.6, 0.7, 0.75, 0.8,
  //    0.85, 0.9, 0.95, 1.0 },
  // false );

  //make_res_plots( "pn", "p_{n}", "sel_CCNp0pi", {2},
  //  { 0., 0.125, 0.225, 0.325, 0.425, 0.525, 0.65, 0.85 },
  //  true, true );

  //make_res_plots( "delta_alphaT", "#delta#alpha_{T}", "sel_CCNp0pi", {1},
  //  { 0, 0.35, 0.85, 1.35, 1.85, 2.3, 2.7, 2.95, M_PI },
  //  true, false );

  //make_res_plots( "delta_pTx", "#deltap_{Tx}", "sel_CCNp0pi", {1},
  // { -0.6, -0.45, -0.35, -0.25, -0.15, -0.075, 0, 0.075, 0.15, 0.25,
  //   0.35, 0.45, 0.6 },
  // false, true );

  //make_res_plots( "delta_pTy", "#deltap_{Ty}", "sel_CCNp0pi", {1},
  // { -0.8, -0.55, -0.39, -0.2125, -0.05, 0.1, 0.225, 0.3375, 0.5 },
  // false, true );

  //make_res_plots( "delta_pL", "#deltap_{L}", "sel_CCNp0pi", {1},
  //  { -0.8, -0.6, -0.475, -0.35, -0.225, -0.115,
  //    -0.0285, 0.0575, 0.145, 0.230, 0.315,
  //    0.4 }, false, true );

  //make_res_plots( "delta_phiT", "#delta#phi_{T}", "sel_CCNp0pi", {1},
  //{ 0., 0.075, 0.2, 0.35, 0.5, 0.7, 0.9, 1.15, 1.4, 1.65, 1.9, 2.35,
  //  2.8, M_PI },
  //false, false );

  //make_res_plots( "delta_pT", "#deltap_{T}", "sel_CCNp0pi", {1},
  //  { 0, 0.1, 0.2, 0.3, 0.4, 0.525, 0.675, 0.9 },
  //  false, true );

  /// Muon-proton opening angle
  //make_res_plots( "TMath::ACos( (p3_mu.X()*p3_lead_p.X() + "
  //  "p3_mu.Y()*p3_lead_p.Y() + p3_mu.Z()*p3_lead_p.Z()) / p3_mu.Mag()"
  //  "/ p3_lead_p.Mag() )", "#theta_{#mu-p}", "sel_CCNp0pi", std::set<int>{1},

  //  { 0, 0.52, 0.78, 1.0, 1.15, 1.35, 1.5, 1.65, 1.8, 1.95, 2.1,
  //    2.35, 2.62, M_PI },

  //  false, false, DEFAULT_TRUE_BINS,

  //  "TMath::ACos( (mc_p3_mu.X()*mc_p3_lead_p.X() + "
  //  "mc_p3_mu.Y()*mc_p3_lead_p.Y() + mc_p3_mu.Z()*mc_p3_lead_p.Z())"
  //  " / mc_p3_mu.Mag() / mc_p3_lead_p.Mag() )"
  //);

  //make_res_plots( "p3_mu.Mag()", "p_{#mu}", "sel_CCNp0pi", {1},
  //  { 0.1, 0.17, 0.24, 0.3, 0.48, 0.75, 1.14, 2.5 },
  //  false );

  //make_res_plots( "p3_mu.CosTheta()", "cos(#theta_{#mu})",
  // "sel_CCNp0pi && p3_mu.Mag() >= 0.1 && p3_mu.Mag() < 0.48"
  //   " && p3_lead_p.Mag() >= 0.25 && p3_lead_p.Mag() < 0.5"
  //   " && p3_lead_p.CosTheta() >= 0.5 && p3_lead_p.CosTheta() < 1.", {1},
  // { -1, -0.1, 0.35, 0.7, 0.85, 1.00 },
  // false );

  //make_res_plots( "p3_mu.CosTheta()", "cos(#theta_{#mu})",
  // "sel_CCNp0pi", {1},
  // { -1, -0.85, -0.775, -0.7, -0.625, -0.55, -0.475, -0.4, -0.325,
  //   -0.25, -0.175, -0.1, -0.025, 0.05, 0.125, 0.2, 0.275, 0.35, 0.425, 0.5,
  //    0.575, 0.65, 0.725, 0.8, 0.85, 0.875, 0.9, 0.925, 0.950, 0.975, 1.00 },
  // false );

  //make_res_plots( "p3_mu.Mag()", "p_{#mu}", "sel_CCNp0pi", {1},
  //  { 0.1, 0.17, 0.21, 0.24, 0.27, 0.3, 0.38, 0.48, 0.75, 1.14, 2.5, 100. },
  //  false );

  ///////////////////////////////////////////////////////
  // 2D muon binning study
  ///////////////////////////////////////////////////////

  //make_res_plots( "p3_mu.CosTheta()", "cos#theta_{#mu}",
  //  "sel_CCNp0pi && p3_mu.Mag() >= 0.1 && p3_mu.Mag() < 0.17", {1},
  //  { -1, 0., 1. },

  // "sel_CCNp0pi && p3_mu.Mag() >= 0.17 && p3_mu.Mag() < 0.21", {1},
  // { -1, -0.2, 0.4, 1. },

  // "sel_CCNp0pi && p3_mu.Mag() >= 0.21 && p3_mu.Mag() < 0.24", {1},
  // { -1, -0.2, 0.4, 1.00 },

  // "sel_CCNp0pi && p3_mu.Mag() >= 0.24 && p3_mu.Mag() < 0.27", {1},
  // { -1, -0.1, 0.5, 1.00 },

  // "sel_CCNp0pi && p3_mu.Mag() >= 0.27 && p3_mu.Mag() < 0.3", {1},
  // { -1, -0.1, 0.35, 0.6, 1.00 },

  // "sel_CCNp0pi && p3_mu.Mag() >= 0.3 && p3_mu.Mag() < 0.38", {1},
  // { -1, -0.4, -0.1, 0.1, 0.35, 0.5, 0.7, 0.85, 1.00 },

  //"sel_CCNp0pi && p3_mu.Mag() >= 0.38 && p3_mu.Mag() < 0.48", {1},
  //{ -1, 0, 0.5, 0.65, 0.8, 0.92, 1.00 },

  //"sel_CCNp0pi && p3_mu.Mag() >= 0.48 && p3_mu.Mag() < 0.75", {1},
  //{ -1, 0.2, 0.5, 0.65, 0.8, 0.875, 0.950, 1.00 },

  //"sel_CCNp0pi && p3_mu.Mag() >= 0.75 && p3_mu.Mag() < 1.14", {1},
  //{ -1, 0.5, 0.8, 0.875, 0.950, 1.00 },

  //"sel_CCNp0pi && p3_mu.Mag() >= 1.14 && p3_mu.Mag() < 2.5", {1},
  //{ -1, 0.75, 0.85, 0.9, 0.950, 1.00 },

  //false );

  //make_res_plots( "myconfig_mcc9_2D_muon.txt", {1} );

  ///////////////////////////////////////////////////////
  // 2D proton binning study
  ///////////////////////////////////////////////////////

  //make_res_plots( "p3_lead_p.Mag()", "p_{lead p}",
  // "sel_CCNp0pi", {1},
  //{ 0.250, 0.325, 0.4, 0.45, 0.5, 0.550, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9,
  //  0.975, 1.20 },
  //false );

  //make_res_plots( "p3_lead_p.CosTheta()", "cos#theta_{p}",
  // "sel_CCNp0pi", {1},
  // { -1, -0.8, -0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.5, 0.6, 0.7, 0.75, 0.8,
  //    0.85, 0.9, 0.95, 1.0 },
  // false );

  //make_res_plots( "p3_lead_p.CosTheta()", "cos#theta_{p}",
  //  "sel_CCNp0pi && p3_lead_p.Mag() >= 0.250 && p3_lead_p.Mag() < 0.325", {1},
  //  { -1, -0.5, 0.1, 0.6, 1.0 },

  //  "sel_CCNp0pi && p3_lead_p.Mag() >= 0.325 && p3_lead_p.Mag() < 0.4", {1},
  //  { -1, -0.7, -0.4, 0, 0.4, 0.6, 0.8, 1.0 },

  //  "sel_CCNp0pi && p3_lead_p.Mag() >= 0.4 && p3_lead_p.Mag() < 0.45", {1},
  //   { -1, -0.5, -0.1, 0.2, 0.5, 0.65, 0.85, 1.0 },

  //  "sel_CCNp0pi && p3_lead_p.Mag() >= 0.45 && p3_lead_p.Mag() < 0.5", {1},
  //   { -1, -0.4, 0, 0.2, 0.4, 0.55, 0.65, 0.8, 0.92, 1.0 },

//CCCC
  //  "sel_CCNp0pi && p3_lead_p.Mag() >= 0.5 && p3_lead_p.Mag() < 0.550", {1},
  //  { -1, -0.2, 0.2, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 },

  //  "sel_CCNp0pi && p3_lead_p.Mag() >= 0.550 && p3_lead_p.Mag() < 0.6", {1},
  //  { -1, 0, 0.3, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 },

//CCCC
  //  "sel_CCNp0pi && p3_lead_p.Mag() >= 0.6 && p3_lead_p.Mag() < 0.65", {1},
  //  { -1, 0.1, 0.37, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 },

  //  "sel_CCNp0pi && p3_lead_p.Mag() >= 0.65 && p3_lead_p.Mag() < 0.7", {1},
  //  { -1, 0.3, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 },

//CCCC
  //  "sel_CCNp0pi && p3_lead_p.Mag() >= 0.7 && p3_lead_p.Mag() < 0.75", {1},
  //  { -1, 0.45, 0.65, 0.75, 0.82, 0.9, 1.0 },

//CCCC
  //  "sel_CCNp0pi && p3_lead_p.Mag() >= 0.75 && p3_lead_p.Mag() < 0.8", {1},
  //  { -1, 0.55, 0.7, 0.8, 0.87, 1.0 },

//CCCC
  //  "sel_CCNp0pi && p3_lead_p.Mag() >= 0.8 && p3_lead_p.Mag() < 0.85", {1},
  // { -1, 0.65, 0.78, 0.89, 1.0 },

//CCCC
  //  "sel_CCNp0pi && p3_lead_p.Mag() >= 0.85 && p3_lead_p.Mag() < 0.9", {1},
  //  { -1, 0.73, 0.86, 1.0 },

//CCCC 150
    //"sel_CCNp0pi && p3_lead_p.Mag() >= 0.9 && p3_lead_p.Mag() < 0.975", {1},
    //{ -1, 0.77, 0.88, 1.0 },

//CCCC 150
  //  "sel_CCNp0pi && p3_lead_p.Mag() >= 0.975 && p3_lead_p.Mag() < 1.20", {1},
  //  { -1, 0.84, 1.0 },

  //false );

  //make_res_plots( "myconfig_mcc9_2D_proton.txt", {1} );
}