Classes | Namespaces | Functions | Variables
datamc.C File Reference
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncCuts.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncExtra.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Vars/XsecTunes.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "TFile.h"
#include "THStack.h"
#include "TH1D.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TPad.h"
#include "TGraph.h"
#include "TArrow.h"
#include "NDAna/muonid/NDXSecMuonPID.h"

Go to the source code of this file.

Classes

struct  Bounds_t
 

Namespaces

 nuebarccinc
 

Functions

void nuebarccinc::One ()
 
void load_libs_muonid ()
 
void DrawArrow (TPad *p, double x, TString arrow_opt, int line_color=kGray+1, int line_style=1, int line_width=3)
 
void DrawBounds (TPad *p, double x, int line_color=kGray+1, int line_style=1, int line_width=3)
 
void datamc (std::string input_file_name="", std::string plot_dump="./plots/")
 

Variables

const ana::Cut nuebarccinc::kMuonIDProd4Cut = ana::kMuonIDProd4 < -0.2
 

Function Documentation

void datamc ( std::string  input_file_name = "",
std::string  plot_dump = "./plots/" 
)

Definition at line 52 of file datamc.C.

References file_size_ana::axes, allInOneTrainingPlots::axis, plot_validation_datamc::c, nuebarccinc::chns, nuebarccinc::chns_style, plot_validation_datamc_2018::color, Cut(), nuebarccinc::SelDef::cut, update_sam_good_runs_metadata::cuts, DrawArrow(), DrawBounds(), APDHVSetting::dummy, HTMLTools::entry(), nuebarccinc::kCVNNueIDCut, ana::kDistAllBottom, ana::kDistAllEast, ana::kDistAllFront, ana::kDistAllTop, ana::kDistAllWest, nuebarccinc::kDQ, findDuplicateFiles::key, nuebarccinc::kFrontPlanes, nuebarccinc::kMuonIDProd4Cut, nuebarccinc::kNHits, ana::kNoShift, nuebarccinc::kNuebarCCIncContainmentLoose, nuebarccinc::kNuebarCCIncFiducialLoose, nuebarccinc::kNumChns, ana::kPPFXFluxCVWgt, ana::kVtxX, ana::kVtxY, ana::kVtxZ, ana::kXSecCVWgt2020, latex(), MECModelEnuComparisons::leg, load_libs_muonid(), loaders, ana::LoadFrom(), nuebarccinc::SelDef::name, nuebarccinc::One(), confusionMatrixTree::out, plot_validation_datamc::p1, plot_validation_datamc::p2, POT, nuebarccinc::PROD5_MC_RHC_NOMINAL, nuebarccinc::PROD5_RHC_DATA, PandAna.reco_validation.prod5_pid_validation::ratio(), PandAna.Demos.tute_pid_validation::specs, string, tmp, ana::UniqueName(), nuebarccinc::vNueCCIncFiducialMax(), and nuebarccinc::vNueCCIncFiducialMin().

Referenced by BinErrors().

