mec_tuning.C

Go to the documentation of this file.

/*
 * mec_tuning.C
 *
 *  Construct numu inclusive plots with Empirical MEC weighted to match Valencia MEC in true q_0
 *
 *  Created on: Mar 29, 2017
 *      Author: J. Wolcott <jwolcott@fnal.gov>
 */



#include "TCanvas.h"
#include "TGraphErrors.h"
#include "TH1D.h"
#include "TH2D.h"
#include "THStack.h"
#include "TLegend.h"

//#include "CAFAna/Analysis/Ana2017Loaders.h"
#include "CAFAna/Fit/MinuitFitter.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2017.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Prediction/PredictionScaleComp.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/PPFXWeights.h"

#include "CAFAna/Core/ISyst.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Systs/XSecSysts.h"

#include "OscLib/OscCalcPMNSOpt.h"
#include "NumuEnergy/func/Numu2017Fits.h"

//#include "Weights.h"
//#include "DummyExtrapStuff.h"
#include "MECTuningUtils.h"

#include <fstream>

using ana::operator+;

namespace jw
{

  const ana::Var GetWeightFromShifts( const ana::SystShifts& sys_shifts )
  {
    const ana::Var kWeightFromShifts( [ sys_shifts ]( const caf::SRProxy *sr )
    {
      caf::SRProxy* face_palm = const_cast<caf::SRProxy*>( sr );
      double weight = 1.0;
      sys_shifts.Shift( face_palm, weight );
      return weight;
    });
    return kWeightFromShifts;
  }

  const ana::Cut GetCutIsFitMEC( const bool isRHC )
  {
    const ana::Cut kIsFitMEC( [ isRHC ]( const caf::SRProxy * sr )
    {
      if ( !ana::kIsDytmanMEC( sr ) || !sr->mc.nu[0].iscc ) return false;
      if      ( !isRHC && sr->mc.nu[0].pdg ==  14 ) return true;
      else if (  isRHC && sr->mc.nu[0].pdg == -14 ) return true;
      else return false;
    });
    return kIsFitMEC;
  }

  //===============================================================================
  // I/O
  //===============================================================================

  std::string GetMECTuningDirectory()
  {
    return Form( "/nova/ana/users/%s/mec_tuning", getenv( "USER" ) );
  }

  std::string GetNeutrinoTag( const bool isRHC )
  {
    return isRHC ? "numubar" : "numu";
  }

  std::string GetFitTag( const bool isRHC, const std::string shift_non_mec )
  {
    std::string fit_tag = GetNeutrinoTag( isRHC );
    if ( !shift_non_mec.empty() ) fit_tag += "_shift_non_mec_" + shift_non_mec;
    return fit_tag;
  }

  std::string GetMECWeightsFileName( const bool isRHC, const std::string shift_non_mec )
  {
    return "mec_weights_" + GetFitTag( isRHC, shift_non_mec ) + ".root";
  }

  //===============================================================================
  // Load numu MEC weights for RHC fit
  //===============================================================================

  TH2D* hist_input_mec_weights_numu = NULL;
  void LoadWeightsTunedNumuMEC( const std::string shift_non_mec )
  {
    std::string numu_mec_weights_file_name = GetMECTuningDirectory() + "/" + GetMECWeightsFileName( false, shift_non_mec );
    if ( gSystem->AccessPathName( numu_mec_weights_file_name.c_str() ) ) throw std::runtime_error( "FHC MEC weights file not found:\n" + numu_mec_weights_file_name );
    std::cout << "Loading FHC MEC weights from file:\n" << numu_mec_weights_file_name << std::endl;
    TFile* numu_mec_weights_file = new TFile( numu_mec_weights_file_name.c_str() );
    hist_input_mec_weights_numu = (TH2D*)numu_mec_weights_file->Get( "numu_mec_weights_smoothed" );
  }

