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 
53 #include <fstream>
54 
55 
56 
58  {
59  return Form( "/nova/ana/users/%s/mec_tuning_2020", getenv( "USER" ) );
60  }
61 
62 
64 (
65  const std::string beam="fhc",
66  const std::string fit="bins", //(doublegauss(B)(B1), gauss, bins)
67  double chisq=1 //scaling factor for chisquare
68 )
69 {
70  const bool isRHC = beam == "rhc";
71 
72 const ana::Cut kSelCuts = ana::kNumu2020ND;
73 
74 // we need to create vector of ISyst before we read predictions (or macro breaks)
75 //------------------------------------ could be more adaptable --------------------------------------------//
76  std::vector<std::unique_ptr<const ana::CompNormSyst>> systs= ana::systsQ3Q0(24, 0.0, 1.2, 24, 0.0, 1.2,
77  ana::GetCutIsFitMEC( isRHC ) && kSelCuts);
78 //std::vector<std::unique_ptr<const ana::CompNormSyst>> systs= ana::systsQ3Q0(20, 0.0, 1.0, 16, 0.0, 0.8,
79  // ana::GetCutIsFitMEC( isRHC ) && kSelCuts);
80  std::vector<const ana::ISyst*> systs_ptrs;
81  if (fit=="bins"){
82  for ( const auto & s : systs ){
83  systs_ptrs.push_back( s.get() );
84  }
85  }
86  if (fit=="gauss"){ systs_ptrs=ana::systs_params_gauss(); std::cout << "systs_ptrs.size() = "<< systs_ptrs.size()<<std::endl;}
87  if (fit=="gausssupp"){ systs_ptrs=ana::systs_params_gauss_supp(); std::cout << "systs_ptrs.size() = "<< systs_ptrs.size()<<std::endl;}
88  if (fit=="doublegauss"){ systs_ptrs=ana::systs_params_double_gauss(); }
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  ana::Spectrum spec_fit_data = *ana::Spectrum::LoadFrom(inFile, "spec_data").release();
103  ana::Spectrum spec_fit_mec_raw=*ana::Spectrum::LoadFrom(inFile, "spec_mc_mec_raw").release();
104  ana::Spectrum spec_fit_non_mec_raw = *ana::Spectrum::LoadFrom(inFile, "spec_mc_nonmec_raw").release();
105  ana::Spectrum spec_fit_mc_raw = *ana::Spectrum::LoadFrom(inFile, "spec_mc_raw").release();
106  inFile->Close();
107 
108  std::cout<< "pot data: "<< spec_fit_data.POT()<<std::endl;
109  std::cout<< "pot mc: "<< spec_fit_mc_raw.POT()<<std::endl;
110 
111  ana::PredictionInterp * pred_interp= ana::LoadFromFile<ana::PredictionInterp>(file_name, "pred_interp").release();
112  ana::PredictionInterp * pred_interpMEC= ana::LoadFromFile<ana::PredictionInterp>(file_name, "pred_interpMEC").release();
113 
114  ana::BigChi2SingleSampleExperiment expt(pred_interp,spec_fit_data, chisq);
115  //ana::SingleSampleExperiment expt( pred_interp, spec_fit_data );
117  std::vector<const ana::ISyst*> systs_ptrs2 = systs_ptrs;
118  //systs_ptrs2.erase(systs_ptrs2.begin()+6, systs_ptrs2.begin() + 12); //only let peak 1 float (higher q0q3)
119  // systs_ptrs2.erase(systs_ptrs2.begin(), systs_ptrs2.begin() + 6); // only let peak 2 float (low q0q3)
120  // systs_ptrs2.erase(systs_ptrs2.begin()+7, systs_ptrs2.begin() + 8); // remove mean_q0_2
121  //systs_ptrs2.erase(systs_ptrs2.begin()+5, systs_ptrs2.begin() + 6); systs_ptrs2.erase(systs_ptrs2.end()-2, systs_ptrs2.end() -1); //remove correlations
122  //systs_ptrs2.erase(systs_ptrs2.end()-1, systs_ptrs2.end()); //rm baseline
123  //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
124 
125  ana::MinuitFitter fitter( &expt, {}, systs_ptrs2);// , ana::MinuitFitter::FitOpts::kCareful );
126 
127  ana::SystShifts shifts;
128  for ( auto & syst : systs_ptrs2 )
129  {
130  shifts.SetShift( syst, 0 );
131  }
132  double chi= fitter.Fit( &calc, shifts, ana::IFitter::kVerbose )->EvalMetricVal();
133  // the results of the fitting operation are now in the shifts object.
134 
135  //==================================================
136  // Output files
137  std::string fit_results_file_name = out_dir + "/mec_fit_results_" + beam + "_" + fit + "_2020.txt";//_bigchisq.txt";
138  std::ofstream fit_results_file( fit_results_file_name.c_str(), std::ofstream::out );
139  for ( auto & syst : systs_ptrs2 )
140  {
141  fit_results_file << syst->ShortName() << " fitted shift = " << shifts.GetShift( syst )<< std::endl;
142  std::cout << syst->ShortName() << " fitted shift = " << shifts.GetShift( syst )<< std::endl;
143  }
144  fit_results_file << "ChiSq scale in fitter =" << chisq << std::endl;
145  fit_results_file << "ChiSq from in fitter = " << chi << std::endl;
146  // Save spectra before and after tune.
147  double pot = spec_fit_data.POT();
148  auto h_data = spec_fit_data.ToTH2( pot );
149  auto h_mc_raw = spec_fit_mc_raw.ToTH2( pot );
150  auto h_non_mec_raw = spec_fit_non_mec_raw.ToTH2( pot );
151  ana::Spectrum spec_fit_mc_tuned( pred_interp->PredictSyst( &calc, shifts ) );
152  auto h_mc_tuned = spec_fit_mc_tuned.ToTH2( pot );
153  ana::Spectrum spec_fit_mc_tunedMEC( pred_interpMEC->PredictSyst( &calc, shifts ) );
154  TH2 * h_mc_tunedMEC = spec_fit_mc_tunedMEC.ToTH2( pot );
155  auto h_mec_raw= spec_fit_mec_raw.ToTH2(pot);
156  TH2 * ratio_mecs = (TH2*)h_mc_tunedMEC->Clone("ratio_mecs");
157  ratio_mecs->Divide(h_mec_raw);
158  new TCanvas;
159  ratio_mecs->SetTitle("ratio raw/tuned pred MEC");
160  ratio_mecs->Draw("colz");
161  std::cout << "here 2"<<std::endl;
162  ana::Spectrum spec_fit_mc_tuned_man( pred_interp->PredictSyst( &calc, shifts ) );
163  auto h_mc_tuned_man = spec_fit_mc_tuned_man.ToTH2( pot );
164 
165  h_mc_tuned->SetTitle("tuned");
166  h_mc_tuned->Draw("colz");
167  new TCanvas;
168  h_mc_raw->SetTitle("raw");
169  h_mc_raw->Draw("colz");
170  new TCanvas; h_data->SetTitle("data");
171  h_data->Draw("colz");
172  new TCanvas; TH2D *ratio = (TH2D*)h_mc_tuned->Clone("ratio");
173  ratio->Divide(h_data);
174  ratio->SetTitle("tuned 2D gauss / data");
175  ratio->Draw("colz");
176  TH1 * datay=(TH1*)h_data->ProjectionY("datay");
177  TH1 * datax=(TH1*)h_data->ProjectionX("datax");
178  TH1 * rawy=(TH1*)h_mc_raw->ProjectionY("rawy"); rawy->SetLineColor(kRed);
179  TH1 * rawx=(TH1*)h_mc_raw->ProjectionX("rawx"); rawx->SetLineColor(kRed);
180  TH1 * tunedy=(TH1*)h_mc_tuned->ProjectionY("tunedy"); tunedy->SetLineColor(kBlue);
181  TH1 * tunedx=(TH1*)h_mc_tuned->ProjectionX("tunedx"); tunedx->SetLineColor(kBlue);
182  TH1 * tunedym=(TH1*)h_mc_tunedMEC->ProjectionY("tunedym"); tunedym->SetLineColor(kGreen+2);
183  TH1 * tunedxm=(TH1*)h_mc_tunedMEC->ProjectionX("tunedxm"); tunedxm->SetLineColor(kGreen+2);
184  TH1 * untunedymfile = (TH1*)h_mec_raw->ProjectionY("untunedymfile");
185  TH1 * untunedxmfile = (TH1*)h_mec_raw->ProjectionX("untunedxmfile");
186 
187  //visaulize dist before and after fit.
188  new TCanvas;
189  datay->Draw("");
190  tunedy->Draw("hist same");
191  tunedym->Draw("hist same");
192  untunedymfile->Draw("hist same");
193  rawy->Draw("hist same");
194  new TCanvas;
195  datax->Draw("");
196  tunedx->Draw("hist same");
197  tunedxm->Draw("hist same");
198  untunedxmfile->Draw("hist same");
199  rawx->Draw("hist same");
200 
201  std::string fit_hists_file_name = out_dir + "/mec_fit_hists_" + beam + "_" +fit+ "_2020_chis"+std::to_string(chisq)+".root";//_bigchisq.root";
202  TFile fit_hists_file( fit_hists_file_name.c_str(), "recreate" );
203  fit_hists_file.cd();
204  h_non_mec_raw->Write("nonmec_raw");
205  h_data->Write( "data" );
206  h_mc_raw->Write( "mc_raw" );
207  h_mc_tuned->Write( "mc_tuned" );
208  h_mc_tunedMEC->Write("mc_tuned_MEC");
209  h_mec_raw->Write("mc_mec_raw");
210  spec_fit_mc_tunedMEC.SaveTo(& fit_hists_file, "spec_tunedMEC") ;
211  new TCanvas;
212  h_mc_tunedMEC->SetTitle("tuned mec");
213  h_mc_tunedMEC->Draw("colz");
214  new TCanvas;
215  h_mec_raw->SetTitle("raw MEC");// h_mec_raw->SetMaximum(h_mc_tunedMEC->GetMaximum());
216  h_mec_raw->Draw("colz");
217 
218  // for the normalization bins of MEC, save a 2D histo with the weights
219  // ugly device to save weights
220  TH2D* hist_tuned_mec_weights = new TH2D( "tuned_mec_weights", ";True |#vec{q}| (GeV);True q_{0} (GeV);Weight",
221  24, 0., 1.2, 24, 0.0, 1.2 );
222  if (fit=="bins"){
223  std::vector<ana::Q3Q0Bin> binsi;
224  for ( int i = 0; i < 24; ++i )
225  {
226  for ( int j = 0; j < 24; ++j )
227  { double binWidthQ3=1.2/24;
228  double binWidthQ0=1.2/24;
229  double q3 = binWidthQ3 * ( i + 0.5 );
230  double q0 = binWidthQ0 * ( j + 0.5 );
231  if ( q3 * q3 < q0 * q0 ) continue;
232  ana::Q3Q0Bin b;
233  b.loQ3 = q3 - binWidthQ3 / 2;
234  b.hiQ3 = q3 + binWidthQ3 / 2;
235  b.loQ0 = q0 - binWidthQ0 / 2;
236  b.hiQ0 = q0 + binWidthQ0 / 2;
237  binsi.push_back( b );
238  }
239  }
240  std::cout << " binsi.size() " << binsi.size()<< std::endl;
241  for ( std::size_t compIdx = 0; compIdx < binsi.size(); compIdx++ )
242  {
243  // 24, 0.0, 1.2, 24, 0.0, 1.2,
244  // 20, 0.0, 1.0, 16, 0.0, 0.
245  auto b = binsi[ compIdx ];
246  float q3 = ( b.loQ3 + b.hiQ3 ) / 2;
247  float q0 = ( b.loQ0 + b.hiQ0 ) / 2;
248  const auto & syst = systs[ compIdx ];
249  double sigma_scale = syst.get()->GetSigmaScale();
250  double scale_factor = 1 + sigma_scale * shifts.GetShift( syst.get() );
251  double the_shift=shifts.GetShift( syst.get() );
252  if (scale_factor < 0) scale_factor=0;
253  hist_tuned_mec_weights->SetBinContent( hist_tuned_mec_weights->FindFixBin( q3, q0 ), scale_factor );
254  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;
255  }
256  new TCanvas;
257  hist_tuned_mec_weights->Draw("colz");
258  hist_tuned_mec_weights->Write( "mec_weights" );
259  }
260 
261  TH2D* hist_tuned_mec_weights_smoothed = (TH2D*)hist_tuned_mec_weights->Clone( "hist_tuned_mec_weights_smoothed" );
262  hist_tuned_mec_weights_smoothed->Smooth( 1, "k3a" );
263  hist_tuned_mec_weights_smoothed->Write( "mec_weights_smoothed" );
264 
265  //fit_hists_file.Close();
266 
267  if (fit=="doublegauss" || fit=="doublegaussB" || fit=="doublegaussB1"){
268  for ( auto & syst : systs_ptrs2 )
269  {
270  std::cout << syst->ShortName() << "old param = " <<
271  ana::CalcMECDoubleGaussEnhShiftedParam( syst->ShortName(), 0 ) <<
272  ", new param = " <<
273  ana::CalcMECDoubleGaussEnhShiftedParam( syst->ShortName(), shifts.GetShift( syst ))<< std::endl;
274  }
275  }
276 
277 
278 }
Implements systematic errors by interpolation between shifted templates.
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:226
std::vector< const ISyst * > systs_params_double_gauss_extra()
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
std::vector< const ISyst * > systs_params_double_gauss()
TH1 * ratio(TH1 *h1, TH1 *h2)
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &shift) const override
osc::OscCalcDumb calc
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
ifstream inFile
Definition: AnaPlotMaker.h:34
const XML_Char * s
Definition: expat.h:262
expt
Definition: demo5.py:34
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))
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:607
std::vector< const ISyst * > systs_params_gauss()
std::string getenv(std::string const &name)
std::string GetMECTuningDirectory()
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
#define pot
double CalcMECDoubleGaussEnhShiftedParam(std::string shift_param, double shift_sigma)
const double j
Definition: BetheBloch.cxx:29
double POT() const
Definition: Spectrum.h:231
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut GetCutIsFitMEC(const bool isRHC)
const hit & b
Definition: hits.cxx:21
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void mec_tuning_fitter_2020(const std::string beam="fhc", const std::string fit="bins", double chisq=1)
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:78
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