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 "TF2.h"
#include "TStyle.h"
#include "TLatex.h"


#include <fstream>



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



  const ana::Cut kNumu2020PIDLoosePTP(    [](const caf::SRProxy* sr)     
    {
        return (sr->sel.remid.pid > 0.7 && sr->sel.cvnloosepreselptp.numuid > 0.82);
    });

// use this to plot the weights

Double_t func0(Double_t *val, Double_t *par)
{
   Float_t x = val[0];
   Float_t y = val[1];
   Double_t f = //par[12]     // baseline: 
              + par[0]*exp( -0.5 *( pow((y-par[1])/par[3],2) + pow((x-par[2])/par[4],2)  - 2*par[5]* (y-par[1])*(x-par[2])/(par[3]*par[4] ))/(1-par[5]*par[5]  ) ) 
              + par[6]*exp( -0.5 *( pow((y-par[7])/par[9],2) + pow((x-par[8])/par[10],2) - 2*par[11]*(y-par[7])*(x-par[8])/(par[9]*par[10]))/(1-par[11]*par[11]) )
              +par[12];
   if (f<0) f=0;
   return f;
}


void mec_tuning_fitter_2020
( 
  const std::string beam="fhc",  // fhc,rhc,
  const std::string fit="doublegaussB", //(doublegauss(B)(B1), gauss, bins) 
  std::string def="mcnp",
  std::string nonmec_shift = "down", //up, down, blank
  double chisq=1 //scaling factor for chisquare 
)
{
  if(def!="mcnp") def = "prod5";
  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"){
     if (shift_non_mec=="up") systs_ptrs=ana::systs_params_double_gauss_extraUP();
       else if (shift_non_mec=="down") systs_ptrs=ana::systs_params_double_gauss_extraDOWN();
         else 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";
  if (def=="mcnp")  file_name = "mec_tuning_predictions_"+beam+"_"+fit+"_"+nonmec_shift+"2020_mcnp_genierw.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::string file_name2="/nova/ana/users/mcasales/mec_tuning_2020_systs/mec_tuning_predictions_"+beam+"_doublegaussB_2020_mcnp_half.root";
 if (def=="mcnp") file_name2="/nova/ana/users/mcasales/mec_tuning_2020_systs/mec_tuning_predictions_"+beam+"_doublegaussB_2020_mcnp_subset.root";
 TFile *inFile2 = new TFile( file_name2.c_str(), "READ" );  
 ana::Spectrum spec_fit_data = *ana::Spectrum::LoadFrom(inFile2,"spec_data").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
  //systs_ptrs2.erase(systs_ptrs2.begin()+6, systs_ptrs2.begin() + 7);  // remove norm 2

  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 + "_"+nonmec_shift+"2020_mcnp_test.txt";//_bigchisq.txt";
  if ( def=="mcnp" ) fit_results_file_name = out_dir + "/mec_fit_results_" + beam + "_" + fit + "_"+nonmec_shift+"2020_mcnp_test.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->SetMarkerStyle(8);
  datay->Draw("");
  tunedy->Draw("hist same");
  tunedym->Draw("hist same");
  untunedymfile->Draw("hist same");
  rawy->Draw("hist same");
  gPad->SaveAs((out_dir +"/mec_tuned_q0_"+beam+"_"+def+"_"+nonmec_shift+".pdf").c_str());

  new TCanvas;
  gPad->SetGridy();
  TH1* ratiorawy = (TH1*)datay->Clone("ratiorawy"); ratiorawy->SetLineColor(kRed);
  ratiorawy->GetYaxis()->SetTitle("Data/MC");
  ratiorawy->SetMinimum(0.8);
  ratiorawy->Divide(rawy);
  TH1* ratiotunedy = (TH1*)datay->Clone("ratiotunedy"); ratiotunedy->SetLineColor(kBlue);
  ratiotunedy->Divide(tunedy);
  ratiorawy->Draw("hist");
  ratiotunedy->Draw("hist same");
  gPad->SaveAs((out_dir +"/mec_tuned_q3_"+beam+"_"+def+"_"+nonmec_shift+"_ratio.pdf").c_str());

  new TCanvas;
  datax->SetMarkerStyle(8);
  datax->Draw("");  
  tunedx->Draw("hist same");
  tunedxm->Draw("hist same");
  untunedxmfile->Draw("hist same");
  rawx->Draw("hist same"); 
  gPad->SaveAs((out_dir +"/mec_tuned_q3_"+beam+"_"+def+"_"+nonmec_shift+".pdf").c_str());

  new TCanvas;
  gPad->SetGridy();
  TH1* ratiorawx = (TH1*)datax->Clone("ratiorawx"); ratiorawx->SetLineColor(kRed);
  ratiorawx->GetYaxis()->SetTitle("Data/MC");
  ratiorawx->SetMinimum(0.8);
  ratiorawx->Divide(rawx);
  TH1* ratiotunedx = (TH1*)datax->Clone("ratiotunedx"); ratiotunedx->SetLineColor(kBlue);
  ratiotunedx->Divide(tunedx);
  ratiorawx->Draw("hist");
  ratiotunedx->Draw("hist same");
  gPad->SaveAs((out_dir +"/mec_tuned_q0_"+beam+"_"+def+"_"+nonmec_shift+"_ratio.pdf").c_str());


  std::string fit_hists_file_name = out_dir + "/mec_fit_hists_" + beam + "_" +fit+nonmec_shift+ "_2020_chis"+std::to_string(chisq)+".root";//_bigchisq.root";
  if (def=="mcnp") fit_hists_file_name = out_dir + "/mec_fit_hists_" + beam + "_" +fit+nonmec_shift+ "_2020_chis"+std::to_string(chisq)+"_mcnp_test.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;
    }
  }

  TLatex latex;
  latex.SetTextSize(0.045);
  latex.SetTextAlign(13);

 TF2 *f_old = new TF2("f_old",func0, 0,1.2,0,0.8,13);

 f_old->SetParameter( 0,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystNorm_1", 0 , ""));
 f_old->SetParameter( 1,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ0_1", 0 , ""));
 f_old->SetParameter( 2,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ3_1", 0 , ""));   
 f_old->SetParameter( 3,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ0_1", 0 , ""));
 f_old->SetParameter( 4,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ3_1", 0 , ""));
 f_old->SetParameter( 5,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystCorr_1", 0 , ""));
 f_old->SetParameter( 6,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystNorm_2", 0 , ""));
 f_old->SetParameter( 7,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ0_2", 0 , ""));
 f_old->SetParameter( 8,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ3_2", 0 , ""));
 f_old->SetParameter( 9,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ0_2", 0 , ""));
 f_old->SetParameter( 10, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ3_2", 0 , ""));
 f_old->SetParameter( 11, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystCorr_2", 0 , ""));
 f_old->SetParameter( 12, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystBaseline", 0 , ""));
// f->SetContour(48);
 f_old->SetNpx(200);
 f_old->SetNpy(200);
 new TCanvas("c_weights2","c_weights2",600,500);
 gPad->SetLeftMargin(0.15); 
 gPad->SetBottomMargin(0.15);
 gPad->SetRightMargin(0.175);
 f_old->SetMaximum(42);
 f_old->SetMinimum(-.0001); 
// f_old->SetMinimum(0); 
 f_old->Draw("colz");
 f_old->GetZaxis()->SetTitle("MEC Weight");
 f_old->Draw("colz");   
 f_old->GetYaxis()->SetTitle("True q_{0} (GeV)");
 f_old->GetXaxis()->SetTitle("True #left|#vec{q}#right| (GeV)");
 f_old->GetZaxis()->SetTitle("Weights");
 f_old->GetYaxis()->SetTitleOffset(1.1);
 f_old->GetXaxis()->SetTitleOffset(1.2);
 f_old->GetZaxis()->SetTitleOffset(1.15);
 f_old->GetYaxis()->CenterTitle();
 f_old->GetXaxis()->CenterTitle();
 f_old->GetZaxis()->CenterTitle();
// f_old->GetXaxis()->SetRangeUser(0,1.2 );
// f_old->GetYaxis()->SetRangeUser(0,0.8);
 latex.DrawLatexNDC(.2,.85,("#color[0]{#nu and #bar{#nu} MEC Weights old}"));
// prelim->Draw();
 gPad->Update();
 gPad->SaveAs((out_dir +"/mec_weights_official.pdf").c_str());


 TF2 *f_new = new TF2("f_new",func0, 0,1.2,0,0.8,13);

 f_new->SetParameter( 0,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystNorm_1", shifts.GetShift(systs_ptrs2[0]) , nonmec_shift));
 f_new->SetParameter( 1,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ0_1", shifts.GetShift(systs_ptrs2[1]) , nonmec_shift));
 f_new->SetParameter( 2,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ3_1", shifts.GetShift(systs_ptrs2[2]) , nonmec_shift));
 f_new->SetParameter( 3,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ0_1", shifts.GetShift(systs_ptrs2[3]) , nonmec_shift));
 f_new->SetParameter( 4,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ3_1", shifts.GetShift(systs_ptrs2[4]) , nonmec_shift));
 f_new->SetParameter( 5,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystCorr_1", shifts.GetShift(systs_ptrs2[5]) , nonmec_shift));
 f_new->SetParameter( 6,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystNorm_2", shifts.GetShift(systs_ptrs2[6]) , nonmec_shift));
 f_new->SetParameter( 7,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ0_2", shifts.GetShift(systs_ptrs2[7]) , nonmec_shift));
 f_new->SetParameter( 8,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ3_2", shifts.GetShift(systs_ptrs2[8]) , nonmec_shift));
 f_new->SetParameter( 9,  ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ0_2", shifts.GetShift(systs_ptrs2[9]) , nonmec_shift));
 f_new->SetParameter( 10, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ3_2", shifts.GetShift(systs_ptrs2[10]) , nonmec_shift));
 f_new->SetParameter( 11, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystCorr_2", shifts.GetShift(systs_ptrs2[11]) , nonmec_shift));
 f_new->SetParameter( 12, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystBaseline", shifts.GetShift(systs_ptrs2[12]) , nonmec_shift));
  
 f_new->SetNpx(200);
 f_new->SetNpy(200);
 new TCanvas("c_weights3","c_weights3",600,500);
 gPad->SetLeftMargin(0.15); 
 gPad->SetBottomMargin(0.15);
 gPad->SetRightMargin(0.175);
 f_new->SetMaximum(42);
// f_new->SetMinimum(-.0001); 
 f_new->SetMinimum(-.0001); 
 f_new->Draw("colz");
 f_new->GetZaxis()->SetTitle("MEC Weight");
 f_new->Draw("colz");   
 f_new->GetYaxis()->SetTitle("True q_{0} (GeV)");
 f_new->GetXaxis()->SetTitle("True #left|#vec{q}#right| (GeV)");
 f_new->GetZaxis()->SetTitle("Weights");
 f_new->GetYaxis()->SetTitleOffset(1.1);
 f_new->GetXaxis()->SetTitleOffset(1.2);
 f_new->GetZaxis()->SetTitleOffset(1.15);
 f_new->GetYaxis()->CenterTitle();
 f_new->GetXaxis()->CenterTitle();
 f_new->GetZaxis()->CenterTitle();
 f_new->GetXaxis()->SetRangeUser(0,1.2 );
 f_new->GetYaxis()->SetRangeUser(0,1.2 );
 latex.DrawLatexNDC(.2,.85,("#color[0]{#nu and #bar{#nu} MEC Weights after this fit}"));
 //prelim->Draw();
 gPad->Update();
 gPad->SaveAs((out_dir +"/mec_weights_"+beam+"_"+def+nonmec_shift+".pdf").c_str());
 

}