43 TString
options=
"fake2017_part1")
46 bool fake2017 =
options.Contains(
"fake2017");
47 bool realData =
options.Contains(
"realData");
48 bool fake2018 =
options.Contains(
"fake2018");
50 bool part1 =
options.Contains(
"part1");
51 bool part2 =
options.Contains(
"part2");
59 std::vector <const IPrediction * > preds;
60 std::vector <std::pair <TH1D*, double > >
cosmics;
61 std::vector <Spectrum * >
data;
62 std::vector <const IExperiment * > expts;
65 double th23, dmsq, dcp;
66 if(fake2017) {th23 = 0.51; dmsq = 2.44e-3; dcp = 1.21;}
82 TFile* histFHCNH =
new TFile(
"/nova/ana/nu_e_ana/Ana2018/Results/FHCOnly/contours/delta_th23/syst/hist_contours_2018_joint_realData_FHCOnly_onlyNHcombo_systs.root",
"read");
83 TFile* histFHCIH =
new TFile(
"/nova/ana/nu_e_ana/Ana2018/Results/FHCOnly/contours/delta_th23/syst/hist_contours_2018_joint_realData_FHCOnly_onlyIHcombo_systs.root",
"read");
85 TH2F* histNH_pdf =
new TH2F(
"histNH_pdf",
"histNH_pdf", 30, 0, 2*
M_PI, 30, 0.3, 0.7);
86 TH2F* histIH_pdf =
new TH2F(
"histIH_pdf",
"histIH_pdf", 30, 0, 2*
M_PI, 30, 0.3, 0.7);
88 auto histNH = (TH2F*)histFHCNH->Get(
"delta_th23_NH");
89 auto histIH = (TH2F*)histFHCIH->Get(
"delta_th23_IH");
91 for(
int i = 1;
i<=30;
i++){
92 for(
int j = 1;
j<=30;
j++){
93 double binNH = histNH->GetBinContent(
i,
j);
94 histNH_pdf->SetBinContent(
i,
j,
exp(-binNH/2.0));
95 double binIH = histIH->GetBinContent(
i,
j);
96 histIH_pdf->SetBinContent(
i,
j,
exp(-binIH/2.0));
101 histNH_pdf->Draw(
"colz");
102 gPad->Print(
"plots/pdf_transform_hist2d_delta_th23_nh.pdf");
105 histIH_pdf->Draw(
"colz");
106 gPad->Print(
"plots/pdf_transform_hist2d_delta_th23_ih.pdf");
109 cout<<
"integrals: "<<histNH_pdf->Integral()<<
" and ih is "<<histIH_pdf->Integral()<<
endl;
110 double sep = histNH_pdf->Integral()/(histNH_pdf->Integral()+histIH_pdf->Integral());
111 cout<<
"probablilty to get the NH is "<<histNH_pdf->Integral()/(histNH_pdf->Integral()+histIH_pdf->Integral())<<
" and to get the ih is "<<histIH_pdf->Integral()/(histNH_pdf->Integral()+histIH_pdf->Integral())<<
endl;
113 TFile* dmsqFHC =
new TFile(
"/nova/ana/nu_e_ana/Ana2018/Results/FHCOnly/slices/syst/hist_slices_2017_joint_realData_FHCOnlycombo_systs_dmsq_noOct.root",
"read");
115 TH1F* histdmsqNH_pdf =
new TH1F(
"histdmsqNH_pdf",
"histdmsqNH_pdf", 60, 2
e-3, 3
e-3);
116 TH1F* histdmsqIH_pdf =
new TH1F(
"histdmsqIH_pdf",
"histdmsqIH_pdf", 60, -3
e-3, -2
e-3);
118 auto histdmsqNH = (TH1F*)dmsqFHC->Get(
"slice_dmsq_NH");
119 auto histdmsqIH = (TH1F*)dmsqFHC->Get(
"slice_dmsq_IH");
121 for(
int i = 1;
i<=60;
i++){
122 double binNH = histdmsqNH->GetBinContent(
i);
123 histdmsqNH_pdf->SetBinContent(
i,
exp(-binNH/2.0));
124 double binIH = histdmsqIH->GetBinContent(
i);
125 histdmsqIH_pdf->SetBinContent(
i,
exp(-binIH/2.0));
129 histdmsqNH_pdf->Draw(
"hist");
130 gPad->Print(
"plots/pdf_transform_dmsq32_nh.pdf");
133 histdmsqIH_pdf->Draw(
"hist");
134 gPad->Print(
"plots/pdf_transform_dmsq32_ih.pdf");
139 for (
int j = 0;
j<
num;
j++){
146 double hie = rnd.Uniform(1);
147 cout<<
"hie is "<<hie<<
" hie separator is "<<sep<<
endl;
149 if(hie<sep) isNH =
true;
151 double thisdcp, thisth23, thisdmsq, thisth13;
153 histNH_pdf->GetRandom2(thisdcp, thisth23);
154 thisdmsq= histdmsqNH_pdf->GetRandom();
157 histIH_pdf->GetRandom2(thisdcp, thisth23);
158 thisdmsq = histdmsqIH_pdf->GetRandom();
160 thisth13 = rnd.Gaus(0.082, 0.004);
162 calc_fake->SetdCP(thisdcp*
M_PI);
163 calc_fake->SetTh23(
asin(
sqrt(thisth23)));
164 calc_fake->SetDmsq32(thisdmsq);
165 calc_fake->SetTh13(
asin(
sqrt(thisth13))/2);
173 cout<<
"osc pars in exp "<<
j<<
" are "<<thisdmsq<<
" "<<thisth23<<
" "<<thisdcp<<
endl;
181 auto sh = rnd.Gaus(0, 1);
182 cout<<
s->ShortName()<<
"("<<sh<<
"), ";
198 cout<<
"mock components: nue "<< tempnue->Integral()<<
" nue bar "<<tempnuebar->Integral()<<
" tau "<<temptau->Integral()<<
" beam "<<tempbeam->Integral()<<
" nc "<<tempnc->Integral()<<
" cc "<<tempcc->Integral()<<
endl;
215 calc_fake->SetdCP(1.20786*
M_PI);
216 calc_fake->SetTh23(
asin(
sqrt(0.513445)));
217 calc_fake->SetDmsq32(2.46343
e-3);
218 calc_fake->SetTh13(
asin(
sqrt(0.0820218))/2);
225 tot2018bkg->Add(beam2018);
227 tot2018bkg->Add(nc2018);
230 auto Nws2018fhc = ws2018fhc->Integral();
231 auto Ntot2018bkg = tot2018bkg->Integral();
233 cout<<
"2018 ws is "<<ws2018fhc->Integral()<<
" and tot bkg is "<<tot2018bkg->Integral()<<
endl;
235 TH1F* ws =
new TH1F(
"ws",
"ws", 100, 0, 10);
236 TH1F*
bkg =
new TH1F(
"bkg",
"bkg", 100, 0, 10);
238 TH1F* ws_norm =
new TH1F(
"wsnorm",
"ws norm", 100, 0, 10);
239 TH1F* bkg_norm =
new TH1F(
"bkgnorm",
"bkg norm", 100, 0, 10);
241 for (
int j = 0;
j<
num;
j++){
248 thisbkg->Add(thiscc);
250 thisbkg->Add(thisnc);
252 thisbkg->Add(thisbeam);
254 auto Nws = thisws->Integral();
255 auto Nbkg = thisbkg->Integral();
260 ws_norm->Fill(Nws/Nws2018fhc);
261 bkg_norm->Fill(
Nbkg/Ntot2018bkg);
267 gPad->Print((
"plots/events_ws_"+
std::to_string(num)+
".pdf").c_str());
271 gPad->Print((
"plots/events_bkg_"+
std::to_string(num)+
".pdf").c_str());
274 ws_norm->Draw(
"hist");
275 gPad->Print((
"plots/events_ws_norm_"+
std::to_string(num)+
".pdf").c_str());
278 bkg_norm->Draw(
"hist");
279 gPad->Print((
"plots/events_bkg_norm_"+
std::to_string(num)+
".pdf").c_str());
290 TFile*
file =
new TFile(
"throw_exp_nuesignif_fake2017_part1_with_syst_10000.root",
"read");
293 calc_fake->SetdCP(1.20786*
M_PI);
294 calc_fake->SetTh23(
asin(
sqrt(0.513445)));
295 calc_fake->SetDmsq32(2.46343
e-3);
296 calc_fake->SetTh13(
asin(
sqrt(0.0820218))/2);
303 tot2018bkg->Add(beam2018);
305 tot2018bkg->Add(nc2018);
307 tot2018bkg->Add(cc2018);
309 auto Nws2018fhc = ws2018fhc->Integral();
310 auto Ntot2018bkg = tot2018bkg->Integral();
312 cout<<
"2018 ws is "<<ws2018fhc->Integral()<<
" and tot bkg is "<<tot2018bkg->Integral()<<
endl;
314 TH1F* ws_norm =
new TH1F(
"wsnorm",
"ws norm", 50, 0, 10);
315 TH1F* bkg_norm =
new TH1F(
"bkgnorm",
"bkg norm", 100, 0, 10);
317 for (
int j = 0;
j<10000;
j++){
323 thisbkg->Add(thiscc);
325 thisbkg->Add(thisnc);
327 thisbkg->Add(thisbeam);
329 auto Nws = thisws->Integral();
330 auto Nbkg = thisbkg->Integral();
332 ws_norm->Fill(Nws/Nws2018fhc);
333 bkg_norm->Fill(
Nbkg/Ntot2018bkg);
337 TFile* writefile =
new TFile(
filename,
"recreate");
343 tempbkg->Add(tempbeam);
345 tempbkg->Add(tempnc);
347 tempbkg->Add(tempcc);
349 for (
int k = 0; k<
num; k++){
351 auto thistempnue = (TH1*)tempnue->Clone(
UniqueName().c_str());
352 auto thistempbkg = (TH1*)tempbkg->Clone(
UniqueName().c_str());
354 auto Nws = ws_norm->GetRandom();
355 auto Nbkg = bkg_norm->GetRandom();
357 cout<<
" Nws "<<Nws<<
" before scale "<<thistempnue->Integral();
358 thistempnue->Scale(Nws);
359 cout<<
" after scale "<<thistempnue->Integral()<<
endl;
361 cout<<
" Nbkg "<<
Nbkg<<
" before scale "<<thistempbkg->Integral();
362 thistempbkg->Scale(
Nbkg);
363 cout<<
" after scale "<<thistempbkg->Integral()<<
endl;
365 cout<<
"cosmic "<<cosmics[
i].first->Integral()<<
endl;
371 auto mockexp =
total.MockData(POT);
372 cout<<
"mock exp "<<mockexp.Integral(POT)<<
endl;
373 mockexp.SaveTo(writefile, (
std::to_string(k)+
"_mockexp").c_str());
384 calc_fake->SetdCP(1.20786*
M_PI);
385 calc_fake->SetTh23(
asin(
sqrt(0.513445)));
386 calc_fake->SetDmsq32(2.46343
e-3);
387 calc_fake->SetTh13(
asin(
sqrt(0.0820218))/2);
394 tot2018bkg->Add(beam2018);
396 tot2018bkg->Add(nc2018);
398 tot2018bkg->Add(cc2018);
399 tot2018bkg->Add(cosmics[i].first);
400 tot2018bkg->Add(ws2018fhc);
401 cout<<
"total 2018 bkg is "<<tot2018bkg->Integral();
405 auto realdata = data[
i]->ToTH1(POT);
407 realdata->SetLineColor(
kGreen+2);
408 realdata->Draw(
"hist");
409 tot2018bkg->Draw(
"hist same");
410 gPad->Print(
"plots/check_the_2018data_bkg.pdf");
412 cout<<
" total 2018 rhc data is "<<realdata->Integral();
414 cout<<
" likelihood btw real data and fhc predicted bkg is "<<LLreal<<
endl;
416 TH1F*
chi2 =
new TH1F(
"chi2",
"chi2", 40, 0, 40);
417 TH1F*
integral =
new TH1F(
"integrals",
"integrals", 40, 0, 40);
419 for (
int j = 0;
j<
num;
j++){
422 auto integr = thisdata->Integral();
423 integral->Fill(integr);
427 TFile* writefile =
new TFile(
"finalhisti_"+
options+
".root",
"recreate");
428 chi2->Write(
"histchi");
431 cout<<
"integral for chiN > LLreal "<<chi2->Integral((
int)LLreal,40);
434 integral->Draw(
"hist");
439 auto line =
new TLine(LLreal, 0, LLreal, 5000);
Cuts and Vars for the 2020 FD DiF Study.
Simple record of shifts applied to systematic parameters.
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
const double kAna2018RHCPOT
const XML_Char const XML_Char * data
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
Charged-current interactions.
fvar< T > exp(const fvar< T > &x)
void nuebar_signif(int num=1, bool throwexp=true, bool corrSysts=false, TString options="fake2017_part1")
std::vector< float > Spectrum
std::vector< double > POT
Neutral-current interactions.
Both neutrinos and antineutrinos.
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
std::string to_string(ModuleType mt)
All neutrinos, any flavor.
void SetShift(const ISyst *syst, double shift, bool force=false)
std::string UniqueName()
Return a different string each time, for creating histograms.