43 TString
options=
"joint_realData_both", TString nametag =
"test_fake2019")
45 bool nueOnly =
options.Contains(
"nueOnly");
46 bool numuOnly =
options.Contains(
"numuOnly");
47 bool joint =
options.Contains(
"joint");
48 assert (nueOnly || numuOnly || joint);
50 bool FHCOnly =
options.Contains(
"FHCOnly");
51 bool RHCOnly =
options.Contains(
"RHCOnly");
52 bool both =
options.Contains(
"both");
53 assert (FHCOnly || RHCOnly || both);
55 bool fake2019 = nametag.Contains(
"fake2019");
68 std::pair <TH1D*, double>
cos;
72 std::vector <predictions> preds;
73 std::vector <predictions> preds_numu_only;
74 std::vector <Spectrum * >
data;
75 std::vector <Spectrum * > data_numu_only;
76 std::vector <const IExperiment * > expts;
77 std::vector <const IExperiment * > expts_numu_only;
80 double th23, dmsq, dcp;
81 if(fake2019) {th23 = 0.565; dmsq = 2.48e-3; dcp = 1.99;}
84 if(FHCOnly || both ) {
87 if(RHCOnly || both ) {
93 if(FHCOnly || both ) {
96 for (
int i = 0;
i< nnumu;
i++ ){
101 if(RHCOnly || both ) {
104 for (
int i = 0;
i< nnumu;
i++ ){
111 TFile* histNH_file =
new TFile(
"/cvmfs/nova.osgstorage.org/analysis/nue/Ana2019/goodness/hist_contours_2019_"+
options+
"_onlyNH_systs.root",
"read");
112 TFile* histIH_file =
new TFile(
"/cvmfs/nova.osgstorage.org/analysis/nue/Ana2019/goodness/hist_contours_2019_"+
options+
"_onlyIH_systs.root",
"read");
114 TH2F* histNH_pdf =
new TH2F(
"histNH_pdf",
"histNH_pdf", 30, 0, 2*
M_PI, 30, 0.3, 0.7);
115 TH2F* histIH_pdf =
new TH2F(
"histIH_pdf",
"histIH_pdf", 30, 0, 2*
M_PI, 30, 0.3, 0.7);
117 auto histNH = (TH2F*)histNH_file->Get(
"delta_th23_NH");
118 auto histIH = (TH2F*)histIH_file->Get(
"delta_th23_IH");
120 for(
int i = 1;
i<=30;
i++){
121 for(
int j = 1;
j<=30;
j++){
122 double binNH = histNH->GetBinContent(
i,
j);
123 histNH_pdf->SetBinContent(
i,
j,
exp(-binNH/2.0));
124 double binIH = histIH->GetBinContent(
i,
j);
125 histIH_pdf->SetBinContent(
i,
j,
exp(-binIH/2.0));
130 histNH_pdf->Draw(
"colz");
131 gPad->Print(outdir+
"pdf_transform_hist2d_delta_th23_nh_"+
options+
"_"+nametag+
".pdf");
134 histIH_pdf->Draw(
"colz");
135 gPad->Print(outdir+
"pdf_transform_hist2d_delta_th23_ih_"+
options+
"_"+nametag+
".pdf");
137 cout<<
"integrals: "<<histNH_pdf->Integral()<<
" and ih is "<<histIH_pdf->Integral()<<
endl;
138 double sep = histNH_pdf->Integral()/(histNH_pdf->Integral()+histIH_pdf->Integral());
140 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;
142 TFile* dmsq_file =
new TFile(
"/cvmfs/nova.osgstorage.org/analysis/nue/Ana2019/goodness/hist_slices_2019_"+
options+
"_systs_dmsq_noOct.root",
"read");
144 TH1F* histdmsqNH_pdf =
new TH1F(
"histdmsqNH_pdf",
"histdmsqNH_pdf", 60, 2
e-3, 3
e-3);
145 TH1F* histdmsqIH_pdf =
new TH1F(
"histdmsqIH_pdf",
"histdmsqIH_pdf", 60, -3
e-3, -2
e-3);
147 auto histdmsqNH = (TH1F*)dmsq_file->Get(
"slice_dmsq_NH");
148 auto histdmsqIH = (TH1F*)dmsq_file->Get(
"slice_dmsq_IH");
150 for(
int i = 1;
i<=60;
i++){
151 double binNH = histdmsqNH->GetBinContent(
i);
152 histdmsqNH_pdf->SetBinContent(
i,
exp(-binNH/2.0));
153 double binIH = histdmsqIH->GetBinContent(
i);
154 histdmsqIH_pdf->SetBinContent(
i,
exp(-binIH/2.0));
158 histdmsqNH_pdf->Draw(
"hist");
159 gPad->Print(outdir+
"pdf_transform_dmsq32_nh_"+
options+
"_"+nametag+
".pdf");
162 histdmsqIH_pdf->Draw(
"hist");
163 gPad->Print(outdir+
"pdf_transform_dmsq32_ih_"+
options+
"_"+nametag+
".pdf");
169 for (
int j = 0;
j<
num;
j++){
178 double hie = rnd.Uniform(1);
179 cout<<
"hie is "<<hie<<
" hie separator is "<<sep<<
endl;
181 if(hie<sep) isNH =
true;
183 double thisdcp, thisth23, thisdmsq, thisth13;
185 histNH_pdf->GetRandom2(thisdcp, thisth23);
186 thisdmsq= histdmsqNH_pdf->GetRandom();
189 histIH_pdf->GetRandom2(thisdcp, thisth23);
190 thisdmsq = histdmsqIH_pdf->GetRandom();
192 thisth13 = rnd.Gaus(0.082, 0.004);
194 calc_fake->SetdCP(thisdcp*
M_PI);
195 calc_fake->SetTh23(
asin(
sqrt(thisth23)));
196 calc_fake->SetDmsq32(thisdmsq);
197 calc_fake->SetTh13(
asin(
sqrt(thisth13))/2);
205 cout<<
"osc pars in exp "<<
j<<
" are "<<thisdmsq<<
" "<<thisth23<<
" "<<thisdcp<<
endl;
211 std::vector <systshifts> shifts;
216 auto sh = rnd.Gaus(0, 1);
217 shifts.push_back({sh,
s->ShortName().c_str()});
218 cout<<
s->ShortName().c_str()<<
"("<<sh<<
"), ";
225 for(
int i = 0;
i < (
int) preds.size(); ++
i){
227 double POT = preds[
i].pot;
237 if(FHCOnly || RHCOnly){
242 for(
auto &sh:shifts){
if(
s->ShortName() == sh.systname) {sshifts.
SetShift(
s, sh.shift);
cout<<
s->ShortName()<<
" "<<sh.shift<<
", ";} }
247 auto tempmock = preds[
i].pred->PredictSyst(calc_fake, sshifts).MockData(POT);
248 auto tempfake = preds[
i].pred->PredictSyst(calc_fake, sshifts).FakeData(POT);
250 tempmock +=
Spectrum(preds[
i].
cos.first, POT,0).MockData(POT);
253 auto thisdatamock =
new Spectrum(tempmock);
254 thisdatamock->SaveTo(d, (
std::to_string(i)+
"_check_mockdata").c_str());
255 auto thisdatafake =
new Spectrum(tempfake);
256 thisdatafake->SaveTo(d, (
std::to_string(i)+
"_check_fakedata").c_str());
257 auto thisdata =thisdatamock;
259 cout<<i<<
" "<<preds[
i].name<<
" POT "<<POT<<
" tot MC "<<preds[
i].pred->Predict(calc_fake).Integral(POT)<<
" cos "<<preds[
i].cos.first->Integral()<<
" cos er "<<preds[
i].cos.second<<
" analyze data "<<thisdata->Integral(POT)<<
endl;
273 std::cout <<
"Adding Dmsq32ConstraintPDG2017\n";
276 std::cout <<
"Creating Multiexperiment with a total of " 277 << expts.size() <<
" experiments\n\n" <<
std::endl;
284 std::cout <<
"Systematics for the fit:\n";
287 std::cout <<
"\n\nSystematic correlations...\n";
288 if(!nueOnly && ! numuOnly && corrSysts){
292 for(
int i =0;
i < nnumu; ++
i) exptThis->SetSystCorrelations(
i+1, notfornumu);
297 for(
int i =0;
i < nnumu; ++
i) exptThis->SetSystCorrelations(
i+1, notfornumu);
303 for(
int i =0;
i < 8; ++
i) exptThis->SetSystCorrelations(
i+2, notfornumu);
307 std::cout <<
"Systematics for the numu only fit:\n";
320 std::vector <SystShifts> seedShifts = {};
328 std::vector<double> seeds;
333 std::vector <th23helper> th23seeds;
339 bool octantSlice =
true;
347 std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
348 std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};
357 double minchi23 = 1E20;
358 for(
int hie:{-1, 1}){
359 for (
auto & thseed:th23seeds){
361 std::cout <<
"\n\nFinding best fit " << (hie > 0 ?
"NH " :
"IH ")
362 << thseed.label <<
", " 370 Fitter fit23(exptThis, fitvars, systs);
371 auto thisminchi = fit23.Fit(calc, auxShifts,
374 { thseed.var, thseed.seeds },
379 minchi23=
min(minchi23, thisminchi);
382 lchi[0] = thisminchi;
384 lchi[2] = thseed.var->GetValue(calc);
386 lchi[4] = kFitSinSq2Theta13.GetValue(calc);
387 lchi.Write((hie>0?
"NH_":
"IH_")+thseed.label+
"_osc_pars");
395 chi.Write(
"best_chi");
398 std::cout <<
"\nFound overall minchi " << minchi23 <<
"\n\n";
406 if(RHCOnly || FHCOnly) {min = 40; max = 140; bins = 100;}
407 TH1F* chisq =
new TH1F(
"chisq",
"best chisq in each mock exp", bins, min, max);
409 if(both) LLreal = 157.146;
410 if(RHCOnly) LLreal = 82.1438;
411 if(FHCOnly) LLreal = 73.2594;
414 outdir =
"/pnfs/nova/scratch/users/lkolupae/2019/goodness/";
416 for (
int j=0;
j<
num;
j++){
419 for (
int k =0; k<5; k++){
421 double minchi = min[0];
422 if(minchi>LLreal) counter++;
428 TLine*
line =
new TLine(LLreal, 0, LLreal, 20);
429 line->SetLineColor(
kRed);
430 cout<<
"chisq intgr "<<chisq->Integral()<<
endl;
432 cout<<
" out of "<<chisq->Integral()<<
endl;
436 gPad->Print(
"exp_"+
options+
".pdf");
Cuts and Vars for the 2020 FD DiF Study.
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
double GetValue(const osc::IOscCalcAdjustable *osc) const override
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Forward to wrapped Var's GetValue()
Simple record of shifts applied to systematic parameters.
const Dmsq32Constraint kDmsq32ConstraintPDG2017(2.45e-3, 0.05e-3, 2.52e-3, 0.05e-3)
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
static SystShifts Nominal()
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
const XML_Char const XML_Char * data
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
fvar< T > exp(const fvar< T > &x)
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
std::vector< float > Spectrum
std::vector< double > POT
static float min(const float a, const float b, const float c)
Combine multiple component experiments.
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
void goodness_of_fit(int num=1, bool throwexp=true, bool corrSysts=false, TString options="joint_realData_both", TString nametag="test_fake2019")
Interface definition for fittable variables.
double GetPOT(bool isFHC=true)
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
Standard interface to all prediction techniques.
const FitSinSq2Theta13 kFitSinSq2Theta13
std::string to_string(ModuleType mt)
std::pair< Spectrum *, double > cos
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::vector< const ISyst * > GetJointFitSystematicList(bool corrSysts, bool nueExclusive=false, bool numuExclusive=false, bool isFHC=true, bool isRHC=true, bool intersection=true, bool ptExtrap=true)
void SetShift(const ISyst *syst, double shift, bool force=false)
Compare a single data spectrum to the MC + cosmics expectation.