29 if (sr->
slc.
nhit == 0)
return float(-1);
47 float eph = (
nhit > 0) ? cale/
nhit : -1;
57 std::vector<Spectrum*>
ret;
68 cuts.push_back(cuts.back() &&
kCVNe>0.75);
78 std::vector<HistAxis>
axes;
93 std::vector<Spectrum*>
ret;
104 std::vector<Spectrum*>
ret;
155 TFile *
out =
new TFile(fname.c_str(),
"recreate");
157 std::vector<std::string> cutname = {
"nocut",
"dq",
"fiducial",
"contain",
158 "presel",
"ptp",
"cvn"};
159 assert(cutname.size() == da_cutflow2.size());
160 for (
int cIdx = 0; cIdx <
int(cutname.size()); cIdx++){
161 mc_cutflow2[cIdx]->SaveTo(out, (
"prod2_mc_"+cutname[cIdx]).c_str());
162 mc_cutflow3[cIdx]->SaveTo(out, (
"prod3_mc_"+cutname[cIdx]).c_str());
164 da_cutflow2[cIdx]->SaveTo(out, (
"prod2_da_"+cutname[cIdx]).c_str());
165 da_cutflow3[cIdx]->SaveTo(out, (
"prod3_da_"+cutname[cIdx]).c_str());
168 std::vector<std::string> varname = {
"calE",
"NHit",
"EPH",
"shwE",
169 "shwEPH",
"HadE",
"pTp"};
170 assert(varname.size() == da_kin2.size());
171 for (
int kIdx = 0; kIdx <
int(varname.size()); kIdx++){
172 mc_kin2[kIdx]->SaveTo(out, (
"prod2_mc_"+varname[kIdx]).c_str());
173 mc_kin3[kIdx]->SaveTo(out, (
"prod3_mc_"+varname[kIdx]).c_str());
174 mc_2D2[kIdx]->SaveTo(out, (
"prod2_mc_runvs"+varname[kIdx]).c_str());
175 mc_2D3[kIdx]->SaveTo(out, (
"prod3_mc_runvs"+varname[kIdx]).c_str());
177 da_kin2[kIdx]->SaveTo(out, (
"prod2_da_"+varname[kIdx]).c_str());
178 da_kin3[kIdx]->SaveTo(out, (
"prod3_da_"+varname[kIdx]).c_str());
179 da_2D2[kIdx]->SaveTo(out, (
"prod2_da_runvs"+varname[kIdx]).c_str());
180 da_2D3[kIdx]->SaveTo(out, (
"prod3_da_runvs"+varname[kIdx]).c_str());
186 h->GetXaxis()->CenterTitle();
187 h->GetYaxis()->CenterTitle();
192 TCanvas *
can =
new TCanvas();
196 double da_pot2 = da_prod2->
POT();
197 double da_pot3 = da_prod3->
POT();
200 double mc_pot2 = mc_prod2->
POT();
201 double mc_pot3 = mc_prod3->
POT();
203 TH1 *da2 = da_prod2->
ToTH1(da_pot3);
204 TH1 *da3 = da_prod3->
ToTH1(da_pot3);
205 TH1 *mc2 = mc_prod2->
ToTH1(da_pot3);
206 TH1 *mc3 = mc_prod3->
ToTH1(da_pot3);
208 TH1 *rat_da2 = da_prod2->
ToTH1(da_pot3);
209 TH1 *rat_da3 = da_prod3->
ToTH1(da_pot3);
210 TH1 *rat_mc3 = mc_prod3->
ToTH1(da_pot3);
213 if (name ==
"EPH" || name ==
"shwEPH") nrebin = 10;
214 if (name ==
"calE" || name ==
"shwE" ||
215 name ==
"HadE" || name ==
"NHit") nrebin = 20;
220 rat_da2->Rebin(nrebin);
221 rat_da3->Rebin(nrebin);
222 rat_mc3->Rebin(nrebin);
224 rat_da2->Divide(mc2);
225 rat_da3->Divide(mc2);
226 rat_mc3->Divide(mc2);
228 TPad *pRat =
new TPad((
"rat"+name) .c_str(),
"",0,0,1,0.35);
229 TPad *pSpec=
new TPad((
"spec"+name).c_str(),
"",0,0.35,1,1);
235 gPad->SetBottomMargin(0);
236 mc2->SetLineColor(
kBlue);
237 mc3->SetLineColor(
kRed);
238 da2->SetLineColor(
kGreen+2);
239 da3->SetLineColor(kBlack);
240 da2->SetMinimum(1
e-3);
242 double max =
std::max(da2->GetMaximum(),da3->GetMaximum());
243 max =
std::max(max,mc2->GetMaximum());
244 max =
std::max(max,mc3->GetMaximum());
246 da2->SetMaximum(1.2*max);
247 da2->GetXaxis()->SetTitle(
"");
248 da2->GetYaxis()->SetTitle(
"Events");
251 da2->GetXaxis()->SetTitleSize(0.06);
256 TLegend *
leg =
new TLegend(0.6,0.65,0.85,0.85);
257 leg->AddEntry(da2,
"Data prod2",
"le");
258 leg->AddEntry(mc2,
"MC prod2",
"le");
259 leg->AddEntry(da3,
"Data prod3",
"le");
260 leg->AddEntry(mc3,
"MC prod3",
"le");
263 rat_mc3->SetLineColor(
kRed);
264 rat_da2->SetLineColor(
kGreen+2);
265 rat_da3->SetLineColor(kBlack);
268 gPad->SetTopMargin(0);
269 gPad->SetBottomMargin(0.2);
270 rat_mc3->GetXaxis()->SetTitle(xtitle.c_str());
271 rat_mc3->GetYaxis()->SetTitle(
"Various Ratios");
272 rat_mc3->SetTitle(
"");
274 rat_mc3->GetXaxis()->SetTitleSize(0.10);
275 rat_mc3->GetXaxis()->SetLabelSize(0.08);
276 rat_mc3->GetXaxis()->SetTitleOffset(0.8);
277 rat_mc3->GetYaxis()->SetTitleSize(0.10);
278 rat_mc3->GetYaxis()->SetLabelSize(0.08);
279 rat_mc3->GetYaxis()->SetTitleOffset(0.5);
280 rat_mc3->GetYaxis()->SetRangeUser(0.5001,1.4999);
282 rat_da2->Draw(
"same");
283 rat_da3->Draw(
"same");
285 can->SaveAs((
"Fig_nddatastability_"+name+
".png").c_str());
290 TCanvas *
can =
new TCanvas();
296 double pot = da_prod3->
POT();
298 TH2F *mc2 = (TH2F*)mc_prod2->
ToTH2(pot);
299 TH2F *mc3 = (TH2F*)mc_prod3->
ToTH2(pot);
300 TH2F *da2 = (TH2F*)da_prod2->
ToTH2(pot);
301 TH2F *da3 = (TH2F*)da_prod3->
ToTH2(pot);
304 if (name ==
"HadE" || name ==
"shwE")
311 TH1D *pmc2 = (TH1D*)mc2->ProfileX((name+
"pmc2").c_str());
312 TH1D *pmc3 = (TH1D*)mc3->ProfileX((name+
"pmc3").c_str());
313 TH1D *pda2 = (TH1D*)da2->ProfileX((name+
"pda2").c_str());
314 TH1D *pda3 = (TH1D*)da3->ProfileX((name+
"pda3").c_str());
316 pmc2->SetLineColor(
kBlue);
317 pmc3->SetLineColor(
kRed);
318 pda2->SetLineColor(
kGreen+2);
319 pda3->SetLineColor(kBlack);
321 pda2->GetXaxis()->SetTitle(
"Run");
322 pda2->GetYaxis()->SetTitle(ytitle.c_str());
327 for (
int i = 1;
i <= pda3->GetNbinsX();
i++){
328 double curcont = pda3->GetBinContent(
i);
335 pda2->GetYaxis()->SetRangeUser(0.5*avgval,1.5*avgval);
342 TLegend *
leg =
new TLegend(0.6,0.65,0.85,0.85);
343 leg->AddEntry(pda2,
"Data prod2",
"le");
344 leg->AddEntry(pmc2,
"MC prod2",
"le");
345 leg->AddEntry(pda3,
"Data prod3",
"le");
346 leg->AddEntry(pmc3,
"MC prod3",
"le");
350 can->SaveAs((
"Fig_nddatastability_"+name+
"_vsrun.png").c_str());
361 Make2DPlot(fname,
"shwEPH",
"ShwE / ShwNHit [GeV]");
365 MakePlot(fname,
"calE",
"CalE [GeV]");
367 MakePlot(fname,
"EPH",
"CalE / NHit [GeV]");
368 MakePlot(fname,
"HadE",
"HadE [GeV]");
369 MakePlot(fname,
"shwE",
"ShwE [GeV]");
370 MakePlot(fname,
"shwEPH",
"ShwE / ShwNHit [GeV]");
372 MakePlot(fname,
"nocut",
"CalE [GeV]");
374 MakePlot(fname,
"fiducial",
"CalE [GeV]");
375 MakePlot(fname,
"contain",
"CalE [GeV]");
376 MakePlot(fname,
"presel",
"CalE [GeV]");
caf::Proxy< unsigned int > nshwlid
T max(const caf::Proxy< T > &a, T b)
caf::Proxy< caf::SRFuzzyK > fuzzyk
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
_HistAxis< Var > HistAxis
Cuts and Vars for the 2020 FD DiF Study.
const Cut kNue2017NDEnergy
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.
Proxy for caf::StandardRecord.
std::vector< Spectrum * > MakeKinSpectra(SpectrumLoader *l, Var wei, Cut sel)
const Cut kNueDQ2017CVN([](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;return true;})
void SetSpillCut(const SpillCut &cut)
void CenterTitles(TH1 *histo)
void Make2DPlot(std::string fname, std::string name, std::string ytitle)
const Cut kNue2017NDProngLength
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
const Var kHadCalE([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return float(sr->slc.calE);if(sr->vtx.elastic.fuzzyk.nshwlid< 1) return float(sr->slc.calE);return((sr->slc.calE- sr->vtx.elastic.fuzzyk.png[0].shwlid.calE));})
Representation of a spectrum in any variable, with associated POT.
caf::Proxy< caf::SRElastic > elastic
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
void CreateAndFillSpectra(std::string fname)
const Var kPtP
Transverse momentum fraction in slice.
caf::Proxy< unsigned int > nhit
const Cut kNue2017NDNHits
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
const Var kShwLIDE([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return float(-1);if(sr->vtx.elastic.fuzzyk.nshwlid< 1) return float(-1);return float(sr->vtx.elastic.fuzzyk.png[0].shwlid.calE);})
void MakePlot(std::string fname, std::string name, std::string xtitle)
virtual void Go() override
Load all the registered spectra.
const Cut kNue2017NDFrontPlanes
No hits in the first 5 planes of the detector, doc 12879.
std::vector< float > Spectrum
std::vector< HistAxis > InitAxes()
const SystShifts kNoShift
const Cut kNue2017NDPresel
caf::Proxy< bool > IsValid
caf::Proxy< caf::SRSlice > slc
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Cut kNue2017NDContain([](const caf::SRProxy *sr){for(unsigned int ix=0;ix< sr->vtx.elastic.fuzzyk.nshwlid;++ix){TVector3 start=sr->vtx.elastic.fuzzyk.png[ix].shwlid.start;TVector3 stop=sr->vtx.elastic.fuzzyk.png[ix].shwlid.stop;if(std::min(start.X(), stop.X())< -170.0) return false;if(std::max(start.X(), stop.X()) > 170.0) return false;if(std::min(start.Y(), stop.Y())< -170.0) return false;if(std::max(start.Y(), stop.Y()) > 170.0) return false;if(std::min(start.Z(), stop.Z())< 100.0) return false;if(std::max(start.Z(), stop.Z()) > 1225.0) return false;}return true;})
Loose containtment on start and end of all showers, docdb-12943.
std::vector< Spectrum * > MakeKin2DSpectra(SpectrumLoader *l, Var wei, Cut sel)
assert(nhit_max >=nhit_nbins)
std::vector< Spectrum * > MakeCutFlowSpectra(SpectrumLoader *l, Var wei)
const SpillCut kStandardSpillCuts
Apply this unless you're doing something special.
caf::Proxy< caf::SRVertexBranch > vtx
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
void nd_datastability(bool fillSpectra, std::string fname)
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
const Var kEPH([](const caf::SRProxy *sr){if(sr->slc.nhit==0) return float(-1);return float(sr->slc.calE/sr->slc.nhit);})
void MakePlots(std::string fname)
const Cut kNue2017NDFiducial([](const caf::SRProxy *sr){assert(sr->vtx.elastic.IsValid &&"Must apply DQ cuts");if(sr->vtx.elastic.vtx.X()< -100.0) return false;if(sr->vtx.elastic.vtx.X() > 160.0) return false;if(sr->vtx.elastic.vtx.Y()< -160.0) return false;if(sr->vtx.elastic.vtx.Y() > 100.0) return false;if(sr->vtx.elastic.vtx.Z()< 150.0) return false;if(sr->vtx.elastic.vtx.Z() > 900.0) return false;return true;})
Cut on ND vertex position.
const Var kShwLIDEPH([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return float(-1);if(sr->vtx.elastic.fuzzyk.nshwlid< 1) return float(-1);float cale=sr->vtx.elastic.fuzzyk.png[0].shwlid.calE;int nhit=sr->vtx.elastic.fuzzyk.png[0].shwlid.nhit;float eph=(nhit > 0)?cale/nhit:-1;return eph;})