19 #include "TLegendEntry.h" 26 #include "NDAna/muonid/NDXSecMuonPID.h" 34 const int ret = gSystem->Load(
"libNDAnamuonid");
48 void DrawArrow(TPad *
p,
double x, TString arrow_opt,
int line_color = kGray+1,
int line_style = 1,
int line_width = 3);
56 std::map<std::string, const HistAxis*>
axes;
57 std::map<std::string, Bounds_t>
bounds;
65 bounds[
"top" ] = {50 };
66 bounds[
"bottom"] = {30 };
67 bounds[
"east" ] = {50 };
68 bounds[
"west" ] = {30 };
69 bounds[
"front" ] = {150};
71 axes[
"vtxx"] =
new HistAxis(
"Vertex in X Direction (cm)",
72 Binning::Simple(40, -200, 200),
74 axes[
"vtxy"] =
new HistAxis(
"Vertex in Y Direction (cm)",
75 Binning::Simple(40, -200, 200),
78 axes[
"vtxz"] =
new HistAxis(
"Vertex in Z Direction (cm)",
79 Binning::Simple(80, 0, 1200),
81 axes[
"top"] =
new HistAxis(
"Min. Prong Distance from Top (cm)",
82 Binning::Simple(40, 0, 400),
84 axes[
"bottom"] =
new HistAxis(
"Min. Prong Distance from Bottom (cm)",
85 Binning::Simple(40, 0, 400),
87 axes[
"east"] =
new HistAxis(
"Min. Prong Distance from East (cm)",
88 Binning::Simple(40, 0, 400),
90 axes[
"west" ] =
new HistAxis(
"Min. Prong Distance from West (cm)",
91 Binning::Simple(40, 0, 400),
93 axes[
"front"] =
new HistAxis(
"Min. Prong Distance from Front (cm)",
94 Binning::Simple(80, 0, 1200),
97 std::map<std::string, const Cut *>
cuts;
105 if(input_file_name ==
"") {
106 std::map<std::string, SpectrumLoader*>
loaders;
113 std::map<std::string, Spectrum*>
specs;
115 for(
auto axis_it = axes.begin(); axis_it != axes.end(); axis_it++) {
116 for(
auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
117 key =
"data_" + cut_it->first +
"_" + axis_it->first;
135 loaders.at(
"data")->Go();
136 loaders.at(
"mc" )->Go();
138 TFile *
out =
new TFile(
"datamc.root",
"recreate");
139 for(
auto spec_it = specs.begin(); spec_it != specs.end(); spec_it++) {
140 spec_it->second->SaveTo(out, spec_it->first.c_str());
148 TFile * input = TFile::Open(input_file_name.c_str());
150 std::map<std::string, TH1D*> hists_data;
151 std::map<std::string, TH1D*> hists_mc;
153 for(
auto axis_it = axes.begin(); axis_it != axes.end(); axis_it++) {
154 for(
auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
155 key =
"data_" + cut_it->first +
"_" + axis_it->first;
158 hists_data[
key] =
tmp->ToTH1(POT);
162 hists_mc[
key] =
tmp->ToTH1(POT);
169 for(
auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
170 key =
"data_" + cut_it->first +
"_" +
axis;
171 hists_data.at(key)->Rebin(4);
174 hists_mc.at(key)->Rebin(4);
179 for(
auto axis_it = axes.begin(); axis_it != axes.end(); axis_it++) {
180 for(
auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
182 TPad *
p1 =
new TPad(
"p1",
"", 0, 0, 1, 1);
183 TPad *
p2 =
new TPad(
"p2",
"", 0, 0, 1, 1);
184 p1->SetLeftMargin(0.15);
185 p2->SetLeftMargin(0.15);
186 p1->SetBottomMargin(0.35);
187 p2->SetTopMargin(0.65);
190 TLegend *
leg =
new TLegend(.65, .89, .89, .65,
"",
"NDC");
192 leg->SetFillStyle(0);
195 key =
"data_" + cut_it->first +
"_" + axis_it->first;
196 auto data = hists_data.at(key);
198 data->SetMarkerStyle(8);
199 data->SetMarkerSize(1);
200 data->GetYaxis()->CenterTitle();
201 data->GetXaxis()->CenterTitle();
202 data->GetYaxis()->SetTitleOffset(1.15);
205 data->GetYaxis()->SetRangeUser(0.0001,
data->GetMaximum() * 1.5);
208 data->SetLabelSize(0);
211 THStack * stack =
new THStack((key +
"_stack").c_str(),
"");
214 hists_mc.at(
"mc_" + cut_it->first +
"_" + axis_it->first +
"_nuebarcc")->Add(hists_mc.at(
"mc_" + cut_it->first +
"_" + axis_it->first +
"_nuebarcc_nonfiducial"));
217 if(ichn == 2)
continue;
222 stack->Add(hists_mc.at(key));
224 leg->AddEntry((TObject*)0,
"",
"");
225 leg->AddEntry(
data,
"Data",
"lp");
226 stack->Draw(
"hist same");
227 data->Draw(
"ep same");
232 ratio->Divide(hists_mc.at(
"mc_" + cut_it->first +
"_" + axis_it->first +
"_mc"));
233 ratio->GetYaxis()->SetRangeUser(.5, 1.5);
234 ratio->GetYaxis()->SetTitle(
"Data / MC");
242 DrawArrow(p1, bounds.at(axis_it->first).min,
">");
243 if(bounds.at(axis_it->first).max)
244 DrawArrow(p1, bounds.at(axis_it->first).max,
"<");
246 DrawBounds(p2, bounds.at(axis_it->first).min);
247 if(bounds.at(axis_it->first).max)
248 DrawBounds(p2, bounds.at(axis_it->first).max);
250 c->SaveAs((plot_dump +
"/datamc_" + cut_it->first +
"_" + axis_it->first +
".pdf").c_str());
260 TGraph*
one =
new TGraph();
261 one->SetPoint(0, -1e6, 1);
262 one->SetPoint(1, 1e6, 1);
263 one->SetLineStyle(7);
264 one->SetLineColor(kGray+3);
273 TGraph*
g =
new TGraph();
274 g->SetPoint(0, x, p->GetY1());
275 g->SetPoint(1, x, p->GetY2());
276 g->SetLineStyle(line_style);
277 g->SetLineColor(line_color);
278 g->SetLineWidth(line_width);
289 bool left = arrow_opt.Contains(
"<");
290 double offset = 0.07 * (p->GetX2() - p->GetX1());
291 if(left) offset *= -1;
293 double y = p->GetY1() + 0.8 * (p->GetY2() - p->GetY1());
297 TLine *
l1 =
new TLine(x, 0, x, y);
298 TArrow * arrow =
new TArrow(x, y, x+offset, y, 0.04, arrow_opt.Data());
299 if(left) arrow =
new TArrow(x+offset, y, x, y, 0.04, arrow_opt.Data());
301 arrow->SetLineColor(line_color);
302 l1 ->SetLineColor(line_color);
304 arrow->SetLineStyle(line_style);
305 l1 ->SetLineStyle(line_style);
307 arrow->SetLineWidth(line_width);
308 l1 ->SetLineWidth(line_width);
const ana::Cut kCVNNueIDCut
::xsd::cxx::tree::bounds< char > bounds
_HistAxis< Var > HistAxis
Cuts and Vars for the 2020 FD DiF Study.
const ana::Var kNHits([](const caf::SRProxy *sr){return sr->slc.nhit;})
void DrawBounds(TPad *p, double x, int line_color=kGray+1, int line_style=1, int line_width=3)
const Var kDistAllBottom([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngbottom)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngbottom);})
Distance of all showers in slice from the bottom edge of detector.
const TVector3 vNueCCIncFiducialMin(-130,-140, 150)
const Var kDistAllWest([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngwest)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngwest);})
Distance of all showers in slice from the west edge of detector.
const Var kDistAllTop([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngtop)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngtop);})
Distance of all showers in slice from the top edge of detector.
void datamc(std::string input_file_name="", std::string plot_dump="./plots/")
const ana::Cut kFrontPlanes
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
void DrawArrow(TPad *p, double x, TString arrow_opt, int line_color=kGray+1, int line_style=1, int line_width=3)
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
const SelDef chns[kNumChns]
const XML_Char const XML_Char * data
const Var kDistAllEast([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngeast)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngeast);})
Distance of all showers in slice from the east edge of detector.
const ana::Cut kNuebarCCIncContainmentLoose([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.nshwlid< 1) return false;TVector3 start=sr->vtx.elastic.fuzzyk.png[0].shwlid.start;TVector3 stop=sr->vtx.elastic.fuzzyk.png[0].shwlid.stop;for(uint ix=0;ix< sr->vtx.elastic.fuzzyk.nshwlid;ix++){TVector3 start_muoncatcher=sr->vtx.elastic.fuzzyk.png[ix].shwlid.start;TVector3 stop_muoncatcher=sr->vtx.elastic.fuzzyk.png[ix].shwlid.stop;if(std::max(start_muoncatcher.Z(), stop_muoncatcher.Z()) > 1250.0) return false;}if(sr->sel.nuecosrej.distallpngtop< 30) return false;if(sr->sel.nuecosrej.distallpngbottom< 10) return false;if(sr->sel.nuecosrej.distallpngeast< 30) return false;if(sr->sel.nuecosrej.distallpngwest< 10) return false;if(sr->sel.nuecosrej.distallpngfront< 100) return false;return true;})
const std::string PROD5_MC_RHC_NOMINAL
const std::string PROD5_RHC_DATA
std::vector< float > Spectrum
const SystShifts kNoShift
std::vector< double > POT
const TVector3 vNueCCIncFiducialMax(150, 140, 800)
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
const Var kDistAllFront([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngfront)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngfront);})
Distance of all showers in slice from the front edge of detector.
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Style chns_style[kNumChns]
const ana::Cut kDQ([](const caf::SRProxy *sr){if(sr->sel.nuecosrej.hitsperplane >=8) return false;if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.fuzzyk.npng==0) return false;if(sr->vtx.elastic.fuzzyk.nshwlid==0) return false;if(sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx<=5||sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity<=5) return false;if(sr->vtx.elastic.fuzzyk.png[0].shwlid.gap > 100) return false; if(sr->vtx.elastic.fuzzyk.nshwlid >=2){if((sr->vtx.elastic.fuzzyk.png[0].shwlid.dir). Dot(sr->vtx.elastic.fuzzyk.png[1].shwlid.dir)< -0.95) return false;}float xyhitdiff=abs(sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx-sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity);float xyhitsum=sr->vtx.elastic.fuzzyk.png[0].shwlid.nhitx+sr->vtx.elastic.fuzzyk.png[0].shwlid.nhity;if(xyhitdiff/xyhitsum >=0.4) return false;float nhitall=0.0;for(unsigned int j=0;j< sr->vtx.elastic.fuzzyk.nshwlid;j++){nhitall+=sr->vtx.elastic.fuzzyk.png[j].shwlid.nhit;}if(nhitall/sr->slc.nhit<=0.7) return false;return true;})
const ana::Cut kMuonIDProd4Cut
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
std::string UniqueName()
Return a different string each time, for creating histograms.
const ana::Cut kNuebarCCIncFiducialLoose([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->vtx.elastic.vtx.X()< -165.0) return false;if(sr->vtx.elastic.vtx.X() > 165.0) return false;if(sr->vtx.elastic.vtx.Y()< -165.0) return false;if(sr->vtx.elastic.vtx.Y() > 165.0) return false;if(sr->vtx.elastic.vtx.Z()< 100.0) return false;if(sr->vtx.elastic.vtx.Z() > 1000.0) return false;return true;})