datamc.C
Go to the documentation of this file.
1 #include "CAFAna/Core/Spectrum.h"
3 
4 // analysis cuts/vars
9 
11 #include "CAFAna/Vars/XsecTunes.h"
13 
14 #include "TFile.h"
15 #include "THStack.h"
16 #include "TH1D.h"
17 #include "TCanvas.h"
18 #include "TLegend.h"
19 #include "TLegendEntry.h"
20 #include "TPad.h"
21 #include "TGraph.h"
22 #include "TArrow.h"
23 
24 
25 // temporary measure to use prod4 trained muonid
26 #include "NDAna/muonid/NDXSecMuonPID.h"
27 namespace nuebarccinc {
28  const ana::Cut kMuonIDProd4Cut = ana::kMuonIDProd4 < -0.2;
29  void One();
30 }
31 
33 {
34  const int ret = gSystem->Load("libNDAnamuonid");
35  if(ret != 0) exit(1);
36 }
37 
38 
39 
40 
41 using namespace ana;
42 
43 struct Bounds_t {
44  double min;
45  double max = 0;
46 };
47 
48 void DrawArrow(TPad * p, double x, TString arrow_opt, int line_color = kGray+1, int line_style = 1, int line_width = 3);
49 
50 void DrawBounds(TPad * p, double x, int line_color = kGray+1, int line_style = 1, int line_width = 3);
51 
52 void datamc(std::string input_file_name = "",
53  std::string plot_dump = "./plots/")
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 }
257 
259 {
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);
265  one->Draw("same");
266 }
267 
268 void DrawBounds(TPad * p, double x, int line_color, int line_style, int line_width)
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 }
284 
285 void DrawArrow(TPad * p, double x, TString arrow_opt, int line_color, int line_style, int line_width)
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 }
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
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
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
void datamc(std::string input_file_name="", std::string plot_dump="./plots/")
Definition: datamc.C:52
const char * p
Definition: xmltok.h:285
TH1 * ratio(TH1 *h1, TH1 *h2)
const ana::Cut kFrontPlanes
const ana::Cut kNHits
double min
Definition: datamc.C:44
Float_t tmp
Definition: plot.C:36
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:570
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
const TVector3 vNueCCIncFiducialMax(150, 140, 800)
void load_libs_muonid()
Definition: datamc.C:32
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:104
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
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
exit(0)
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
auto one()
Definition: PMNS.cxx:49
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
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:106
void One()
Definition: datamc.C:258