10 #include "TGraphErrors.h" 26 #include "CAFAna/Core/Var.h" 63 return Form(
"/nova/ana/users/%s/mec_tuning_2020_systs",
getenv(
"USER" ) );
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]) )
97 if(def!=
"mcnp") def =
"prod5";
104 std::vector<std::unique_ptr<const ana::CompNormSyst>>
systs=
ana::systsQ3Q0(24, 0.0, 1.2, 24, 0.0, 1.2,
108 std::vector<const ana::ISyst*> systs_ptrs;
110 for (
const auto &
s : systs ){
111 systs_ptrs.push_back(
s.get() );
117 if (fit==
"doublegaussB" || fit==
"doublegaussB1"){
127 if (def==
"mcnp") file_name =
"mec_tuning_predictions_"+
beam+
"_"+fit+
"_"+nonmec_shift+
"2020_mcnp_genierw.root";
133 TFile *
inFile =
new TFile( file_name.c_str(),
"READ" );
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" );
155 std::vector<const ana::ISyst*> systs_ptrs2 = systs_ptrs;
167 for (
auto & syst : systs_ptrs2 )
176 std::string fit_results_file_name = out_dir +
"/mec_fit_results_" +
beam +
"_" + fit +
"_"+nonmec_shift+
"2020_mcnp_test.txt";
177 if ( def==
"mcnp" ) fit_results_file_name = out_dir +
"/mec_fit_results_" +
beam +
"_" + fit +
"_"+nonmec_shift+
"2020_mcnp_test.txt";
179 for (
auto & syst : systs_ptrs2 )
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;
184 fit_results_file <<
"ChiSq scale in fitter =" << chisq <<
std::endl;
185 fit_results_file <<
"ChiSq from in fitter = " << chi <<
std::endl;
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 );
192 auto h_mc_tuned = spec_fit_mc_tuned.
ToTH2( pot );
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);
199 ratio_mecs->SetTitle(
"ratio raw/tuned pred MEC");
200 ratio_mecs->Draw(
"colz");
203 auto h_mc_tuned_man = spec_fit_mc_tuned_man.
ToTH2( pot );
205 h_mc_tuned->SetTitle(
"tuned");
206 h_mc_tuned->Draw(
"colz");
208 h_mc_raw->SetTitle(
"raw");
209 h_mc_raw->Draw(
"colz");
210 new TCanvas;
h_data->SetTitle(
"data");
212 new TCanvas; TH2D *
ratio = (TH2D*)h_mc_tuned->Clone(
"ratio");
214 ratio->SetTitle(
"tuned 2D gauss / data");
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");
229 datay->SetMarkerStyle(8);
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());
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());
250 datax->SetMarkerStyle(8);
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());
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());
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";
274 TFile fit_hists_file( fit_hists_file_name.c_str(),
"recreate" );
276 h_non_mec_raw->Write(
"nonmec_raw");
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") ;
284 h_mc_tunedMEC->SetTitle(
"tuned mec");
285 h_mc_tunedMEC->Draw(
"colz");
287 h_mec_raw->SetTitle(
"raw MEC");
288 h_mec_raw->Draw(
"colz");
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 );
295 std::vector<ana::Q3Q0Bin> binsi;
296 for (
int i = 0;
i < 24; ++
i )
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;
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 );
313 for ( std::size_t compIdx = 0; compIdx < binsi.size(); compIdx++ )
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;
329 hist_tuned_mec_weights->Draw(
"colz");
330 hist_tuned_mec_weights->Write(
"mec_weights" );
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" );
338 if (fit==
"doublegauss" || fit==
"doublegaussB" || fit==
"doublegaussB1"){
339 for (
auto & syst : systs_ptrs2 )
341 std::cout << syst->ShortName() <<
" old param = " <<
349 latex.SetTextSize(0.045);
350 latex.SetTextAlign(13);
352 TF2 *f_old =
new TF2(
"f_old",
func0, 0,1.2,0,0.8,13);
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);
378 f_old->GetZaxis()->SetTitle(
"MEC Weight");
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();
391 latex.DrawLatexNDC(.2,.85,(
"#color[0]{#nu and #bar{#nu} MEC Weights old}"));
394 gPad->SaveAs((out_dir +
"/mec_weights_official.pdf").c_str());
397 TF2 *f_new =
new TF2(
"f_new",
func0, 0,1.2,0,0.8,13);
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);
421 f_new->SetMinimum(-.0001);
423 f_new->GetZaxis()->SetTitle(
"MEC Weight");
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}"));
439 gPad->SaveAs((out_dir +
"/mec_weights_"+
beam+
"_"+def+nonmec_shift+
".pdf").c_str());
caf::Proxy< caf::SRCVNResult > cvnloosepreselptp
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.
double CalcMECDoubleGaussEnhShiftedParam(std::string shift_param, double shift_sigma, std::string shifted)
Simple record of shifts applied to systematic parameters.
Proxy for caf::StandardRecord.
Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &shift) const override
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
std::vector< const ISyst * > systs_params_double_gauss_extraDOWN()
Representation of a spectrum in any variable, with associated POT.
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
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)
Optimized version of OscCalcPMNS.
caf::Proxy< caf::SRRemid > remid
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()
std::vector< const ISyst * > systs_params_double_gauss_extraUP()
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 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
std::string to_string(ModuleType mt)
void SetShift(const ISyst *syst, double shift, bool force=false)
caf::Proxy< float > numuid
std::vector< const ISyst * > systs_params_gauss_supp()
float h_data[xbins][ybins]
Perform MINUIT fits in one or two dimensions.