Functions
sensitivity_plot.C File Reference
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Analysis/Plots.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Experiment/Dmsq32Constraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/ReactorExperiment.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/FC/FCSurface.h"
#include "3FlavorAna/Prediction/PredictionSystJoint2018.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "CAFAna/Vars/FitVars.h"
#include "Utilities/rootlogon.C"
#include "OscLib/IOscCalc.h"
#include "DDTUtilities/func/MathUtil.h"
#include "3FlavorAna/Ana2020/joint_fit_2020_loader_tools.h"
#include "TCanvas.h"
#include "TMultiGraph.h"
#include "TMarker.h"
#include "TBox.h"
#include "TLatex.h"
#include "TColor.h"
#include "TGraph.h"
#include "TVectorD.h"
#include "TF1.h"
#include "TLegend.h"
#include "TLine.h"
#include "TStyle.h"
#include "TSystem.h"
#include "TSpline.h"
#include "TGraphSmooth.h"
#include "TGraphAsymmErrors.h"
#include "TGaxis.h"
#include <algorithm>
#include <vector>
#include <string>
#include <fstream>
#include <sstream>

Go to the source code of this file.

Functions

bool is_file_exist (const char *fileName)
 
double sigma (TH1F *hist, double percentile)
 
void sensitivity_plot (string hie="IH", string tag="test")
 

Function Documentation

bool is_file_exist ( const char *  fileName)

Definition at line 49 of file sensitivity_plot.C.

References om::cout, allTimeWatchdog::endl, and infile.

Referenced by sensitivity_plot().

50 {
51  std::ifstream infile(fileName);
52  std::cout<<fileName<<std::endl;
53  return infile.good();
54 }
plotROC.fileName
fileName
Definition: plotROC.py:78
infile
string infile
Definition: INukeNucleonCorr.cxx:60
om::cout
OStream cout
Definition: OStream.cxx:6
void sensitivity_plot ( string  hie = "IH",
string  tag = "test" 
)

Definition at line 63 of file sensitivity_plot.C.

References allInOneTrainingPlots::axis, band, plot_validation_datamc::c, demo5::c1, demo5::c2, om::cout, delta, APDHVSetting::dummy, allTimeWatchdog::endl, stan::math::fabs(), shutoffs::filename, make_syst_table_plots::h, MECModelEnuComparisons::i, infile, makeTrainCVSamples::int, is_file_exist(), calib::j, ana::kBlessedLabelFontSize, ana::kBlessedTitleFontSize, kBlue, kViolet, MECModelEnuComparisons::leg, make_syst_table_plots::line, extractScale::mean, extractScale::means, min(), stan::math::minus(), getGoodRuns4SAM::n, num, Preliminary(), sigma(), std::sqrt(), fillBadChanDBTables::step, getGoodRuns4SAM::tag, art::to_string(), ana::UniqueName(), and ana::XAxisDeltaCPLabels().

