mrbrem_plots.C
Go to the documentation of this file.
1 #include <fstream>
4 #include "CAFAna/Core/Spectrum.h"
7 #include "CAFAna/Core/Ratio.h"
8 #include "CAFAna/Vars/Vars.h"
10 #include "CAFAna/Analysis/Calcs.h"
11 #include "CAFAna/Cuts/NueCutsSecondAna.h"
16 #include "CAFAna/Cuts/SpillCuts.h"
17 #include "CAFAna/Vars/XsecTunes.h"
18 #include "3FlavorAna/Vars/MRVars.h"
20 
22 
23 using namespace ana;
24 
25 #include "OscLib/IOscCalc.h"
26 #include "OscLib/OscCalcDumb.h"
27 #include "OscLib/OscCalcPMNSOpt.h"
28 
29 #include "TCanvas.h"
30 #include "TObject.h"
31 #include "TArrow.h"
32 #include "TColor.h"
33 #include "TFile.h"
34 #include "TLine.h"
35 #include "TH1.h"
36 #include "TText.h"
37 #include "TLegend.h"
38 #include "TLatex.h"
39 #include "TPad.h"
40 #include "TGraphAsymmErrors.h"
41 
42 void Prelim()
43 {
44  TLatex* prelim = new TLatex(.93, .9, "NOvA Preliminary");
45  prelim->SetTextColor(kBlue);
46  prelim->SetNDC();
47  prelim->SetTextSize(2/30.);
48  prelim->SetTextAngle(270);
49  prelim->SetTextAlign(12);
50  prelim->Draw();
51 }
52 
53 TCanvas *MakeMePlot(std::vector<TH1*> Histos, std::vector<TString> Names, bool isFHC, string yname = "none");
54 
56 
57  std::string filename = "MRbrem_spectra_with_weights_"+mode+".root";
58 
59  struct HistDef
60  {
62  HistAxis axis;
63  int lowedge;
64  int highedge;
65  int rebin;
66  };
67 
68  struct CutDef
69  {
71  Cut cut;
72  };
73 
74  const int kNumSels = 10;
75  const CutDef sels[kNumSels] = {
76  {"core_presel", kNue2020CorePresel_MRBrem},
77  {"core_full", kNue2020FD_MRBrem},
78  {"periph_presel", kNue2020PeripheralBasicEventCut_MRBrem},
79  {"periph_full", kNue2020FDPeripheral_MRBrem},
80  {"periph_preselbdt", kNue2020PeripheralBasicEventCut_MRBrem},
81  {"periph_fullbdt", kNue2020FDPeripheral_MRBrem},
82  {"periph_preselcvn", kNue2020PeripheralBasicEventCut_MRBrem},
83  {"periph_fullcvn", kNue2020FDPeripheral_MRBrem},
84  {"periph_preselen", kNue2020PeripheralBasicEventCut_MRBrem},
85  {"periph_fullen", kNue2020FDPeripheral_MRBrem},
86  };
87  const int kNumVars = 14;
88  const HistDef defs[kNumVars] = {
89  {"shwangle", {"Cos(#theta_{Z})", Binning::Simple(60,0,1.0), kShwAngle}, 0, 60, 2},
90  {"cvne", {"CVNe", Binning::Simple(120, 0.0, 1.0), kCVNe_looseptp}, 0, 120, 2},
91  {"corebdt", {"Core BDT", Binning::Simple(120, 0.0, 1.0), kCosPIDCoreBDT}, 35, 90, 1},
92  {"peribdt", {"Peripheral BDT", Binning::Simple(120, 0.0, 1.0), kCosPIDPeriBDT}, 40, 90, 1},
93  {"shwe", {"Shower Energy (GeV)", Binning::Simple(120, 0.0, 5.0), kShwE}, 0, 120, 4},
94  {"shwidth", {"Shower Width (cm)", Binning::Simple(100, 0.0, 10.0), kShwWidth}, 0, 120, 4},
95  {"shwlen", {"Shower Length (cm)", Binning::Simple(100, 0.0, 1000.0), kShwLen}, 0, 100, 2},
96  {"shwstx", {"Shower Start position X (cm)", Binning::Simple(64, -800, 800), kShwStartX}, 0, 64, 2},
97  {"shwsty", {"Shower Start position Y (cm)", Binning::Simple(64, -800, 800), kShwStartY}, 0, 64, 2},
98  {"shwstz", {"Shower Start position Z (cm)", Binning::Simple(104,0,5962), kShwStartZ}, 0, 104, 2},
99  {"nhit", {"NHit", Binning::Simple(120,0.0,150), kNHit}, 0, 120, 4},
100  {"ptp", {"pTp", Binning::Simple(120,0.0,1.0), kPtP}, 0, 120, 4},
101  {"asym", {"SparsenessAsymm", Binning::Simple(60,-1.0,1.0), kSparsenessAsymm}, 0, 60, 1},
102  {"dist", {"DistNearestWall (cm)", Binning::Simple(8,0,800), kDistNearestWall}, 0, 8, 1},
103  };
104 
105  Spectrum* spec_mrbrem_data[kNumSels][kNumVars];
106  Spectrum* spec_mrbrem_cry[kNumSels][kNumVars];
107  Spectrum* spec_genie_fd[kNumSels][kNumVars];
108  Spectrum* spec_mrbrem_data_w_weights[kNumSels][kNumVars];
109  Spectrum* spec_mrbrem_cry_w_weights[kNumSels][kNumVars];
110 
111  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
112  const char* selName = sels[selIdx].name.c_str();
113  for (int varIdx=0; varIdx<kNumVars;varIdx++){
114  const char* varName = defs[varIdx].name.c_str();
115  spec_mrbrem_data[selIdx][varIdx] = LoadFromFile<Spectrum>(filename, TString::Format("%s/mrbrem_data_%s", selName, varName).Data()).release();
116  spec_mrbrem_cry[selIdx][varIdx] = LoadFromFile<Spectrum>(filename, TString::Format("%s/mrbrem_mc_%s", selName, varName).Data()).release();
117  spec_genie_fd[selIdx][varIdx] = LoadFromFile<Spectrum>(filename, TString::Format("%s/genie_fd_%s", selName, varName).Data()).release();
118  spec_mrbrem_data_w_weights[selIdx][varIdx] = LoadFromFile<Spectrum>(filename, TString::Format("%s/mrbrem_data_w_weight_%s", selName, varName).Data()).release();
119  spec_mrbrem_cry_w_weights[selIdx][varIdx] = LoadFromFile<Spectrum>(filename, TString::Format("%s/mrbrem_mc_w_weight_%s", selName, varName).Data()).release();
120  }
121  }
122 
123  double livetime = spec_mrbrem_data[0][0]->Livetime();
124  double pot = spec_genie_fd[0][0]->POT();
125 
126  TH1D* hspec_mrbrem_data[kNumSels][kNumVars];
127  TH1D* hspec_mrbrem_cry[kNumSels][kNumVars];
128  TH1D* hspec_genie_fd[kNumSels][kNumVars];
129  TH1D* hspec_mrbrem_data_w_weights[kNumSels][kNumVars];
130  TH1D* hspec_mrbrem_cry_w_weights[kNumSels][kNumVars];
131  TH1D* hspec_genie_fd_w_weights[kNumSels][kNumVars];
132 
133  double cry_factor_w_weight[kNumSels][kNumVars];
134  double genie_factor_w_weight[kNumSels][kNumVars];
135  double cry_factor[kNumSels][kNumVars];
136  double genie_factor[kNumSels][kNumVars];
137 
138  for (int selIdx = 0; selIdx < kNumSels; ++selIdx){
139  string selName = sels[selIdx].name;
140  for (int varIdx=0; varIdx<kNumVars;varIdx++){
141  string varName = defs[varIdx].name;
142  hspec_mrbrem_data[selIdx][varIdx] = spec_mrbrem_data[selIdx][varIdx]->ToTH1(livetime, kBlack, kSolid, kLivetime);
143  hspec_mrbrem_cry[selIdx][varIdx] = spec_mrbrem_cry[selIdx][varIdx]->ToTH1(livetime, kRed, kSolid, kLivetime);
144  hspec_genie_fd[selIdx][varIdx] = spec_genie_fd[selIdx][varIdx]->ToTH1(pot, kBlue);
145  hspec_mrbrem_data_w_weights[selIdx][varIdx] = spec_mrbrem_data_w_weights[selIdx][varIdx]->ToTH1(livetime, kBlack, kSolid, kLivetime);
146  hspec_mrbrem_cry_w_weights[selIdx][varIdx] = spec_mrbrem_cry_w_weights[selIdx][varIdx]->ToTH1(livetime,kRed, kSolid, kLivetime);
147  hspec_genie_fd_w_weights[selIdx][varIdx] = spec_genie_fd[selIdx][varIdx]->ToTH1(pot, kBlue);
148 
149  cry_factor[selIdx][varIdx] = hspec_mrbrem_data[selIdx][varIdx]->Integral()/hspec_mrbrem_cry[selIdx][varIdx]->Integral();
150  cry_factor_w_weight[selIdx][varIdx] =hspec_mrbrem_data_w_weights[selIdx][varIdx]->Integral()/hspec_mrbrem_cry_w_weights[selIdx][varIdx]->Integral();
151 
152  genie_factor[selIdx][varIdx] = hspec_mrbrem_data[selIdx][varIdx]->Integral()/hspec_genie_fd[selIdx][varIdx]->Integral();
153  genie_factor_w_weight[selIdx][varIdx] = hspec_mrbrem_data_w_weights[selIdx][varIdx]->Integral()/hspec_genie_fd[selIdx][varIdx]->Integral();
154 
155  if (selIdx%2!=0){
156  cry_factor[selIdx][varIdx] = cry_factor[selIdx-1][varIdx];
157  cry_factor_w_weight[selIdx][varIdx] = cry_factor_w_weight[selIdx-1][varIdx];
158  genie_factor[selIdx][varIdx] = genie_factor[selIdx-1][varIdx];
159  genie_factor_w_weight[selIdx][varIdx] = genie_factor_w_weight[selIdx-1][varIdx];
160  }
161 
162  std::cout<<cry_factor[selIdx][varIdx]<<" "<<genie_factor[selIdx][varIdx]<<endl;
163  hspec_mrbrem_cry[selIdx][varIdx]->Scale(cry_factor[selIdx][varIdx]);
164  hspec_mrbrem_cry_w_weights[selIdx][varIdx]->Scale(cry_factor_w_weight[selIdx][varIdx]);
165  hspec_genie_fd[selIdx][varIdx]->Scale(genie_factor[selIdx][varIdx]);
166  hspec_genie_fd_w_weights[selIdx][varIdx]->Scale(genie_factor_w_weight[selIdx][varIdx]);
167 
168  for(auto &h:{hspec_mrbrem_data[selIdx][varIdx], hspec_mrbrem_cry[selIdx][varIdx], hspec_genie_fd[selIdx][varIdx], hspec_mrbrem_data_w_weights[selIdx][varIdx], hspec_mrbrem_cry_w_weights[selIdx][varIdx], hspec_genie_fd_w_weights[selIdx][varIdx]}){
169  if(defs[varIdx].rebin>1) h->Rebin(defs[varIdx].rebin);
170  h->GetXaxis()->SetRange(defs[varIdx].lowedge, defs[varIdx].highedge/defs[varIdx].rebin);
171  }
172 
173  TCanvas* c1 = MakeMePlot({hspec_genie_fd[selIdx][varIdx], hspec_mrbrem_cry[selIdx][varIdx], hspec_mrbrem_data[selIdx][varIdx]}, {"FD MC", "CRY MC", "Data"}, true);
174  c1->SaveAs(("plots/mrbrem-"+mode+"-plot-"+varName+"-"+selName+"-before-reweighting.pdf").c_str());
175  TCanvas* c2 = MakeMePlot({hspec_genie_fd_w_weights[selIdx][varIdx], hspec_mrbrem_cry_w_weights[selIdx][varIdx], hspec_mrbrem_data_w_weights[selIdx][varIdx]}, {"FD MC", "CRY MC", "Data"}, true);
176  c2->SaveAs(("plots/mrbrem-"+mode+"-plot-"+varName+"-"+selName+"-after-reweighting.pdf").c_str());
177 
178  if(selIdx>0 && (selIdx%2!=0)){
179  //std::string additionaltag = "";
180  //if (selIdx ==5) additionaltag = "bdt";
181  TH1D* hEfficData = (TH1D*)hspec_mrbrem_data[selIdx][varIdx]->Clone();
182  hEfficData->Divide(hspec_mrbrem_data[selIdx-1][varIdx]);
183 
184  TH1D* hEfficMC = (TH1D*)hspec_mrbrem_cry[selIdx][varIdx]->Clone();
185  hEfficMC->Divide(hspec_mrbrem_cry[selIdx-1][varIdx]);
186 
187  TCanvas* c3 = MakeMePlot({hEfficMC, hEfficData}, {"CRY MC", "Data"}, true, "Efficiency");
188  c3->SaveAs(("plots/mrbrem-"+mode+"-plot-"+varName+"-"+selName+"-before-reweighting-effic.pdf").c_str());
189 
190  TH1D* hEfficData_w_weights = (TH1D*)hspec_mrbrem_data_w_weights[selIdx][varIdx]->Clone();
191  hEfficData_w_weights->Divide(hspec_mrbrem_data_w_weights[selIdx-1][varIdx]);
192 
193  TH1D* hEfficMC_w_weights = (TH1D*)hspec_mrbrem_cry_w_weights[selIdx][varIdx]->Clone();
194  hEfficMC_w_weights->Divide(hspec_mrbrem_cry_w_weights[selIdx-1][varIdx]);
195 
196  TCanvas* c4 = MakeMePlot({hEfficMC_w_weights, hEfficData_w_weights}, {"CRY MC", "Data"}, true, "Efficiency");
197  c4->SaveAs(("plots/mrbrem-"+mode+"-plot-"+varName+"-"+selName+"-after-reweighting-effic.pdf").c_str());
198  }
199  }
200  }
201 
202 }
203 
204 
205 TCanvas *MakeMePlot(std::vector<TH1*> Histos, std::vector<TString> Names, bool isFHC, string yname){
206 
207  TCanvas *c = new TCanvas("c", "canvas", 800, 700);
208  TPad* p1 = new TPad("", "", 0, 0, 1, 1);
209  TPad* p2 = new TPad("", "", 0, 0, 1, 1);
210  p1->SetBottomMargin(.3);
211  p2->SetTopMargin(.7);
212 
213  for(auto p:{p1,p2}){
214  p->SetLeftMargin(0.1);
215  p->SetFillStyle(0);
216  p->Draw();
217  }
218  p1->cd();
219 
220  auto denom = (TH1*)Histos[0]->Clone();
221  double max = denom->GetMaximum();
222  for (int i=1; i<Histos.size(); i++){
223  max = std::max(Histos[i]->GetMaximum(), max);
224  }
225  denom->SetMaximum(max*1.4);
226  denom->SetMinimum(0);
227  denom->GetYaxis()->CenterTitle();
228  denom->GetXaxis()->CenterTitle();
229  denom->GetXaxis()->SetLabelOffset(999);
230  if(yname !="none") denom->GetYaxis()->SetTitle(yname.c_str());
231 
232  denom->Draw("hist E1");
233  for(int i=1; i<Histos.size(); i++){
234  Histos[i]->Draw("hist E1 same");
235  }
236  denom->Draw("hist E1 same");
237  p1->RedrawAxis();
238  gPad->Update();
239  TText *tt;
240  if(isFHC){
241  auto *tt = new TText(gPad->GetUxmax()*0.7, 1.02*gPad->GetUymax(), "FHC beam");
242  tt->SetTextColor(kBlue);
243  tt->SetTextSize(0.05);
244  tt->Draw();}
245  else{
246  auto *tt = new TText(gPad->GetUxmax()*0.7, 1.02*gPad->GetUymax(), "RHC beam");
247  tt->SetTextColor(kBlue);
248  tt->SetTextSize(0.05);
249  tt->Draw();}
250 
251  TLegend *leg = new TLegend(0.52,0.76,0.78,0.87);
252  leg->SetTextSize(0.025);
253  for(int i=0; i<Histos.size(); i++) leg->AddEntry(Histos[i], Names[i], "l");
254  leg->Draw();
255 
256  gPad->RedrawAxis();
257  p2->cd();
258  auto temp = (TH1*)Histos[1]->Clone();
259  temp->Divide(denom);
260  temp->Draw();
261  for(int i=2; i<Histos.size(); i++){
262  auto temp = (TH1*)Histos[i]->Clone();
263  temp->Divide(denom);
264  temp->Draw("same");
265  }
266  temp->GetYaxis()->SetRangeUser(0.8, 1.25);
267  temp->GetYaxis()->SetTitle("Ratio");
268  p2->RedrawAxis();
269 
270  return c;
271 
272 }
273 
const Var kShwE([](const caf::SRProxy *sr){double maxe=-99.0;if(!sr->vtx.elastic.IsValid) return-99999.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-99999.0;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.npng;i++){if(maxe< sr->vtx.elastic.fuzzyk.png[i].shwlid.shwE){maxe=sr->vtx.elastic.fuzzyk.png[i].shwlid.shwE;}}return maxe;})
Definition: MRVars.h:9
T max(const caf::Proxy< T > &a, T b)
const int kNumVars
Definition: vars.h:14
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Cut kNue2020PeripheralBasicEventCut_MRBrem
Definition: NueCuts2020.h:147
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.
Definition: Spectrum.cxx:148
const Cut kNue2020FD_MRBrem
Definition: NueCuts2020.h:144
const TString name
const Cut kNue2020CorePresel_MRBrem
Definition: NueCuts2020.h:143
const char * p
Definition: xmltok.h:285
std::string name
Definition: NuePlotLists.h:12
const Var kShwStartY
string filename
Definition: shutoffs.py:106
const Var kDistNearestWall([](const caf::SRProxy *sr){float disttonearestwall=std::min(float(kDistAllTop(sr)), float(kDistNotTop(sr)));return disttonearestwall;})
Definition: MRVars.h:10
const Var kSparsenessAsymm
Definition: NueVars.cxx:92
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Var kShwWidth
c2
Definition: demo5.py:33
Definition: type_traits.h:56
const Var kPtP
Transverse momentum fraction in slice.
Definition: NueVars.cxx:90
void Prelim()
Definition: mrbrem_plots.C:42
const Var kShwStartX
const Cut sels[kNumSels]
Definition: vars.h:44
TCanvas * MakeMePlot(std::vector< TH1 * > Histos, std::vector< TString > Names, bool isFHC, string yname="none")
Definition: mrbrem_plots.C:205
const Var kNHit
Definition: Vars.cxx:71
#define pot
Struct to hold cut information.
const HistDef defs[kNumVars]
Definition: vars.h:15
const Cut kNue2020FDPeripheral_MRBrem(kNue2020FDPeripheralFunc_MRBrem)
const int kNumSels
Definition: vars.h:43
TLatex * prelim
Definition: Xsec_final.C:133
double POT() const
Definition: Spectrum.h:227
static bool isFHC
HistAxis axis
Definition: NuePlotLists.h:13
OStream cout
Definition: OStream.cxx:6
void mrbrem_plots(std::string mode="fhc")
Definition: mrbrem_plots.C:55
const Cut cut
Definition: exporter_fd.C:30
double livetime
Definition: saveFDMCHists.C:21
c1
Definition: demo5.py:24
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const Var kShwAngle([](const caf::SRProxy *sr){double maxe=-99.0;double cosz=-1.0;if(!sr->vtx.elastic.IsValid) return-99999.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-99999.0;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.npng;i++){if(maxe< sr->vtx.elastic.fuzzyk.png[i].shwlid.shwE){maxe=sr->vtx.elastic.fuzzyk.png[i].shwlid.shwE;cosz=sr->vtx.elastic.fuzzyk.png[i].shwlid.dir.z;}}return cosz;})
Definition: MRVars.h:8
const Var kShwStartZ
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
const Var kShwLen
enum BeamMode kBlue
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Var kCVNe_looseptp
Definition: Vars.cxx:36
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
const Var kCosPIDCoreBDT
2020 nue cosmic rejection BDT variable - core
Definition: NueVars.cxx:112
const Var kCosPIDPeriBDT
2020 nue cosmic rejection BDT variable - peripheral
Definition: NueVars.cxx:113
enum BeamMode string