mec_tuning_fitter_2020.C
Go to the documentation of this file.
1 
2 /*
3  * run this after mec_tuning_preds_2020.C
4  * Reads the filled predictions and data spectrum
5  * to fit MEC either with 2D gaussian (MINERvA-style)
6  * or 200 normalization bin dhifts (NOvA2018-style)
7  */
8 
9 #include "TCanvas.h"
10 #include "TGraphErrors.h"
11 #include "TH1D.h"
12 #include "TH2D.h"
13 #include "THStack.h"
14 #include "TLegend.h"
15 
16 //#include "CAFAna/Analysis/Ana2017Loaders.h"
18 #include "CAFAna/Core/Loaders.h"
20 #include "CAFAna/Core/Spectrum.h"
26 #include "CAFAna/Core/Var.h"
27 //#include "CAFAna/Cuts/SpillCuts.h"
28 #include "CAFAna/Cuts/TruthCuts.h"
33 //#include "CAFAna/Vars/GenieWeights.h"
34 //#include "CAFAna/Vars/NumuVars.h"
35 //#include "CAFAna/Vars/PPFXWeights.h"
36 
37 #include "CAFAna/Core/ISyst.h"
38 #include "CAFAna/Systs/Systs.h"
39 #include "CAFAna/Systs/XSecSysts.h"
40 #include "CAFAna/Systs/RPASysts.h"
41 
42 #include "OscLib/OscCalcPMNSOpt.h"
44 
46 
47 //#include "Weights.h"
48 //#include "DummyExtrapStuff.h"
49 #include "MECTuningUtils.h"
50 #include "TuningSysts2020.h"
51 
52 #include "TF2.h"
53 #include "TStyle.h"
54 #include "TLatex.h"
55 
56 
57 #include <fstream>
58 
59 
60 
62  {
63  return Form( "/nova/ana/users/%s/mec_tuning_2020_systs", getenv( "USER" ) );
64  }
65 
66 
67 
68  const ana::Cut kNumu2020PIDLoosePTP( [](const caf::SRProxy* sr)
69  {
70  return (sr->sel.remid.pid > 0.7 && sr->sel.cvnloosepreselptp.numuid > 0.82);
71  });
72 
73 // use this to plot the weights
74 
75 Double_t func0(Double_t *val, Double_t *par)
76 {
77  Float_t x = val[0];
78  Float_t y = val[1];
79  Double_t f = //par[12] // baseline:
80  + 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] ) )
81  + 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]) )
82  +par[12];
83  if (f<0) f=0;
84  return f;
85 }
86 
87 
89 (
90  const std::string beam="fhc", // fhc,rhc,
91  const std::string fit="doublegaussB", //(doublegauss(B)(B1), gauss, bins)
92  std::string def="mcnp",
93  std::string nonmec_shift = "down", //up, down, blank
94  double chisq=1 //scaling factor for chisquare
95 )
96 {
97  if(def!="mcnp") def = "prod5";
98  const bool isRHC = beam == "rhc";
99 
100 const ana::Cut kSelCuts = ana::kNumu2020ND;
101 
102 // we need to create vector of ISyst before we read predictions (or macro breaks)
103 //------------------------------------ could be more adaptable --------------------------------------------//
104  std::vector<std::unique_ptr<const ana::CompNormSyst>> systs= ana::systsQ3Q0(24, 0.0, 1.2, 24, 0.0, 1.2,
105  ana::GetCutIsFitMEC( isRHC ) && kSelCuts);
106 //std::vector<std::unique_ptr<const ana::CompNormSyst>> systs= ana::systsQ3Q0(20, 0.0, 1.0, 16, 0.0, 0.8,
107  // ana::GetCutIsFitMEC( isRHC ) && kSelCuts);
108  std::vector<const ana::ISyst*> systs_ptrs;
109  if (fit=="bins"){
110  for ( const auto & s : systs ){
111  systs_ptrs.push_back( s.get() );
112  }
113  }
114  if (fit=="gauss"){ systs_ptrs=ana::systs_params_gauss(); std::cout << "systs_ptrs.size() = "<< systs_ptrs.size()<<std::endl;}
115  if (fit=="gausssupp"){ systs_ptrs=ana::systs_params_gauss_supp(); std::cout << "systs_ptrs.size() = "<< systs_ptrs.size()<<std::endl;}
116  if (fit=="doublegauss"){ systs_ptrs=ana::systs_params_double_gauss(); }
117  if (fit=="doublegaussB" || fit=="doublegaussB1"){
118  if (shift_non_mec=="up") systs_ptrs=ana::systs_params_double_gauss_extraUP();
119  else if (shift_non_mec=="down") systs_ptrs=ana::systs_params_double_gauss_extraDOWN();
120  else systs_ptrs=ana::systs_params_double_gauss_extra();
121  }
122 // --------------------------------------------------------------------------------------------------------//
123  // read hists and predictions
124  const std::string out_dir = GetMECTuningDirectory();
125  // If ran on grid, you should place this file in your directory (ana::GetMECTuningDirectory())
126  std::string file_name = "mec_tuning_predictions_"+beam+"_"+fit+"_2020.root";
127  if (def=="mcnp") file_name = "mec_tuning_predictions_"+beam+"_"+fit+"_"+nonmec_shift+"2020_mcnp_genierw.root";
128 
129  file_name=out_dir +"/"+ file_name;
130 
131  std::cout << "reading file from "<<file_name << std::endl;
132 
133  TFile *inFile = new TFile( file_name.c_str(), "READ" );
134 // ana::Spectrum spec_fit_data = *ana::Spectrum::LoadFrom(inFile, "spec_data").release();
135  ana::Spectrum spec_fit_mec_raw=*ana::Spectrum::LoadFrom(inFile, "spec_mc_mec_raw").release();
136  ana::Spectrum spec_fit_non_mec_raw = *ana::Spectrum::LoadFrom(inFile, "spec_mc_nonmec_raw").release();
137  ana::Spectrum spec_fit_mc_raw = *ana::Spectrum::LoadFrom(inFile, "spec_mc_raw").release();
138  inFile->Close();
139 
140  std::string file_name2="/nova/ana/users/mcasales/mec_tuning_2020_systs/mec_tuning_predictions_"+beam+"_doublegaussB_2020_mcnp_half.root";
141  if (def=="mcnp") file_name2="/nova/ana/users/mcasales/mec_tuning_2020_systs/mec_tuning_predictions_"+beam+"_doublegaussB_2020_mcnp_subset.root";
142  TFile *inFile2 = new TFile( file_name2.c_str(), "READ" );
143  ana::Spectrum spec_fit_data = *ana::Spectrum::LoadFrom(inFile2,"spec_data").release();
144  inFile->Close();
145 
146  std::cout<< "pot data: "<< spec_fit_data.POT()<<std::endl;
147  std::cout<< "pot mc: "<< spec_fit_mc_raw.POT()<<std::endl;
148 
149  ana::PredictionInterp * pred_interp= ana::LoadFromFile<ana::PredictionInterp>(file_name, "pred_interp").release();
150  ana::PredictionInterp * pred_interpMEC= ana::LoadFromFile<ana::PredictionInterp>(file_name, "pred_interpMEC").release();
151 
152  ana::BigChi2SingleSampleExperiment expt(pred_interp,spec_fit_data, chisq);
153  //ana::SingleSampleExperiment expt( pred_interp, spec_fit_data );
155  std::vector<const ana::ISyst*> systs_ptrs2 = systs_ptrs;
156  //systs_ptrs2.erase(systs_ptrs2.begin()+6, systs_ptrs2.begin() + 12); //only let peak 1 float (higher q0q3)
157  //systs_ptrs2.erase(systs_ptrs2.begin(), systs_ptrs2.begin() + 6); // only let peak 2 float (low q0q3)
158  //systs_ptrs2.erase(systs_ptrs2.begin()+7, systs_ptrs2.begin() + 8); // remove mean_q0_2
159  //systs_ptrs2.erase(systs_ptrs2.begin()+5, systs_ptrs2.begin() + 6); systs_ptrs2.erase(systs_ptrs2.end()-2, systs_ptrs2.end() -1); //remove correlations
160  //systs_ptrs2.erase(systs_ptrs2.end()-1, systs_ptrs2.end()); //rm baseline
161  //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
162  //systs_ptrs2.erase(systs_ptrs2.begin()+6, systs_ptrs2.begin() + 7); // remove norm 2
163 
164  ana::MinuitFitter fitter( &expt, {}, systs_ptrs2);// , ana::MinuitFitter::FitOpts::kCareful );
165 
166  ana::SystShifts shifts;
167  for ( auto & syst : systs_ptrs2 )
168  {
169  shifts.SetShift( syst, 0 );
170  }
171  double chi= fitter.Fit( &calc, shifts, ana::IFitter::kVerbose )->EvalMetricVal();
172  // the results of the fitting operation are now in the shifts object.
173 
174  //==================================================
175  // Output files
176  std::string fit_results_file_name = out_dir + "/mec_fit_results_" + beam + "_" + fit + "_"+nonmec_shift+"2020_mcnp_test.txt";//_bigchisq.txt";
177  if ( def=="mcnp" ) fit_results_file_name = out_dir + "/mec_fit_results_" + beam + "_" + fit + "_"+nonmec_shift+"2020_mcnp_test.txt";//_bigchisq.txt";
178  std::ofstream fit_results_file( fit_results_file_name.c_str(), std::ofstream::out );
179  for ( auto & syst : systs_ptrs2 )
180  {
181  fit_results_file << syst->ShortName() << " fitted shift = " << shifts.GetShift( syst )<< std::endl;
182  std::cout << syst->ShortName() << " fitted shift = " << shifts.GetShift( syst )<< std::endl;
183  }
184  fit_results_file << "ChiSq scale in fitter =" << chisq << std::endl;
185  fit_results_file << "ChiSq from in fitter = " << chi << std::endl;
186  // Save spectra before and after tune.
187  double pot = spec_fit_data.POT();
188  auto h_data = spec_fit_data.ToTH2( pot );
189  auto h_mc_raw = spec_fit_mc_raw.ToTH2( pot );
190  auto h_non_mec_raw = spec_fit_non_mec_raw.ToTH2( pot );
191  ana::Spectrum spec_fit_mc_tuned( pred_interp->PredictSyst( &calc, shifts ) );
192  auto h_mc_tuned = spec_fit_mc_tuned.ToTH2( pot );
193  ana::Spectrum spec_fit_mc_tunedMEC( pred_interpMEC->PredictSyst( &calc, shifts ) );
194  TH2 * h_mc_tunedMEC = spec_fit_mc_tunedMEC.ToTH2( pot );
195  auto h_mec_raw= spec_fit_mec_raw.ToTH2(pot);
196  TH2 * ratio_mecs = (TH2*)h_mc_tunedMEC->Clone("ratio_mecs");
197  ratio_mecs->Divide(h_mec_raw);
198  new TCanvas;
199  ratio_mecs->SetTitle("ratio raw/tuned pred MEC");
200  ratio_mecs->Draw("colz");
201  std::cout << "here 2"<<std::endl;
202  ana::Spectrum spec_fit_mc_tuned_man( pred_interp->PredictSyst( &calc, shifts ) );
203  auto h_mc_tuned_man = spec_fit_mc_tuned_man.ToTH2( pot );
204 
205  h_mc_tuned->SetTitle("tuned");
206  h_mc_tuned->Draw("colz");
207  new TCanvas;
208  h_mc_raw->SetTitle("raw");
209  h_mc_raw->Draw("colz");
210  new TCanvas; h_data->SetTitle("data");
211  h_data->Draw("colz");
212  new TCanvas; TH2D *ratio = (TH2D*)h_mc_tuned->Clone("ratio");
213  ratio->Divide(h_data);
214  ratio->SetTitle("tuned 2D gauss / data");
215  ratio->Draw("colz");
216  TH1 * datay=(TH1*)h_data->ProjectionY("datay");
217  TH1 * datax=(TH1*)h_data->ProjectionX("datax");
218  TH1 * rawy=(TH1*)h_mc_raw->ProjectionY("rawy"); rawy->SetLineColor(kRed);
219  TH1 * rawx=(TH1*)h_mc_raw->ProjectionX("rawx"); rawx->SetLineColor(kRed);
220  TH1 * tunedy=(TH1*)h_mc_tuned->ProjectionY("tunedy"); tunedy->SetLineColor(kBlue);
221  TH1 * tunedx=(TH1*)h_mc_tuned->ProjectionX("tunedx"); tunedx->SetLineColor(kBlue);
222  TH1 * tunedym=(TH1*)h_mc_tunedMEC->ProjectionY("tunedym"); tunedym->SetLineColor(kGreen+2);
223  TH1 * tunedxm=(TH1*)h_mc_tunedMEC->ProjectionX("tunedxm"); tunedxm->SetLineColor(kGreen+2);
224  TH1 * untunedymfile = (TH1*)h_mec_raw->ProjectionY("untunedymfile");
225  TH1 * untunedxmfile = (TH1*)h_mec_raw->ProjectionX("untunedxmfile");
226 
227  //visaulize dist before and after fit.
228  new TCanvas;
229  datay->SetMarkerStyle(8);
230  datay->Draw("");
231  tunedy->Draw("hist same");
232  tunedym->Draw("hist same");
233  untunedymfile->Draw("hist same");
234  rawy->Draw("hist same");
235  gPad->SaveAs((out_dir +"/mec_tuned_q0_"+beam+"_"+def+"_"+nonmec_shift+".pdf").c_str());
236 
237  new TCanvas;
238  gPad->SetGridy();
239  TH1* ratiorawy = (TH1*)datay->Clone("ratiorawy"); ratiorawy->SetLineColor(kRed);
240  ratiorawy->GetYaxis()->SetTitle("Data/MC");
241  ratiorawy->SetMinimum(0.8);
242  ratiorawy->Divide(rawy);
243  TH1* ratiotunedy = (TH1*)datay->Clone("ratiotunedy"); ratiotunedy->SetLineColor(kBlue);
244  ratiotunedy->Divide(tunedy);
245  ratiorawy->Draw("hist");
246  ratiotunedy->Draw("hist same");
247  gPad->SaveAs((out_dir +"/mec_tuned_q3_"+beam+"_"+def+"_"+nonmec_shift+"_ratio.pdf").c_str());
248 
249  new TCanvas;
250  datax->SetMarkerStyle(8);
251  datax->Draw("");
252  tunedx->Draw("hist same");
253  tunedxm->Draw("hist same");
254  untunedxmfile->Draw("hist same");
255  rawx->Draw("hist same");
256  gPad->SaveAs((out_dir +"/mec_tuned_q3_"+beam+"_"+def+"_"+nonmec_shift+".pdf").c_str());
257 
258  new TCanvas;
259  gPad->SetGridy();
260  TH1* ratiorawx = (TH1*)datax->Clone("ratiorawx"); ratiorawx->SetLineColor(kRed);
261  ratiorawx->GetYaxis()->SetTitle("Data/MC");
262  ratiorawx->SetMinimum(0.8);
263  ratiorawx->Divide(rawx);
264  TH1* ratiotunedx = (TH1*)datax->Clone("ratiotunedx"); ratiotunedx->SetLineColor(kBlue);
265  ratiotunedx->Divide(tunedx);
266  ratiorawx->Draw("hist");
267  ratiotunedx->Draw("hist same");
268  gPad->SaveAs((out_dir +"/mec_tuned_q0_"+beam+"_"+def+"_"+nonmec_shift+"_ratio.pdf").c_str());
269 
270 
271  std::string fit_hists_file_name = out_dir + "/mec_fit_hists_" + beam + "_" +fit+nonmec_shift+ "_2020_chis"+std::to_string(chisq)+".root";//_bigchisq.root";
272  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";
273 
274  TFile fit_hists_file( fit_hists_file_name.c_str(), "recreate" );
275  fit_hists_file.cd();
276  h_non_mec_raw->Write("nonmec_raw");
277  h_data->Write( "data" );
278  h_mc_raw->Write( "mc_raw" );
279  h_mc_tuned->Write( "mc_tuned" );
280  h_mc_tunedMEC->Write("mc_tuned_MEC");
281  h_mec_raw->Write("mc_mec_raw");
282  spec_fit_mc_tunedMEC.SaveTo(& fit_hists_file, "spec_tunedMEC") ;
283  new TCanvas;
284  h_mc_tunedMEC->SetTitle("tuned mec");
285  h_mc_tunedMEC->Draw("colz");
286  new TCanvas;
287  h_mec_raw->SetTitle("raw MEC");// h_mec_raw->SetMaximum(h_mc_tunedMEC->GetMaximum());
288  h_mec_raw->Draw("colz");
289 
290  // for the normalization bins of MEC, save a 2D histo with the weights
291  // ugly device to save weights
292  TH2D* hist_tuned_mec_weights = new TH2D( "tuned_mec_weights", ";True |#vec{q}| (GeV);True q_{0} (GeV);Weight",
293  24, 0., 1.2, 24, 0.0, 1.2 );
294  if (fit=="bins"){
295  std::vector<ana::Q3Q0Bin> binsi;
296  for ( int i = 0; i < 24; ++i )
297  {
298  for ( int j = 0; j < 24; ++j )
299  { double binWidthQ3=1.2/24;
300  double binWidthQ0=1.2/24;
301  double q3 = binWidthQ3 * ( i + 0.5 );
302  double q0 = binWidthQ0 * ( j + 0.5 );
303  if ( q3 * q3 < q0 * q0 ) continue;
304  ana::Q3Q0Bin b;
305  b.loQ3 = q3 - binWidthQ3 / 2;
306  b.hiQ3 = q3 + binWidthQ3 / 2;
307  b.loQ0 = q0 - binWidthQ0 / 2;
308  b.hiQ0 = q0 + binWidthQ0 / 2;
309  binsi.push_back( b );
310  }
311  }
312  std::cout << " binsi.size() " << binsi.size()<< std::endl;
313  for ( std::size_t compIdx = 0; compIdx < binsi.size(); compIdx++ )
314  {
315  // 24, 0.0, 1.2, 24, 0.0, 1.2,
316  // 20, 0.0, 1.0, 16, 0.0, 0.
317  auto b = binsi[ compIdx ];
318  float q3 = ( b.loQ3 + b.hiQ3 ) / 2;
319  float q0 = ( b.loQ0 + b.hiQ0 ) / 2;
320  const auto & syst = systs[ compIdx ];
321  double sigma_scale = syst.get()->GetSigmaScale();
322  double scale_factor = 1 + sigma_scale * shifts.GetShift( syst.get() );
323  double the_shift=shifts.GetShift( syst.get() );
324  if (scale_factor < 0) scale_factor=0;
325  hist_tuned_mec_weights->SetBinContent( hist_tuned_mec_weights->FindFixBin( q3, q0 ), scale_factor );
326  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;
327  }
328  new TCanvas;
329  hist_tuned_mec_weights->Draw("colz");
330  hist_tuned_mec_weights->Write( "mec_weights" );
331  }
332 
333  TH2D* hist_tuned_mec_weights_smoothed = (TH2D*)hist_tuned_mec_weights->Clone( "hist_tuned_mec_weights_smoothed" );
334  hist_tuned_mec_weights_smoothed->Smooth( 1, "k3a" );
335  hist_tuned_mec_weights_smoothed->Write( "mec_weights_smoothed" );
336 
337  //fit_hists_file.Close();
338  if (fit=="doublegauss" || fit=="doublegaussB" || fit=="doublegaussB1"){
339  for ( auto & syst : systs_ptrs2 )
340  {
341  std::cout << syst->ShortName() << " old param = " <<
342  ana::CalcMECDoubleGaussEnhShiftedParam( syst->ShortName(), 0 , "") <<
343  ", new param = " <<
344  ana::CalcMECDoubleGaussEnhShiftedParam( syst->ShortName(), shifts.GetShift( syst ), "")<< std::endl;
345  }
346  }
347 
348  TLatex latex;
349  latex.SetTextSize(0.045);
350  latex.SetTextAlign(13);
351 
352  TF2 *f_old = new TF2("f_old",func0, 0,1.2,0,0.8,13);
353 
354  f_old->SetParameter( 0, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystNorm_1", 0 , ""));
355  f_old->SetParameter( 1, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ0_1", 0 , ""));
356  f_old->SetParameter( 2, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ3_1", 0 , ""));
357  f_old->SetParameter( 3, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ0_1", 0 , ""));
358  f_old->SetParameter( 4, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ3_1", 0 , ""));
359  f_old->SetParameter( 5, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystCorr_1", 0 , ""));
360  f_old->SetParameter( 6, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystNorm_2", 0 , ""));
361  f_old->SetParameter( 7, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ0_2", 0 , ""));
362  f_old->SetParameter( 8, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ3_2", 0 , ""));
363  f_old->SetParameter( 9, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ0_2", 0 , ""));
364  f_old->SetParameter( 10, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ3_2", 0 , ""));
365  f_old->SetParameter( 11, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystCorr_2", 0 , ""));
366  f_old->SetParameter( 12, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystBaseline", 0 , ""));
367 // f->SetContour(48);
368  f_old->SetNpx(200);
369  f_old->SetNpy(200);
370  new TCanvas("c_weights2","c_weights2",600,500);
371  gPad->SetLeftMargin(0.15);
372  gPad->SetBottomMargin(0.15);
373  gPad->SetRightMargin(0.175);
374  f_old->SetMaximum(42);
375  f_old->SetMinimum(-.0001);
376 // f_old->SetMinimum(0);
377  f_old->Draw("colz");
378  f_old->GetZaxis()->SetTitle("MEC Weight");
379  f_old->Draw("colz");
380  f_old->GetYaxis()->SetTitle("True q_{0} (GeV)");
381  f_old->GetXaxis()->SetTitle("True #left|#vec{q}#right| (GeV)");
382  f_old->GetZaxis()->SetTitle("Weights");
383  f_old->GetYaxis()->SetTitleOffset(1.1);
384  f_old->GetXaxis()->SetTitleOffset(1.2);
385  f_old->GetZaxis()->SetTitleOffset(1.15);
386  f_old->GetYaxis()->CenterTitle();
387  f_old->GetXaxis()->CenterTitle();
388  f_old->GetZaxis()->CenterTitle();
389 // f_old->GetXaxis()->SetRangeUser(0,1.2 );
390 // f_old->GetYaxis()->SetRangeUser(0,0.8);
391  latex.DrawLatexNDC(.2,.85,("#color[0]{#nu and #bar{#nu} MEC Weights old}"));
392 // prelim->Draw();
393  gPad->Update();
394  gPad->SaveAs((out_dir +"/mec_weights_official.pdf").c_str());
395 
396 
397  TF2 *f_new = new TF2("f_new",func0, 0,1.2,0,0.8,13);
398 
399  f_new->SetParameter( 0, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystNorm_1", shifts.GetShift(systs_ptrs2[0]) , nonmec_shift));
400  f_new->SetParameter( 1, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ0_1", shifts.GetShift(systs_ptrs2[1]) , nonmec_shift));
401  f_new->SetParameter( 2, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ3_1", shifts.GetShift(systs_ptrs2[2]) , nonmec_shift));
402  f_new->SetParameter( 3, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ0_1", shifts.GetShift(systs_ptrs2[3]) , nonmec_shift));
403  f_new->SetParameter( 4, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ3_1", shifts.GetShift(systs_ptrs2[4]) , nonmec_shift));
404  f_new->SetParameter( 5, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystCorr_1", shifts.GetShift(systs_ptrs2[5]) , nonmec_shift));
405  f_new->SetParameter( 6, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystNorm_2", shifts.GetShift(systs_ptrs2[6]) , nonmec_shift));
406  f_new->SetParameter( 7, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ0_2", shifts.GetShift(systs_ptrs2[7]) , nonmec_shift));
407  f_new->SetParameter( 8, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystMeanQ3_2", shifts.GetShift(systs_ptrs2[8]) , nonmec_shift));
408  f_new->SetParameter( 9, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ0_2", shifts.GetShift(systs_ptrs2[9]) , nonmec_shift));
409  f_new->SetParameter( 10, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystSigmaQ3_2", shifts.GetShift(systs_ptrs2[10]) , nonmec_shift));
410  f_new->SetParameter( 11, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystCorr_2", shifts.GetShift(systs_ptrs2[11]) , nonmec_shift));
411  f_new->SetParameter( 12, ana::CalcMECDoubleGaussEnhShiftedParam( "MECDoubleGaussEnhSystBaseline", shifts.GetShift(systs_ptrs2[12]) , nonmec_shift));
412 
413  f_new->SetNpx(200);
414  f_new->SetNpy(200);
415  new TCanvas("c_weights3","c_weights3",600,500);
416  gPad->SetLeftMargin(0.15);
417  gPad->SetBottomMargin(0.15);
418  gPad->SetRightMargin(0.175);
419  f_new->SetMaximum(42);
420 // f_new->SetMinimum(-.0001);
421  f_new->SetMinimum(-.0001);
422  f_new->Draw("colz");
423  f_new->GetZaxis()->SetTitle("MEC Weight");
424  f_new->Draw("colz");
425  f_new->GetYaxis()->SetTitle("True q_{0} (GeV)");
426  f_new->GetXaxis()->SetTitle("True #left|#vec{q}#right| (GeV)");
427  f_new->GetZaxis()->SetTitle("Weights");
428  f_new->GetYaxis()->SetTitleOffset(1.1);
429  f_new->GetXaxis()->SetTitleOffset(1.2);
430  f_new->GetZaxis()->SetTitleOffset(1.15);
431  f_new->GetYaxis()->CenterTitle();
432  f_new->GetXaxis()->CenterTitle();
433  f_new->GetZaxis()->CenterTitle();
434  f_new->GetXaxis()->SetRangeUser(0,1.2 );
435  f_new->GetYaxis()->SetRangeUser(0,1.2 );
436  latex.DrawLatexNDC(.2,.85,("#color[0]{#nu and #bar{#nu} MEC Weights after this fit}"));
437  //prelim->Draw();
438  gPad->Update();
439  gPad->SaveAs((out_dir +"/mec_weights_"+beam+"_"+def+nonmec_shift+".pdf").c_str());
440 
441 
442 }
caf::Proxy< caf::SRCVNResult > cvnloosepreselptp
Definition: SRProxy.h:1254
void latex(double x, double y, std::string txt, double ang=0, int align=22, double size=0.042)
Implements systematic errors by interpolation between shifted templates.
enum BeamMode kRed
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
double CalcMECDoubleGaussEnhShiftedParam(std::string shift_param, double shift_sigma, std::string shifted)
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
constexpr T pow(T x)
Definition: pow.h:75
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &shift) const override
caf::Proxy< float > pid
Definition: SRProxy.h:1136
Int_t par
Definition: SimpleIterate.C:24
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
std::vector< const ISyst * > systs_params_double_gauss_extraDOWN()
osc::OscCalcDumb calc
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
ifstream inFile
Definition: AnaPlotMaker.h:34
const XML_Char * s
Definition: expat.h:262
expt
Definition: demo5.py:34
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
std::vector< const ISyst * > systs_params_gauss()
std::string getenv(std::string const &name)
std::vector< const ISyst * > systs_params_double_gauss_extra()
std::string GetMECTuningDirectory()
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
#define pot
caf::StandardRecord * sr
const double j
Definition: BetheBloch.cxx:29
caf::Proxy< caf::SRRemid > remid
Definition: SRProxy.h:1269
std::vector< std::unique_ptr< const CompNormSyst > > systsQ3Q0(int nbinsQ3=20, float fitMinQ3=0.0, float fitMaxQ3=1.0, int nbinsQ0=16, float fitMinQ0=0.0, float fitMaxQ0=0.8, Cut cut=GetCutIsFitMEC(false))
double POT() const
Definition: Spectrum.h:227
std::vector< const ISyst * > systs_params_double_gauss()
std::vector< const ISyst * > systs_params_double_gauss_extraUP()
OStream cout
Definition: OStream.cxx:6
void mec_tuning_fitter_2020(const std::string beam="fhc", const std::string fit="doublegaussB", std::string def="mcnp", std::string nonmec_shift="down", double chisq=1)
const Cut GetCutIsFitMEC(const bool isRHC)
Double_t func0(Double_t *val, Double_t *par)
const hit & b
Definition: hits.cxx:21
const ana::Cut kNumu2020PIDLoosePTP([](const caf::SRProxy *sr){return(sr->sel.remid.pid > 0.7 &&sr->sel.cvnloosepreselptp.numuid > 0.82);})
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
enum BeamMode kGreen
enum BeamMode kBlue
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:81
caf::Proxy< float > numuid
Definition: SRProxy.h:907
std::vector< const ISyst * > systs_params_gauss_supp()
float h_data[xbins][ybins]
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string