test_caf_validation_plots.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Spectrum.h"
8 
9 #include "TH1.h"
10 #include "TFile.h"
11 #include "TLegend.h"
12 #include "TCanvas.h"
13 #include "Rtypes.h"
15 
16 using namespace ana;
17 
18 struct HistDef {
21 };
22 struct CutDef {
25  Color_t color;
26 };
27 
28 /////////////////////////////////////////////////////////////////////////////
29 std::vector<HistDef> InitHistDefs()
30 {
31  std::vector<HistDef> ret;
32 
33  ret.push_back({"trkl",{"Trk L (cm)",Binning::Simple(100,0,20),kTrkLength}});
34  ret.push_back({"numue",{"E_{#nu} (reco)",Binning::Simple(50,0,5),kCCE}});
35  ret.push_back({"hade",{"E_{had} (reco)",Binning::Simple(50,0,2),kHadE}});
36 
37  return ret;
38 }
39 
40 /////////////////////////////////////////////////////////////////////////////
41 std::vector<CutDef> InitCutDefs()
42 {
43  std::vector<CutDef> ret;
44 
45  ret.push_back({"#nu_{#mu} CC", kIsNumuCC && !kIsAntiNu, kBlue});
46  ret.push_back({"NC", kIsNC && !kIsAntiNu, kMagenta});
47  ret.push_back({"#nu_{e} CC", kIsBeamNue && !kIsAntiNu, kGreen+2});
48  ret.push_back({"anti-#nu", kIsAntiNu, kRed});
49  return ret;
50 }
51 
52 /////////////////////////////////////////////////////////////////////////////
53 std::vector< std::vector<Spectrum*> > InitSpectra(SpectrumLoader &load,
54  Cut sel,
55  std::vector<HistDef> defs,
56  std::vector<CutDef> types)
57 {
58  std::vector< std::vector<Spectrum*> > ret;
59  for (int dIdx = 0; dIdx < int(defs.size()); dIdx++){
60  std::vector<Spectrum*> cur;
61  for (int tIdx = 0; tIdx < int(types.size()); tIdx++){
62  cur.push_back(new Spectrum(load,defs[dIdx].axis,sel&&types[tIdx].cut));
63  }
64  ret.push_back(cur);
65  }
66  return ret;
67 }
68 
69 /////////////////////////////////////////////////////////////////////////////
70 void MakePlot(std::vector<Spectrum*> spect1, std::vector<Spectrum*> spect2,
71  std::string canname, std::vector<CutDef> cuts)
72 {
73  assert(spect1.size()>0 && spect2.size()==spects1.size() && cuts.size()==spects1.size());
74 
75  // Dataset 1 hists are solid, 2 are dashed
76  Style_t kStyle1 = 1;
77  Style_t kStyle2 = 7;
78  // Take POT from second dataset - assume it's the smaller set
79  double pot = spect2[0]->POT();
80 
81  std::vector<TH1*> h1, h2;
82  for (int i = 0; i < int(cuts.size()); i++){
83  h1.push_back(spect1[i]->ToTH1(pot,cuts[i].color,kStyle1));
84  h2.push_back(spect2[i]->ToTH1(pot,cuts[i].color,kStyle2));
85  }
86 
87  TCanvas *can = new TCanvas();
88 
89  for (int i = 0; i < int(cuts.size()); i++){
90  h1[i]->Draw("hist,same");
91  h2[i]->Draw("hist,same");
92  }
93 
94  TLegend* leg = AutoPlaceLegend(0.3,0.2);
95  for (int i = 0; i < int(cuts.size()); i++){
96  leg->AddEntry(h1[i],(cuts[i].name+" dataset 1").c_str(),"l");
97  leg->AddEntry(h2[i],(cuts[i].name+" dataset 2").c_str(),"l");
98  }
99  leg->Draw();
100  can->SaveAs(("/nusoft/app/web/htdoc/nova/users/dpershey/crude_caf_validation/SA_FD_NonSwap_MC_BirksBVsBirksC/"+canname+".png").c_str());
101 }
102 
103 
104 /////////////////////////////////////////////////////////////////////////////
106  bool runSmallSample=false)
107 {
108  std::string valid_dir = "/nusoft/app/web/htdoc/nova/users/dpershey/crude_caf_validation/"+dataset1+" vs "+dataset2;
109 
110  std::string name1 =
111  runSmallSample ? "defname: "+dataset1+" with stride 20" : dataset1;
112  std::string name2 =
113  runSmallSample ? "defname: "+dataset2+" with stride 20" : dataset2;
114 
115  SpectrumLoader load1(name1);
116  SpectrumLoader load2(name2);
117 
118  // Set up cuts and axes
119  Cut sel = kNumuFD;
120  std::vector<HistDef> defs = InitHistDefs();
121  std::vector<CutDef> cuts = InitCutDefs();
122  // Make some spectra
123  std::vector< std::vector<Spectrum*> > spect1 = InitSpectra(load1, sel,
124  defs, cuts);
125  std::vector< std::vector<Spectrum*> > spect2 = InitSpectra(load2, sel,
126  defs, cuts);
127  assert(spect.size() > 0);
128 
129  // Fill your spectra
130  load1.Go();
131  load2.Go();
132 
133  // Assume that the second dataset is smaller, scale hists to that POT
134  for (int dIdx = 0; dIdx < int(defs.size()); dIdx++)
135  MakePlot(spect1[dIdx],spect2[dIdx],defs[dIdx].name,cuts);
136 
137 }
const Var kHadE
Definition: NumuVars.h:23
const XML_Char * name
Definition: expat.h:151
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< CutDef > InitCutDefs()
void load(std::string lib)
Definition: load_libs.C:3
const Cut kIsAntiNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return sr->mc.nu[0].pdg< 0;})
Is this truly an antineutrino?
Definition: TruthCuts.h:53
const Cut kIsBeamNue(CCFlavSel(12, 12))
Select CC .
std::vector< HistDef > InitHistDefs()
const Cut kNumuFD
Definition: NumuCuts.h:53
const Var kTrkLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].len/100;})
Definition: NumuVars.h:65
const Cut kIsNC([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return!sr->mc.nu[0].iscc;})
Is this a Neutral Current event?
Definition: TruthCuts.h:8
void MakePlot(std::vector< Spectrum * > spect1, std::vector< Spectrum * > spect2, std::string canname, std::vector< CutDef > cuts)
const Var kCCE
Definition: NumuVars.h:21
#define pot
Struct to hold cut information.
virtual void Go() override
Load all the registered spectra.
const HistDef defs[kNumVars]
Definition: vars.h:15
std::vector< float > Spectrum
Definition: Constants.h:527
HistAxis axis
Definition: NuePlotLists.h:13
TH1F * h2
Definition: plot.C:45
std::vector< std::vector< Spectrum * > > InitSpectra(SpectrumLoader &load, Cut sel, std::vector< HistDef > defs, std::vector< CutDef > types)
TH1F * h1
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
assert(nhit_max >=nhit_nbins)
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:715
void test_caf_validation_plots(std::string dataset1, std::string dataset2, bool runSmallSample=false)
std::string name
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38