64 {
65  int n = 62;
66  int k = n-2;
67 
68  double step = (int) (2.0/k * 10000.0 )/10000.0;
69  cout<<"step "<<step<<endl;
70 
71  TGraph *ihres = new TGraph(1);
72  ihres->SetMarkerStyle(20);
73  ihres->SetMarkerSize(1.5);
74  ihres->SetMarkerColor(kPink-2);
75  ihres->SetPoint(0, 1.99, 1.66);
76 
77  TGraph *lores = new TGraph(1);
78  lores->SetMarkerStyle(20);
79  lores->SetMarkerSize(1.5);
80  lores->SetMarkerColor(kPink-2);
81  lores->SetPoint(0, 1.99, 1.17);
82 
83  for (auto &hname:{"ihrej", "lorej"}){
84  auto c = new TCanvas(ana::UniqueName().c_str());
85  c->SetFillStyle(0);
86  c->SetLeftMargin(0.14);
87  c->SetBottomMargin(0.15);
88  double up = 5.55;
89  if (tag == "delta2025NHUO") up = 9.55;
90  string htitle = ";#delta_{CP} / #pi;#sqrt{LL_{IH}-LL_{NH}}";
91  if (hname == "lorej") htitle = ";#delta_{CP} / #pi;#sqrt{LL_{LO}-LL_{UO}}";
92  TH2F* h = new TH2F("", htitle.c_str(), n, -.0166,step*(k+1)+0.0166, 100, 0, up);
93  h->GetXaxis()->CenterTitle();
94  h->GetYaxis()->CenterTitle();
95  h->GetXaxis()->SetTitleSize(kBlessedTitleFontSize);
96  h->GetYaxis()->SetTitleSize(kBlessedTitleFontSize);
97  h->GetXaxis()->SetLabelSize(kBlessedLabelFontSize);
98  h->GetYaxis()->SetLabelSize(kBlessedLabelFontSize);
99  h->GetXaxis()->SetTitleOffset(0.8);
100  h->GetYaxis()->SetTitleOffset(0.8);
101 
102  TGraph *means = new TGraph(n);
103  means->SetMarkerColor(kPink-2);
104  means->SetMarkerStyle(20);
105 
106  TGraph *minsigma = new TGraph(n);
107  minsigma->SetMarkerColor(kPink-2);
108  minsigma->SetMarkerStyle(20);
109  TGraph *min2sigma = new TGraph(n);
110  min2sigma->SetMarkerColor(kPink-2);
111  min2sigma->SetMarkerStyle(20);
112 
113  TGraph *plsigma = new TGraph(n);
114  plsigma->SetMarkerColor(kPink-2);
115  plsigma->SetMarkerStyle(20);
116  TGraph *pl2sigma = new TGraph(n);
117  pl2sigma->SetMarkerColor(kPink-2);
118  pl2sigma->SetMarkerStyle(20);
119 
120  for (int k=0; k<n; k++){
121  int i = k;
122  double delta = step*i;
123  cout<<"delta is "<<delta<<endl;
124  TH1F* tmp_for_delta = new TH1F("tmp", "tmp", 100, 0, up);
125  int num = 5;
126  if (tag == "delta2025NHUO") num = 50;
127  if(k==61) i=0;
128  for (int j = 1; j<=20; j++){
129  TString filename ("/pnfs/nova/scratch/users/lkolupae/sens/prod5/"+tag+"/stat_bestfit_for_joint_realData_both_"+hie+"_deltacp_" + std::to_string(i*step) + "_"+std::to_string(num)+"_" + "thr"+std::to_string(j)+".root");
130  if (is_file_exist(filename)){
131  TFile * infile = new TFile(filename,"read");
132  for (int j =0; j<num; j++){
133  auto ihlo =* (TVectorD*)infile->Get(("universe_"+std::to_string(j)+"/chi_"+std::to_string(j)+"_IH_LO").c_str());
134  auto ihuo =* (TVectorD*)infile->Get(("universe_"+std::to_string(j)+"/chi_"+std::to_string(j)+"_IH_UO").c_str());
135  auto nhlo =* (TVectorD*)infile->Get(("universe_"+std::to_string(j)+"/chi_"+std::to_string(j)+"_NH_LO").c_str());
136  auto nhuo =* (TVectorD*)infile->Get(("universe_"+std::to_string(j)+"/chi_"+std::to_string(j)+"_NH_UO").c_str());
137  double chi_ihlo = ihlo[0];
138  double chi_ihuo = ihuo[0];
139  double chi_nhlo = nhlo[0];
140  double chi_nhuo = nhuo[0];
141  double lobf = min(chi_nhlo, chi_ihlo);
142  double ihbf = min(chi_ihlo, chi_ihuo);
143 
144  double uobf = min(chi_nhuo, chi_ihuo);
145  double nhbf = min(chi_nhlo, chi_nhuo);
146  double ihrej = sqrt(fabs(ihbf-nhbf));
147  double lorej = sqrt(fabs(lobf-uobf));
148  if (hname=="ihrej") {
149  h->Fill(delta, ihrej);
150  tmp_for_delta->Fill(ihrej);
151  }
152  if (hname=="lorej") {
153  h->Fill(delta, lorej);
154  tmp_for_delta->Fill(lorej);}
155  }
156  infile->Close();
157  }
158  else cout<<"no such file"<<endl;
159  }
160  cout<<" for delta n. bin is "<<h->GetXaxis()->FindBin(delta)<<endl;
161  double mean = sigma(tmp_for_delta, 0.5);
162  means->SetPoint(k, delta, mean);
163  double minus = sigma(tmp_for_delta, 0.16);
164  double plus = sigma(tmp_for_delta, 0.84);
165  double minus2 = sigma(tmp_for_delta, 0.024);
166  double plus2 = sigma(tmp_for_delta, 0.976);
167  minsigma->SetPoint(k, delta, minus);
168  plsigma->SetPoint(k, delta, plus);
169  min2sigma->SetPoint(k, delta, minus2);
170  pl2sigma->SetPoint(k, delta, plus2);
171  cout<<"entries in the hist "<<tmp_for_delta->GetEntries()<<" mean "<<mean<<" th1 integral "<<tmp_for_delta->Integral()<<" minus sigma "<<minus<<endl;
172  cout<<"entries in the th2 hist"<<h->GetEntries()<<endl;
173  }
174  means->SetLineColor(kPink-2);
175  minsigma->SetLineColor(kViolet+2);
176  plsigma->SetLineColor(kViolet+2);
177 
178  h->Draw("colz");
179  means->Draw("C same");
180  minsigma->Draw("C same");
181  plsigma->Draw("C same");
182  if (tag == "delta2025NHUO") min2sigma->Draw("C same");
183  if (tag == "delta2025NHUO") pl2sigma->Draw("C same");
184  gPad->Print(("sens_exp_"+tag+"_"+hname+"_"+hie+".pdf").c_str());
185 
186  auto c2 = new TCanvas(ana::UniqueName().c_str());
187  c2->SetFillStyle(0);
188  c2->SetLeftMargin(0.14);
189  c2->SetBottomMargin(0.15);
190  string htitle3 = ";best fit #delta_{CP};#sqrt{LL_{IH}-LL_{NH}}";
191  if (hname == "lorej") htitle3 = ";best fit #delta_{CP};#sqrt{LL_{LO}-LL_{UO}}";
192  TH2F* axis2 = new TH2F("3", htitle3.c_str(), n, 0, 2, 100, 0, up-1);
193  axis2->GetXaxis()->CenterTitle();
194  axis2->GetYaxis()->CenterTitle();
195  axis2->GetXaxis()->SetTitleSize(kBlessedTitleFontSize);
196  axis2->GetYaxis()->SetTitleSize(kBlessedTitleFontSize);
197  axis2->GetXaxis()->SetLabelSize(kBlessedLabelFontSize);
198  axis2->GetYaxis()->SetLabelSize(kBlessedLabelFontSize);
199  axis2->GetXaxis()->SetTitleOffset(0.8);
200  axis2->GetYaxis()->SetTitleOffset(0.8);
201  axis2->Draw();
202  XAxisDeltaCPLabels(axis2);
203  means->Draw("C same");
204  minsigma->Draw("C same");
205  plsigma->Draw("C same");
206  if (tag == "delta2025NHUO") min2sigma->Draw("C same");
207  if (tag == "delta2025NHUO") pl2sigma->Draw("C same");
208  gPad->Print(("sens_beforesmooth_"+tag+"_"+hname+"_"+hie+".pdf").c_str());
209 
210  auto c1 = new TCanvas(ana::UniqueName().c_str());
211  c1->SetFillStyle(0);
212  c1->SetLeftMargin(0.14);
213  c1->SetBottomMargin(0.15);
214  string htitle2 = ";best fit #delta_{CP};#sqrt{LL_{IH}-LL_{NH}}";
215  if (hname == "lorej") htitle2 = ";best fit #delta_{CP};#sqrt{LL_{LO}-LL_{UO}}";
216  TH2F* axis = new TH2F("2", htitle2.c_str(), n, 0, 2, 100, 0, up-1);
217  axis->GetXaxis()->CenterTitle();
218  axis->GetYaxis()->CenterTitle();
219  axis->GetXaxis()->SetTitleSize(kBlessedTitleFontSize);
220  axis->GetYaxis()->SetTitleSize(kBlessedTitleFontSize);
221  axis->GetXaxis()->SetLabelSize(kBlessedLabelFontSize);
222  axis->GetYaxis()->SetLabelSize(kBlessedLabelFontSize);
223  axis->GetXaxis()->SetTitleOffset(0.8);
224  axis->GetYaxis()->SetTitleOffset(0.8);
225  axis->Draw();
226  XAxisDeltaCPLabels(axis);
227  TColor *c_coral = new TColor(100, 255/255.0, 228/255.0, 218/255.0);
228  TColor *c_coral2 = new TColor(102, 255/255.0, 169/255.0, 155/255.0);
229  TColor *c_wave = new TColor(101, 0/255.0, 167/255.0, 200/255.0);
230  Color_t kFillTestColor = ((hie=="NH")? 101:102);
231  Int_t kColor1 = TColor::GetColorTransparent(kFillTestColor, 0.3);
232  Int_t kColor2 = TColor::GetColorTransparent(kBlue-10, 0.6);
233  if (tag == "delta2025NHUO") kColor1 = TColor::GetColorTransparent(kFillTestColor, 0.7);
234  if (hie=="IH") kColor1 = 100;
235 
236  TGraphSmooth *minsigma_smooth = new TGraphSmooth("normal");
237  TGraph * m_sm = minsigma_smooth->SmoothKern(minsigma,"normal",0.2);
238  m_sm->SetLineColor(kColor1);
239  m_sm->DrawClone("l");
240  TGraphSmooth *min2sigma_smooth = new TGraphSmooth("normal");
241  TGraph * m2_sm = min2sigma_smooth->SmoothKern(min2sigma,"normal",0.2);
242  m2_sm->SetLineColor(kColor1);
243  TGraphSmooth *plsigma_smooth = new TGraphSmooth("normal");
244  TGraph * p_sm = plsigma_smooth->SmoothKern(plsigma,"normal",0.2);
245  p_sm->SetLineColor(kColor1);
246  p_sm->DrawClone("l");
247  TGraphSmooth *pl2sigma_smooth = new TGraphSmooth("normal");
248  TGraph * p2_sm = pl2sigma_smooth->SmoothKern(pl2sigma,"normal",0.2);
249  p2_sm->SetLineColor(kColor1);
250  cout<<" N bins after smooth "<<p_sm->GetN()<<" "<<m_sm->GetN()<<endl;
251  double diff_array[p_sm->GetN()];
252  double diff_array2[p2_sm->GetN()];
253  for (int i=0;i<p_sm->GetN();i++) diff_array[i] = p_sm->GetY()[i] - m_sm->GetY()[i];
254  for (int i=0;i<p2_sm->GetN();i++) diff_array2[i] = p2_sm->GetY()[i] - m2_sm->GetY()[i];
255  TGraphAsymmErrors* band = new TGraphAsymmErrors(p_sm->GetN(), p_sm->GetX(), p_sm->GetY(), 0, 0, diff_array, 0 );
256  TGraphAsymmErrors* band2 = new TGraphAsymmErrors(p2_sm->GetN(), p2_sm->GetX(), p2_sm->GetY(), 0, 0, diff_array2, 0 );
257  band->SetFillColor(kColor1);
258  band->Draw("E3 same");
259 
260  TGraphSmooth *mean_smooth = new TGraphSmooth("normal");
261  TGraph * mean_sm = mean_smooth->SmoothKern(means,"normal",0.2);
262  mean_sm->SetLineColor(kBlack);
263  mean_sm->DrawClone("L");
264 
265  double val2018 = 1.57;
266  if(hname == "lorej") val2018 = 1.88;
267 
268  if (hname == "ihrej") ihres->Draw("lp+ same");
269  if (hname == "lorej") lores->Draw("lp+ same");
270 
271  auto line = new TLine(0, val2018, 2, val2018);
272  line->SetLineColor(kColor2);
273  line->SetLineWidth(2);
274  mean_sm->SetLineWidth(2);
275  mean_sm->DrawClone("L same");
276 
277  TLegend * leg = new TLegend(0.16,0.8,0.88,0.88);
278  leg->SetTextSize(kBlessedLabelFontSize*0.75);
279  leg->SetFillStyle(0);
280  leg->SetMargin(1.3*leg->GetMargin());
281  leg->SetNColumns(4);
282  leg->SetMargin(1.4*leg->GetMargin());
283  TGraph * dummy = new TGraph;
284  dummy->SetLineWidth(3);
285  dummy->SetLineColor(kColor2);
286  dummy->SetLineColor(kBlack);
287  leg->AddEntry(dummy->Clone(),"median","l");
288  dummy->SetFillColor(kColor1);
289  leg->AddEntry(dummy->Clone(),"68\% exp","f");
290  if (tag == "delta2025NHUO") {
291  dummy->SetFillColor(kColor2);
292  leg->AddEntry(dummy->Clone(),"95\% exp","f");
293  }
294  leg->AddEntry(ihres,"2019 fit","p");
295  leg->Draw();
296 
297  Preliminary();
298  gPad->RedrawAxis();
299 
300  gPad->Print(("sens_smooth_"+tag+"_"+hname+"_"+hie+".pdf").c_str());
301 
302  }
303 }
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
delta
double delta
Definition: runWimpSim.h:98
make_syst_table_plots.line
line
Definition: make_syst_table_plots.py:268
band
TGraphAsymmErrors * band
Definition: Xsec_final.C:100
ana::XAxisDeltaCPLabels
void XAxisDeltaCPLabels(TH1 *axes, bool t2kunits)
Label the x-axis with fractions of pi.
Definition: NuePlotStyle.cxx:203
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
shutoffs.filename
string filename
Definition: shutoffs.py:106
sigma
double sigma(TH1F *hist, double percentile)
Definition: sensitivity_plot.C:56
demo5.c2
c2
Definition: demo5.py:33
infile
string infile
Definition: INukeNucleonCorr.cxx:60
make_syst_table_plots.h
h
Definition: make_syst_table_plots.py:240
plot_validation_datamc.c
c
Definition: plot_validation_datamc.py:198
calib::j
const double j
Definition: BetheBloch.cxx:29
Preliminary
void Preliminary()
Definition: draw_decomp_plots.C:55
om::cout
OStream cout
Definition: OStream.cxx:6
min
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
is_file_exist
bool is_file_exist(const char *fileName)
Definition: sensitivity_plot.C:49
num
int num
Definition: f2_nu.C:119
kViolet
enum BeamMode kViolet
demo5.c1
c1
Definition: demo5.py:24
ana::kBlessedLabelFontSize
const Float_t kBlessedLabelFontSize
Definition: Style.h:90
art::to_string
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
ana::kBlessedTitleFontSize
const Float_t kBlessedTitleFontSize
Definition: Style.h:89
kBlue
enum BeamMode kBlue
ana::UniqueName
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
stan::math::minus
T minus(const T &x)
Definition: minus.hpp:15
double sigma ( TH1F *  hist,
double  percentile 
)

Definition at line 56 of file sensitivity_plot.C.

Referenced by sensitivity_plot().

56  {
57  double xq[1] = {percentile};
58  double yq[1];
59  hist->GetQuantiles(1, yq, xq);
60  return yq[0];
61 }
Generated on Thu Jan 21 2021 14:16:19 for NOvA by   doxygen 1.8.11