joint_fit_2017_fc_submit.C
Go to the documentation of this file.
1 #include <stdio.h>
2 #include <string>
3 #include <iostream>
4 
5 #include "TH1.h"
6 #include "TH2.h"
7 #include "TF1.h"
8 #include "TGraph.h"
9 #include "TFile.h"
10 #include "TMath.h"
11 
12 // Pulls our results and gets a dchisq for every analysis bin.
13 // Lets you only as many points as you need to get a reasonable estimate for
14 // the critical value given this chisq.
15 // We'll then submit an FC job for each bin of the slice / surf.
17 {
18  double sliceN = 0;
19  double totN = 0;
20  for (std::string plot : {"delta", "ssth23", "dmsq32"}){
21 
22  TFile *in1=TFile::Open(("slices/hist_slices_2017_"+plot+".root").c_str());
23  TH1 *slice_IH;
24  TH1 *slice_NH;
25  if (plot=="delta"){
26  slice_IH = (TH1F*)in1->Get("slice_delta_IH");
27  slice_NH = (TH1F*)in1->Get("slice_delta_IH");
28  }
29  else if (plot=="ssth23"){
30  slice_IH = (TH1F*)in1->Get("slice_th23_IH");
31  slice_NH = (TH1F*)in1->Get("slice_th23_IH");
32  }
33  else if (plot=="dmsq32"){
34  slice_IH = (TH1F*)in1->Get("slice_dmsq_IH");
35  slice_NH = (TH1F*)in1->Get("slice_dmsq_IH");
36  }
37  else{
38  std::cout << "plot " << plot << " not found " << std::endl;
39  break;
40  }
41  std::vector<TH1*> slices = {slice_IH,slice_NH};
42 
43  for (int hier = 0; hier <= 1; hier++){
44  for (int i = 1; i <= slice_IH->GetNbinsX(); i++){
45 
46  double gaus = slices[hier]->GetBinContent(i); // Gaussian significance
47  if (gaus > 4 && (plot=="ssth23"||plot=="dmsq32"))
48  continue; // Don't even try for these plots at high sigma
49  double pval = TMath::Erfc(gaus/sqrt(2)); // Transform to pvalue
50  int N = 10 * (1/pval); // Number of trials to get 4 above gaus
51  N /= 100; // But we do 100 points / grid job
52  N = std::max(19,N);
53 
54  std::string cmd = "sh fc_submit_slice.sh ";
55  cmd += std::to_string(i-1)+" ";
56  cmd += std::to_string(hier)+" ";
57  cmd += std::to_string(N+2) + " ";
58  cmd += plot;
59  std::cout << cmd << std::endl;
60  system(cmd.c_str());
61 
62  totN += N;
63  sliceN += N;
64  }
65  }
66  delete slice_IH; delete slice_NH;
67  in1->Close();
68  }
69  std::cout << "For the slice plots " << 100*sliceN << std::endl;
70 
71 
72  // And now do the contour
73 
74  // Harder to write down closed form expression for 2D critical values
75  // Calculate the up value, round up to nearest
76  double surfN = 0;
77  int surfNPts = 0;
78  for (std::string plot : {"deltassth23", "ssth23dmsq32"}){
79  TFile *in2 = TFile::Open("contours4/hist_contours_2017_jointrealDatacombo_dmsq_systs.root");
80  TH2 *surf_IH = (TH2F*)in2->Get("th23_dm32_IH");
81  TH2 *surf_NH = (TH2F*)in2->Get("th23_dm32_NH");
82 
83  std::vector<TH2*> surfs = {surf_IH,surf_NH};
84  for (int hier = 1; hier <= 1; hier++){
85  //for (int hier = 0; hier <= 1; hier++){
86  for (int j = 1; j <= surfs[hier]->GetNbinsY(); j++){
87  for (int i = 1; i <= surfs[hier]->GetNbinsX(); i++){
88 
89  double dchisq = surfs[hier]->GetBinContent(i,j); // Gaussian sig
90  // 3 sig 2D up value is 11.83, if dchisq is much higher, don't even
91  // bother running FC, becomes increasingly hard to believe anyway
92  if (dchisq > 14) continue;
93  double pval = exp(-dchisq/2); // Analytic for for 2D
94  int N = 10 * (1/pval); // Number of trials to get 10 above gaus
95  N /= 100; // But we do 100 points / grid job
96  N = std::max(5,N); // Throw at least 400 expts, enough for 3 sigma
97 
98 
99  std::string cmd = "sh fc_submit_surf.sh ";
100  int bin = (i-1) + 20*(j-1); // FC bin to pass to FC macro
101  cmd += std::to_string(bin)+" ";
102  cmd += std::to_string(hier)+" ";
103  cmd += std::to_string(N) + " ";
104  cmd += plot;
105  std::cout << cmd << std::endl;
106  system(cmd.c_str());
107 
108  totN += N;
109  surfN += N;
110  surfNPts++;
111  }
112  }
113  }
114  in2->Close();
115  }
116 
117  std::cout << "Total FCPoints needed " << 100*totN << std::endl;
118  std::cout << "For the slice plots " << 100*sliceN << std::endl;
119  std::cout << "For the surf plots " << 100*surfN << std::endl;
120  std::cout << "Number of surf points FC'd: " << surfNPts << std::endl;
121 
122 }
T max(const caf::Proxy< T > &a, T b)
system("rm -rf microbeam.root")
T sqrt(T number)
Definition: d0nt_math.hpp:156
string cmd
Definition: run_hadd.py:52
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
const double j
Definition: BetheBloch.cxx:29
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void joint_fit_2017_fc_submit()
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)