39 std::string fileout = folder +
"/dataMCNDAnaOut.root";
41 TFile* spSaved =
new TFile(filein.c_str());
42 TFile* plotsOut =
new TFile(fileout.c_str(),
"RECREATE");
56 TIter
next(SavedFile->GetListOfKeys());
58 while ((key=(TKey*)
next()))
60 const char* fDatKey = key->GetName();
63 if(strDatKey.find(
"MC") != std::string::npos) {
continue; }
66 strMCKey.replace(strMCKey.size()-4, 4,
"MC");
69 fPlotLabel.erase(fPlotLabel.end()-4, fPlotLabel.end());
103 TH1* hDat = data .
ToTH1(pot);
104 TH1* hMCTot = mcspecs[0].
ToTH1(pot);
105 TH1* hMCNC = mcspecs[1].
ToTH1(pot);
106 TH1* hMCNuMu = mcspecs[2].
ToTH1(pot);
107 TH1* hMCNuE = mcspecs[3].
ToTH1(pot);
110 if(xLabel.find(
"PT/P") != std::string::npos) {
111 xLabel.replace(xLabel.find(
"PT/P"), 4,
"p_{T}/p");
113 if(xLabel.find(
"RecoE") != std::string::npos) {
114 xLabel =
"Visible Energy (GeV)";
116 if(xLabel.find(
"RecoELong") != std::string::npos) {
117 xLabel =
"Visible Energy (GeV)";
121 if(plotlabel.find(
"RecoE") != std::string::npos) { yLabel =
"10^{3} Events / 0.1 GeV / 8.04 #times10^{20} POT"; }
122 if(plotlabel.find(
"RecoELong") != std::string::npos) { yLabel =
"10^{3} Events / 0.1 GeV / 8.04 #times10^{20} POT"; }
123 else yLabel =
"10^{3} Events / 8.04 #times10^{20} POT";
124 std::string histTitle =
";" + xLabel +
";" + yLabel;
127 TH1* hNumuStack = (TH1*)hMCNuMu->Clone();
130 TH1* hNueStack = (TH1*)hNumuStack->Clone();
131 hNueStack->Add(hMCNuE);
134 TH1* hNCStack = (TH1*)hNueStack->Clone();
135 hNCStack->Add(hMCNC);
138 double maxDat = hDat ->GetBinContent(hDat ->GetMaximumBin());
139 double maxMCTot = hMCTot ->GetBinContent(hMCTot ->GetMaximumBin());
140 double maxMCNC = hMCNC ->GetBinContent(hMCNC ->GetMaximumBin());
141 double maxMCNuE = hMCNuE ->GetBinContent(hMCNuE ->GetMaximumBin());
142 double maxMCNuMu = hMCNuMu->GetBinContent(hMCNuMu->GetMaximumBin());
149 int nbins = hDat->GetNbinsX();
151 hDat->SetBinError(
i,
sqrt(hDat->GetBinContent(
i)));
155 double ndScale = 1./1000.;
156 hDat ->Scale(ndScale);
157 hMCTot ->Scale(ndScale);
158 hMCNC ->Scale(ndScale);
159 hMCNuE ->Scale(ndScale);
160 hMCNuMu->Scale(ndScale);
161 hNCStack ->Scale(ndScale);
162 hNueStack ->Scale(ndScale);
163 hNumuStack->Scale(ndScale);
172 hDat->SetMarkerStyle(kFullCircle);
179 double xL = 0.5, xR = 0.7, yB = 0.64, yT = 0.84;
180 TLegend*
leg =
new TLegend(xL, yB, xR, yT);
182 leg->AddEntry(hDat,
"ND Data",
"lep");
183 leg->AddEntry(hMCTot,
"Total Prediction",
"l");
184 leg->AddEntry(hMCNC,
"NC Prediction",
"l");
185 leg->AddEntry(hMCNuE,
"#nu_{e} CC Background",
"l");
186 leg->AddEntry(hMCNuMu,
"#nu_{#mu} CC Background",
"l");
187 leg->SetY1(leg->GetY2() - leg->GetNRows()*0.05);
189 TH1* hRat = (TH1*)hDat->Clone(
"hRat");
190 hRat->Divide(hMCTot);
191 TH1* hLine = (TH1*)hMCTot->Clone(
"");
193 double maxvalRat = 1., minvalRat = 1.;
195 if(hDat->GetBinContent(
i) > 0) {
196 double errpct = hDat->GetBinError(
i)/hDat->GetBinContent(
i);
197 hRat->SetBinError(
i, errpct);
200 double binval = hRat->GetBinContent(
i);
201 double errval = hRat->GetBinError(
i);
203 maxvalRat =
std::max(maxvalRat, binval + errval);
204 minvalRat =
std::min(minvalRat, binval - errval);
207 hLine->SetBinContent(
i, 1.);
210 hRat->GetYaxis()->SetTitle(
"Data / Total MC");
211 hRat->SetLineColor(kBlack);
212 hLine->GetYaxis()->SetTitle(
"Data / Total MC");
215 hRat ->SetMaximum(1.05*maxvalRat);
216 hLine->SetMaximum(1.05*maxvalRat);
217 hRat ->SetMinimum(0.95*minvalRat);
218 hLine->SetMinimum(0.95*minvalRat);
222 if(plotlabel.find(
"NM1") == std::string::npos &&
223 plotlabel.find(
"ECalLongAllE") == std::string::npos) {
224 TCanvas* cnorat =
new TCanvas(plotlabel.c_str(), plotlabel.c_str(), 800, 500);
225 gPad->SetFillStyle(0);
226 hDat->SetMinimum(0.001);
227 hMCTot->SetMinimum(0.001);
228 hMCNC->SetMinimum(0.001);
229 hMCNuE->SetMinimum(0.001);
230 hMCNuMu->SetMinimum(0.001);
233 hMCTot->Draw(
"hist same");
234 hMCNC->Draw(
"hist same");
235 hMCNuE->Draw(
"hist same");
236 hMCNuMu->Draw(
"hist same");
241 plotsOut->WriteTObject(cnorat);
244 TCanvas* cStack =
new TCanvas((plotlabel +
"Stack").c_str(),
245 (plotlabel +
"Stack").c_str(), 800, 500);
246 gPad->SetFillStyle(0);
248 TLegend* legStack =
new TLegend(xL, yB, xR, yT);
250 legStack->AddEntry(hDat,
"ND Data",
"lep");
251 legStack->AddEntry(hNCStack,
"NC 3 Flavor Prediction",
"f");
252 legStack->AddEntry(hNueStack,
"#nu_{e} CC Background",
"f");
253 legStack->AddEntry(hNumuStack,
"#nu_{#mu} CC Background",
"f");
254 legStack->SetY1(legStack->GetY2() - legStack->GetNRows()*0.05);
256 hNCStack ->Draw(
"hist");
257 hNueStack ->Draw(
"hist same");
258 hNumuStack->Draw(
"hist same");
263 if(plotlabel.find(
"RecoELong") != std::string::npos) {
265 textFO = fopen(
"out2.txt",
"a+");
267 fprintf(textFO,
"Event numbers at the %2.2s for %s %.2fe20 POT:\n",
"ND", plotlabel.c_str(), 8.04);
268 if(hDat) { fprintf(textFO,
"%12.12s, ",
"Data"); }
269 else { fprintf(textFO,
"%14.14s",
""); }
270 fprintf(textFO,
"%12.12s, %12.12s, %12.12s, %12.12s, %12.12s\n",
271 "All",
"NC (3 Flav)",
"Numu",
"Nue",
"FOM");
273 double nAll = hMCTot ->Integral();
274 double nNC3F = hMCNC ->Integral();
275 double nNumu = hMCNuMu->Integral();
276 double nNue = hMCNuE ->Integral();
277 double fom = nNC3F/
sqrt(nAll);
279 double nData = hDat->Integral(1, hDat->GetNbinsX());
280 fprintf(textFO,
"%10E, ", nData);
283 fprintf(textFO,
"%14.14s",
"");
285 fprintf(textFO,
"%10E, %10E, %10E, %10E, %10E\n\n",
286 nAll, nNC3F, nNumu, nNue, fom);
290 plotsOut->WriteTObject(cStack);
294 textF = fopen(
"out.txt",
"a+");
298 PlotStack(mcspecs, plotsOut, textF, plotlabel,
"ND",
303 TCanvas*
canvas =
new TCanvas((plotlabel +
"Rat").c_str(),
304 (plotlabel +
"Rat").c_str(), 800, 800);
305 gPad->SetFillStyle(0);
306 TPad* pSpecs =
new TPad(
"pSpecs",
"", 0., 0.375, 1., 1.);
307 TPad* pRatio =
new TPad(
"pRatio",
"", 0., 0., 1., 0.375);
313 hDat->GetXaxis()->SetRangeUser(0,5);
314 hDat->GetXaxis()->SetNdivisions(-405);
315 hMCTot->Draw(
"hist same");
316 hMCNC->Draw(
"hist same");
317 hMCNuE->Draw(
"hist same");
318 hMCNuMu->Draw(
"hist same");
324 hLine->Draw(
"hist ][");
325 hRat->Draw(
"E1 same");
327 plotsOut->WriteTObject(canvas);
337 double xL = 0.5, xR = 0.7, yB = 0.62, yT = 0.87;
338 double texX = xL + .01;
339 double texY = yB - 0.045;
343 if(name.find(
"RecoELong") != std::string::npos) {
351 if(name.find(
"VtxX") != std::string::npos) {
359 if(name.find(
"VtxY") != std::string::npos) {
368 if(name.find(
"VtxZ") != std::string::npos) {
376 if(name.find(
"CVN") != std::string::npos) {
383 if(name.find(
"NHit") != std::string::npos) {
388 if(name.find(
"PTP") != std::string::npos) {
396 if(name.find(
"EpH") != std::string::npos) {
400 if(name.find(
"DistTop") != std::string::npos) {
void PlotStack(Spectrum spectra[], TDirectory *rootOut, FILE *textOFS, std::string name, std::string det, double POT, double potEquiv, PlotOptions opt, Spectrum *dataspec)
void SetHistOptions(TH1 *h, double max, std::string title, int ndiv, Color_t col, bool fill)
Set common options for a TLegend.
T max(const caf::Proxy< T > &a, T b)
Cuts and Vars for the 2020 FD DiF Study.
Spectrum NCTotalComponent() const override
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.
void CenterTitles(TH1 *histo)
TCanvas * canvas(const char *nm, const char *ti, const char *a)
const Color_t kNumuBackgroundColor
void make_plots(const ana::Spectrum &data, const ana::CheatDecomp &MC, std::string plotlabel, TFile *plotsOut)
Representation of a spectrum in any variable, with associated POT.
const XML_Char const XML_Char * data
static std::unique_ptr< CheatDecomp > LoadFrom(TDirectory *dir, const std::string &name)
Spectrum NueComponent() const override
Spectrum AntiNumuComponent() const override
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
const Color_t kBeamNueBackgroundColor
Spectrum NumuComponent() const override
void make_dataMC(TFile *SavedFile, TFile *OutFile)
void SetLegendOptions(TLegend *leg)
Set common options for a TLegend.
void DataMCNDAna_nus17(std::string folder)
const Color_t kTotalMCColor
T min(const caf::Proxy< T > &a, T b)
const Color_t kNCBackgroundColor
Just return the ND truth spectra as the decomposition.
PlotOptions GetPlotOptions(std::string name)
Spectrum AntiNueComponent() const override