18 #include "CAFAna/nue/Ana2018/joint_fit_2018_tools.h" 43 void POTtag(
bool both,
bool FHCOnly,
bool RHCOnly)
51 ltx1->SetTextAlign(11);
55 ltx->SetTextAlign(11);
61 void MakeMaps ( std::vector <const IFitVar*> profvars,
62 std::map<const IFitVar*, TGraph*> &profMap,
63 std::vector <const ISyst* >
systs,
64 std::map<const ISyst*, TGraph*> &systMap)
67 profMap.insert(std::pair<const IFitVar*, TGraph*> (
var,
new TGraph()));
68 for (
const ISyst* syst : systs)
69 systMap.insert(std::pair<const ISyst*, TGraph*> (syst,
new TGraph()));
73 std::vector <const IFitVar*> profvars,
74 std::vector <TString > profnames,
75 std::map<const IFitVar*, TGraph*> &profMap,
76 std::vector <const ISyst* >
systs,
77 std::map<const ISyst*, TGraph*> &systMap)
79 TDirectory *
tmp = gDirectory;
81 for (
int i = 0;
i < (
int) profvars.size(); ++
i){
82 profMap[profvars[
i]]->Write((prefix +
"_" + profnames[
i]));
84 for (
const ISyst*
s : systs)
85 systMap[
s]->Write((prefix +
"_" +
s->ShortName()));
91 bool corrSysts =
true,
93 bool th23slice =
false,
94 bool octantSlice =
true,
95 bool dmsqSlice =
false,
96 bool th13Slice =
false)
102 bool nueOnly =
options.Contains(
"nueOnly");
103 bool numuOnly =
options.Contains(
"numuOnly");
104 bool joint =
options.Contains(
"joint");
105 assert (nueOnly || numuOnly || joint);
107 bool FHCOnly =
options.Contains(
"FHCOnly");
108 bool RHCOnly =
options.Contains(
"RHCOnly");
109 bool both =
options.Contains(
"both");
110 assert (FHCOnly || RHCOnly || both);
113 if(corrSysts) suffix+=
"_systs";
114 else suffix+=
"_stats";
121 if(th23slice) suffix+=
"_th23";
122 if(dmsqSlice) suffix+=
"_dmsq";
123 if(th13Slice) suffix+=
"_th13";
124 if(!octantSlice)suffix+=
"_noOct";
127 TString
outfilename (outdir +
"hist_slices_2017_" + suffix);
128 TString outFCfilename (outdir +
"slices_FCInput_2018_" + suffix);
129 if(joint && corrSysts && !th23slice && ! dmsqSlice && !th13Slice)
130 outFCfilename = outdir +
"slices_FCInput_2018_"+ suffix;
139 std::vector <const IPrediction * > preds;
140 std::vector <const IPrediction * > preds_numu_only;
141 std::vector <std::pair <TH1D*, double > >
cosmics;
142 std::vector <std::pair <TH1D*, double > > cosmics_numu_only;
143 std::vector <Spectrum * >
data;
144 std::vector <Spectrum * > data_numu_only;
145 std::vector <const IExperiment * > expts;
146 std::vector <const IExperiment * > expts_numu_only;
152 if(FHCOnly || both ) {
156 if(RHCOnly || both ) {
164 if(FHCOnly || both ) {
167 preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
168 preds_numu_only.insert(preds_numu_only.end(),numu_preds.begin(), numu_preds.end());
171 cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
172 cosmics_numu_only.insert(cosmics_numu_only.end(),numu_cosmics.begin(), numu_cosmics.end());
174 if(RHCOnly || both ) {
177 preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
178 preds_numu_only.insert(preds_numu_only.end(),numu_preds.begin(), numu_preds.end());
181 cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
182 cosmics_numu_only.insert(cosmics_numu_only.end(),numu_cosmics.begin(), numu_cosmics.end());
186 cout<<
"\n\n Predictions used for the fit:\n";
188 for(
int i = 0;
i < (
int) preds.size(); ++
i){
193 if(nueOnly || both) {
206 auto thisdata =
GetFakeData (preds[
i],calc_fake, POT, cosmics[i].first);
207 cout<<i<<
" POT "<<POT<<
" tot MC "<<preds[
i]->Predict(calc_fake).Integral(POT)<<
" cos "<<cosmics[
i].first->Integral()<<
" cos er "<<cosmics[
i].second<<
" analyze data "<<thisdata->Integral(POT)<<
endl;
211 cout<<
"Make numu only experiment to get the seeds"<<
endl;
213 for(
int i = 0;
i < (
int) preds_numu_only.size(); ++
i){
221 auto thisdata =
GetFakeData (preds_numu_only[
i],calc_fake, POT, cosmics_numu_only[i].first);
222 cout<<i<<
" POT "<<POT<<
" tot MC "<<preds_numu_only[
i]->Predict(calc_fake).Integral(POT)<<
" cos "<<cosmics_numu_only[
i].first->Integral()<<
" cos er "<<cosmics_numu_only[
i].second<<
" analyze data "<<thisdata->Integral(POT)<<
endl;
223 expts_numu_only.push_back(
new SingleSampleExperiment(preds_numu_only[i],*thisdata, cosmics_numu_only[i].first,cosmics_numu_only[i].
second));
232 else {
std::cout <<
"No reactor constraint, that's th13 slice\n";}
234 std::cout <<
"Adding Dmsq32ConstraintPDG2017\n";
237 std::cout <<
"Creating Multiexperiment with a total of " 238 << expts.size() <<
" experiments\n\n" <<
std::endl;
241 std::cout <<
"Creating Multiexperiment of numu only SimpleExp with a total of " 242 << expts_numu_only.size() <<
" experiments\n\n" <<
std::endl;
248 std::cout <<
"Systematics for the fit:\n";
250 std::vector< const ISyst * >
systs = {};
251 std::vector< const ISyst * > systs_numu_only = {};
270 std::vector <SystShifts> seedShifts = {};
292 cout<<
"------------------- Start prestage seeds --------------------------"<<
endl;
294 double minchi_numu = 1E20;
295 double pre_seed_th23;
296 double pre_seed_dmsq;
300 double maxmix = 0.514;
301 double numu_pre_seedLONH, numu_pre_seedUONH, numu_pre_seedLOIH, numu_pre_seedUOIH, dmsq_numu_pre_seedNH, dmsq_numu_pre_seedIH;
304 std::cout <<
"\n\nFinding test best fit " << (hie>0?
"NH " :
"IH ") <<
"\n\n";
307 auto thisminchi = fitnumu_only.Fit(calc, auxShifts ,
312 if(thisminchi < minchi_numu) minchi_numu = thisminchi;
315 pre_seed_dmsq = kFitDmSq32.GetValue(calc);
318 numu_pre_seedLONH = ((pre_seed_th23>maxmix)?(2*maxmix-pre_seed_th23):pre_seed_th23);
319 numu_pre_seedUONH = ((pre_seed_th23>maxmix)?pre_seed_th23:(2*maxmix-pre_seed_th23));
324 cout<<
"------------------- End prestage seeds --------------------------"<<
endl;
328 std::vector<double> seeds;
334 std::vector <th23helper> th23seeds;
340 else th23seeds.push_back({ {numu_pre_seedLONH, 0.5, numu_pre_seedUONH}, &
kFitSinSqTheta23,
""});
342 std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
343 std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};
350 double minchi23 = 1E20;
351 for(
int hie:{-1, 1}){
352 for (
auto & thseed:th23seeds){
354 std::cout <<
"\n\nFinding best fit " << (hie > 0 ?
"NH " :
"IH ")
355 << thseed.label <<
", " 365 Fitter fit23(exptThis, fitvars, systs);
366 cout<<
" pre dmsq seed is "<<pre_seed_dmsq<<
endl;
367 auto thisminchi = fit23.Fit(calc, auxShifts,
370 { thseed.var, thseed.seeds },
374 minchi23=
min(minchi23, thisminchi);
379 TString
str =
"Best fit " ;
380 for (
auto &
v:fitvars){
384 shifts->SetTitle(str);
386 TLine *
l=
new TLine(gPad->GetUxmin(),0,gPad->GetUxmax(),0);
388 gPad->Print(outdir +
"debug_slice_shifts_bestfit_" + suffix +
389 (hie>0?
"NH":
"IH") + thseed.label +
".pdf");
396 std::cout <<
"\nFound overall minchi " << minchi23 <<
"\n\n";
402 TFile *
outfile =
new TFile(outfilename+
".root",
"recreate");
403 TFile * outFCfile =
new TFile(outFCfilename+
".root",
"recreate");
408 v.Write(
"overall_minchi");
413 for(
int hie: {-1, +1}){
418 std::cout <<
"Starting profile " << (hie>0?
"NH ":
"IH") <<
"\n\n";
424 std::vector <const IFitVar * > profvars;
425 std::vector <TString> profvarnames;
426 std::map <const IFitVar *, std::vector <double> >profseeds;
431 for(
auto const & thseed:th23seeds){
433 const IFitVar* kCurTheta23 = thseed.var;
436 if(!th23slice && !dmsqSlice && !th13Slice){
441 profdef.profvarnames = {
"DmSq32",
"SinSq2Theta13",
"SinSqTheta23"};
442 profdef.profseeds = {{kCurTheta23, thseed.seeds}, { &
kFitDmSq32, {hie*pre_seed_dmsq}} };
443 profdef.shortname=
"delta";
449 profdef.minX = (nueOnly ? 0:0.3);
450 profdef.maxX = (nueOnly ? 1:0.7);
452 profdef.profvarnames = {
"DmSq32",
"SinSq2Theta13",
"DeltaCPpi"};
454 profdef.shortname=
"th23";
462 profdef.profvarnames = {
"DmSq32",
"SinSqTheta23",
"DeltaCPpi"};
464 profdef.shortname=
"th13";
469 profdef.minX = hie > 0 ? 2
e-3 : -3
e-3;
470 profdef.maxX = hie > 0 ? 3
e-3 : -2
e-3;
472 profdef.profvarnames = {
"SinSqTheta23",
"SinSq2Theta13",
"DeltaCPpi"};
473 profdef.profseeds = { {kCurTheta23, thseed.seeds},
475 profdef.shortname=
"dmsq";
478 std::map<const IFitVar*, TGraph*> profVarsMap;
479 std::map<const ISyst*, TGraph*> systMap;
480 MakeMaps (profdef.profvars, profVarsMap, systs, systMap);
483 std::cout <<
" Profile "<< thseed.label <<
"\n";
486 profdef.fitvar, steps, profdef.minX, profdef.maxX,
492 profVarsMap, systMap);
494 TString hieroctstr = (hie > 0 ?
"NH" :
"IH") + thseed.label;
497 slice->Write((TString )
"slice_" + profdef.shortname +
"_" + hieroctstr);
499 profdef.profvars, profdef.profvarnames, profVarsMap, systs, systMap);
518 TFile *
infile =
new TFile (outfilename+
".root",
"read");
521 auto mins =* (
TVectorD*)infile->Get(
"overall_minchi");
522 double minchi23 = mins[0];
524 std::vector <TString> sliceNames = {};
525 if (!th23slice && ! dmsqSlice && !th13Slice) sliceNames = {
"slice_delta_"};
526 if (th23slice) sliceNames = {
"slice_th23_"};
527 if (dmsqSlice) sliceNames = {
"slice_dmsq_"};
528 if (th13Slice) sliceNames = {
"slice_th13_"};
533 std::vector<double> seeds;
537 std::vector <th23helper> th23seeds;
547 for(TString sliceName:sliceNames){
548 std::vector < std::pair <TGraph*, TString> >
graphs;
550 if(sliceName.Contains(
"delta")) {
551 for (TString hie:{
"NH",
"IH"}){
552 for(
auto const & thseed:th23seeds){
553 auto h = (TH1D*) infile ->
Get( sliceName + hie + thseed.label);
555 std::cerr <<
"Problem " << sliceName + hie + thseed.label <<
"\n";
558 graphs.push_back({
new TGraph(
h), hie +
" " + thseed.label});
565 for (TString hie:{
"NH",
"IH"}) {
566 for(
auto const & thseed:th23seeds){
567 auto h = (TH1D*) infile ->
Get( sliceName + hie + thseed.label);
568 if(!
h) {
std::cerr <<
"Problem " << sliceName + hie + thseed.label <<
"\n"; }
569 graphs.push_back({
new TGraph(
h),hie +
" " + thseed.label});
573 if(sliceName.Contains(
"th23")) {
575 (nueOnly ? 0 : 0.32), (nueOnly ? 1 : 0.72));
578 if(sliceName.Contains(
"dmsq")) {
581 if(sliceName.Contains(
"th13")) {
587 for (
auto &gr:graphs) {
588 if(gr.second.Contains(
"IH")) gr.first->SetLineColor(
kPrimColorIH);
590 if(gr.second.Contains(
"LO")) gr.first->SetLineStyle(7);
591 gr.first->SetLineWidth(3);
593 for (
int i =0;
i < gr.first->GetN(); ++
i){
594 gr.first->SetPoint(
i,1000*
fabs(gr.first->GetX()[
i]),gr.first->GetY()[
i]);
597 gr.first->Draw(
"l same");
600 POTtag(both, FHCOnly, RHCOnly);
604 auto leg =
SliceLegend(graphs, !dmsqSlice && !th23slice && !th13Slice, 0.7, dmsqSlice);
609 for(TString
ext: {
".pdf",
".root"}) {
610 gPad->Print(outdir +
"slice_" + suffix +
621 TFile * inFCfile =
new TFile(outFCfilename+
".root",
"read");
624 std::vector <TString> strprof = {
"DmSq32",
"SinSq2Theta13",
"SinSqTheta23"};
625 if(th23slice) strprof = {
"DmSq32",
"SinSq2Theta13",
"DeltaCPpi"};
626 if(th13Slice) strprof = {
"DmSq32",
"DeltaCPpi",
"SinSqTheta23"};
627 if(dmsqSlice) strprof = {
"SinSqTheta23",
"SinSq2Theta13",
"DeltaCPpi"};
629 std::vector< const ISyst * >
systs = {};
634 for (
const ISyst*
s : systs){
635 strprof.push_back(
s->ShortName());
638 for(
auto const &
str:strprof){
640 double miny=-1.5,
maxy=1.5;
641 if(
str==
"DeltaCPpi") {miny=0;
maxy=2;}
642 if(
str==
"DmSq32") {miny=2.2;
maxy=2.9;}
643 if(
str==
"SinSq2Theta13") {miny=0.080;
maxy=0.086;}
644 if(
str==
"SinSqTheta23") {miny=0.3;
maxy=0.8;}
646 auto axes =
new TH2F(
"",
";#delta_{CP}/#pi;" +
str, 100, 0, 2, 100, miny,
maxy);
647 if(dmsqSlice)
axes =
new TH2F(
"",
";|#Deltam^{2}_{32}| (10^{-3} eV^{2});" +
str, 100, 2.2, 2.8, 100, miny,
maxy);
648 if(th23slice)
axes =
new TH2F(
"",
";sin^{2}#theta_{23};" +
str, 100, 0.3, 0.7, 100, miny,
maxy);
649 if(th13Slice)
axes =
new TH2F(
"",
";sin^{2}#theta_{13};" +
str, 100, 0.01, 0.2, 100, miny,
maxy);
652 for (TString hie:{
"NH",
"IH"}){
653 for(
auto const & thseed:th23seeds){
654 auto gr = (TGraph*) inFCfile->Get(hie + thseed.label +
"_" +
str);
656 std::cerr <<
"Problem " << hie + thseed.label +
"_" +
str <<
"\n";
660 for (
int i = 0;
i < (
int) gr->GetN(); ++
i){
661 gr->SetPoint(
i, gr->GetX()[
i], 1000*
fabs(gr->GetY()[
i]));
665 for (
int i = 0;
i < (
int) gr->GetN(); ++
i){
666 gr->SetPoint(
i, 1000*
fabs(gr->GetX()[
i]), gr->GetY()[
i]);
670 if(thseed.label.Contains(
"LO")) gr->SetLineStyle(7);
672 gr->SetMarkerColor(gr->GetLineColor());
677 gPad->Print(outdir +
"debug_slice_prof_" + suffix +
"_"+
str +
".pdf");
const Color_t kPrimColorIH
void DrawSliceCanvas(TString slicename, const double ymax, const double xmin=0, const double xmax=2.)
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)
fvar< T > fabs(const fvar< T > &x)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
const FitDmSq32 kFitDmSq32
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
void SaveMaps(TDirectory *out, std::map< std::string, IDecomp * > decomps_nominal, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > decomps_shifted, std::map< std::string, PredictionNoExtrap * > predsNE_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionNoExtrap * > > > predsNE_shifted, std::map< std::string, PredictionSterile * > predsSt_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionSterile * > > > predsSt_shifted)
Save all of the objects in the input maps to the out directory/file.
void MakeMaps(std::vector< const IFitVar * > profvars, std::map< const IFitVar *, TGraph * > &profMap, std::vector< const ISyst * > systs, std::map< const ISyst *, TGraph * > &systMap)
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
static SystShifts Nominal()
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Encapsulate code to systematically shift a caf::SRProxy.
virtual void SetDmsq32(const T &dmsq32)=0
const double kAna2018RHCPOT
const XML_Char const XML_Char * data
void demoFitSlices(bool createFile=false, bool corrSysts=true, TString options="joint_both", bool th23slice=false, bool octantSlice=true, bool dmsqSlice=false, bool th13Slice=false)
string outfilename
knobs that need extra care
const DummyAnaSyst kAnaCalibrationSyst("Calibration","AbsCalib")
const Color_t k3SigmaConfidenceColorNH
void POTtag(bool both, bool FHCOnly, bool RHCOnly)
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
virtual T GetDmsq32() const
TLegend * SliceLegend(std::vector< std::pair< TGraph *, TString > > &graphs, bool isDelta)
std::vector< double > POT
static float min(const float a, const float b, const float c)
Combine multiple component experiments.
const double kAna2017FullDetEquivPOT
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
const Float_t kBlessedLabelFontSize
const FitSinSq2Theta13 kFitSinSq2Theta13
const double kAna2018FHCPOT
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
Compare a single data spectrum to the MC + cosmics expectation.
TGraph * SqrtProfile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double minchi, std::vector< const IFitVar * > profVars, std::vector< const ISyst * > profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &systsMap, MinuitFitter::FitOpts opts)
Forward to Profile but sqrt the result for a crude significance.