  const ana::Var kWeightTunedNumuMEC( []( const caf::SRProxy * sr )
  {
    if ( !ana::kIsDytmanMEC( sr ) || !sr->mc.nu[0].iscc ) return 1.0;
    if ( sr->mc.nu[0].pdg != 14 ) return 1.0;
    if ( hist_input_mec_weights_numu == NULL ) throw std::runtime_error( "FHC MEC true (q0,q3) weights histogram is NULL" );

    double q0 = ana::kTrueQ0( sr );
    double q3 = ana::kTrueQ3( sr );

    if ( q0 < hist_input_mec_weights_numu->GetYaxis()->GetXmin() ||
         q0 > hist_input_mec_weights_numu->GetYaxis()->GetXmax() ||
         q3 < hist_input_mec_weights_numu->GetXaxis()->GetXmin() ||
         q3 > hist_input_mec_weights_numu->GetXaxis()->GetXmax() ) return 1.0;

    double weight = hist_input_mec_weights_numu->Interpolate( q3, q0 );
    //double weight = hist_input_mec_weights_numu->GetBinContent( hist_input_mec_weights_numu->FindFixBin( q3, q0 ) );
    if ( weight < 0 ) weight = 0.0;
    return weight;
  });

  //===============================================================================
  // q0 & q3 bins
  //===============================================================================

  struct Q3Q0Bin
  {
    float loQ3;
    float hiQ3;
    float loQ0;
    float hiQ0;
  };

  ana::Cut Q3Q0CutFactory( float loQ3, float hiQ3, float loQ0, float hiQ0 )
  {
    return std::move( ana::Cut ( [ loQ3, hiQ3, loQ0, hiQ0 ]( const caf::SRProxy * sr )
             {
               double q3 = ana::kTrueQ3( sr );
               double q0 = ana::kTrueQ0( sr );
               return q3 > loQ3 && q3 < hiQ3 && q0 > loQ0 && q0 < hiQ0;
             }
    ));
  }

  /*
  // we don't want penalties used in fitting
  class PenaltyFreeSystShifts : public ana::SystShifts
  {
    public:
      virtual ~PenaltyFreeSystShifts() {};

      //double Penalty() const override;
      double Penalty() const;
  };
  double PenaltyFreeSystShifts::Penalty() const
  {
    std::cout << "Penalty is 0!" << std::endl;
    return 0;
  }
  */

  class IRescaledSigmaSyst : public ana::ISyst
  {
    public:
      IRescaledSigmaSyst(const std::string & shortName, const std::string & latexName, double sigmaScale=1.0)
        : ana::ISyst(shortName, latexName), fSigmaScale(sigmaScale) {};
      double GetSigmaScale() const    { return fSigmaScale; };
      void   SetSigmaScale(double sc) { fSigmaScale = sc; };
    protected:
      double fSigmaScale;
  };

  class CompNormSyst : public IRescaledSigmaSyst
  {
    public:
      CompNormSyst(const ana::Cut & selCut, double sigmaScale=1.0)
        : fID(sInstanceCount++),
          IRescaledSigmaSyst("CompNormShift_" + std::to_string(sInstanceCount), "Component Normalization Shift " + std::to_string(sInstanceCount), sigmaScale),
          fSelCut(selCut)
      {}
      void Shift(double sigma, caf::SRProxy* sr, double& weight) const override;

    private:
      const ana::Cut fSelCut;
      unsigned int fID;

      static unsigned int sInstanceCount;
  };
  unsigned int CompNormSyst::sInstanceCount = 0;
  void CompNormSyst::Shift(double sigma, caf::SRProxy* sr, double& weight) const
  {
    if (!fSelCut(sr))
      return;

    double wgt = 1 + fSigmaScale * sigma;
    if (wgt < 0)
      wgt = 0;
    weight *= wgt;
  }

  class MECInitStateNPFracShift : public IRescaledSigmaSyst
  {
    public:
      MECInitStateNPFracShift( double sigmaScale = 1.0 )
        : IRescaledSigmaSyst( "MECInitStateCompShift", "MEC initial state np fraction", sigmaScale )
      {};

      // this guy is a bit hacksy since it doesn't make sense
      // for the fraction of np pairs to be negative or greater than 1
      void Shift( double sigma, caf::SRProxy* sr, double& weight ) const override
      {
        if ( !ana::kIsDytmanMEC( sr ) ) return;

        // Correct nominal np fraction for GENIE Empirical MEC
        const double nominalNPfrac = 0.89;
        // 10% nominal uncertainty
        double newNPfrac = nominalNPfrac + 0.1 * fSigmaScale * sigma;
        // pin to being between 0 and 1
        if      ( newNPfrac > 1 ) newNPfrac = 1;
        else if ( newNPfrac < 0 ) newNPfrac = 0;

        if ( ana::kHitNuc( sr ) == 2000000201 ) // np
          weight *= newNPfrac / nominalNPfrac;
        else
          weight *= ( 1 - newNPfrac ) / ( 1 - nominalNPfrac );
      } // MECInitStateNPFracShift::Shift()
  };  // class MECInitStateNPFracShift

} // namespace jw

