mec_nux_tester_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 #include "TCanvas.h"
9 #include "TGraphErrors.h"
10 #include "TH1D.h"
11 #include "TH2D.h"
12 #include "THStack.h"
13 #include "TLegend.h"
14 
15 //#include "CAFAna/Analysis/Ana2017Loaders.h"
16 #include "CAFAna/Fit/MinuitFitter.h"
17 #include "CAFAna/Core/Loaders.h"
18 #include "CAFAna/Core/LoadFromFile.h"
19 #include "CAFAna/Core/Spectrum.h"
20 #include "CAFAna/Core/Ratio.h"
21 #include "CAFAna/Core/SpectrumLoader.h"
22 #include "3FlavorAna/Cuts/NumuCuts.h"
23 #include "3FlavorAna/Cuts/NumuCuts2017.h"
24 #include "3FlavorAna/Cuts/NumuCuts2018.h"
25 #include "3FlavorAna/Cuts/NumuCuts2020.h"
26 #include "CAFAna/Core/Var.h"
27 //#include "CAFAna/Cuts/SpillCuts.h"
28 #include "CAFAna/Cuts/TruthCuts.h"
29 #include "CAFAna/Experiment/SingleSampleExperiment.h"
30 #include "CAFAna/Prediction/PredictionScaleComp.h"
31 #include "CAFAna/Prediction/PredictionInterp.h"
32 #include "CAFAna/Prediction/PredictionNoOsc.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"
43 #include "NumuEnergy/func/Numu2017Fits.h"
44 
45 #include "StandardRecord/Proxy/SRProxy.h"
46 
47 //#include "Weights.h"
48 //#include "DummyExtrapStuff.h"
49 #include "MECTuningUtils.h"
50 #include "TuningSysts2020.h"
51 //#include "DoubleGaussParams.h"
52 
53 #include "TF2.h"
54 #include "TStyle.h"
55 #include "TLatex.h"
56 
57 
58 #include <fstream>
59 
60 std::string GetMECTuningDirectory()
61 {
62  return Form( "/nova/ana/users/%s/mec_tuning_2020", getenv( "USER" ) );
63 }
64 
65 
66 void mec_nux_tester_2020
67  (
68  const std::string beam="fhc",
69  const std::string fit="bins", //(doublegauss(B)(B1), gauss, bins)
70  double chisq=1 //scaling factor for chisquare
71 )
72 {
73  const bool isRHC = beam == "rhc";
74 
75  const ana::Cut kSelCuts = ana::kNumu2020ND;//LoosePTP;
76 
77  // we need to create vector of ISyst before we read predictions (or macro breaks)
78  //------------------------------------ could be more adaptable --------------------------------------------//
79  std::vector<std::unique_ptr<const ana::CompNormSyst>> systs= ana::systsQ3Q0(24, 0.0, 1.2, 24, 0.0, 1.2,
80  ana::GetCutIsFitMEC( isRHC ) && kSelCuts);
81  //std::vector<std::unique_ptr<const ana::CompNormSyst>> systs= ana::systsQ3Q0(20, 0.0, 1.0, 16, 0.0, 0.8,
82  // ana::GetCutIsFitMEC( isRHC ) && kSelCuts);
83  std::vector<const ana::ISyst*> systs_ptrs;
84  if (fit=="bins"){
85  for ( const auto & s : systs ){
86  systs_ptrs.push_back( s.get() );
87  }
88  }
89  if (fit=="doublegaussB" || fit=="doublegaussB1"){systs_ptrs=ana::systs_params_double_gauss_extra();}
90 
91  // --------------------------------------------------------------------------------------------------------//
92  // read hists and predictions
93  const std::string out_dir = GetMECTuningDirectory();
94  // If ran on grid, you should place this file in your directory (ana::GetMECTuningDirectory())
95  std::string file_name = "mec_tuning_predictions_"+beam+"_"+fit+"_2020.root";
96 
97  file_name=out_dir +"/"+ file_name;
98 
99  std::cout << "reading file from "<<file_name << std::endl;
100 
101  TFile *inFile = new TFile( file_name.c_str(), "READ" );
102 
103  // ana::Spectrum spec_fit_data = *ana::Spectrum::LoadFrom(inFile, "spec_mc_mec_raw").release();//*ana::Spectrum::LoadFrom(inFile2,"spec_data").release();
104  ana::Spectrum spec_fit_data = *ana::Spectrum::LoadFrom(inFile, "spec_mc_raw").release();//*ana::Spectrum::LoadFrom(inFile2,"spec_data").release();
105  ana::Spectrum spec_fit_mec_raw=*ana::Spectrum::LoadFrom(inFile, "spec_mc_mec_raw").release();
106  ana::Spectrum spec_fit_non_mec_raw = *ana::Spectrum::LoadFrom(inFile, "spec_mc_nonmec_raw").release();
107  ana::Spectrum spec_fit_mc_raw = *ana::Spectrum::LoadFrom(inFile, "spec_mc_raw").release();
108 
109 
110  inFile->Close();
111 
112  ana::Ratio spec_fit_flat = ana::Ratio(spec_fit_mc_raw, spec_fit_mc_raw);
113 
114  std::cout<< "pot data: "<< spec_fit_data.POT()<<std::endl;
115  std::cout<< "pot mc: "<< spec_fit_mc_raw.POT()<<std::endl;
116 
117  ana::PredictionInterp * pred_interp= ana::LoadFromFile<ana::PredictionInterp>(file_name, "pred_interp").release();
118  ana::PredictionInterp * pred_interpMEC= ana::LoadFromFile<ana::PredictionInterp>(file_name, "pred_interpMEC").release();
119 
120  ana::BigChi2SingleSampleExperiment expt(pred_interp,spec_fit_data, chisq);
121  osc::OscCalcPMNSOpt calc;
122  std::vector<const ana::ISyst*> systs_ptrs2 = systs_ptrs;
123 
124  std::vector<int> bins_vec = {30,10,10,10,10,10,10,10,10,10,10,10,10};
125 
126  //just one floating parameter
127  int w=2000;
128  int h=1000;
129  TCanvas *c1 = new TCanvas("c1","c1",w,h);
130  c1->SetWindowSize( w+(w-c1->GetWw()), h+(h-c1->GetWh()) );
131  c1->Divide(2, 1);
132  double pot = spec_fit_data.POT();
133  auto h_data = spec_fit_data.ToTH2( pot );
134  auto h_mc_raw = spec_fit_mc_raw.ToTH2( pot );
135  auto h_non_mec_raw = spec_fit_non_mec_raw.ToTH2( pot );
136 
137  TVirtualPad *c1_1 = c1->cd(1);
138  c1_1->SetRightMargin(0.15);
139  c1_1->SetLeftMargin(0.15);
140  c1_1->SetBottomMargin(0.15);
141  h_data->SetTitle("Data");
142  h_data->Draw("colz");
143  gPad->Update();
144 
145 
146  TVirtualPad *c1_2 = c1->cd(2);
147  c1_2->SetRightMargin(0.15);
148  c1_2->SetLeftMargin(0.15);
149  c1_2->SetBottomMargin(0.15);
150  h_mc_raw->SetTitle("Nominal MC");
151  h_mc_raw->Draw("colz");
152  gPad->Update();
153 
154  c1->Print("graphs_onefloat.pdf(");
155 
156  TH2 * ratio_nom = (TH2*)h_mc_raw->Clone("ratio_nom");
157  ratio_nom->Divide(h_data);
158 // ratio_nom->SetMinimum(0.99);
159 // ratio_nom->SetMaximum(1.01);
160 
161  int param_limit = 13;
162  std::vector<double> shift_x, chisq_y;
163 
164  for(int param=0;param < param_limit; param++){
165  std::cout<<"Ok, starting in again on param "<<param<<std::endl;
166  TCanvas *c2 = new TCanvas("c2","c2",w,h);
167  c2->SetRightMargin(0.15);
168  c2->SetLeftMargin(0.15);
169  c2->SetBottomMargin(0.15);
170 
171  c2->SetWindowSize( (w+(w-c1->GetWw()))-100, (h+(h-c1->GetWh()))-100 );
172  c2->Divide(2, 1);
173  std::vector<const ana::ISyst*> systs_ptrs_temp = systs_ptrs;
174  systs_ptrs_temp.erase(systs_ptrs_temp.begin()+param+1, systs_ptrs_temp.end());
175  systs_ptrs_temp.erase(systs_ptrs_temp.begin(), systs_ptrs_temp.begin()+param);
176  // systs_ptrs_temp.erase(systs_ptrs_temp.begin(), systs_ptrs_temp.end());
177  // systs_ptrs_temp.erase(systs_ptrs_temp.begin(), systs_ptrs_temp.begin()+1);
178  // systs_ptrs_temp.erase(systs_ptrs_temp.begin()+1, systs_ptrs_temp.end());
179  for (int j=0;j<systs_ptrs_temp.size();j++)
180  std::cout<<"Systs "<<systs_ptrs_temp[j]->ShortName().c_str()<<std::endl;
181  ana::SystShifts shifts_temp;
182  ana::MinuitFitter fitter_temp( &expt, {}, systs_ptrs_temp);
183 
184 
185  //Get the nominal shift function
186  shifts_temp.SetShift(systs_ptrs[param],0);
187  double chitemp = fitter_temp.Fit( &calc, shifts_temp, ana::IFitter::kVerbose )->EvalMetricVal();
188  ana::Spectrum spec_fit_mc_zerotuned( pred_interp->PredictSyst( &calc, shifts_temp ) );
189  auto h_mc_zerotuned = spec_fit_mc_zerotuned.ToTH2( pot );
190 
191 
192  for (int i = -bins_vec[0];i<bins_vec[0];i++){
193 // shifts_temp.SetShift(systs_ptrs_temp[param],i);
194  shifts_temp.SetShift(systs_ptrs[param],i);
195  double chitemp = fitter_temp.Fit( &calc, shifts_temp, ana::IFitter::kVerbose )->EvalMetricVal();
196 
197  chisq_y.push_back(chitemp);
198  shift_x.push_back(i);
199 
200  std::cout<<"I finished the fit for "<<systs_ptrs_temp[param]->ShortName().c_str()<<" shift "<<i<<std::endl;
201 
202  for ( auto & syst : systs_ptrs_temp )
203  {
204  std::cout << syst->ShortName() << "old param = " <<ana::CalcMECDoubleGaussEnhShiftedParam( syst->ShortName(), 0 , "") <<", new param = " <<ana::CalcMECDoubleGaussEnhShiftedParam( syst->ShortName(), shifts_temp.GetShift( syst ), "")<< std::endl;
205  }
206  std::cout<<"I made it through the systematics"<<std::endl;
207 
208  ana::Spectrum spec_fit_mc_tuned( pred_interp->PredictSyst( &calc, shifts_temp ) );
209  auto h_mc_tuned = spec_fit_mc_tuned.ToTH2( pot );
210 
211  TVirtualPad *c2_1 = c2->cd(1);
212  c2_1->SetRightMargin(0.15);
213  c2_1->SetLeftMargin(0.15);
214  c2_1->SetBottomMargin(0.15);
215 
216 
217  h_mc_tuned->SetTitle(Form("Tuned Distribution Shift %s %i",systs_ptrs_temp[param]->ShortName().c_str(),i));
218  gStyle->SetTitleFontSize(0.05);
219  h_mc_tuned->Draw("colz");
220  gPad->Update();
221  std::cout<<"Plotted the Tuned Distribution for "<<systs_ptrs_temp[param]->ShortName().c_str()<<" "<<i<<std::endl;
222 
223  TH2 * ratio = (TH2*)h_mc_tuned->Clone("ratio");
224  ratio->Divide(h_mc_zerotuned);
225 // ratio->SetMinimum(0.99);
226 // ratio->SetMaximum(1.01);
227 
228  TVirtualPad *c2_2 = c2->cd(2);
229  c2_2->SetRightMargin(0.15);
230  c2_2->SetLeftMargin(0.15);
231  c2_2->SetBottomMargin(0.15);
232  ratio->SetTitle(Form("Ratio Tuned to Nominal Shift %i",i));
233  ratio->Draw("colz");
234  gPad->Update();
235 
236  c2->Print("graphs_onefloat.pdf(");
237 
238 // TH2 * double_ratio = (TH2*)ratio->Clone("double_ratio");
239 // double_ratio->Divide(ratio_nom);
240 // double_ratio->SetTitle(Form("Double Ratio Tuned to Nominal Shift %i",i));
241 // double_ratio->Draw("colz");
242 // gPad->Update();
243 //
244 // c2->Print("graphs_onefloat.pdf(");
245  std::cout<<"Plotted the ratio for "<<i<<std::endl;
246  shifts_temp.SetShift(systs_ptrs_temp[param], 0);
247 
248  std::cout<<"reset the shifts to zero "<<i<<std::endl;
249 
250  }
251 
252  TCanvas *c3 = new TCanvas("c3","c3",w/2,h);
253  c3->SetRightMargin(0.15);
254  c3->SetLeftMargin(0.15);
255  c3->SetBottomMargin(0.15);
256  double *x = &shift_x[0];
257  double *y = &chisq_y[0];
258  TGraph* chi2temp = new TGraph(shift_x.size(),x,y);
259  gStyle->SetTitleFontSize(0.08);
260  chi2temp->SetTitle(systs_ptrs_temp[param]->ShortName().c_str());
261  chi2temp->Draw("AC*");
262  c3->Print("graphs_onefloat.pdf(");
263 
264  shift_x.clear();
265  chisq_y.clear();
266  std::cout<<"cleared vectors ready to start again"<<std::endl;
267 
268  }
269  c1->Clear();
270  c1->Print("graphs_onefloat.pdf)");
271 }
272 
273 
274 
275 
ana::PredictionInterp
Implements systematic errors by interpolation between shifted templates.
Definition: PredictionInterp.h:23
ana::Spectrum::ToTH2
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
systs
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
make_mec_shifts_plots.isRHC
string isRHC
Definition: make_mec_shifts_plots.py:100
ana::SystShifts
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
ana::kNumu2020ND
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
ana::PredictionInterp::PredictSyst
Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &shift) const override
Definition: PredictionInterp.cxx:335
ana::CalcMECDoubleGaussEnhShiftedParam
double CalcMECDoubleGaussEnhShiftedParam(std::string shift_param, double shift_sigma, const GaussParams initialParams)
Definition: MECTuningSysts.cxx:162
PandAna.reco_validation.prod5_pid_validation.ratio
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
Definition: prod5_pid_validation.py:554
calc
osc::OscCalcDumb calc
Definition: make_nueFDprediction_kinematics.C:70
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
submit_syst.y
y
Definition: submit_syst.py:249
demo5.c2
c2
Definition: demo5.py:33
inFile
ifstream inFile
Definition: AnaPlotMaker.h:34
s
const XML_Char * s
Definition: expat.h:262
demo5.expt
expt
Definition: demo5.py:34
submit_syst.x
x
Definition: submit_syst.py:249
ana::Spectrum::LoadFrom
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
mec_nux_tester_2020
void mec_nux_tester_2020(const std::string beam="fhc", const std::string fit="bins", double chisq=1)
Definition: mec_nux_tester_2020.C:67
cet::getenv
std::string getenv(std::string const &name)
make_syst_table_plots.h
h
Definition: make_syst_table_plots.py:240
osc::_OscCalcPMNSOpt
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
pot
#define pot
Definition: template_basic.C:49
ana::SystShifts::GetShift
T GetShift(const ISyst *syst) const
ana::systs_params_double_gauss_extra
std::vector< const ISyst * > systs_params_double_gauss_extra(const GaussParams init, std::string suffix)
Definition: MECTuningSysts.cxx:284
calib::j
const double j
Definition: BetheBloch.cxx:29
GetMECTuningDirectory
std::string GetMECTuningDirectory()
Definition: mec_nux_tester_2020.C:60
ana::systsQ3Q0
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))
Definition: TuningSysts2020.h:21
ana::Spectrum::POT
double POT() const
Definition: Spectrum.h:227
ana::Ratio
Represent the ratio between two spectra.
Definition: Ratio.h:8
om::cout
OStream cout
Definition: OStream.cxx:6
MECTuningUtils.h
ana::GetCutIsFitMEC
const Cut GetCutIsFitMEC(const bool isRHC)
Definition: MECTuningUtils.h:211
demo5.c1
c1
Definition: demo5.py:24
TuningSysts2020.h
ana::_Cut< caf::SRProxy >
w
Float_t w
Definition: plot.C:20
ana::SystShifts::SetShift
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:81
h_data
float h_data[xbins][ybins]
Definition: Fittingsidebandfittest.C:43
ana::MinuitFitter
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
string
enum BeamMode string
Generated on Sat Mar 6 2021 14:13:31 for NOvA by   doxygen 1.8.11