mec_tuning_fitter_2020.C

Go to the documentation of this file.

/*
 *  run this after mec_tuning_preds_2020.C 
 *  Reads the filled predictions and data spectrum
 *    to fit MEC either with 2D gaussian (MINERvA-style)
 *        or 200 normalization bin dhifts (NOvA2018-style)
 */

#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/LoadFromFile.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 "3FlavorAna/Cuts/NumuCuts2020.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/Prediction/PredictionNoOsc.h"
//#include "CAFAna/Vars/GenieWeights.h"
//#include "CAFAna/Vars/NumuVars.h"
//#include "CAFAna/Vars/PPFXWeights.h"

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

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

#include "StandardRecord/Proxy/SRProxy.h"

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


#include <fstream>



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


void mec_tuning_fitter_2020
( 
  const std::string beam="fhc",
  const std::string fit="bins", //(doublegauss(B)(B1), gauss, bins) 
  double chisq=1 //scaling factor for chisquare 
)
{
  const bool isRHC = beam == "rhc";

const ana::Cut kSelCuts = ana::kNumu2020ND; 

// we need to create vector of ISyst before we read predictions (or macro breaks)
//------------------------------------ could be more adaptable --------------------------------------------//
  std::vector<std::unique_ptr<const ana::CompNormSyst>> systs= ana::systsQ3Q0(24, 0.0, 1.2, 24, 0.0, 1.2,
                                                               ana::GetCutIsFitMEC( isRHC ) && kSelCuts); 
//std::vector<std::unique_ptr<const ana::CompNormSyst>> systs= ana::systsQ3Q0(20, 0.0, 1.0, 16, 0.0, 0.8, 
  //                                                            ana::GetCutIsFitMEC( isRHC ) && kSelCuts);
  std::vector<const ana::ISyst*> systs_ptrs;
  if (fit=="bins"){
     for ( const auto & s : systs ){
      systs_ptrs.push_back( s.get() );
    }
  }
  if (fit=="gauss"){ systs_ptrs=ana::systs_params_gauss(); std::cout << "systs_ptrs.size() = "<< systs_ptrs.size()<<std::endl;}
  if (fit=="gausssupp"){ systs_ptrs=ana::systs_params_gauss_supp(); std::cout << "systs_ptrs.size() = "<< systs_ptrs.size()<<std::endl;}
  if (fit=="doublegauss"){ systs_ptrs=ana::systs_params_double_gauss(); }
  if (fit=="doublegaussB" || fit=="doublegaussB1"){systs_ptrs=ana::systs_params_double_gauss_extra();}
  
// --------------------------------------------------------------------------------------------------------//
 // read hists and predictions
  const std::string out_dir = GetMECTuningDirectory();
  // If ran on grid, you should place this file in your directory (ana::GetMECTuningDirectory())
  std::string file_name = "mec_tuning_predictions_"+beam+"_"+fit+"_2020.root";

  file_name=out_dir +"/"+ file_name;

  std::cout << "reading file from "<<file_name << std::endl;

  TFile *inFile = new TFile( file_name.c_str(), "READ" );
  ana::Spectrum spec_fit_data = *ana::Spectrum::LoadFrom(inFile, "spec_data").release();
  ana::Spectrum spec_fit_mec_raw=*ana::Spectrum::LoadFrom(inFile, "spec_mc_mec_raw").release();
  ana::Spectrum spec_fit_non_mec_raw = *ana::Spectrum::LoadFrom(inFile, "spec_mc_nonmec_raw").release();
  ana::Spectrum spec_fit_mc_raw = *ana::Spectrum::LoadFrom(inFile, "spec_mc_raw").release();
  inFile->Close();

  std::cout<< "pot data: "<< spec_fit_data.POT()<<std::endl;
  std::cout<< "pot mc: "<< spec_fit_mc_raw.POT()<<std::endl;

  ana::PredictionInterp * pred_interp= ana::LoadFromFile<ana::PredictionInterp>(file_name, "pred_interp").release();
  ana::PredictionInterp * pred_interpMEC= ana::LoadFromFile<ana::PredictionInterp>(file_name, "pred_interpMEC").release();

  ana::BigChi2SingleSampleExperiment expt(pred_interp,spec_fit_data, chisq);
  //ana::SingleSampleExperiment expt( pred_interp, spec_fit_data );
  osc::OscCalcPMNSOpt calc;
  std::vector<const ana::ISyst*> systs_ptrs2 = systs_ptrs;
  //systs_ptrs2.erase(systs_ptrs2.begin()+6, systs_ptrs2.begin() + 12);  //only let peak 1 float (higher q0q3)
  //  systs_ptrs2.erase(systs_ptrs2.begin(), systs_ptrs2.begin() + 6);   // only let peak 2 float (low q0q3)
  // systs_ptrs2.erase(systs_ptrs2.begin()+7, systs_ptrs2.begin() + 8);  // remove mean_q0_2
  //systs_ptrs2.erase(systs_ptrs2.begin()+5, systs_ptrs2.begin() + 6);  systs_ptrs2.erase(systs_ptrs2.end()-2, systs_ptrs2.end() -1); //remove correlations
  //systs_ptrs2.erase(systs_ptrs2.end()-1, systs_ptrs2.end()); //rm baseline
  //systs_ptrs2.erase(systs_ptrs2.begin()+1, systs_ptrs2.begin()+6); systs_ptrs2.erase(systs_ptrs2.begin()+2, systs_ptrs2.begin()+7); // only leave normalizations

  ana::MinuitFitter fitter( &expt, {}, systs_ptrs2);// ,  ana::MinuitFitter::FitOpts::kCareful );

  ana::SystShifts shifts;
  for ( auto & syst : systs_ptrs2 )
  {
    shifts.SetShift( syst, 0 );
  }
 double chi= fitter.Fit( &calc, shifts, ana::IFitter::kVerbose )->EvalMetricVal();
  // the results of the fitting operation are now in the shifts object.
  
  //==================================================
  // Output files
  std::string fit_results_file_name = out_dir + "/mec_fit_results_" + beam + "_" + fit + "_2020.txt";//_bigchisq.txt";
  std::ofstream fit_results_file( fit_results_file_name.c_str(), std::ofstream::out );
  for ( auto & syst : systs_ptrs2 )
  {
      fit_results_file << syst->ShortName() << " fitted shift = " << shifts.GetShift( syst )<< std::endl;
      std::cout << syst->ShortName() << " fitted shift = " << shifts.GetShift( syst )<< std::endl;
  }
  fit_results_file << "ChiSq scale in fitter =" << chisq << std::endl;
  fit_results_file << "ChiSq from in fitter = " << chi << std::endl;
  // Save spectra before and after tune. 
  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 );
  ana::Spectrum spec_fit_mc_tunedMEC( pred_interpMEC->PredictSyst( &calc, shifts ) );
  TH2  * h_mc_tunedMEC = spec_fit_mc_tunedMEC.ToTH2( pot );
  auto h_mec_raw= spec_fit_mec_raw.ToTH2(pot);
  TH2 * ratio_mecs = (TH2*)h_mc_tunedMEC->Clone("ratio_mecs");
  ratio_mecs->Divide(h_mec_raw);
  new TCanvas;
  ratio_mecs->SetTitle("ratio raw/tuned pred MEC");
  ratio_mecs->Draw("colz");
  std::cout << "here 2"<<std::endl;
  ana::Spectrum spec_fit_mc_tuned_man( pred_interp->PredictSyst( &calc, shifts ) );
  auto h_mc_tuned_man = spec_fit_mc_tuned_man.ToTH2( pot );

  h_mc_tuned->SetTitle("tuned");
  h_mc_tuned->Draw("colz");
  new TCanvas;
  h_mc_raw->SetTitle("raw");
  h_mc_raw->Draw("colz");
  new TCanvas; h_data->SetTitle("data");
  h_data->Draw("colz"); 
  new TCanvas; TH2D *ratio = (TH2D*)h_mc_tuned->Clone("ratio");
  ratio->Divide(h_data);
  ratio->SetTitle("tuned 2D gauss / data");
  ratio->Draw("colz");
  TH1 * datay=(TH1*)h_data->ProjectionY("datay");
  TH1 * datax=(TH1*)h_data->ProjectionX("datax");
  TH1 * rawy=(TH1*)h_mc_raw->ProjectionY("rawy"); rawy->SetLineColor(kRed);
  TH1 * rawx=(TH1*)h_mc_raw->ProjectionX("rawx"); rawx->SetLineColor(kRed);
  TH1 * tunedy=(TH1*)h_mc_tuned->ProjectionY("tunedy"); tunedy->SetLineColor(kBlue);
  TH1 * tunedx=(TH1*)h_mc_tuned->ProjectionX("tunedx"); tunedx->SetLineColor(kBlue);
  TH1 * tunedym=(TH1*)h_mc_tunedMEC->ProjectionY("tunedym"); tunedym->SetLineColor(kGreen+2);
  TH1 * tunedxm=(TH1*)h_mc_tunedMEC->ProjectionX("tunedxm"); tunedxm->SetLineColor(kGreen+2);
  TH1 * untunedymfile = (TH1*)h_mec_raw->ProjectionY("untunedymfile"); 
  TH1 * untunedxmfile = (TH1*)h_mec_raw->ProjectionX("untunedxmfile"); 

  //visaulize dist before and after fit.
  new TCanvas;
  datay->Draw("");
  tunedy->Draw("hist same");
  tunedym->Draw("hist same");
  untunedymfile->Draw("hist same");
  rawy->Draw("hist same");
  new TCanvas;
  datax->Draw("");  
  tunedx->Draw("hist same");
  tunedxm->Draw("hist same");
  untunedxmfile->Draw("hist same");
  rawx->Draw("hist same"); 

  std::string fit_hists_file_name = out_dir + "/mec_fit_hists_" + beam + "_" +fit+ "_2020_chis"+std::to_string(chisq)+".root";//_bigchisq.root";
  TFile fit_hists_file( fit_hists_file_name.c_str(), "recreate" );
  fit_hists_file.cd();
  h_non_mec_raw->Write("nonmec_raw");
  h_data->Write( "data" );
  h_mc_raw->Write( "mc_raw" );
  h_mc_tuned->Write( "mc_tuned" );
  h_mc_tunedMEC->Write("mc_tuned_MEC");
  h_mec_raw->Write("mc_mec_raw");
  spec_fit_mc_tunedMEC.SaveTo(& fit_hists_file, "spec_tunedMEC") ;
  new TCanvas;
  h_mc_tunedMEC->SetTitle("tuned mec");
  h_mc_tunedMEC->Draw("colz");
  new TCanvas;
  h_mec_raw->SetTitle("raw MEC");// h_mec_raw->SetMaximum(h_mc_tunedMEC->GetMaximum());
  h_mec_raw->Draw("colz");

 // for the normalization bins of MEC, save a 2D histo with the weights
 // ugly device to save weights 
  TH2D* hist_tuned_mec_weights = new TH2D( "tuned_mec_weights", ";True |#vec{q}| (GeV);True q_{0} (GeV);Weight",
                                           24, 0., 1.2, 24, 0.0, 1.2 );   
 if (fit=="bins"){
   std::vector<ana::Q3Q0Bin> binsi;
   for ( int i = 0; i < 24; ++i )
   {
     for ( int j = 0; j < 24; ++j )
     { double binWidthQ3=1.2/24;  
       double binWidthQ0=1.2/24;
       double q3 =  binWidthQ3 * ( i + 0.5 );
       double q0 =  binWidthQ0 * ( j + 0.5 );
       if ( q3 * q3 < q0 * q0 ) continue;
       ana::Q3Q0Bin b;
       b.loQ3 = q3 - binWidthQ3 / 2;
       b.hiQ3 = q3 + binWidthQ3 / 2;
       b.loQ0 = q0 - binWidthQ0 / 2;
       b.hiQ0 = q0 + binWidthQ0 / 2;
       binsi.push_back( b );
    }
  }  
  std::cout << " binsi.size() " << binsi.size()<< std::endl;
  for ( std::size_t compIdx = 0; compIdx < binsi.size(); compIdx++ )
  {
   // 24, 0.0, 1.2, 24, 0.0, 1.2,
   // 20, 0.0, 1.0, 16, 0.0, 0.
    auto b = binsi[ compIdx ];
    float q3 = ( b.loQ3 + b.hiQ3 ) / 2;
    float q0 = ( b.loQ0 + b.hiQ0 ) / 2;         
    const auto & syst = systs[ compIdx ];
    double sigma_scale = syst.get()->GetSigmaScale();
    double scale_factor = 1 + sigma_scale * shifts.GetShift( syst.get() );
    double the_shift=shifts.GetShift( syst.get() );
    if (scale_factor < 0) scale_factor=0;
    hist_tuned_mec_weights->SetBinContent( hist_tuned_mec_weights->FindFixBin( q3, q0 ), scale_factor );
    std::cout << Form( "q3 = %.2f, q0 = %.2f, shift = %.2f , sigma_scale = %.1f, scale_factor = %.2f", q3, q0, the_shift, sigma_scale, scale_factor ) << std::endl;
  }
  new TCanvas;
  hist_tuned_mec_weights->Draw("colz");
  hist_tuned_mec_weights->Write(  "mec_weights" );
  }

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

  //fit_hists_file.Close();

  if (fit=="doublegauss" || fit=="doublegaussB" || fit=="doublegaussB1"){
    for ( auto & syst : systs_ptrs2 )
    {
      std::cout << syst->ShortName() << "old param = " << 
      ana::CalcMECDoubleGaussEnhShiftedParam( syst->ShortName(), 0 ) << 
      ",  new param = " << 
      ana::CalcMECDoubleGaussEnhShiftedParam( syst->ShortName(), shifts.GetShift( syst ))<< std::endl;
    }
  }


}