std::vector<std::string> GetFiles( const std::string playlist, const int n_files_max = -1 )
{
  std::ifstream input( playlist.c_str() );
  std::vector<std::string> files;
  std::string s;
  int count = 0;
  while( getline( input, s ) )
  {
    count++;
    if ( n_files_max > -1 && count > n_files_max ) break;
    files.push_back( s );
  }
  return files;
}

////////////////////////

void mec_tuning
( 
  const std::string beam,
  const std::string shift_non_mec = ""
)
{

  if ( beam != "fhc" && beam != "rhc" ) throw std::runtime_error( "Beam must be \"fhc\" or \"rhc\"" );
  const bool isRHC = beam == "rhc";

  if ( isRHC ) jw::LoadWeightsTunedNumuMEC( shift_non_mec );

  ana::SystShifts kShiftsNonMEC;
  if ( shift_non_mec.empty() )
  {
    kShiftsNonMEC = ana::kNoShift;
  }
  else if ( shift_non_mec == "down" )
  {
    kShiftsNonMEC.SetShift( &ana::kMAQEGenieReducedSyst2018, +1 );
    kShiftsNonMEC.SetShift( &ana::kRPACCQESuppSyst2018, +1 );
    kShiftsNonMEC.SetShift( &ana::kRPACCQEEnhSyst2018, +1 );
    kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightCCQEPauliSupViaKF ), -1 );
    kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMaCCRES ), -1 );
    kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMvCCRES ), -1 );
    kShiftsNonMEC.SetShift( &ana::kRPARESSyst2018, -1 );
    //kShiftsNonMEC.SetShift( &ana::kDirectRelHadEScaleSyst2017, -1 );
  }
  else if ( shift_non_mec == "up"   )
  {
    kShiftsNonMEC.SetShift( &ana::kMAQEGenieReducedSyst2018, -1 );
    kShiftsNonMEC.SetShift( &ana::kRPACCQESuppSyst2018, -1 );
    kShiftsNonMEC.SetShift( &ana::kRPACCQEEnhSyst2018, -1 );
    kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightCCQEPauliSupViaKF ), +1 );
    kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMaCCRES ), +1 );
    kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMvCCRES ), +1 );
    kShiftsNonMEC.SetShift( &ana::kRPARESSyst2018, +1 );
    //kShiftsNonMEC.SetShift( &ana::kDirectRelHadEScaleSyst2017, +1 );
  }
  else throw std::runtime_error( Form( "Unrecognized shift: \"%s\"", shift_non_mec.c_str() ) );

  const ana::Var wgt_non_mec       = ana::kRescaleMAQE * ana::kRPAWeightCCQE2018 * ana::kRPAWeightRES2017 * ana::kFixNonres1Pi * ana::kRescaleHighWDIS;
  const ana::Var wgt_shift_non_mec = !shift_non_mec.empty() ? jw::GetWeightFromShifts( kShiftsNonMEC ) : ana::kUnweighted;
  const ana::Var wgt_mec_numu      = isRHC ? jw::kWeightTunedNumuMEC : ana::kUnweighted;
  const ana::Var wgt_fit           = ana::kPPFXFluxCVWgt * wgt_non_mec * wgt_shift_non_mec * wgt_mec_numu;

  const ana::Cut kSelCuts = ana::kNumuCutND2018;

  const int nbinsQ3 = 20; 
  const float fitMinQ3 = 0.0;
  const float fitMaxQ3 = 1.0;
  const float binWidthQ3 = ( fitMaxQ3 - fitMinQ3 ) / nbinsQ3;

  const int nbinsQ0 = 16;
  const float fitMinQ0 = 0.0;
  const float fitMaxQ0 = 0.8;
  const float binWidthQ0 = ( fitMaxQ0 - fitMinQ0 ) / nbinsQ0;

  std::vector<jw::Q3Q0Bin> binsQ3Q0;
  //std::vector<std::unique_ptr<const ana::ISyst>> systsQ3Q0;
  std::vector<std::unique_ptr<const jw::CompNormSyst>> systsQ3Q0;

  for ( int i = 0; i < nbinsQ3; ++i )
  {
    for ( int j = 0; j < nbinsQ0; ++j )
    {
      double q3 = fitMinQ3 + binWidthQ3 * ( i + 0.5 );
      double q0 = fitMinQ0 + binWidthQ0 * ( j + 0.5 );
      if ( q3 * q3 < q0 * q0 ) continue;
      jw::Q3Q0Bin b;
      b.loQ3 = q3 - binWidthQ3 / 2;
      b.hiQ3 = q3 + binWidthQ3 / 2;
      b.loQ0 = q0 - binWidthQ0 / 2;
      b.hiQ0 = q0 + binWidthQ0 / 2;
      binsQ3Q0.push_back( b );

      double sig_scale = 1.0;
      //double sig_scale = ( q3 < 0.5 && q0 < 0.2 ) ? 2.0 : 1.0;

      systsQ3Q0.emplace_back( std::make_unique<jw::CompNormSyst>( jw::Q3Q0CutFactory( b.loQ3, b.hiQ3, b.loQ0, b.hiQ0 ) && jw::GetCutIsFitMEC( isRHC ) && kSelCuts, sig_scale ) );
    }
  }

  std::cout << "systsQ3Q0.size() = " << systsQ3Q0.size() << std::endl;

  std::vector<const ana::ISyst*> systs_ptrs;
  for ( const auto & s : systsQ3Q0 )
  {
    systs_ptrs.push_back( s.get() );
  }

  std::string defn_data = "";
  std::string defn_mc   = "";
  if ( beam == "fhc" )
  {
    defn_data = "defname: prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_fhc_period2_v1_addShortSimpleCVN_goodruns_numu2018 with limit 5";
    defn_mc   = "defname: prod_sumdecaf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_period2_v1_numu2018 with limit 10";
  }
  else if ( beam == "rhc" )
  {
    defn_data = "prod_sumdecaf_R17-09-05-prod4recopreview.f_nd_numi_rhc_epoch6a_v1_addShortSimpleCVN_goodruns_numu2018 with limit 5";
    defn_mc   = "prod_sumdecaf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_epoch6a_v1_numu2018 with limit 10";
  }
  if ( defn_data.empty() ) throw std::runtime_error( "Data SAM Definition is empty" );
  if ( defn_mc  .empty() ) throw std::runtime_error( "MC SAM Definition is empty" );

  ana::Loaders loaders;
  loaders.SetLoaderPath( defn_data, caf::kNEARDET, ana::Loaders::kData );
  loaders.SetLoaderPath( defn_mc  , caf::kNEARDET, ana::Loaders::kMC   );
  loaders.SetSpillCut( ana::kStandardSpillCuts );

  ana::HistAxis axis_fit( "Reco |#vec{q}| (GeV)"   , ana::Binning::Simple( 20, 0.0, 1.0 ), ana::kRecoQmag,
                          "Reco E_{had, vis} (GeV)", ana::Binning::Simple( 20, 0.0, 0.4 ), ana::kNumuHadVisE );

  ana::Spectrum spec_fit_data
  (
    loaders.GetLoader( caf::kNEARDET, ana::Loaders::kData ),
    axis_fit,
    kSelCuts
  );

  ana::Spectrum spec_fit_mc_raw
  (
    loaders.GetLoader( caf::kNEARDET, ana::Loaders::kMC ),
    axis_fit,
    kSelCuts,
    ana::kNoShift,
    wgt_fit
  );

  ana::Spectrum spec_fit_non_mec_raw
  ( 
    loaders.GetLoader( caf::kNEARDET, ana::Loaders::kMC ), 
    axis_fit, 
    kSelCuts && !ana::kIsDytmanMEC,
    ana::kNoShift, 
    wgt_fit 
  );

  osc::NoOscillations no_osc_calc;
  jw::NDPredGenerator pred_gen( axis_fit, kSelCuts, wgt_fit );
  ana::PredictionInterp pred_interp( systs_ptrs, &no_osc_calc, pred_gen, loaders );

  loaders.Go();

  ana::SingleSampleExperiment expt( &pred_interp, spec_fit_data );
  osc::OscCalcPMNSOpt calc;
  ana::MinuitFitter fitter( &expt, {}, systs_ptrs );

  ana::SystShifts shifts;
  //jw::PenaltyFreeSystShifts shifts;
  for ( auto & syst : systs_ptrs )
  {
    shifts.SetShift( syst, 0 );
  }

  fitter.Fit( &calc, shifts, ana::IFitter::kVerbose );
  // the results of the fitting operation are now in the shifts object.

  //==================================================
  // Output files
  //==================================================

  const std::string out_dir = jw::GetMECTuningDirectory();
  if ( gSystem->AccessPathName( out_dir.c_str() ) ) gSystem->mkdir( out_dir.c_str(), true );

  const std::string fit_tag = jw::GetFitTag( isRHC, shift_non_mec );

  std::string fit_results_file_name = out_dir + "/mec_fit_results_" + fit_tag + ".txt";
  std::ofstream fit_results_file( fit_results_file_name.c_str(), std::ofstream::out );

  for ( auto & syst : systs_ptrs )
  {
    fit_results_file << syst->ShortName() << " fitted shift = " << shifts.GetShift( syst ) << std::endl;
  }

  TH2D* hist_tuned_mec_weights = new TH2D( "tuned_mec_weights", ";True |#vec{q}| (GeV);True q_{0} (GeV);Weight",
                                           nbinsQ3, fitMinQ3, fitMaxQ3, nbinsQ0, fitMinQ0, fitMaxQ0 );

  for ( std::size_t compIdx = 0; compIdx < binsQ3Q0.size(); compIdx++ )
  {
    auto b = binsQ3Q0[ compIdx ];
    float q3 = ( b.loQ3 + b.hiQ3 ) / 2;
    float q0 = ( b.loQ0 + b.hiQ0 ) / 2;
    const auto & syst = systsQ3Q0[ compIdx ];
    double sigma_scale = syst.get()->GetSigmaScale();
    double scale_factor = 1 + sigma_scale * shifts.GetShift( syst.get() );
    hist_tuned_mec_weights->SetBinContent( hist_tuned_mec_weights->FindFixBin( q3, q0 ), scale_factor );
    std::cout << Form( "q3 = %.2f, q0 = %.2f, sigma_scale = %.1f, scale_factor = %.2f", q3, q0, sigma_scale, scale_factor ) << std::endl;
  }

  TH2D* hist_tuned_mec_weights_smoothed = (TH2D*)hist_tuned_mec_weights->Clone( "tuned_mec_weights_smoothed" );
  hist_tuned_mec_weights_smoothed->Smooth( 1, "k3a" );

  std::string weights_file_name = out_dir + "/" + jw::GetMECWeightsFileName( isRHC, shift_non_mec );
  TFile weights_file( weights_file_name.c_str(), "recreate" );

  const std::string neutrino_tag = jw::GetNeutrinoTag( isRHC );
  hist_tuned_mec_weights         ->Write( Form( "%s_mec_weights", neutrino_tag.c_str() ) );
  hist_tuned_mec_weights_smoothed->Write( Form( "%s_mec_weights_smoothed", neutrino_tag.c_str() ) );
  if ( isRHC && jw::hist_input_mec_weights_numu != NULL ) jw::hist_input_mec_weights_numu->Write();
  weights_file.Close();

  double pot = spec_fit_data.POT();
  auto h_data = spec_fit_data.ToTH2( pot );
  auto h_mc_raw = spec_fit_mc_raw.ToTH2( pot );
  auto h_non_mec_raw = spec_fit_non_mec_raw.ToTH2( pot );
  ana::Spectrum spec_fit_mc_tuned( pred_interp.PredictSyst( &calc, shifts ) );
  auto h_mc_tuned = spec_fit_mc_tuned.ToTH2( pot );

  std::string fit_hists_file_name = out_dir + "/mec_fit_hists_" + fit_tag + ".root";
  TFile fit_hists_file( fit_hists_file_name.c_str(), "recreate" );
  fit_hists_file.cd();

  spec_fit_data.SaveTo(& fit_hists_file,  "spec_data" ) ;
  h_data->Write( "data" );
  h_mc_raw->Write( "mc_raw" );
  h_non_mec_raw->Write( "non_mec_raw" );
  h_mc_tuned->Write( "mc_tuned" );

  fit_hists_file.Close();

}