compare_fits.C
Go to the documentation of this file.
6 #include "CAFAna/Core/Spectrum.h"
7 #include "CAFAna/FC/FCSurface.h"
12 #include "CAFAna/Core/Spectrum.h"
17 #include "CAFAna/Systs/Systs.h"
19 #include "CAFAna/Vars/FitVars.h"
21 #include "OscLib/IOscCalc.h"
22 #include "Utilities/rootlogon.C"
23 
24 #include "TFile.h"
25 #include "TStyle.h"
26 #include "TGraph.h"
27 #include "TCanvas.h"
28 #include "TH2.h"
29 
30 using namespace ana;
31 void compare_fits(bool isTH23dcp = false) {
32 
33  struct fitcase{
34  string filename;
35  string nickname;
36  int color;
37  int style;
38  };
39 
40  std::string TagName = "compare_joint_";
41 
42  std::string BaseLoc = "/nova/app/users/karlwarb/Workspace/Analysis2020/Sensitivity/fit_results_Both_";
43  std::string StName = "hist_contours_numuOnly_fake2018_both_";
44  std::string EndName = "_dmsq_stats_combo_looseptp.root";
45  std::string sDipSep = "OptDip_SepQuant/" + StName + "OptDip_SepQuant" + EndName;
46  std::string sDipComb = "OptDip_CombQuant/" + StName + "OptDip_CombQuant" + EndName;
47  std::string sFullSep = "OptFull_SepQuant/" + StName + "OptFull_SepQuant" + EndName;
48  std::string sFullComb = "OptFull_CombQuant/"+ StName + "OptFull_CombQuant" + EndName;
49 
50  int Col1 = kAzure+7;
51  int Col2 = kMagenta-3;
52  int Col3 = kRed-10;
53  int Col4 = kGreen+2;
54 
55  std::vector <fitcase> cases;
56  cases.push_back( {BaseLoc+sDipSep , "OptDip_SepQuant" , Col1, 1 } );
57  cases.push_back( {BaseLoc+sDipComb , "OptDip_CombQuant" , Col2, 10} );
58  cases.push_back( {BaseLoc+sFullSep , "OptFull_SepQuant" , Col3, 3 } );
59  cases.push_back( {BaseLoc+sFullComb, "OptFull_CombQuant", Col4, 4 } );
60 
61  std::string sForceBad = "ForceBad/" + StName + "ForceBad" + EndName;
62  int Col5 = kBlack;
63  cases.push_back( {BaseLoc+sForceBad, "ForceBad", Col5, 1 } );
64 
65  for ( int hie:{+1, -1} ) {
66  std::string hietag= (hie>0?"NH":"IH");
67 
68  if( isTH23dcp ) {
69  //TFile* infile = new TFile (("save_prod4_cosmic_results_joint_"+hietag+"_th23delta.root").c_str());
70 
71  std::cout << "\n Before the loop" << std::endl;
72  for (int j= 1; j< 4; j++){
73  std::cout << "\t Looking at loop j = " << j << std::endl;
74 
75  //DrawContourCanvas("delta_th23", 0., 2., 0., 1);
76  DrawContourCanvas("delta_th23", 0., 2., 0.2, 0.75);
77  TagName = "compare_joint_";
78  auto legend = new TLegend(0.2,0.2,0.5,0.4);
79 
80  for (int i = 0; i < int(cases.size()); i++){
81  std::cout << "\t\t Looking at loop i = " << i << std::endl;
82 
83  TH1D* dummy = new TH1D("tmp", "tmp", 5, 0, 5);
84  TagName+=cases[i].nickname+"_";
85  TFile * infile = new TFile (cases[i].filename.c_str(),"read");
86  auto mins =* (TVectorD*)infile->Get("overall_minchi");
87  double minchi23 = mins[0];
88  std::cout << "\t\t Min chi in file " << cases[i].nickname << " : " << minchi23 << std::endl;
89 
90  auto surf_th23_dcp = *FrequentistSurface::LoadFrom(infile, (TString)"surf_delta_th23_"+(hie>0?"NH":"IH"));
91  TH2 * surf_th23_dcp_xSigma;
92  if (j==1) surf_th23_dcp_xSigma = Gaussian68Percent2D(surf_th23_dcp);
93  if (j==2) surf_th23_dcp_xSigma = Gaussian2Sigma2D(surf_th23_dcp);
94  if (j==3) surf_th23_dcp_xSigma = Gaussian3Sigma2D(surf_th23_dcp);
95  std::cout << "\t\t Before color" << std::endl;
96  dummy->SetLineColor(cases[i].color);
97  std::cout << "\t\t Before leg" << std::endl;
98  legend->AddEntry(dummy,cases[i].nickname.c_str(),"l");
99  surf_th23_dcp.DrawContour(surf_th23_dcp_xSigma, cases[i].style, cases[i].color, minchi23);
100 
101  std::cout << "\t Get keys for delta_th23" << std::endl;
102  if (infile->GetListOfKeys()->Contains((TString)"delta_th23_"+std::to_string(j)+"_sigma0") && hie<0){
103  TGraph * gr = (TGraph*)infile->Get((TString)"delta_th23_"+std::to_string(j)+"_sigma0");
104  gr->SetLineColor(kGray+2);
105  legend->AddEntry(gr, "prod4", "l");
106  gr->Draw("l same");
107  }
108 
109  std::cout << "\t Enter if (hie > 0)" << std::endl;
110  if (hie > 0){
111  int i =2;
112  for ( int k = 0; k < i; k++ ){
113  TGraph * gr = (TGraph*)infile->Get((TString)"delta_th23_"+std::to_string(j)+"_sigma"+std::to_string(k));
114  gr->SetLineColor(kGray+2);
115  if (k==0) legend->AddEntry(gr, "prod4","l");
116  gr->Draw("l same");
117  }
118  }
119  }
120 
121  std::cout << "\t And draw!" << std::endl;
122  Simulation();
123  legend->Draw();
124  gPad->Print((TString)"compare_surf_th23_dcp_"+TagName+(hie>0 ? "nh_" : "ih_")+std::to_string(j)+"contour_wcosmic.pdf");
125  }
126  }
127  // If !isTH23dcp
128  else {
129  std::cout << "Looking at !isTH23dcp" << std::endl;
130  std::string surf= "surf_th23_dm32_" + hietag;
131 
132  for (int j= 1; j< 4; j++){
133  std::cout << "\t Looking at loop j = " << j << std::endl;
134 
135  if (hie>0) DrawContourCanvas("th23_dm32", 0.35, 0.65, 2.20, 2.70);
136  else DrawContourCanvas("th23_dm32", 0.35, 0.65, -2.90, -2.30);
137 
138  TagName = "joint_";
139  auto legend = new TLegend(0.15,0.18,0.45,0.35);
140  legend -> SetFillColor(0);
141 
142  for (int i = 0; i < int(cases.size()); i++){
143  std::cout << "\t\t Looking at loop i = " << i << " --> Nickname " << cases[i].nickname << std::endl;
144 
145  TH1D* dummy = new TH1D("tmp", "tmp", 5, 0, 5);
146  TagName += cases[i].nickname+"_";
147 
148  TFile * infile = new TFile (cases[i].filename.c_str(),"read");
149  auto mins =* (TVectorD*)infile->Get("overall_minchi");
150  double minchi23 = mins[0];
151 
152  std::cerr << "\t\t\t Min chi in file " << minchi23 << std::endl;
153  auto surf_th23_dm32 = *FrequentistSurface::LoadFrom(infile, surf.c_str() ) ;
154  TH2 * surf_th23_dm32_xSigma;
155  if (j==1) {
156  std::cerr << "\t\t\t Make Gaussian68Percent2D" << std::endl;
157  surf_th23_dm32_xSigma = Gaussian68Percent2D(surf_th23_dm32);
158  }
159  else if (j==2) {
160  std::cerr << "\t\t\t Make Gaussian2Sigma2D" << std::endl;
161  surf_th23_dm32_xSigma = Gaussian2Sigma2D(surf_th23_dm32);
162  }
163  else if (j==3) {
164  std::cerr << "\t\t\t Make Gaussian3Sigma2D" << std::endl;
165  surf_th23_dm32_xSigma = Gaussian3Sigma2D(surf_th23_dm32);
166  }
167  std::cout << "\t\t\t Adding Colour " << cases[i].color << ", Style " << cases[i].style << ", Name " << cases[i].nickname << std::endl;
168  dummy->SetLineColor(cases[i].color);
169  dummy->SetLineStyle(cases[i].style);
170  legend->AddEntry(dummy,cases[i].nickname.c_str(),"l");
171 
172  surf_th23_dm32.DrawContour(surf_th23_dm32_xSigma, cases[i].style, cases[i].color, minchi23);
173  //}
174 
175  if (infile->GetListOfKeys()->Contains((TString)"dmsq_"+std::to_string(j)+"_sigma0")){
176  std::cout << "\t\t Adding Sigma0 graph" << std::endl;
177  TGraph * gr = (TGraph*)infile->Get((TString)"dmsq_"+std::to_string(j)+"_sigma0");
178  gr->SetLineColor(kGray+2);
179  legend->AddEntry(gr, "prod4", "l");
180  gr->Draw("l same");
181  }
182  }
183  Simulation();
184  legend->Draw();
185  std::string PadName = "Compare_"+hietag+"_"+TagName+std::to_string(j)+"contour_wcosmic.pdf";
186  gPad->Print( PadName.c_str() );
187  }
188  }
189  } // for (int hie:{+1} )
190 
191 }
void Simulation()
Definition: tools.h:16
bin1_2sigma SetFillColor(3)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void compare_fits(bool isTH23dcp=false)
Definition: compare_fits.C:31
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
OStream cerr
Definition: OStream.cxx:7
TH2 * Gaussian2Sigma2D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 2D in gaussian approximation.
string filename
Definition: shutoffs.py:106
string infile
TH2 * Gaussian3Sigma2D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 2D in gaussian approximation.
const double j
Definition: BetheBloch.cxx:29
void DrawContourCanvas(TString surfName, double minx=0, double maxx=2, double miny=0, double maxy=1)
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
static std::unique_ptr< FrequentistSurface > LoadFrom(TDirectory *dir, const std::string &name)
surf
Definition: demo5.py:35