54 {
56  std::map<std::string, const HistAxis*> axes;
57  std::map<std::string, Bounds_t> bounds;
58 
59  bounds["vtxx"] = {nuebarccinc::vNueCCIncFiducialMin.X(),
61  bounds["vtxy"] = {nuebarccinc::vNueCCIncFiducialMin.Y(),
63  bounds["vtxz"] = {nuebarccinc::vNueCCIncFiducialMin.Z(),
65  bounds["top" ] = {50 };
66  bounds["bottom"] = {30 };
67  bounds["east" ] = {50 };
68  bounds["west" ] = {30 };
69  bounds["front" ] = {150};
70 
71  axes["vtxx"] = new HistAxis("Vertex in X Direction (cm)",
72  Binning::Simple(40, -200, 200),
73  kVtxX);
74  axes["vtxy"] = new HistAxis("Vertex in Y Direction (cm)",
75  Binning::Simple(40, -200, 200),
76  kVtxY);
77 
78  axes["vtxz"] = new HistAxis("Vertex in Z Direction (cm)",
79  Binning::Simple(80, 0, 1200),
80  kVtxZ);
81  axes["top"] = new HistAxis("Min. Prong Distance from Top (cm)",
82  Binning::Simple(40, 0, 400),
83  kDistAllTop);
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),
89  kDistAllEast);
90  axes["west" ] = new HistAxis("Min. Prong Distance from West (cm)",
91  Binning::Simple(40, 0, 400),
92  kDistAllWest);
93  axes["front"] = new HistAxis("Min. Prong Distance from Front (cm)",
94  Binning::Simple(80, 0, 1200),
96 
97  std::map<std::string, const Cut *> cuts;
98  cuts["selection"] = new Cut(nuebarccinc::kDQ &&
104  cuts["selection_and_cvn"] = new Cut(*cuts.at("selection") && nuebarccinc::kCVNNueIDCut);
105  if(input_file_name == "") {
106  std::map<std::string, SpectrumLoader*> loaders;
108  loaders["data"] = new SpectrumLoader(nuebarccinc::PROD5_RHC_DATA);
109 
110  // central value weight to be applied to all samples
111  const Var * kCVWeight = new Var(kPPFXFluxCVWgt * kXSecCVWgt2020);
112 
113  std::map<std::string, Spectrum*> specs;
114  std::string key = "";
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;
118  specs[key] = new Spectrum(*loaders["data"],
119  *axis_it->second,
120  *cut_it->second,
121  kNoShift,
122  *kCVWeight);
123  for(auto ichn = 0; ichn < nuebarccinc::kNumChns; ichn++) {
124  key = "mc_" + cut_it->first + "_" + axis_it->first + "_" + nuebarccinc::chns[ichn].name;
125  specs[key] = new Spectrum(*loaders["mc"],
126  *axis_it->second,
127  *cut_it->second && nuebarccinc::chns[ichn].cut,
128  kNoShift,
129  *kCVWeight);
130 
131  }
132  }
133  }
134 
135  loaders.at("data")->Go();
136  loaders.at("mc" )->Go();
137 
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());
141  }
142 
143  out->Close();
144  }
145 
146  else {
147 
148  TFile * input = TFile::Open(input_file_name.c_str());
149 
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;
156  auto tmp = Spectrum::LoadFrom(input, key.c_str());
157  double POT = tmp->POT();
158  hists_data[key] = tmp->ToTH1(POT);
159  for(auto ichn = 0; ichn < nuebarccinc::kNumChns; ichn++) {
160  key = "mc_" + cut_it->first + "_" + axis_it->first + "_" + nuebarccinc::chns[ichn].name;
161  tmp = Spectrum::LoadFrom(input, key.c_str());
162  hists_mc[key] = tmp->ToTH1(POT);
163  }
164  }
165  }
166 
167  // rebin vtx z and min distance from front
168  for(std::string axis : {"front", "vtxz"}) {
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);
172  for(auto ichn = 0; ichn < nuebarccinc::kNumChns; ichn++) {
173  key = "mc_" + cut_it->first + "_" + axis + "_" + nuebarccinc::chns[ichn].name;
174  hists_mc.at(key)->Rebin(4);
175  }
176  }
177  }
178 
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++) {
181  TCanvas * c = new TCanvas(UniqueName().c_str());
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);
188  p1->SetFillStyle(0);
189  p2->SetFillStyle(0);
190  TLegend * leg = new TLegend(.65, .89, .89, .65, "", "NDC");
191  leg->SetNColumns(2);
192  leg->SetFillStyle(0);
193 
194  p1->cd();
195  key = "data_" + cut_it->first + "_" + axis_it->first;
196  auto data = hists_data.at(key);
197 
198  data->SetMarkerStyle(8);
199  data->SetMarkerSize(1);
200  data->GetYaxis()->CenterTitle();
201  data->GetXaxis()->CenterTitle();
202  data->GetYaxis()->SetTitleOffset(1.15);
203  auto ratio = (TH1D*) data->Clone(UniqueName().c_str());
204 
205  data->GetYaxis()->SetRangeUser(0.0001, data->GetMaximum() * 1.5);
206  data->Draw("e0 p");
207  // data->SetTickLength(0);
208  data->SetLabelSize(0);
209 
210  auto dummy = (TH1D*) data->Clone(UniqueName().c_str());
211  THStack * stack = new THStack((key + "_stack").c_str(), "");
212 
213  // add nonfiducial to nuebar
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"));
215 
216  for(auto ichn = 1; ichn < nuebarccinc::kNumChns; ichn++) {
217  if(ichn == 2) continue; // don't fill MC total and skip nonfiducial
218  key = "mc_" + cut_it->first + "_" + axis_it->first + "_" + nuebarccinc::chns[ichn].name;
219  hists_mc.at(key)->SetFillColor(nuebarccinc::chns_style[ichn].color);
220  auto entry = leg->AddEntry((TObject*)0, nuebarccinc::chns_style[ichn].latex.c_str(), "");
221  entry->SetTextColor(nuebarccinc::chns_style[ichn].color);
222  stack->Add(hists_mc.at(key));
223  }
224  leg->AddEntry((TObject*)0, "", "");
225  leg->AddEntry(data, "Data", "lp");
226  stack->Draw("hist same");
227  data->Draw("ep same");
228 
229  leg->Draw();
230 
231  p2->cd();
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");
235  ratio->Draw("ep");
237 
238  c->cd();
239  p1->Draw();
240  p2->Draw("same");
241 
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, "<");
245 
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);
249 
250  c->SaveAs((plot_dump + "/datamc_" + cut_it->first + "_" + axis_it->first + ".pdf").c_str());
251  }
252  }
253 
254 
255  }
256 }
void latex(double x, double y, std::string txt, double ang=0, int align=22, double size=0.042)
const ana::Cut kCVNNueIDCut
::xsd::cxx::tree::bounds< char > bounds
Definition: Database.h:226
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
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)
Definition: datamc.C:268
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.
Definition: NueVars.h:33
const TVector3 vNueCCIncFiducialMin(-130,-140, 150)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
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.
Definition: NueVars.h:36
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.
Definition: NueVars.h:30
const ana::Cut kFrontPlanes
Float_t tmp
Definition: plot.C:36
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
const Var kVtxX
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
void DrawArrow(TPad *p, double x, TString arrow_opt, int line_color=kGray+1, int line_style=1, int line_width=3)
Definition: datamc.C:285
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:74
const SelDef chns[kNumChns]
const XML_Char const XML_Char * data
Definition: expat.h:268
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.
Definition: NueVars.h:39
void Cut(double x)
Definition: plot_outliers.C:1
const int kNumChns
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
Definition: Constants.h:739
const SystShifts kNoShift
Definition: SystShifts.cxx:22
std::vector< double > POT
const TVector3 vNueCCIncFiducialMax(150, 140, 800)
void load_libs_muonid()
Definition: datamc.C:32
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.
Definition: NueVars.h:45
const Var kVtxY
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Style chns_style[kNumChns]
std::vector< Loaders * > loaders
Definition: syst_header.h:386
const Var kVtxZ
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
Definition: datamc.C:28
def entry(str)
Definition: HTMLTools.py:26
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
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;})
const Var kXSecCVWgt2020
Definition: XsecTunes.h:105
void One()
Definition: datamc.C:258
enum BeamMode string
void DrawArrow ( TPad *  p,
double  x,
TString  arrow_opt,
int  line_color = kGray+1,
int  line_style = 1,
int  line_width = 3 
)

