mec_tuning_preds_2020.C
Go to the documentation of this file.
1 /*
2  * Fill the useful spectra and prediction for mec tuning fit
3  * either with 2D gaussian (MINERvA-style),
4  * two 2D gaussians with / without overall normalization
5  * or lots of normalization bin shifts (like 2018/9 analyses )
6  *
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/Core/Loaders.h"
17 #include "CAFAna/Core/Spectrum.h"
23 #include "CAFAna/Core/Var.h"
24 #include "CAFAna/Cuts/SpillCuts.h"
25 #include "CAFAna/Cuts/TruthCuts.h"
33 
34 #include "CAFAna/Core/ISyst.h"
35 #include "CAFAna/Systs/Systs.h"
36 #include "CAFAna/Systs/XSecSysts.h"
37 #include "CAFAna/Systs/RPASysts.h"
38 #include "CAFAna/Systs/RESSysts.h"
39 
40 #include "OscLib/OscCalcPMNSOpt.h"
42 
44 
47 
48 #include <fstream>
49 
51  {
52  return Form( "/nova/ana/users/%s/mec_tuning_2020_systs", getenv( "USER" ) );
53  }
54 
57  {
58  std::string numu_mec_weights_file_name = "/pnfs/nova/persistent/users/mcasales/prod5/mec_fit_hists_fhc_bins_2020_chis100_ext.root";
59  //miniprod5/mec_fit_hists_fhc_bins_mp5.root" ;
60  std::cout << "Loading FHC MEC weights from file:\n" << numu_mec_weights_file_name << std::endl;
61  TFile *numu_mec_weights_file = TFile::Open( ana::pnfs2xrootd(numu_mec_weights_file_name).c_str());
62  hist_input_mec_weights_numu = (TH2D*)numu_mec_weights_file->Get( "mec_weights_smoothed" );
63  }
64 
65  const ana::Var kWeightTunedNumuMEC( []( const caf::SRProxy * sr )
66  {
67  if ( !ana::kIsDytmanMEC( sr ) || !sr->mc.nu[0].iscc ) return 1.0;
68  if ( sr->mc.nu[0].pdg != 14 ) return 1.0;
69  if ( hist_input_mec_weights_numu == NULL ) throw std::runtime_error( "FHC MEC true (q0,q3) weights histogram is NULL" );
70 
71  double q0 = ana::kTrueQ0( sr );
72  double q3 = ana::kTrueQ3( sr );
73 
74  if ( q0 < hist_input_mec_weights_numu->GetYaxis()->GetXmin() ||
75  q0 > hist_input_mec_weights_numu->GetYaxis()->GetXmax() ||
76  q3 < hist_input_mec_weights_numu->GetXaxis()->GetXmin() ||
77  q3 > hist_input_mec_weights_numu->GetXaxis()->GetXmax() ) return 1.0;
78 
79  double weight = hist_input_mec_weights_numu->Interpolate( q3, q0 );
80  //double weight = hist_input_mec_weights_numu->GetBinContent( hist_input_mec_weights_numu->FindFixBin( q3, q0 ) );
81  if ( weight < 0 ) weight = 0.0;
82  return weight;
83  });
84 
85 
87 (
88  const std::string beam = "fhc", // so far only fhc
89  const std::string shift_non_mec = "", //up, down
90  const std::string fit="doublegaussB", // gauss (for 2D gaussian only), gausssupp (2D gaussian + RPA/QE suppression) , bins (2018-style 200-norm bins)
91  const int test = 0
92 )
93 {
94  if ( beam != "fhc" && beam != "rhc" ) throw std::runtime_error( "Beam must be \"fhc\" or \"rhc\"" );
95  if ( fit!="gausssupp" && fit != "gauss" && fit != "bins" && fit != "doublegauss" && fit != "doublegaussB"
96  && fit != "doublegaussB1" && fit!= "doublegaussSupp" && fit!= "doublegaussSupp2" ) throw std::runtime_error( "Fit must be \"gauss\", \"bins\" or \"doublegauss\" " );
97  const bool isRHC = beam == "rhc";
98  const bool MECbins= isRHC && fit=="bins";
99  if ( MECbins ) LoadWeightsTunedNumuMEC();
100  const bool Gauss2 = ( fit == "doublegauss" || fit== "doublegaussSupp" || fit=="doublegaussSupp2" || fit=="doublegaussB" || fit=="doublegaussB1" );
101  const bool Gauss = fit =="gauss";
102 
103  ana::SystShifts kShiftsNonMEC;
104  if ( shift_non_mec != "down" && shift_non_mec != "up" )
105  {
106  kShiftsNonMEC = ana::kNoShift;
107  }
108  else if ( shift_non_mec == "down" )
109  {
110  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZNormCCQE), +1);
111  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA1CCQE), +1);
112  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA2CCQE), +1);
113  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA3CCQE), +1);
114  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA4CCQE), +1);
115  kShiftsNonMEC.SetShift( &ana::kRPACCQESuppSyst2020, +1 );
116  kShiftsNonMEC.SetShift( &ana::kRPACCQEEnhSyst2020, +1 );
117  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMaCCRES ), -1 );
118  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMvCCRES ), -1 );
119  kShiftsNonMEC.SetShift( &ana::kRESLowQ2SuppressionSyst2020, +1 );
120  }
121  else if ( shift_non_mec == "up" )
122  {
123  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZNormCCQE), -1);
124  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA1CCQE), -1);
125  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA2CCQE), -1);
126  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA3CCQE), -1);
127  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst(rwgt::fReweightZExpA4CCQE), -1);
128  kShiftsNonMEC.SetShift( &ana::kRPACCQESuppSyst2020, -1 );
129  kShiftsNonMEC.SetShift( &ana::kRPACCQEEnhSyst2020, -1 );
130  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMaCCRES ), +1 );
131  kShiftsNonMEC.SetShift( ana::GetGenieKnobSyst( rwgt::fReweightMvCCRES ), +1 );
132  kShiftsNonMEC.SetShift( &ana::kRESLowQ2SuppressionSyst2020, -1 );
133  }
134 
135  else throw std::runtime_error( Form( "Unrecognized shift: \"%s\"", shift_non_mec.c_str() ) );
136 
137  const ana::Var wgt_non_mec = ana::kUnweighted; // * ana::kMinosResSupp;
138  const ana::Var wgt_shift_non_mec = !shift_non_mec.empty() ? ana::GetWeightFromShifts( kShiftsNonMEC ) : ana::kUnweighted;
139  const ana::Var wgt_mec_numu = MECbins ? kWeightTunedNumuMEC : ana::kUnweighted;
140  const ana::Var wgt_mec_2gauss = Gauss2 ? ana::kMECDoubleGaussEnh : ana::kUnweighted;
141  const ana::Var wgt_mec_gauss = Gauss ? ana::kMECGaussEnh : ana::kUnweighted;
142  const ana::Var wgt_fit = ana::kPPFXFluxCVWgt * wgt_non_mec * wgt_shift_non_mec * wgt_mec_numu
143  * wgt_mec_gauss * wgt_mec_2gauss;// * wgt_mec_2gaussB *wgt_mec_2gaussB1;
144 
145  const ana::Cut kSelCuts = ana::kNumu2020ND;
146 
147 // --------------------------- could be more adaptable -----------------------------------------
148  std::vector<const ana::ISyst*> systs_ptrs;
149  // need a less awkward way to declare this // extended range
150  std::vector<std::unique_ptr<const ana::CompNormSyst>> systs= ana::systsQ3Q0(24, 0.0, 1.2, 24, 0.0, 1.2,
151  ana::GetCutIsFitMEC( isRHC ) && kSelCuts);
152  if (fit=="bins"){
153  for ( const auto & s : systs ){
154  systs_ptrs.push_back( s.get() );
155  }
156  }
157  if (fit=="gauss"){ systs_ptrs=ana::systs_params_gauss();}
158  if (fit=="gausssupp"){ systs_ptrs=ana::systs_params_gauss_supp(); }
159  if (fit=="doublegauss"){ systs_ptrs=ana::systs_params_double_gauss(); }
160  if (fit=="doublegaussB" || fit=="doublegaussB1"){systs_ptrs=ana::systs_params_double_gauss_extra();}
161  if (fit== "doublegaussSupp") { systs_ptrs= ana::systs_params_double_gauss_supp();}
162  if (fit== "doublegaussSupp2"){ systs_ptrs= ana::systs_params_double_gauss_supp();}
163 
164  std::cout << "Fitting (systs_ptrs.size()) "<< systs_ptrs.size()<< " systematic parameters" << std::endl;
165 // -------------------------------------------------------------------------------------------
166 // datasets
167  int strided=0;
168  int stridem=0;
169  std::string defn_data = "";
170  std::string defn_mc = "";
171 if ( beam == "fhc" )
172  {
173  defn_data = "defname: prod_caf_R19-11-18-prod5reco.d_nd_numi_fhc_full_v1";
174  defn_mc = "defname: prod_caf_R19-11-18-prod5reco.d_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1";
175  if(fit != "bins" ) strided=0; stridem=5;
176  }
177  else if ( beam == "rhc" )
178  {
179  defn_data = "defname: prod_caf_R19-11-18-prod5reco.g_nd_numi_rhc_full_v1_goodruns";
180  defn_mc = "defname: prod_caf_R19-11-18-prod5reco.d_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1";
181  if (fit !="bins" ) strided=0; stridem=5;
182  }
183 
184  defn_data += test ? " with limit 150" : " with stride "+ std::to_string(strided);
185  defn_mc += test ? " with limit 200" : " with stride "+std::to_string(stridem);
186 
187  if ( defn_data.empty() ) throw std::runtime_error( "Data SAM Definition is empty" );
188  if ( defn_mc .empty() ) throw std::runtime_error( "MC SAM Definition is empty" );
189 
191  loaders.SetLoaderPath( defn_data, caf::kNEARDET, ana::Loaders::kData );
192  loaders.SetLoaderPath( defn_mc , caf::kNEARDET, ana::Loaders::kMC );
194 
195  //Define range for data-mc agreement comparison;
196  double q0_max=0;
197  double q3_max=0;
198  double q0_min=0;
199  double q3_min=0;
200  int binsq0=0;
201  int binsq3=0;
202  if (fit=="bins" ){
203  q3_max=1.4; binsq3=28;
204  q0_max=0.4; binsq0=20;
205  }
206  if (fit=="gauss" || fit=="gausssupp" || fit=="doublegauss" || fit=="doublegaussB" || fit=="doublegaussSupp" || fit=="doublegaussSupp2"){
207  q3_max=1.4; binsq3=28;
208  q0_max=0.4; binsq0=20;//-1;
209  // q0_min=0.02;
210  }
211  // removing first bin of ehad which never agrees >:(
212  if (fit=="doublegaussB1"){
213  q3_max=1.4; binsq3=28;
214  q0_max=0.4; binsq0=20-2;
215  q0_min=0.02*2;
216  }
217  ana::HistAxis axis_fit( "Reco |#vec{q}| (GeV)" , ana::Binning::Simple( binsq3, q3_min, q3_max ), ana::kRecoQmag,
218  "Reco E_{had, vis} (GeV)", ana::Binning::Simple( binsq0, q0_min, q0_max ), ana::kNumuHadVisE );
219 
220  ana::Spectrum spec_fit_data
221  (
223  axis_fit,
224  kSelCuts
225  );
226 
227  ana::Spectrum spec_fit_mc_raw
228  (
230  axis_fit,
231  kSelCuts,
233  wgt_fit
234  );
235 
236  ana::Spectrum spec_fit_non_mec_raw
237  (
239  axis_fit,
240  kSelCuts && !ana::kIsDytmanMEC,
241  ana::kNoShift,
242  wgt_fit
243  );
244 
245  ana::Spectrum spec_fit_mec_raw
246  (
248  axis_fit,
249  kSelCuts && ana::kIsDytmanMEC,
250  ana::kNoShift,
251  wgt_fit
252  );
253 
254  osc::NoOscillations no_osc_calc;
255  ana::NoOscPredictionGenerator pred_gen( loaders.GetLoader(caf::kNEARDET,ana::Loaders::kMC), axis_fit, kSelCuts , wgt_fit);
256  // create a prediction that only has MEC events to see the how MEC only changes
257  ana::NoOscPredictionGenerator pred_genMEC( loaders.GetLoader(caf::kNEARDET,ana::Loaders::kMC), axis_fit, kSelCuts && ana::kIsDytmanMEC , wgt_fit);
258 
259  ana::PredictionInterp pred_interp( systs_ptrs, &no_osc_calc, pred_gen, loaders );
260  ana::PredictionInterp pred_interpMEC( systs_ptrs, &no_osc_calc, pred_genMEC, loaders );
261 
262  loaders.Go();
263 
264  const std::string out_dir = GetMECTuningDirectory();
265  if ( gSystem->AccessPathName( out_dir.c_str() ) ) gSystem->mkdir( out_dir.c_str(), true );
266 
267  std::string file_name = "mec_tuning_predictions_"+beam+"_"+fit+"_"+shift_non_mec+"2020";
268  if(test==1) file_name=out_dir +"/"+ file_name + "_test";
269  file_name=file_name + ".root";
270 
271  TFile file( file_name.c_str(), "recreate" );
272  file.cd();
273  spec_fit_data.SaveTo(& file, "spec_data" ) ;
274  spec_fit_non_mec_raw.SaveTo(& file, "spec_mc_nonmec_raw") ;
275  spec_fit_mec_raw.SaveTo(& file, "spec_mc_mec_raw") ;
276  spec_fit_mc_raw.SaveTo(& file, "spec_mc_raw") ;
277  pred_interp.SaveTo(& file, "pred_interp") ;
278  pred_interpMEC.SaveTo(& file, "pred_interpMEC") ;
279 
280  // store your parameters
281 TH1D * init_h = new TH1D( "initial pars", "Gaussian parameters",
282  systs_ptrs.size() , 0., systs_ptrs.size());
283 
284 if (fit=="doublegaussB1" || fit=="doublegaussB" || fit=="doublegauss" || fit=="doublegaussSupp" || fit=="doublegaussSupp2" ){
285  int bin=1;
286  std::string initial_pars_file_name = out_dir + "/mec_fit_initial_" + beam + "_" + fit + "_"+shift_non_mec+"2020.txt";//_bigchisq.txt";
287  std::ofstream fit_initial_file( initial_pars_file_name.c_str(), std::ofstream::out );
288  for ( auto & syst : systs_ptrs )
289  {
290  fit_initial_file << syst->ShortName() << "old param = " <<
291  ana::CalcMECDoubleGaussEnhShiftedParam( syst->ShortName(), 0 )<< std::endl;
292  init_h->SetBinContent(bin, ana::CalcMECDoubleGaussEnhShiftedParam( syst->ShortName(),0));
293  init_h->GetXaxis()->SetBinLabel(bin,(syst->ShortName()).c_str());
294  bin++;
295  }
296  }
297 
298  init_h->Write("initial_gauss");
299 
300  file.Close();
301  std::cout << "Wrote file: "<< file_name << std::endl;
302 
303 
304 }
Near Detector underground.
Definition: SREnums.h:10
const ana::Var kWeightTunedNumuMEC([](const caf::SRProxy *sr){if(!ana::kIsDytmanMEC(sr)||!sr->mc.nu[0].iscc) return 1.0;if(sr->mc.nu[0].pdg!=14) return 1.0;if(hist_input_mec_weights_numu==NULL) throw std::runtime_error("FHC MEC true (q0,q3) weights histogram is NULL");double q0=ana::kTrueQ0(sr);double q3=ana::kTrueQ3(sr);if(q0< hist_input_mec_weights_numu->GetYaxis() ->GetXmin()|| q0 > hist_input_mec_weights_numu->GetYaxis() ->GetXmax()|| q3< hist_input_mec_weights_numu->GetXaxis() ->GetXmin()|| q3 > hist_input_mec_weights_numu->GetXaxis() ->GetXmax()) return 1.0;double weight=hist_input_mec_weights_numu->Interpolate(q3, q0);if(weight< 0) weight=0.0;return weight;})
Pass neutrinos through unchanged.
Definition: IOscCalc.h:40
Implements systematic errors by interpolation between shifted templates.
std::vector< const ISyst * > systs_params_double_gauss_extra()
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
correl_xv GetYaxis() -> SetDecimals()
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var weight
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const Var kRecoQmag([](const caf::SRProxy *sr){return sqrt(kRecoQ2(sr)+util::sqr(kHadE(sr)));})
Reconstructed three-momentum transfer.
Definition: NumuVars.h:149
std::vector< const ISyst * > systs_params_double_gauss()
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2043
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:570
const Var kMECDoubleGaussEnh([](const caf::SRProxy *sr){if(!kIsNumuCC(sr)||sr->mc.nu[0].mode!=caf::kMEC) return 1.0;double q0=kTrueQ0(sr);double q3=kTrueQ3(sr);return CalcMECDoubleGaussEnh(q0, q3, kGauss2DNorm_1, 0);})
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: Utilities.cxx:620
const Var kTrueQ0
Definition: TruthVars.h:32
const Var kMECGaussEnh([](const caf::SRProxy *sr){if(!kIsNumuCC(sr)||sr->mc.nu[0].mode!=caf::kMEC) return 1.0;double q0=kTrueQ0(sr);double q3=kTrueQ3(sr);return CalcMECGaussEnh(q0, q3, kGauss2DNorm, 0);})
const Var kNumuHadVisE([](const caf::SRProxy *sr){return kNumuHadCalE(sr)+kNumuHadTrkE(sr);})
Definition: NumuVars.h:124
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
std::string GetMECTuningDirectory()
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
const XML_Char * s
Definition: expat.h:262
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
const Var kTrueQ3
Definition: TruthVars.h:38
if(dump)
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_gauss()
void LoadWeightsTunedNumuMEC()
SpectrumLoaderBase & GetLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Retrieve a specific loader.
Definition: Loaders.cxx:129
std::string getenv(std::string const &name)
double CalcMECDoubleGaussEnhShiftedParam(std::string shift_param, double shift_sigma)
caf::StandardRecord * sr
float bin[41]
Definition: plottest35.C:14
const NOvARwgtSyst kRPACCQEEnhSyst2020("RPAShapeenh2020","RPA shape: higher-Q^{2} enhancement (2020)", novarwgt::kRPACCQEEnhSyst2020)
Definition: RPASysts.h:17
const NOvARwgtSyst kRESLowQ2SuppressionSyst2020("LowQ2RESSupp2020","RES low-Q^2 suppression", novarwgt::kRESLowQ2SuppressionSyst2020)
Definition: RESSysts.h:10
const SystShifts kNoShift
Definition: SystShifts.h:115
OStream cout
Definition: OStream.cxx:6
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2055
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void mec_tuning_preds_2020(const std::string beam="fhc", const std::string shift_non_mec="", const std::string fit="doublegaussB", const int test=0)
const Cut GetCutIsFitMEC(const bool isRHC)
const NOvARwgtSyst kRPACCQESuppSyst2020("RPAShapesupp2020","RPA shape: low-Q^{2} suppression (2020)", novarwgt::kRPACCQESuppSyst2020)
Definition: RPASysts.h:18
std::vector< Loaders * > loaders
Definition: syst_header.h:386
TFile * file
Definition: cellShifts.C:17
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:46
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const NOvARwgtSyst * GetGenieKnobSyst(rwgt::ReweightLabel_t knobIdx, std::string altName, std::string altLabel)
Convenience function to get a GENIE knob syst. (Allows using the GENIE knob name & description as the...
Definition: XSecSysts.cxx:119
void SetLoaderPath(const std::string &path, caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Configure loader via wildcard path.
Definition: Loaders.cxx:25
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:100
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:78
std::vector< const ISyst * > systs_params_gauss_supp()
TH2D * hist_input_mec_weights_numu
const Var GetWeightFromShifts(const SystShifts &sys_shifts)