21 const int N = plot->GetN();
22 low = plot->GetX()[0];
23 hi = plot->GetX()[N-1];
24 float min = plot->GetY()[0];
25 float bf = plot->GetX()[0];
28 for (
int i = 0;
i < N; ++
i )
30 float y = plot->GetY()[
i];
31 if ( y < min ) { min =
y; bf = plot->GetX()[
i]; bfp =
i; }
35 const float chi2th = 1.0;
37 for (
int i = bfp;
i > 0; --
i )
39 float y1 = plot->GetY()[
i-1];
40 float y2 = plot->GetY()[
i];
41 if ( y1 > chi2th && y2 < chi2th )
43 float x1 = plot->GetX()[
i-1];
44 float x2 = plot->GetX()[
i];
45 low = x1 + ( chi2th -
y1 ) * ( x2 - x1 ) / ( y2 -
y1 );
46 std::cout <<
"Crossed chisq of " << chi2th <<
", descending, with x value of: " << low <<
std::endl;
51 for (
int i = bfp;
i < N-1; ++
i )
53 float y1 = plot->GetY()[
i];
54 float y2 = plot->GetY()[
i+1];
55 if ( y1 < chi2th && y2 > chi2th )
57 float x1 = plot->GetX()[
i];
58 float x2 = plot->GetX()[
i+1];
59 hi = x1 + ( chi2th -
y1 ) * ( x2 - x1 ) / ( y2 -
y1 );
60 std::cout <<
"Crossed chisq of " << chi2th <<
", ascending, with x value of: " << hi <<
std::endl;
69 TLatex*
prelim =
new TLatex(.9, .93,
"NOvA Simulation");
70 prelim->SetTextColor(kGray+1);
72 prelim->SetTextSize(1.5/30.);
73 prelim->SetTextAlign(32);
79 txt<<
" & +" << plus <<
" / -" <<
minus;
83 typedef std::vector<std::pair<const ISyst*,const ISyst*> >
Correlations;
84 typedef std::pair<SingleSampleExperiment*,Correlations>
ExptCorr;
88 bool isFHC = beam ==
"fhc";
89 bool corrSysts =
true;
93 if ( !ptExtrap ) decomp +=
"_noPt";
99 exptCorrs.emplace_back( expt,
GetCorrelations( isNue, isFHC, ptExtrap ) );
104 bool isFHC = beam ==
"fhc";
105 bool useSysts =
true;
115 for (
int i = 0;
i < nEhadQuants; ++
i )
119 exptCorrs.emplace_back( expt, corrs );
126 const bool full = option.Contains(
"full" );
127 const bool ptExtrap = option.Contains(
"ptExtrap" );
128 const bool fakeNDData = option.Contains(
"fakeND" );
129 const bool seeds = option.Contains(
"seeds" );
132 if ( ptExtrap ) outStr +=
"ptExtrap_";
133 if ( fakeNDData ) outStr +=
"fakeND_";
134 if ( seeds ) outStr +=
"withSeeds_";
135 outStr += full ?
"full" :
"short";
140 if ( gSystem->AccessPathName( outDir.c_str() ) ) gSystem->mkdir( outDir.c_str() );
158 std::string calcfile =
"/nova/ana/3flavor/Ana2020/Fits/RealData/WithSeeds/bestfits_joint_realData_both_systs.root";
160 TFile *calcf =
new TFile(calcfile.c_str());
161 auto calctemp = LoadFrom<osc::IOscCalc>(calcf,
"NH_UO_calc").
release();
175 std::vector<ExptCorr> exptCorrs;
182 for (
unsigned int i = 0;
i < exptCorrs.size(); ++
i )
185 exptThis->Add( exptCorrs[
i].first );
186 exptThis->SetSystCorrelations(
i, exptCorrs[
i].
second );
191 bool smallgenie =
true;
196 std::vector<const ISyst*> xsecSysts;
198 std::vector<const ISyst*> fluxSysts;
200 std::vector<const ISyst*> calibSysts;
202 std::vector<const ISyst*> lightSysts;
204 std::vector<const ISyst*> leptonRecoSysts;
207 std::vector<const ISyst*> neutronSysts;
209 std::vector<const ISyst*> nearFarSysts;
215 std::map< std::string, std::vector<const ISyst*> > mapSystVecs;
216 mapSystVecs[
"Neutrino Cross Sections" ] = xsecSysts;
217 mapSystVecs[
"Beam Flux" ] = fluxSysts;
218 mapSystVecs[
"Detector Calibration" ] = calibSysts;
219 mapSystVecs[
"Detector Response" ] = lightSysts;
220 mapSystVecs[
"Lepton Reconstruction" ] = leptonRecoSysts;
221 mapSystVecs[
"Neutron Uncertainty" ] = neutronSysts;
222 mapSystVecs[
"Near-Far Uncor." ] = nearFarSysts;
235 const int steps = 60;
240 std::vector<const IFitVar *> profvars;
241 std::map <const IFitVar *, std::vector <double> >profseeds;
250 std::vector<ProfDef> profiles = {
252 ssth23, 0,
"$\\mathrm{sin}^2\\theta_{23}$",
"Uncertainty in sin^{2}#theta_{23}",
"th23"},
254 dCP, 0,
"$\\delta_{CP}/\\pi$",
"Uncertainty in #delta_{CP}/#pi",
"dcp"},
256 dmsq32, 3,
"$\\vert\\Delta m^2_{32}\\vert$ ($\\times 10^{-3}~\\mathrm{eV}^2$)",
"Uncertainty in #Deltam^{2}_{32} (#times10^{-3} eV^{2})",
"dmsq"}
260 TFile*
fOut =
new TFile( ( outDir +
"/syst_table_profiles_" + outStr +
".root" ).c_str(),
"RECREATE" );
263 txt.open( ( outDir +
"/syst_table_" + outStr +
".txt" ).c_str() );
265 txt<<std::setprecision(2);
273 std::map<std::string, double> statup;
274 std::map<std::string, double> statdn;
275 std::map<std::string,std::map<std::string,std::pair<double,double>>> bars;
277 txt<<
"Source of Uncertainty";
278 for(
auto prof: profiles) txt<<
" & " << prof.latex;
283 for(
auto prof: profiles){
289 steps, prof.minx, prof.maxx, -1,
294 plot ->
Write((prof.suffix+
"_stat").c_str());
298 double up = hi-prof.bf;
299 double dn = prof.bf-low;
301 statup[prof.suffix] = up;
302 statdn[prof.suffix] = dn;
307 for(
auto const &systVec: mapSystVecs) {
308 if ( systVec.second.empty() )
313 txt<< systVec.first ;
315 for(
auto prof: profiles){
319 TGraph* systplot = (TGraph*)
Profile(exptThis,
calc,
321 steps, prof.minx, prof.maxx, -1,
325 systplot ->
Write((prof.suffix+
"_" + systVec.first).c_str());
329 double plus =
sqrt(
fabs( (
hi-prof.bf)*(
hi-prof.bf) - (statup[prof.suffix]*statup[prof.suffix]) ) ) *
pow(10,prof.oom);
330 double minus =
sqrt(
fabs( (prof.bf-low)*(prof.bf-low) - (statdn[prof.suffix]*statdn[prof.suffix]) ) ) *
pow(10,prof.oom);
334 bars[prof.suffix][systVec.first] = {
minus,plus};
341 for(
const ISyst* syst: allSysts) {
342 txt<< syst->LatexName() ;
343 for(
auto prof: profiles){
347 TGraph* systplot = (TGraph*)
Profile(exptThis,
calc,
349 steps, prof.minx, prof.maxx, -1,
354 systplot ->
Write((prof.suffix+
"_" + syst->ShortName()).c_str());
358 double plus =
sqrt(
fabs( (
hi-prof.bf)*(
hi-prof.bf) - (statup[prof.suffix]*statup[prof.suffix]) ) ) *
pow(10,prof.oom);
359 double minus =
sqrt(
fabs( (prof.bf-low)*(prof.bf-low) - (statdn[prof.suffix]*statdn[prof.suffix]) ) ) *
pow(10,prof.oom);
363 bars[prof.suffix][syst->ShortName()] = {
minus, plus};
372 txt<<
"Systematic Uncertainty" ;
373 for(
auto prof: profiles){
378 steps, prof.minx, prof.maxx, -1,
382 totplot ->
Write((prof.suffix+
"_totError").c_str());
386 double plus =
sqrt(
fabs( (
hi-prof.bf)*(
hi-prof.bf) - (statup[prof.suffix]*statup[prof.suffix]) ) ) *
pow(10,prof.oom);
387 double minus =
sqrt(
fabs( (prof.bf-low)*(prof.bf-low) - (statdn[prof.suffix]*statdn[prof.suffix]) ) ) *
pow(10,prof.oom);
394 txt<<
"Statistical Uncertainty" ;
395 for(
auto prof: profiles){
397 double plus = statup[prof.suffix] *
pow(10,prof.oom);
398 double minus = statdn[prof.suffix] *
pow(10,prof.oom);
406 for(
auto prof: profiles){
408 ErrorBarChart(bars[prof.suffix],{statdn[prof.suffix]*pow(10,prof.oom),statup[prof.suffix]*pow(10,prof.oom)},prof.root);
409 gPad->SetLeftMargin(0.35);
410 gPad->SetBottomMargin(0.15);
412 gPad->Print( (
outDir +
"/syst_table_" + prof.suffix +
"_" + outStr +
".pdf" ).c_str() );
413 gPad->Print( (
outDir +
"/syst_table_" + prof.suffix +
"_" + outStr +
".png" ).c_str() );
414 gPad->Print( (
outDir +
"/syst_table_" + prof.suffix +
"_" + outStr +
".C" ).c_str() );
Cuts and Vars for the 2020 FD DiF Study.
void AddExptCorrNumu(std::vector< ExptCorr > &exptCorrs, osc::IOscCalc *calc, std::string beam, bool ptExtrap, bool fakeNDData)
fvar< T > fabs(const fvar< T > &x)
Float_t y1[n_points_granero]
double GetValue(const osc::IOscCalcAdjustable *osc) const override
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Forward to wrapped Var's GetValue()
const FitDmSq32 kFitDmSq32
Float_t x1[n_points_granero]
double GetValue(const osc::IOscCalcAdjustable *osc) const override
General interface to oscillation calculators.
float GetLimits(TH1F *plot, float &hi, float &low)
TGraph * Profile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double input_minchi, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &profSystsMap, MinuitFitter::FitOpts opts)
scan in one variable, profiling over all others
void AddNonLoadable2020Systs(std::vector< const ISyst * > &systs, const EAnaType2020 ana)
void Add3FlavorAna2020NeutronSysts(std::vector< const ISyst * > &systs)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::vector< const IPrediction * > GetNumuPredictions2020(const int nq=4, std::string decomp="noPt", osc::IOscCalc *calc=DefaultOscCalc(), bool useSysts=true, std::string beam="fhc", bool isFakeData=false, bool GetFromUPS=false, bool minimizeMemory=true, bool NERSC=false)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Encapsulate code to systematically shift a caf::SRProxy.
Representation of a spectrum in any variable, with associated POT.
std::pair< Spectrum *, double > GetNueCosmics2020(std::string beam, bool GetFromUPS=false, bool NERSC=false)
void Add3FlavorAna2020MichelTagSysts(std::vector< const ISyst * > &systs, const EAnaType2020 ana, const BeamType2020 beam)
TSpline3 hi("hi", xhi, yhi, 18,"0")
std::pair< SingleSampleExperiment *, Correlations > ExptCorr
void syst_table_fit_new(const TString option)
std::vector< std::pair< const ISyst *, const ISyst * > > Correlations
root
Link up the nodes tree #####.
const IPrediction * GetNuePrediction2020(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=true, bool NERSC=false)
void Add3FlavorAna2020LeptonAngleSysts(std::vector< const ISyst * > &systs, const bool ptExtrap)
const double kAna2020FHCPOT
const double kAna2020RHCPOT
void AddExptCorrNue(std::vector< ExptCorr > &exptCorrs, osc::IOscCalc *calc, std::string beam, bool ptExtrap, bool fakeNDData)
std::vector< double > POT
std::vector< std::pair< Spectrum *, double > > GetNumuCosmics2020(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
static float min(const float a, const float b, const float c)
Combine multiple component experiments.
void Add3FlavorAna2020LightSysts(std::vector< const ISyst * > &systs)
void Add3FlavorAna2020NueAcceptSysts(std::vector< const ISyst * > &systs, const EAnaType2020 ana, const BeamType2020 beam, const bool ptExtrap)
void Add3FlavorAna2020NormSysts(std::vector< const ISyst * > &systs, const BeamType2020 beam)
void Add3FlavorAna2020MuEnergySysts(std::vector< const ISyst * > &systs, const EAnaType2020 ana)
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
void WriteUncertainty(ofstream &txt, double plus, double minus)
void Add3FlavorAna2020XSecSysts(std::vector< const ISyst * > &systs, const EAnaType2020 ana, bool smallgenie)
std::vector< const ISyst * > get3FlavorAna2020AllSysts(const EAnaType2020 ana, const bool smallgenie, const BeamType2020 beam, const bool isFit, const bool ptExtrap)
void ErrorBarChart(const std::map< std::string, std::pair< double, double >> &systErrors, const std::pair< double, double > &statErr, const std::string &label)
Make a simple plot of relative size of different errors.
Interface definition for fittable variables.
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
const ReactorExperiment * WorldReactorConstraint2019()
Reactor constraint from PDG 2019.
const FitSinSq2Theta13 kFitSinSq2Theta13
void Add3FlavorAna2020BeamSysts(std::vector< const ISyst * > &systs, const EAnaType2020 ana)
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
void Add3FlavorAna2020CalibSysts(std::vector< const ISyst * > &systs)
virtual IOscCalc * Copy() const override
Compare a single data spectrum to the MC + cosmics expectation.