Definition at line 285 of file datamc.C.

References om::cout, allTimeWatchdog::endl, central_limit::l1, art::left(), PandAna.reco_validation.add_data::offset, and submit_syst::y.

Referenced by datamc().

286 {
287  p->cd();
288  p->Update();
289  bool left = arrow_opt.Contains("<");
290  double offset = 0.07 * (p->GetX2() - p->GetX1());
291  if(left) offset *= -1;
292 
293  double y = p->GetY1() + 0.8 * (p->GetY2() - p->GetY1());
294 
295  std::cout << p->GetY1() << "\t" << p->GetY2() << std::endl;
296 
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());
300 
301  arrow->SetLineColor(line_color);
302  l1 ->SetLineColor(line_color);
303 
304  arrow->SetLineStyle(line_style);
305  l1 ->SetLineStyle(line_style);
306 
307  arrow->SetLineWidth(line_width);
308  l1 ->SetLineWidth(line_width);
309 
310  arrow->Draw();
311  l1->Draw();
312  p->Update();
313 }
const char * p
Definition: xmltok.h:285
OStream cout
Definition: OStream.cxx:6
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:104
void DrawBounds ( TPad *  p,
double  x,
int  line_color = kGray+1,
int  line_style = 1,
int  line_width = 3 
)

Definition at line 268 of file datamc.C.

References MECModelEnuComparisons::g.

Referenced by datamc().

269 {
270  p->cd();
271  p->Update();
272 
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);
279  g->Draw("same");
280 
281  g->Draw();
282  p->Update();
283 }
const char * p
Definition: xmltok.h:285
void load_libs_muonid ( )

Definition at line 32 of file datamc.C.

References exit(), and runNovaSAM::ret.

Referenced by datamc().

33 {
34  const int ret = gSystem->Load("libNDAnamuonid");
35  if(ret != 0) exit(1);
36 }
exit(0)