3 #include "CAFAna/Analysis/Fit.h" 5 #include "CAFAna/Analysis/NuePlotStyle.h" 6 #include "CAFAna/Analysis/Surface.h" 13 #include "CAFAna/Prediction/PredictionSystJoint2018.h" 15 #include "CAFAna/Systs/NumuSysts.h" 20 #include "CAFAna/nue/Ana2019/joint_fit_2019_loader_tools.h" 44 std::vector <Spectrum * >
GetMockData(TString ana_options);
46 void POTtag(
bool both,
bool FHCOnly,
bool RHCOnly);
48 void MakeMaps ( std::vector <const IFitVar*> profvars,
49 std::map<const IFitVar*, TGraph*> &profMap,
50 std::vector <const ISyst* >
systs,
51 std::map<const ISyst*, TGraph*> &systMap);
55 std::vector <const IFitVar*> profvars,
56 std::vector <TString > profnames,
57 std::map<const IFitVar*, TGraph*> &profMap,
58 std::vector <const ISyst* >
systs,
59 std::map<const ISyst*, TGraph*> &systMap);
63 bool corrSysts =
true,
64 TString
options=
"joint_realData_both",
65 bool th23slice =
false,
66 bool octantSlice =
true,
67 bool dmsqSlice =
false,
68 bool th13Slice =
false)
73 bool nueOnly =
options.Contains(
"nueOnly");
74 bool numuOnly =
options.Contains(
"numuOnly");
75 bool joint =
options.Contains(
"joint");
76 assert (nueOnly || numuOnly || joint);
78 bool FHCOnly =
options.Contains(
"FHCOnly");
79 bool RHCOnly =
options.Contains(
"RHCOnly");
80 bool both =
options.Contains(
"both");
81 assert (FHCOnly || RHCOnly || both);
83 bool fake2019 =
options.Contains(
"fake2019");
84 bool realData =
options.Contains(
"realData");
87 if(corrSysts) suffix+=
"_systs";
88 else suffix+=
"_stats";
95 if(th23slice) suffix+=
"_th23";
96 if(dmsqSlice) suffix+=
"_dmsq";
97 if(th13Slice) suffix+=
"_th13";
98 if(!octantSlice)suffix+=
"_noOct";
101 TString
outfilename (outdir +
"hist_slices_2019_" + suffix);
102 TString outFCfilename (outdir +
"slices_FCInput_2019_" + suffix);
114 std::pair <TH1D*, double>
cos;
118 std::vector <predictions> preds;
119 std::vector <Spectrum*>
data;
120 std::vector <const IExperiment * > expts;
123 if(fake2019)
SetFakeCalc(calc_fake, 0.565, 2.48
e-3, 1.99);
126 if(FHCOnly || both ) {
129 if(RHCOnly || both ) {
136 if(FHCOnly || both ) {
139 for (
int i = 0;
i< nnumu;
i++ ){
143 if(RHCOnly || both ) {
146 for (
int i = 0;
i< nnumu;
i++ ){
152 cout<<
"\n\n Predictions used for the fit:\n";
156 if(FHCOnly || both ) {
158 data.insert(data.end(),
temp.begin(),
temp.end());
160 if(RHCOnly || both ) {
162 data.insert(data.end(),
temp.begin(),
temp.end());
166 if(FHCOnly || both ){
168 data.insert(data.end(),
temp.begin(),
temp.end());}
169 if(RHCOnly || both ){
171 data.insert(data.end(),
temp.begin(),
temp.end());}
176 for(
int i = 0;
i <
int(preds.size()); ++
i){
177 double POT = preds[
i].pot;
179 if(realData) thisdata = data[
i];
180 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;
191 else {
std::cout <<
"No reactor constraint, that's th13 slice\n";}
193 std::cout <<
"Adding Dmsq32ConstraintPDG2017\n";
196 std::cout <<
"Creating Multiexperiment with a total of " 197 << expts.size() <<
" experiments\n\n" <<
std::endl;
204 std::cout <<
"Systematics for the fit:\n";
207 std::vector< const ISyst * >
systs = {};
223 std::vector <SystShifts> seedShifts = {};
246 std::vector<double> seeds;
251 std::vector <th23helper> th23seeds;
259 std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
260 std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};
261 double seed_dmsq = 2.5e-3;
268 double minchi23 = 1E20;
269 for(
int hie:{-1, 1}){
270 for (
auto & thseed:th23seeds){
272 std::cout <<
"\n\nFinding best fit " << (hie > 0 ?
"NH " :
"IH ")
273 << thseed.label <<
", " 283 Fitter fit23(exptThis, fitvars, systs);
284 auto thisminchi = fit23.Fit(calc, auxShifts,
287 { thseed.var, thseed.seeds },
291 minchi23=
min(minchi23, thisminchi);
297 std::cout <<
"\nFound overall minchi " << minchi23 <<
"\n\n";
304 TFile *
outfile =
new TFile(outfilename+
".root",
"recreate");
305 TFile * outFCfile =
new TFile(outFCfilename+
".root",
"recreate");
310 v.Write(
"overall_minchi");
315 for(
int hie: {-1, +1}){
320 std::cout <<
"Starting profile " << (hie>0?
"NH ":
"IH") <<
"\n\n";
326 std::vector <const IFitVar * > profvars;
327 std::vector <TString> profvarnames;
328 std::map <const IFitVar *, std::vector <double> >profseeds;
333 for(
auto const & thseed:th23seeds){
335 const IFitVar* kCurTheta23 = thseed.var;
338 if(!th23slice && !dmsqSlice && !th13Slice){
343 profdef.profvarnames = {
"DmSq32",
"SinSq2Theta13",
"SinSqTheta23"};
344 profdef.profseeds = {{kCurTheta23, thseed.seeds}, { &
kFitDmSq32, {hie*seed_dmsq}} };
345 profdef.shortname=
"delta";
351 profdef.minX = (nueOnly ? 0:0.3);
352 profdef.maxX = (nueOnly ? 1:0.7);
354 profdef.profvarnames = {
"DmSq32",
"SinSq2Theta13",
"DeltaCPpi"};
356 profdef.shortname=
"th23";
364 profdef.profvarnames = {
"DmSq32",
"SinSqTheta23",
"DeltaCPpi"};
366 profdef.shortname=
"th13";
371 profdef.minX = hie > 0 ? 2
e-3 : -3
e-3;
372 profdef.maxX = hie > 0 ? 3
e-3 : -2
e-3;
374 profdef.profvarnames = {
"SinSqTheta23",
"SinSq2Theta13",
"DeltaCPpi"};
375 profdef.profseeds = { {kCurTheta23, thseed.seeds},
377 profdef.shortname=
"dmsq";
380 std::map<const IFitVar*, TGraph*> profVarsMap;
381 std::map<const ISyst*, TGraph*> systMap;
382 MakeMaps (profdef.profvars, profVarsMap, systs, systMap);
385 std::cout <<
" Profile "<< thseed.label <<
"\n";
388 profdef.fitvar, steps, profdef.minX, profdef.maxX,
394 profVarsMap, systMap);
396 TString hieroctstr = (hie > 0 ?
"NH" :
"IH") + thseed.label;
399 slice->Write((TString )
"slice_" + profdef.shortname +
"_" + hieroctstr);
401 profdef.profvars, profdef.profvarnames, profVarsMap, systs, systMap);
420 TFile *
infile =
new TFile (outfilename+
".root",
"read");
423 auto mins =* (
TVectorD*)infile->Get(
"overall_minchi");
424 double minchi23 = mins[0];
426 std::vector <TString> sliceNames = {};
427 if (!th23slice && ! dmsqSlice && !th13Slice) sliceNames = {
"slice_delta_"};
428 if (th23slice) sliceNames = {
"slice_th23_"};
429 if (dmsqSlice) sliceNames = {
"slice_dmsq_"};
430 if (th13Slice) sliceNames = {
"slice_th13_"};
435 std::vector<double> seeds;
439 std::vector <th23helper> th23seeds;
449 for(TString sliceName:sliceNames){
450 std::vector < std::pair <TGraph*, TString> >
graphs;
452 if(sliceName.Contains(
"delta")) {
453 for (TString hie:{
"NH",
"IH"}){
454 for(
auto const & thseed:th23seeds){
455 auto h = (TGraph*) infile ->
Get( sliceName + hie + thseed.label);
457 std::cerr <<
"Problem " << sliceName + hie + thseed.label <<
"\n";
460 graphs.push_back({
h, hie +
" " + thseed.label});
466 for (TString hie:{
"NH",
"IH"}) {
467 for(
auto const & thseed:th23seeds){
468 auto h = (TGraph*) infile ->
Get( sliceName + hie + thseed.label);
469 if(!
h) {
std::cerr <<
"Problem " << sliceName + hie + thseed.label <<
"\n"; }
470 graphs.push_back({
h,hie +
" " + thseed.label});
474 if(sliceName.Contains(
"th23")) {
476 (nueOnly ? 0 : 0.32), (nueOnly ? 1 : 0.72));
479 if(sliceName.Contains(
"dmsq")) {
482 if(sliceName.Contains(
"th13")) {
488 for (
auto &gr:graphs) {
489 if(gr.second.Contains(
"IH")) gr.first->SetLineColor(
kPrimColorIH);
491 if(gr.second.Contains(
"LO")) gr.first->SetLineStyle(7);
492 gr.first->SetLineWidth(3);
494 for (
int i =0;
i < gr.first->GetN(); ++
i){
495 gr.first->SetPoint(
i,1000*
fabs(gr.first->GetX()[
i]),gr.first->GetY()[
i]);
498 gr.first->Draw(
"l same");
505 auto leg =
SliceLegend(graphs, !dmsqSlice && !th23slice && !th13Slice, 0.7, dmsqSlice);
510 for(TString
ext: {
".pdf",
".root"}) {
511 gPad->Print(outdir +
"slice_" + suffix +
522 TFile * inFCfile =
new TFile(outFCfilename+
".root",
"read");
525 std::vector <TString> strprof = {
"DmSq32",
"SinSq2Theta13",
"SinSqTheta23"};
526 if(th23slice) strprof = {
"DmSq32",
"SinSq2Theta13",
"DeltaCPpi"};
527 if(th13Slice) strprof = {
"DmSq32",
"DeltaCPpi",
"SinSqTheta23"};
528 if(dmsqSlice) strprof = {
"SinSqTheta23",
"SinSq2Theta13",
"DeltaCPpi"};
530 std::vector< const ISyst * >
systs = {};
535 for (
const ISyst*
s : systs){
536 strprof.push_back(
s->ShortName());
539 for(
auto const &
str:strprof){
541 double miny=-1.5,
maxy=1.5;
542 if(
str==
"DeltaCPpi") {miny=0;
maxy=2;}
543 if(
str==
"DmSq32") {miny=2.2;
maxy=2.9;}
544 if(
str==
"SinSq2Theta13") {miny=0.080;
maxy=0.086;}
545 if(
str==
"SinSqTheta23") {miny=0.3;
maxy=0.8;}
547 auto axes =
new TH2F(
"",
";#delta_{CP}/#pi;" +
str, 100, 0, 2, 100, miny,
maxy);
548 if(dmsqSlice)
axes =
new TH2F(
"",
";|#Deltam^{2}_{32}| (10^{-3} eV^{2});" +
str, 100, 2.2, 2.8, 100, miny,
maxy);
549 if(th23slice)
axes =
new TH2F(
"",
";sin^{2}#theta_{23};" +
str, 100, 0.3, 0.7, 100, miny,
maxy);
550 if(th13Slice)
axes =
new TH2F(
"",
";sin^{2}#theta_{13};" +
str, 100, 0.01, 0.2, 100, miny,
maxy);
553 for (TString hie:{
"NH",
"IH"}){
554 for(
auto const & thseed:th23seeds){
555 auto gr = (TGraph*) inFCfile->Get(hie + thseed.label +
"_" +
str);
557 std::cerr <<
"Problem " << hie + thseed.label +
"_" +
str <<
"\n";
561 for (
int i = 0;
i < (
int) gr->GetN(); ++
i){
562 gr->SetPoint(
i, gr->GetX()[
i], 1000*
fabs(gr->GetY()[
i]));
566 for (
int i = 0;
i < (
int) gr->GetN(); ++
i){
567 gr->SetPoint(
i, 1000*
fabs(gr->GetX()[
i]), gr->GetY()[
i]);
571 if(thseed.label.Contains(
"LO")) gr->SetLineStyle(7);
573 gr->SetMarkerColor(gr->GetLineColor());
578 gPad->Print(outdir +
"debug_slice_prof_" + suffix +
"_"+
str +
".pdf");
586 void MakeMaps ( std::vector <const IFitVar*> profvars,
587 std::map<const IFitVar*, TGraph*> &profMap,
588 std::vector <const ISyst* >
systs,
589 std::map<const ISyst*, TGraph*> &systMap)
592 profMap.insert(std::pair<const IFitVar*, TGraph*> (
var,
new TGraph()));
593 for (
const ISyst* syst : systs)
594 systMap.insert(std::pair<const ISyst*, TGraph*> (syst,
new TGraph()));
598 std::vector <const IFitVar*> profvars,
599 std::vector <TString > profnames,
600 std::map<const IFitVar*, TGraph*> &profMap,
601 std::vector <const ISyst* >
systs,
602 std::map<const ISyst*, TGraph*> &systMap)
604 TDirectory *
tmp = gDirectory;
606 for (
int i = 0;
i < (
int) profvars.size(); ++
i){
607 profMap[profvars[
i]]->Write((prefix +
"_" + profnames[
i]));
609 for (
const ISyst*
s : systs)
610 systMap[
s]->Write((prefix +
"_" +
s->ShortName()));
617 cout<<
"\n\n ============= You're loading Mock data from the file ============="<<
endl;
618 std::vector <Spectrum *> spec;
620 bool RHCmode = ana_options.Contains(
"rhc");
621 bool FHCmode = ana_options.Contains(
"fhc");
622 bool NueData = ana_options.Contains(
"nue");
623 bool NumuData = ana_options.Contains(
"numu");
625 if (FHCmode) spec.push_back(LoadFromFile<Spectrum>(Dir+
"mock_data_for_yn_tutorial.root",
"data_exp_0").
release());
626 if (RHCmode) spec.push_back(LoadFromFile<Spectrum>(Dir+
"mock_data_for_yn_tutorial.root",
"data_exp_1").
release());
630 for (
int i = 2;
i < 6; ++
i) {
632 spec.push_back(
temp);
636 for (
int i = 6;
i < 10; ++
i) {
638 spec.push_back(
temp);
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.
fvar< T > fabs(const fvar< T > &x)
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.
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 XML_Char const XML_Char * data
string outfilename
knobs that need extra care
const DummyAnaSyst kAnaCalibrationSyst("Calibration","AbsCalib")
const Color_t k3SigmaConfidenceColorNH
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
virtual T GetDmsq32() const
TLegend * SliceLegend(std::vector< std::pair< TGraph *, TString > > &graphs, bool isDelta)
void Tutorial2019FitSlices(bool createFile=false, bool corrSysts=true, TString options="joint_realData_both", bool th23slice=false, bool octantSlice=true, bool dmsqSlice=false, bool th13Slice=false)
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)
Interface definition for fittable variables.
double GetPOT(bool isFHC=true)
void MakeMaps(std::vector< const IFitVar * > profvars, std::map< const IFitVar *, TGraph * > &profMap, std::vector< const ISyst * > systs, std::map< const ISyst *, TGraph * > &systMap)
Standard interface to all prediction techniques.
const FitSinSq2Theta13 kFitSinSq2Theta13
std::string to_string(ModuleType mt)
std::pair< Spectrum *, double > cos
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::vector< Spectrum * > GetMockData(TString ana_options)
void POTtag(bool both, bool FHCOnly, bool RHCOnly)
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.