10 #include "TGraphErrors.h" 26 #include "CAFAna/Core/Var.h" 61 return Form(
"/nova/ana/users/%s/mec_tuning_2020_systs",
getenv(
"USER" ) );
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]) )
95 if(def!=
"mcnp") def =
"prod5";
106 std::vector<std::unique_ptr<const ana::CompNormSyst>>
systs=
ana::systsQ3Q0(24, 0.0, 1.2, 24, 0.0, 1.2,
110 std::vector<const ana::ISyst*> systs_ptrs;
112 for (
const auto &
s : systs ){
113 systs_ptrs.push_back(
s.get() );
119 if (fit==
"doublegaussB" || fit==
"doublegaussB1"){
132 if (def==
"mcnp") file_name =
"mec_tuning_predictions_"+
beam+
"_"+fit+
"_"+nonmec_shift+
"2020_mcnp_genierw.root";
138 TFile *
inFile =
new TFile( file_name.c_str(),
"READ" );
154 std::vector<const ana::ISyst*> systs_ptrs2 = systs_ptrs;
166 for (
auto & syst : systs_ptrs2 )
175 std::string fit_results_file_name = out_dir +
"/mec_fit_results_" +
beam +
"_" + fit +
"_"+nonmec_shift+
"2020_mcnp_test.txt";
176 if ( def==
"mcnp" ) fit_results_file_name = out_dir +
"/mec_fit_results_" +
beam +
"_" + fit +
"_"+nonmec_shift+
"2020_mcnp_test.txt";
178 for (
auto & syst : systs_ptrs2 )
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;
183 fit_results_file <<
"ChiSq scale in fitter =" << chisq <<
std::endl;
184 fit_results_file <<
"ChiSq from in fitter = " << chi <<
std::endl;
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 );
191 auto h_mc_tuned = spec_fit_mc_tuned.
ToTH2( pot );
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);
198 ratio_mecs->SetTitle(
"ratio raw/tuned pred MEC");
199 ratio_mecs->Draw(
"colz");
202 auto h_mc_tuned_man = spec_fit_mc_tuned_man.
ToTH2( pot );
204 h_mc_tuned->SetTitle(
"tuned");
205 h_mc_tuned->Draw(
"colz");
207 h_mc_raw->SetTitle(
"raw");
208 h_mc_raw->Draw(
"colz");
209 new TCanvas;
h_data->SetTitle(
"data");
211 new TCanvas; TH2D *
ratio = (TH2D*)h_mc_tuned->Clone(
"ratio");
213 ratio->SetTitle(
"tuned 2D gauss / data");
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");
228 datay->SetMarkerStyle(8);
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());
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());
249 datax->SetMarkerStyle(8);
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());
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());
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";
273 TFile fit_hists_file( fit_hists_file_name.c_str(),
"recreate" );
275 h_non_mec_raw->Write(
"nonmec_raw");
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") ;
283 h_mc_tunedMEC->SetTitle(
"tuned mec");
284 h_mc_tunedMEC->Draw(
"colz");
286 h_mec_raw->SetTitle(
"raw MEC");
287 h_mec_raw->Draw(
"colz");
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 );
294 std::vector<ana::Q3Q0Bin> binsi;
295 for (
int i = 0;
i < 24; ++
i )
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;
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 );
312 for ( std::size_t compIdx = 0; compIdx < binsi.size(); compIdx++ )
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;
328 hist_tuned_mec_weights->Draw(
"colz");
329 hist_tuned_mec_weights->Write(
"mec_weights" );
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" );
337 if (fit==
"doublegauss" || fit==
"doublegaussB" || fit==
"doublegaussB1"){
338 for (
auto & syst : systs_ptrs2 )
340 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);
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);
364 f_old->GetZaxis()->SetTitle(
"MEC Weight");
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();
377 latex.DrawLatexNDC(.2,.85,(
"#color[0]{#nu and #bar{#nu} MEC Weights old}"));
380 gPad->SaveAs((out_dir +
"/mec_weights_initial.pdf").c_str());
383 TF2 *f_new =
new TF2(
"f_new",
func0, 0,1.2,0,0.8,13);
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);
407 f_new->SetMinimum(-.0001);
409 f_new->GetZaxis()->SetTitle(
"MEC Weight");
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}"));
425 gPad->SaveAs((out_dir +
"/mec_weights_"+
beam+
"_"+def+nonmec_shift+
".pdf").c_str());
caf::Proxy< caf::SRCVNResult > cvnloosepreselptp
GaussParams InitialProd5UP
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.
Simple record of shifts applied to systematic parameters.
Proxy for caf::StandardRecord.
Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &shift) const override
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))
Representation of a spectrum in any variable, with associated POT.
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
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.
std::vector< const ISyst * > systs_params_double_gauss_extra(const GaussParams init, std::string suffix)
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_supp(const GaussParams init, std::string suffix="")
GaussParams InitialProd5nominal
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);})
GaussParams InitialProd5DOWN
caf::Proxy< caf::SRIDBranch > sel
void SetShift(const ISyst *syst, double shift, bool force=false)
std::string to_string(ModuleType const mt)
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.