make_ratios.C
Go to the documentation of this file.
1 #include "TH1.h"
2 #include "TString.h"
3 #include "TFile.h"
4 
5 #include <string>
6 #include <iostream>
7 
8 void make_ratios(TString opt, TString filename){
9 
10  TFile* f = new TFile(filename, "READ");
11 
12  std::vector<std::string> dets = {"NearDet", "FarDet"};
13  std::vector<std::string> beams = {"FHC", "RHC"};
14  std::vector<std::string> sels =
15  {"NuMuSel",
16  "NuESel_LowPID",
17  "NuESel_HighPID",
18  "NuESel_Peripheral",
19  "NCSel"
20  };
21 
22  std::vector<std::string> systs;
23  std::string outname;
24  // for calib systs, we have 5x the apparent number of events
25  // for the nominal sample, becasue of the way the grid scripts
26  // work, so want to correct for that
27  float scale=1;
28  if (opt.Contains("eshift", TString::kIgnoreCase)){
29  systs =
30  {
31  "CorrMuEScaleSyst2020",
32  "NeutronEvisPrimariesSyst2018",
33  "PileupMuESyst2020",
34  "UnCorrMuCatMuESyst2020",
35  "UnCorrNDMuEScaleSyst2020",
36  "michel_tagging2020"
37  };
38  outname="eshiftsystratios.root";
39  }
40  else if (opt.Contains("calib", TString::kIgnoreCase)){
41  systs = {
42  "ll",
43  "calib",
44  "calibshape",
45  "cherenkov",
46  "detectoraging"
47  };
48  outname="calibsystratios.root";
49  scale=1./systs.size();
50  }
51 
52  TFile* outFile = new TFile(outname.c_str(), "RECREATE");
53  outFile->cd();
54 
55  for (auto& det : dets){
56 
57  for (auto& beam : beams){
58 
59  for (auto& sel : sels){
60 
61  // no RHC NCSel
62  if (sel == "NCSel" && beam == "RHC") continue;
63 
64  // no peripheral in ND
65  if (sel == "NuESel_Peripheral" && det == "NearDet") continue;
66 
67  std::string nomName = "Total" + beam + det + sel;
68  std::cout << nomName << std::endl;
69  TH1D* hNom = (TH1D*)f->Get((nomName+"Nominal").c_str());
70  hNom->Scale(scale);
71 
72  for (auto& syst : systs){
73 
74  std::string histName = syst + nomName;
75  std::string upName = histName + "Plus1Sigma";
76  std::string dnName = histName + "Minus1Sigma";
77 
78  std::cout << upName << std::endl;
79 
80  TH1D* hUp = (TH1D*)f->Get(upName.c_str());
81  TH1D* hDn = (TH1D*)f->Get(dnName.c_str());
82 
83  hUp->Sumw2();
84  hDn->Sumw2();
85 
86  hUp->Divide(hNom);
87  hDn->Divide(hNom);
88 
89  hUp->Write(upName.c_str());
90  hDn->Write(dnName.c_str());
91 
92  }
93 
94  }
95 
96  }
97 
98  }
99 
100 }
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void make_ratios(TString opt, TString filename)
Definition: make_ratios.C:8
string filename
Definition: shutoffs.py:106
Double_t scale
Definition: plot.C:25
TFile * outFile
Definition: PlotXSec.C:135
const Cut sels[kNumSels]
Definition: vars.h:44
OStream cout
Definition: OStream.cxx:6
enum BeamMode string