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