generate_fc_binlists.C
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <string>
3 #include <iostream>
4 #include <fstream>
5 
6 #include "TH1.h"
7 #include "TH2.h"
8 #include "TF1.h"
9 #include "TGraph.h"
10 #include "TFile.h"
11 #include "TMath.h"
12 
13 // Pulls our results and gets a dchisq for every analysis bin.
14 // Lets you only as many points as you need to get a reasonable estimate for
15 // the critical value given this chisq.
16 // We'll then submit an FC job for each bin of the slice / surf.
18 {
19  std::string binlistdir = std::getenv("FCHELPERANA2018_LIB_PATH");
20  binlistdir += "/binlists/";
21  std::ofstream binlist;
22  std::ofstream summaryFile;
23  summaryFile.open(binlistdir+"binlists_summary.txt");
24 
25  int bincount = 0;
26  int totNbins = 0;
27  int slicesNPts = 0;
28  std::string slicedir = std::string(getenv("FCHELPERANA2018_LIB_PATH"));
29  slicedir += "/slices/";
30 
31  for (std::string plot : {"delta", "th23", "dmsq"}){
32  for(std::string oct : {"LO", "UO"}){
33  for(std::string hier : {"IH", "NH"}){
34  // binlist.open(plot+"_binlist.txt");
35  std::string slicename = "hist_slices_2017_joint_realData_bothcombo_systs";
36  if(plot == "delta") slicename += "";
37  if(plot == "th23") slicename += "_th23_noOct";
38  if(plot == "dmsq") slicename += "_dmsq";
39  slicename += ".root";
40 
41  TFile *in1=TFile::Open((slicedir+slicename).c_str());
42  TH1 *slice;
43  if(plot == "th23" && oct == "LO") oct = "";
44  if(plot == "th23" && oct == "UO") continue;
45 
46  std::string sname = "slice_" + plot + "_" + hier + oct;
47 
48  slice = (TH1F*)in1->Get((sname).c_str());
49  totNbins = slice->GetNbinsX();
50  binlist.open(binlistdir+plot+"_"+oct+hier+"_binlist.txt");
51  bincount = 0;
52  double dchisqcut = 20;
53  for (int i = 1; i <= slice->GetNbinsX(); i++){
54 
55  // need delta chi^2, not sqrt(delta chi^2)
56  double gaus = slice->GetBinContent(i); // Gaussian significance
57  gaus = gaus*gaus;
58  if (gaus > dchisqcut )
59  continue; // Don't even try for these plots at high sigma
60 
61  binlist << std::to_string(i-1) << endl;
62  bincount++;
63  slicesNPts++;
64  }
65  binlist.close();
66  std::string hierstr = (hier == 1) ? "_NH: " : "_IH: ";
67  std::string summ = plot + hier+oct+": " + std::to_string(bincount)
68  +"/"+std::to_string(totNbins)
69  + " bins below delta chisq of " + std::to_string(dchisqcut);
70  std::cout << summ << endl;
71  summaryFile << summ << endl;
72  delete slice;
73  in1->Close();
74  }
75  }
76  }
77  std::cout << "Total NPts for the slice plots: " << slicesNPts << std::endl;
78 
79 
80  // And now do the contour
81 
82  // Harder to write down closed form expression for 2D critical values
83  // Calculate the up value, round up to nearest
84 
85  int surfsNPts = 0;
86  bincount = 0;
87  totNbins = 0;
88  std::string dir = std::string(getenv("FCHELPERANA2018_LIB_PATH"));
89  dir += "/contours/";
90 
91  for (std::string plot : {"deltassth23", "ssth23dmsq32"}){
92  for(std::string hier : {"IH", "NH"}){
93  std::string fname = "/hist_contours_2018_joint_realData_both_only" + hier+"combo_" + plot+"_systs.root";
94  TFile *in2 = TFile::Open((dir+fname).c_str());
95 
96 
97  std::string shortname;
98  TH2 *surf;
99  if(plot == "ssth23dmsq32") {
100  shortname = "th23_dm32_";
101  }
102  else if (plot == "deltassth23") {
103  shortname = "delta_th23_";
104  }
105  else {
106  std::cout << "plot " << plot << " not found " << std::endl;
107  break;
108  }
109  surf = (TH2F*)in2->Get((shortname+hier).c_str());
110 
111  bincount = 0;
112  binlist.open(binlistdir+plot+"_"+hier+"_binlist.txt");
113  int nx = surf->GetNbinsX();
114  int ny = surf->GetNbinsY();
115  totNbins = nx*ny;
116 
117  double dchisqcut = 12;
118  cout << plot+" "+hier+": "+ std::to_string(nx) + " x bins\n";
119  cout << plot+" "+hier+": "+ std::to_string(ny) + " y bins\n";
120  for (int j = 1; j <= surf->GetNbinsY(); j++){
121  for (int i = 1; i <= surf->GetNbinsX(); i++){
122  // Surfaces are saved as up value so no need to square
123 
124  double dchisq = surf->GetBinContent(i,j); // Gaussian sig
125  // 3 sig 2D up value is 11.83, if dchisq is much higher, don't even
126  // bother running FC, becomes increasingly hard to believe anyway
127  if (dchisq > dchisqcut) continue;
128 
129  int bin = (i-1) + nx*(j-1); // FC bin to pass to FC macro
130  binlist << std::to_string(bin) << endl;
131  bincount++;
132  surfsNPts++;
133  }
134  }
135 
136  binlist.close();
137  std::string summ = plot + "_"+hier+": " + std::to_string(bincount)+"/"
138  + std::to_string(totNbins)
139  + " bins below dchisq = "
140  + std::to_string(dchisqcut);
141  std::cout << summ << endl;
142  summaryFile << summ << endl;
143  in2->Close();
144  }
145  }
146 
147  summaryFile << endl;
148  summaryFile << "Number of surf bins FC'd: " << surfsNPts << std::endl;
149  summaryFile << "Number of slice bins FC'd: " << slicesNPts;
150  summaryFile.close();
151 }
std::string getenv(std::string const &name)
const double j
Definition: BetheBloch.cxx:29
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
void generate_fc_binlists()
TDirectory * dir
Definition: macro.C:5
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void plot(std::string label, std::map< std::string, std::map< std::string, Spectrum * >> &plots, bool log)
surf
Definition: demo5.py:35
enum BeamMode string