25 return norm *
exp(-1*
pow(x-mean,2)/(2*var));
35 std::vector<std::string>
one_d = {
"AbsCalibShape",
"Cherenkov",
"XSecTune" };
46 gSystem->MakeDirectory(
"tex_pulls");
47 gSystem->MakeDirectory(
"tex_pulls/plots");
48 gSystem->CopyFile(Form(
"%s/nus/Nus20/Plots/header_nus.tex",
52 std::ofstream maintex;
53 maintex.open(
"tex_pulls/main.tex");
54 maintex <<
"\n\\input{header_nus.tex}\n\n";
55 maintex <<
"\\newpage\n\n";
71 vND(
i) = hND->GetBinContent(
i+1);
73 vFD(
i) = hFD->GetBinContent(
i+1);
74 HistCache::Delete(hND);
75 HistCache::Delete(hFD);
78 std::map<std::string, CovarianceMatrix*> mxND;
80 mxND[itMx->GetName()] = itMx;
81 std::map<std::string, CovarianceMatrix*> mxFD;
83 mxFD[itMx->GetName()] = itMx;
86 maintex <<
"\\section{Near Detector}\n\n";
95 std::replace(systname.begin(), systname.end(),
'_',
' ');
96 maintex <<
"\\subsection{" << systname <<
"}\n\n";
103 maintex <<
"\\section{Far Detector}\n\n";
112 std::replace(systname.begin(), systname.end(),
'_',
' ');
113 maintex <<
"\\subsection{" << systname <<
"}\n\n";
119 maintex <<
"\\end{document}\n";
128 f <<
"\\begin{figure}[h!]\n \\centering\n";
136 std::vector<TH1D*>
bins(nbins);
137 bool allNoEffect =
true;
138 bool allOneDirectional =
true;
139 std::vector<bool> oneDirectional(nbins,
false);
141 bins[
i] =
new TH1D(
"",
"", 60, 1-(3*
sqrt(mx(
i,
i))), 1+(3*
sqrt(mx(
i,
i))));
147 for (
size_t i = 0;
i <= 1e5; ++
i) {
151 else shifts.
SetShift(syst, rnd.Gaus());
155 HistCache::Delete(h);
164 for (
size_t i = 0;
i <
nbins; ++
i) {
168 oneDirectional[
i] =
true;
169 else allOneDirectional =
false;
175 std::vector<float> nonGaussianScore(nbins, 0);
176 for (
size_t i = 0;
i <
nbins; ++
i) {
177 float xrange = 1 - bins[
i]->GetXaxis()->GetBinLowEdge(1);
178 float norm = 0.1*bins[
i]->Integral()*
sqrt(mx(
i,
i));
179 for (
int j = -999;
j <= +999; ++
j) {
180 float x = (0.001 *
float(
j) * xrange)+1;
181 float yHist = bins[
i]->GetBinContent(bins[
i]->
GetXaxis()->FindBin(x));
182 float yGaus = norm *
Gaus(1, mx(
i,
i), x);
183 if (oneDirectional[
i]) yGaus *= 2;
184 float yMin = yHist < yGaus ? yHist : yGaus;
185 nonGaussianScore[
i] += (yMin * 0.1 * xrange) / norm;
190 TCanvas
c(
"",
"", 900, 1600);
191 std::vector<TPad*> pads;
192 std::vector<TGraph*>
graphs;
193 pads.push_back(
new TPad(
"classification",
"", 0, 0.95, 1, 0.999));
196 if (allNoEffect) classification =
"No effect";
197 else if (allOneDirectional) classification =
"One-directional";
198 else classification =
"Two-directional";
199 float meanScore = 0.;
200 for (
size_t i = 0;
i <
nbins; ++
i) meanScore += nonGaussianScore[
i] /
float(nbins);
201 TText*
t =
new TText(0.5, 0.5, Form(
"%s, score = %.1f%%", classification.c_str(), meanScore));
205 for (
size_t i = 0;
i <
nbins; ++
i) {
208 pads.push_back(
new TPad(Form(
"bin%zu",
i),
"", x1, y1, x1+0.25, y1+0.18));
213 float norm = 0.1*bins[
i]->Integral()*
sqrt(mx(
i,
i));
215 x[
j] = 1 + (-3.0*
sqrt(mx(
i,
i))) + (
float(
j)*binWidth);
217 if (oneDirectional[
i]) y[
j] *= 2;
219 graphs.push_back(
new TGraph(nPoints, &x[0], &y[0]));
220 bins[
i]->Draw(
"hist");
221 bins[
i]->SetTitle(Form(
"Bin %zu - %.2f < E (GeV) < %.2f",
i+1,
222 hNom->GetXaxis()->GetBinLowEdge(
i+1), hNom->GetXaxis()->GetBinUpEdge(
i+1)));
223 graphs.back()->Draw(
"c same");
224 graphs.back()->SetLineColor(
kRed-7);
225 graphs.back()->SetLineWidth(3);
226 TText*
tmp =
new TText(0.85, 0.85, Form(
"%.1f%%", nonGaussianScore[
i]));
227 tmp->SetTextAlign(33);
228 tmp->SetTextSize(0.1);
235 for (
auto pad : pads)
pad->Draw();
237 c.SaveAs(Form(
"tex_pulls/plots/gaus_%s_%s.pdf",
238 det.c_str(), systname.c_str()));
239 f <<
" \\includegraphics[height=0.8\\textheight]{" 240 << Form(
"plots/gaus_%s_%s.pdf", det.c_str(), systname.c_str())
241 <<
"}\n\\end{figure}\n\n\\clearpage\n\n";
243 HistCache::Delete(hNom);
250 f <<
"\\begin{figure}[h!]\n \\centering\n";
259 std::vector<std::vector<float>>
y(nbins, std::vector<float>(
nPoints));
265 x[
i] = (6./
float(nPoints-1))*(
float(
i)-(0.5*(nPoints-1)));
269 y[
j][i] = h->GetBinContent(
j+1);
270 HistCache::Delete(h);
274 TCanvas
c(
"",
"", 900, 1600);
275 std::vector<TPad*> pads;
276 std::vector<TGraph*>
graphs;
277 for (
size_t i = 0;
i <
nbins; ++
i) {
280 pads.push_back(
new TPad(Form(
"bin%zu",
i),
"", x1, y1, x1+0.25, y1+0.1999));
282 graphs.push_back(
new TGraph(nPoints, &x[0], &y[
i][0]));
283 graphs.back()->Draw(
"al");
284 graphs.back()->SetTitle(Form(
"Bin %zu - %.2f < E (GeV) < %.2f", i+1,
285 hNom->GetXaxis()->GetBinLowEdge(i+1), hNom->GetXaxis()->GetBinUpEdge(i+1)));
286 graphs.back()->SetLineWidth(3);
290 for (
auto pad : pads)
pad->Draw();
292 c.SaveAs(Form(
"tex_pulls/plots/pull_%s_%s.pdf",
293 det.c_str(), systname.c_str()));
294 f <<
" \\includegraphics[height=0.8\\textheight]{" 295 << Form(
"plots/pull_%s_%s.pdf", det.c_str(), systname.c_str())
296 <<
"}\n\\end{figure}\n\n\\clearpage\n\n";
Cuts and Vars for the 2020 FD DiF Study.
std::string GetSystBaseName(std::string name)
Get the last part of a systematic's ShortName.
std::vector< const ISyst * > getAllXsecSysts_2018()
Get master XSec syst list for 2018 analyses.
fvar< T > fabs(const fvar< T > &x)
Float_t y1[n_points_granero]
void SetNus20Params(osc::OscCalcSterile *calc, std::string type="3flav")
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
virtual const std::string & ShortName() const final
The name printed out to the screen.
double Integral(const Spectrum &s, const double pot, int cvnregion)
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
Simple record of shifts applied to systematic parameters.
Float_t x1[n_points_granero]
const double kAna2018SensitivityFHCNDPOT
std::string FindCAFAnaDir()
Encapsulate code to systematically shift a caf::SRProxy.
void DrawBins(IPrediction *pred, const ISyst *syst, TMatrixD mx, TVectorD vec, size_t nbins, std::ofstream &f, std::string det)
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
std::vector< std::string > one_d
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
correl_xv GetXaxis() -> SetDecimals()
fvar< T > exp(const fvar< T > &x)
virtual Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &syst) const
const covmx::Sample kNusFHCFarDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kFarDet)
std::vector< CovarianceMatrix * > GetMatrices(std::vector< covmx::Sample > samples, bool cvmfs=true)
float Gaus(float mean, float var, float x)
const Binning kNCNDBinning
NC custom binning.
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
fvar< T > floor(const fvar< T > &x)
const Binning kNCFDBinning
Standard interface to all prediction techniques.
void DrawPulls(IPrediction *pred, const ISyst *syst, size_t nbins, std::ofstream &f, std::string det)
const double kAna2018FHCPOT
Prevent histograms being added to the current directory.
const covmx::Sample kNusFHCNearDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kNearDet)
void SetShift(const ISyst *syst, double shift, bool force=false)