makeCalibSystRatioPlots.C
Go to the documentation of this file.
1 /**
2  * \brief This code creates ratio plots of calibration samples
3  *
4  * This code takes the output of CMFCalibSystHistMaker_module.cc
5  * and creates ratio plots using the spectra pulled from those
6  *
7  * \author Adam Lister
8  *
9  * \date 2020/01/09
10  *
11  * Contact: alister1@fnal.gov
12  */
13 
14 enum HistType {
18 };
19 
23  std::vector<std::string> const fileshort;
24 
26  : name(std::string())
27  , latexname(std::string())
28  , fileshort({})
29  {}
30 
32  std::string l,
33  std::vector<std::string> const& fs = std::vector<std::string>())
34  : name(n)
35  , latexname(l)
36  , fileshort(fs)
37  {}
38 
39 };
40 
41 void setHistogramText(TH1D* h){
42  gStyle->SetTitleX(0.5);
43  gStyle->SetTitleAlign(23);
44  h->GetXaxis()->SetTitleFont(43);
45  h->GetYaxis()->SetTitleFont(43);
46  h->GetXaxis()->SetTitleSize(16);
47  h->GetYaxis()->SetTitleSize(16);
48  h->GetXaxis()->SetLabelFont(43);
49  h->GetYaxis()->SetLabelFont(43);
50  h->GetXaxis()->SetLabelSize(16);
51  h->GetYaxis()->SetLabelSize(16);
52  h->GetYaxis()->SetTitleOffset(1.4);
53  h->GetXaxis()->SetTitleOffset(3.0);
54  h->GetXaxis()->CenterTitle();
55  h->GetYaxis()->CenterTitle();
56 }
57 
58 void styleHistogram(TH1D* h, HistType ht){
60 
61  h->SetLineWidth(2);
62  if (ht == kNominal){
63  h->SetLineColor(kBlack);
64  }
65  if (ht == kPlus1Sigma){
66  h->SetLineColor(kGreen+1);
67  h->SetLineStyle(7);
68  }
69  if (ht == kMinus1Sigma){
70  h->SetLineColor(kAzure+1);
71  h->SetLineStyle(7);
72  }
73 
74 }
75 
77 
78  // SETUP
79 
80  std::string location = "./"; // location of syst files
81  std::string analysisTag = "2020";
82 
83  std::vector<CalibrationType> calibrationTypes = {
84  CalibrationType("lightmodel" , "Light Level" , std::vector<std::string>({"lldown" , "llup" })) ,
85  CalibrationType("calibration" , "Calibration" , std::vector<std::string>({"calibup" , "calibdown"})) ,
86  CalibrationType("cherenkov" , "Cherenkov" , std::vector<std::string>({"cherenkov" })) ,
87  CalibrationType("calibshape" , "Calib. Shape", std::vector<std::string>({"calibshape" })) ,
88  CalibrationType("detectoraging" , "Det. Aging" , std::vector<std::string>({"detectoraging" }))
89  } ;
90 
91  std::vector< std::string > detectors = {"Near", "Far"};
92 
93  std::vector< std::string > beams = {"FHC", "RHC"};
94 
95  std::vector< std::string > selections = {
96  "NuESel_LowPID",
97  "NuESel_HighPID",
98  "NuESel_Peripheral",
99  "NuMuSelQ1",
100  "NuMuSelQ2",
101  "NuMuSelQ3",
102  "NuMuSelQ4",
103  "NCSel"
104  };
105 
106  std::vector<std::string> ints = {
107  "Total",
108  //"NuMuCC",
109  //"NuMuBarCC",
110  //"NuECC",
111  //"NuEBarCC",
112  //"NuTauCC",
113  //"NuTauBarCC",
114  //"NC",
115  };
116 
117  // make output file
118  std::string outFileName = "CalibrationSystematics"+analysisTag+".root";
119  TFile* outputFile = new TFile(outFileName.c_str(), "recreate");
120 
121  // get nominal file
122  TFile* nominalSpectraFile =
123  new TFile((location+"syst_shift_nominal.root").c_str(), "read");
124 
125  for (int icalib = 0; icalib < calibrationTypes.size(); icalib++){
126 
127  // get systematically shifted files
128  TFile *upSystFile = new TFile((location+"syst_shift_"+calibrationTypes.at(icalib).fileshort.at(0)+".root").c_str(), "read");
129  TFile *dnSystFile = nullptr;
130  if (calibrationTypes.at(icalib).fileshort.size() == 2)
131  dnSystFile = new TFile((location+"syst_shift_"+calibrationTypes.at(icalib).fileshort.at(1)+".root").c_str(), "read");
132 
133  if (upSystFile->IsZombie()) {
134  std::cout
135  << "File "
136  << upSystFile->GetName()
137  << " is zombie. Does it exist?"
138  << std::endl;
139  continue;
140  }
141 
142  for (int idet = 0; idet < detectors.size(); idet++){
143  for (int ibeam = 0; ibeam < beams.size(); ibeam++){
144  for (int isel = 0; isel < selections.size(); isel++){
145  for (int iint = 0; iint < ints.size(); iint++){
146 
147  // there is no peripheral sample in the near detector
148  if (selections.at(isel) == "NuESel_Peripheral" && detectors.at(idet) == "Near")
149  continue;
150 
151  std::string baseHistName = beams.at(ibeam) +
152  detectors.at(idet) +
153  selections.at(isel);
154 
155  std::string histName = ints.at(iint) + baseHistName;
156 
157  TH1D* nomHist = (TH1D*)nominalSpectraFile->Get(("hists/Ana2020/"+histName).c_str());
158  TH1D* upHist = (TH1D*)upSystFile ->Get(("hists/Ana2020/"+histName).c_str());
159 
160  if (nomHist->Integral() == 0) continue;
161  std::string plotTitle = detectors.at(idet) + " " +
162  beams.at(ibeam) + " " +
163  selections.at(isel) + ", " +
164  calibrationTypes.at(icalib).latexname + ", " +
165  ints.at(iint);
166 
167  nomHist->SetTitle(plotTitle.c_str());
168 
169  styleHistogram(nomHist, kNominal);
170  styleHistogram(upHist , kPlus1Sigma);
171 
172  TH1D* upHistR = (TH1D*)upHist->Clone("upHistR");
173  upHistR->Divide(nomHist);
174  upHistR->SetTitle("");
175  upHistR->GetYaxis()->SetTitle("#frac{Variation}{Nominal}");
176 
177  // setup histograms for plotting
178  TCanvas *c1 = new TCanvas("c1", "c1", 500, 500);
179  TPad *topPad = new TPad("topPad", "", 0.005, 0.3, 0.990, 0.995);
180  TPad *bottomPad = new TPad("bottomPad", "", 0.005, 0.005, 0.995, 0.3);
181 
182  topPad ->SetTopMargin(0.1);
183  topPad ->SetBottomMargin(0.005);
184  bottomPad ->SetTopMargin(0.005);
185  bottomPad ->SetBottomMargin(0.25);
186  topPad ->SetLeftMargin(0.1);
187  bottomPad ->SetLeftMargin(0.1);
188  bottomPad ->SetGridy();
189  topPad ->Draw();
190  bottomPad ->Draw();
191 
192  nomHist->GetYaxis()->SetRangeUser(1e-3, nomHist->GetBinContent(nomHist->GetMaximumBin())*1.25);
193  std::string ytitle = std::string(nomHist->GetYaxis()->GetTitle());
194  if (ytitle.find("GeV") == std::string::npos){
195  if (selections.at(isel).find("NuMu") != std::string::npos){
196  ytitle += " / 0.1 GeV";
197  nomHist->GetYaxis()->SetTitle(ytitle.c_str());
198  }
199  if (selections.at(isel).find("NuE") != std::string::npos){
200  ytitle += " / 0.5 GeV";
201  nomHist->GetYaxis()->SetTitle(ytitle.c_str());
202  }
203  if (selections.at(isel).find("NC") != std::string::npos){
204  ytitle += " / 1.0 GeV";
205  nomHist->GetYaxis()->SetTitle(ytitle.c_str());
206  }
207  }
208  upHistR->GetYaxis()->SetRangeUser(0.7, 1.3);
209  upHistR->GetYaxis()->SetNdivisions(505);
210 
211  topPad ->cd();
212  nomHist ->Draw("hist");
213  upHist ->Draw("hist same");
214  bottomPad ->cd();
215  upHistR ->Draw("hist");
216 
217  outputFile->cd();
218  std::string upRatioName = calibrationTypes.at(icalib).name + histName + "Plus1Sigma";
219  upHistR->SetName(upRatioName.c_str());
220  upHistR->Write();
221 
222  if (dnSystFile){
223  TH1D* dnHist = (TH1D*)dnSystFile->Get(("hists/Ana2020/"+histName).c_str());
224  styleHistogram(dnHist, kMinus1Sigma);
225 
226  TH1D* dnHistR = (TH1D*)dnHist->Clone("dnHistR");
227  dnHistR->Divide(nomHist);
228  dnHistR->SetName((histName+"Ratio").c_str());
229 
230  topPad ->cd();
231 
232  dnHist ->Draw("hist same");
233  bottomPad ->cd();
234  dnHistR ->Draw("hist same");
235 
236  std::string dnRatioName = calibrationTypes.at(icalib).name + histName + "Minus1Sigma";
237  dnHistR->SetName(dnRatioName.c_str());
238  dnHistR->Write();
239  }
240 
241  c1->SetName((baseHistName + calibrationTypes.at(icalib).fileshort.at(0)).c_str());
242  c1->SaveAs((baseHistName + calibrationTypes.at(icalib).name + ints.at(iint) +".png").c_str());
243  c1->SaveAs((baseHistName + calibrationTypes.at(icalib).name + ints.at(iint) +".pdf").c_str());
244  }
245  }
246  }
247  }
248  }
249 }
#define location
void makeCalibSystRatioPlots()
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
TPad * topPad
Definition: PlotSpectra.h:74
std::string const name
std::string const latexname
std::vector< std::string > const fileshort
void styleHistogram(TH1D *h, HistType ht)
OStream cout
Definition: OStream.cxx:6
c1
Definition: demo5.py:24
CalibrationType(std::string n, std::string l, std::vector< std::string > const &fs=std::vector< std::string >())
HistType
This code creates ratio plots of calibration samples.
dictionary ints
Float_t e
Definition: plot.C:35
enum BeamMode kGreen
void setHistogramText(TH1D *h)
enum BeamMode string