nc_bkgd_by_interaction_mode.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Spectrum.h"
3 
8 
11 
12 #include "TH1D.h"
13 #include "TLegend.h"
14 #include "TCanvas.h"
15 #include "THStack.h"
16 
17 using namespace ana;
18 void Plot(std::vector<TH1D*> hists, std::vector<std::string> labels, std::string xtitle, std::string name, bool show_integral = false, bool logy = false);
19 
21 {
22  if(infile == "") {
24 
25  const Cut kSelection = kcDQ && kcNHits && kcFrontPlanes && kcNueCCIncContainment && kcMuonIDCut;
27  const HistAxis * vtx_axis = new HistAxis("Vertex in X Direction (cm)",
28  Binning::Simple(50, -200, 200),
29  kVtxX);
30  const HistAxis * cvn_axis = new HistAxis("CVN2017 nueid",
31  Binning::Simple(50, 0, 1),
32  kvCVN2017);
33  std::vector<Spectrum*> specs_mode_cvn;
34  std::vector<Spectrum*> specs_mode_vtx;
35  std::vector<Spectrum*> specs_prong_cvn;
36  std::vector<Spectrum*> specs_prong_vtx;
37  for(auto icut = 2; icut < kcNumIntChns; icut++) {
38  specs_mode_vtx.push_back(new Spectrum(loader,
39  *vtx_axis,
40  kSelection && interaction_chns[icut].cut && kIsNC,
41  kNoShift,
42  kCVWeight));
43  specs_mode_cvn.push_back(new Spectrum(loader,
44  *cvn_axis,
45  kSelection && interaction_chns[icut].cut && kIsNC,
46  kNoShift,
47  kCVWeight));
48  }
49 
50 
51  for(auto iprong = 0; iprong < kcNumChargedParticles; iprong++) {
52  specs_prong_vtx.push_back(new Spectrum(loader,
53  *vtx_axis,
54  kSelection && charged_particles[iprong].cut && kIsNC,
55  kNoShift,
56  kCVWeight));
57  specs_prong_cvn.push_back(new Spectrum(loader,
58  *cvn_axis,
59  kSelection && charged_particles[iprong].cut && kIsNC,
60  kNoShift,
61  kCVWeight));
62  }
63  loader.Go();
64 
65  TFile * output = new TFile("nc_bkgd_by_interaction_mode.root", "recreate");
66  for(auto icut = 2; icut < kcNumIntChns; icut++) {
67  specs_mode_vtx[icut-2]->SaveTo(output, ("vtx_mode_" + interaction_chns[icut].name).c_str());
68  specs_mode_cvn[icut-2]->SaveTo(output, ("cvn_mode_" + interaction_chns[icut].name).c_str());
69  }
70  for(auto iprong = 0; iprong < kcNumChargedParticles; iprong++) {
71  specs_mode_vtx[iprong]->SaveTo(output, ("vtx_prong_" + charged_particles[iprong].name).c_str());
72  specs_mode_cvn[iprong]->SaveTo(output, ("cvn_prong_" + charged_particles[iprong].name).c_str());
73  }
74  output->Close();
75  }
76  else {
77  TFile * input = TFile::Open(infile.c_str(), "read");
78 
79  std::vector<Spectrum*> specs_mode_cvn;
80  std::vector<Spectrum*> specs_mode_vtx;
81  std::vector<Spectrum*> specs_prong_cvn;
82  std::vector<Spectrum*> specs_prong_vtx;
83  for(auto icut = 2; icut < kcNumIntChns; icut++) {
84  specs_mode_vtx.push_back(Spectrum::LoadFrom(input, ("vtx_mode_" + interaction_chns[icut].name).c_str())).release();
85  specs_mode_cvn.push_back(Spectrum::LoadFrom(input, ("cvn_mode_" + interaction_chns[icut].name).c_str())).release();
86  }
87  for(auto iprong = 0; iprong < kcNumChargedParticles; iprong++) {
88  specs_prong_vtx.push_back(Spectrum::LoadFrom(input, ("vtx_prong_" + charged_particles[iprong].name).c_str())).release();
89  specs_prong_cvn.push_back(Spectrum::LoadFrom(input, ("cvn_prong_" + charged_particles[iprong].name).c_str())).release();
90  }
91 
92  int colors[] = {kRed, kMagenta+2, kBlue-4, kCyan, kGreen, kOrange, kGray+2, kYellow+2};
93 
94  std::vector<TH1D*> hists_mode_cvn;
95  std::vector<TH1D*> hists_mode_vtx;
96  std::vector<TH1D*> hists_prong_cvn;
97  std::vector<TH1D*> hists_prong_vtx;
98  for(auto icut = 2; icut < kcNumIntChns; icut++) {
99  hists_mode_cvn.push_back(specs_mode_cvn[icut-2]->ToTH1(kAna2019RHCPOT));
100  hists_mode_vtx.push_back(specs_mode_vtx[icut-2]->ToTH1(kAna2019RHCPOT));
101  hists_mode_cvn.back()->SetLineColor(colors[icut-2]);
102  hists_mode_vtx.back()->SetLineColor(colors[icut-2]);
103  hists_mode_cvn.back()->SetFillColor(colors[icut-2]);
104  hists_mode_vtx.back()->SetFillColor(colors[icut-2]);
105  hists_mode_cvn.back()->SetTitle("Selected NC");
106  hists_mode_vtx.back()->SetTitle("Selected NC");
107  }
108  for(auto iprong = 0; iprong < kcNumChargedParticles; iprong++) {
109  hists_prong_cvn.push_back(specs_prong_cvn[iprong]->ToTH1(kAna2019RHCPOT));
110  hists_prong_vtx.push_back(specs_prong_vtx[iprong]->ToTH1(kAna2019RHCPOT));
111  hists_prong_cvn.back()->SetLineColor(colors[iprong]);
112  hists_prong_vtx.back()->SetLineColor(colors[iprong]);
113  hists_prong_cvn.back()->SetFillColor(colors[iprong]);
114  hists_prong_vtx.back()->SetFillColor(colors[iprong]);
115  hists_prong_cvn.back()->SetTitle("Selected NC");
116  hists_prong_vtx.back()->SetTitle("Selected NC");
117  }
118 
119  std::vector<std::string> chn_labels;
120  std::vector<std::string> png_labels;
121  for(auto icut = 2; icut < kcNumIntChns; icut++) {
122  chn_labels.push_back(interaction_chns[icut].name);
123  }
124  for(auto iprong = 0; iprong < kcNumChargedParticles; iprong++) {
125  png_labels.push_back(charged_particles[iprong].name);
126  }
127 
128  Plot(hists_prong_cvn, png_labels, "CVN nueid" , plot_dump + "/nc_bkgd_png_cvn_axis.pdf", true, true);
129  Plot(hists_prong_vtx, png_labels, "Vertex X (cm)", plot_dump + "/nc_bkgd_png_vtx_axis.pdf", true);
130  Plot(hists_mode_cvn , chn_labels, "CVN nueid" , plot_dump + "/nc_bkgd_chn_cvn_axis.pdf", true, true);
131  Plot(hists_mode_vtx , chn_labels, "Vertex X (cm)", plot_dump + "/nc_bkgd_chn_vtx_axis.pdf", true);
132  }
133 
134 }
135 
136 void Plot(std::vector<TH1D*> hists, std::vector<std::string> labels, std::string xtitle, std::string name, bool show_integral, bool logy)
137 {
138  THStack * stack = new THStack("hs", "");
139  TCanvas * c = new TCanvas();
140  TLegend * leg = new TLegend(0.4, 0.25, "", "NDC");
141  double maxval = 0;
142 
143  TH1D * sum = (TH1D*) hists[0]->Clone();
144  sum->SetLineColor(kBlack);
145  for(auto ihist = 1u; ihist < hists.size(); ihist++) sum->Add(hists[ihist]);
146  if(!logy) sum->GetYaxis()->SetRangeUser(0.00001, sum->GetMaximum() * 1.5);
147  leg->AddEntry(sum, "Total", "l");
148 
149  for(auto ihist = 0u; ihist < labels.size(); ihist++) {
150  maxval = std::max(maxval, hists[ihist]->GetMaximum());
151  }
152  for(auto ihist = 0u; ihist < labels.size(); ihist++) {
153  double integral = hists[ihist]->Integral();
154  stack->Add(hists[ihist]);
155  if(show_integral)
156  leg->AddEntry(hists[ihist], TString::Format("%.2f: %s", integral, labels[ihist].c_str()), "l");
157  else
158  leg->AddEntry(hists[ihist], TString::Format("%s", labels[ihist].c_str()), "l");
159  if(ihist == 0) {
160  hists[ihist]->SetMaximum(maxval);
161  // hists[ihist]->Draw("hist");
162  }
163  else {
164  // hists[ihist]->Draw("hist same");
165  }
166  }
167 
168  sum->Draw("hist");
169  stack->Draw("hist same");
170 // gPad->Modified();
171 // gPad->Update();
172 // stack->GetYaxis()->SetTitle("Events");
173 // stack->GetXaxis()->SetTitle(xtitle.c_str());
174 // stack->GetXaxis()->CenterTitle();
175 // stack->GetYaxis()->CenterTitle();
176 // stack->Draw("hist");
177 
178  leg->SetTextSize(0.03);
179  leg->Draw();
180  if(logy) c->SetLogy();
181  c->Print(name.c_str());
182 }
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
const Cut kcNueCCIncContainment([](const caf::SRProxy *sr){if(sr->vtx.elastic[0].fuzzyk.nshwlid< 1) return false;TVector3 start=sr->vtx.elastic[0].fuzzyk.png[0].shwlid.start;TVector3 stop=sr->vtx.elastic[0].fuzzyk.png[0].shwlid.stop;for(uint ix=0;ix< sr->vtx.elastic[0].fuzzyk.nshwlid;ix++){TVector3 start_muoncatcher=sr->vtx.elastic[0].fuzzyk.png[ix].shwlid.start;TVector3 stop_muoncatcher=sr->vtx.elastic[0].fuzzyk.png[ix].shwlid.stop;if(std::max(start_muoncatcher.Z(), stop_muoncatcher.Z()) > 1250.0) return false;}if(sr->sel.nuecosrej.distallpngtop< 50) return false;if(sr->sel.nuecosrej.distallpngbottom< 30) return false;if(sr->sel.nuecosrej.distallpngeast< 50) return false;if(sr->sel.nuecosrej.distallpngwest< 30) return false;if(sr->sel.nuecosrej.distallpngfront< 150) return false;return true;})
enum BeamMode kOrange
ofstream output
const std::string NominalMC_entire_RHC_prod4
enum BeamMode kRed
const Cut kcFrontPlanes
Definition: NueCCIncCuts.h:117
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Cut kcNHits
Definition: NueCCIncCuts.h:118
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const ana::Var kCVWeight
void nc_bkgd_by_interaction_mode(std::string infile="", std::string plot_dump="")
std::vector< double > Spectrum
Definition: Constants.h:746
const int kcNumIntChns
Definition: NueCCIncCuts.h:398
TString hists[nhists]
Definition: bdt_com.C:3
const Var kVtxX
bool logy
void Plot(std::vector< TH1D * > hists, std::vector< std::string > labels, std::string xtitle, std::string name, bool show_integral=false, bool logy=false)
const Cut kcMuonIDCut
Definition: NueCCIncCuts.h:123
const Var kvCVN2017([](const caf::SRProxy *sr){return(sr->sel.cvn2017.nueid);})
const Cut kcDQ
Definition: NueCCIncCuts.h:114
int colors[6]
Definition: tools.h:1
string infile
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:563
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
virtual void Go() override
Load all the registered spectra.
loader
Definition: demo0.py:10
const std::map< std::string, ana::Cut > interaction_chns
const SystShifts kNoShift
Definition: SystShifts.cxx:22
const double kAna2019RHCPOT
Definition: Exposures.h:224
const Cut cut
Definition: exporter_fd.C:30
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
Double_t sum
Definition: plot.C:31
enum BeamMode kGreen
enum BeamMode kBlue
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
_HistAxis< Var > HistAxis
Definition: HistAxis.h:10
enum BeamMode string