52 bool fake2017 =
options.Contains(
"fake2017");
53 bool fake2018 =
options.Contains(
"fake2018");
62 std::vector <const IPrediction * > preds;
63 std::vector <std::pair <TH1D*, double > >
cosmics;
64 std::vector <const IExperiment * > expts;
68 double th23, dmsq, dcp;
69 if(fake2017) {th23 = 0.51; dmsq = 2.44e-3; dcp = 1.21;}
71 if(fake2018) {th23 = 0.58; dmsq = 2.51e-3; dcp = 0.17;}
83 preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
85 cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
88 preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
90 cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
103 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");
104 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");
106 TH2F* histNH_pdf =
new TH2F(
"histNH_pdf",
"histNH_pdf", 30, 0, 2*
M_PI, 30, 0.3, 0.7);
107 TH2F* histIH_pdf =
new TH2F(
"histIH_pdf",
"histIH_pdf", 30, 0, 2*
M_PI, 30, 0.3, 0.7);
109 auto histNH = (TH2F*)histFHCNH->Get(
"delta_th23_NH");
110 auto histIH = (TH2F*)histFHCIH->Get(
"delta_th23_IH");
114 for(
int i = 1;
i<=30;
i++){
115 for(
int j = 1;
j<=30;
j++){
116 double binNH = histNH->GetBinContent(
i,
j);
117 histNH_pdf->SetBinContent(
i,
j,
exp(-binNH/2.0));
118 double binIH = histIH->GetBinContent(
i,
j);
119 histIH_pdf->SetBinContent(
i,
j,
exp(-binIH/2.0));
124 histNH_pdf->Draw(
"colz");
125 gPad->Print(
"plots/pdf_transform_hist2d_delta_th23_nh.pdf");
128 histIH_pdf->Draw(
"colz");
129 gPad->Print(
"plots/pdf_transform_hist2d_delta_th23_ih.pdf");
132 cout<<
"integrals: "<<histNH_pdf->Integral()<<
" and ih is "<<histIH_pdf->Integral()<<
endl;
133 double sep = histNH_pdf->Integral()/(histNH_pdf->Integral()+histIH_pdf->Integral());
137 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;
139 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");
141 TH1F* histdmsqNH_pdf =
new TH1F(
"histdmsqNH_pdf",
"histdmsqNH_pdf", 60, 2
e-3, 3
e-3);
142 TH1F* histdmsqIH_pdf =
new TH1F(
"histdmsqIH_pdf",
"histdmsqIH_pdf", 60, -3
e-3, -2
e-3);
144 auto histdmsqNH = (TH1F*)dmsqFHC->Get(
"slice_dmsq_NH");
145 auto histdmsqIH = (TH1F*)dmsqFHC->Get(
"slice_dmsq_IH");
149 for(
int i = 1;
i<=60;
i++){
150 double binNH = histdmsqNH->GetBinContent(
i);
151 histdmsqNH_pdf->SetBinContent(
i,
exp(-binNH/2.0));
152 double binIH = histdmsqIH->GetBinContent(
i);
153 histdmsqIH_pdf->SetBinContent(
i,
exp(-binIH/2.0));
157 histdmsqNH_pdf->Draw(
"hist");
158 gPad->Print(
"plots/pdf_transform_dmsq32_nh.pdf");
161 histdmsqIH_pdf->Draw(
"hist");
162 gPad->Print(
"plots/pdf_transform_dmsq32_ih.pdf");
166 for (
int j = 0;
j<
num;
j++){
175 double hie = rnd.Uniform(1);
176 cout<<
"hie is "<<hie<<
" hie separator is "<<sep<<
endl;
178 if(hie<sep) isNH =
true;
180 double thisdcp, thisth23, thisdmsq, thisth13;
183 histNH_pdf->GetRandom2(thisdcp, thisth23);
184 thisdmsq= histdmsqNH_pdf->GetRandom();
187 histIH_pdf->GetRandom2(thisdcp, thisth23);
188 thisdmsq = histdmsqIH_pdf->GetRandom();
191 thisth13 = rnd.Gaus(0.082, 0.004);
195 calc_fake->SetdCP(thisdcp*
M_PI);
196 calc_fake->SetTh23(
asin(
sqrt(thisth23)));
197 calc_fake->SetDmsq32(thisdmsq);
198 calc_fake->SetTh13(
asin(
sqrt(thisth13))/2);
206 cout<<
"osc pars in exp "<<
j<<
" are "<<thisdmsq<<
" "<<thisth23<<
" "<<thisdcp<<
endl;
215 std::vector <systshifts> shifts;
220 auto sh = rnd.Gaus(0, 1);
221 shifts.push_back({sh,
s->ShortName().c_str()});
222 cout<<
s->ShortName().c_str()<<
"("<<sh<<
"), ";
232 for(
int i=0;
i < 10;
i++){
240 auto check_wo_syst = preds[
i]->Predict(calc_fake).ToTH1(POT);
241 cout<<
" before random syst shifts the integral was "<<check_wo_syst->Integral()<<
endl;
255 for(
auto &sh:shifts){
if(
s->ShortName() == sh.systname) {sshifts.
SetShift(
s, sh.shift);
cout<<
s->ShortName()<<
" "<<sh.shift<<
", ";} }
259 auto temphist = preds[
i]->PredictSyst(calc_fake, sshifts).ToTH1(POT);
261 auto tempcosm = (TH1D*)cosmics[
i].first->Clone(
UniqueName().c_str());
266 cout<<
"numu file was made with "<<denom<<
" livetime, need to scale to the "<<Livetime<<
endl;
267 cout<<
"before scaling "<<tempcosm->Integral()<<
endl;
268 tempcosm->Scale(Livetime/denom);
269 cout<<
"after scaling "<<tempcosm->Integral()<<
endl;
275 cout<<
"mock experiment integral is "<<temphist->Integral()<<
" and cosmic is "<<tempcosm->Integral()<<
endl;
277 auto mockexp =
total.MockData(POT);
278 cout<<
"mock experiment after Pois. fluct "<<mockexp.Integral(POT)<<
endl;
294 for (
int j = 0;
j<
num;
j++){
296 for(
int i = 0;
i < 10; ++
i){
304 auto tempcosm = (TH1D*)cosmics[
i].first->Clone(
UniqueName().c_str());
309 tempcosm->Scale(Livetime/denom);
311 cout<<
i<<
" POT "<<POT<<
" tot MC "<<preds[
i]->Predict(calc_fake).Integral(POT)<<
" cos "<<tempcosm->Integral()<<
" cos er "<<1/
sqrt(tempcosm->Integral())<<
" analyze data "<<thisdata->Integral(POT)<<
endl;
324 std::cout <<
"Creating Multiexperiment with a total of " 325 << expts.size() <<
" experiments\n\n" <<
std::endl;
332 std::cout <<
"Systematics for the fit:\n";
335 std::cout <<
"\n\nSystematic correlations...\n";
340 for(
int i =0;
i < 8; ++
i) exptThis->SetSystCorrelations(
i+2, notfornumu);
347 std::vector <SystShifts> seedShifts = {};
352 std::vector<double> seeds;
357 std::vector <th23helper> th23seeds;
359 bool octantSlice =
true;
367 std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
371 double minchi23 = 1E20;
372 for(
int hie:{-1, 1}){
373 for (
auto & thseed:th23seeds){
375 std::cout <<
"\n\nFinding best fit " << (hie > 0 ?
"NH " :
"IH ")
376 << thseed.label <<
", " 385 auto thisminchi = fit23.
Fit(calc, auxShifts,
388 { thseed.var, thseed.seeds },
392 minchi23=
min(minchi23, thisminchi);
397 TString
str =
"Best fit " ;
398 for (
auto &
v:fitvars){
402 shifts->SetTitle(str);
404 TLine *
l=
new TLine(gPad->GetUxmin(),0,gPad->GetUxmax(),0);
407 (hie>0?
"_NH":
"_IH") + thseed.label +
".pdf");
414 std::cout <<
"\nFound overall minchi " << minchi23 <<
"\n\n";
Cuts and Vars for the 2020 FD DiF Study.
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
TH1 * PlotSystShifts(const SystShifts &shifts, bool sortName)
Simple record of shifts applied to systematic parameters.
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
const double kFutureRHCPOT
const double kFutureFHCPOT
static SystShifts Nominal()
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
const double kFutureFHCLivetime
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
fvar< T > exp(const fvar< T > &x)
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
std::vector< float > Spectrum
const double kFutureRHCLivetime
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
Interface definition for fittable variables.
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
const FitSinSq2Theta13 kFitSinSq2Theta13
std::string to_string(ModuleType mt)
void fake_future_data(int num=1, bool throwexp=true, bool corrSysts=true, TString options="fake2018")
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)
std::string UniqueName()
Return a different string each time, for creating histograms.
const double kAna2018FHCLivetime
const double kAna2018RHCLivetime
Compare a single data spectrum to the MC + cosmics expectation.
Perform MINUIT fits in one or two dimensions.