CAFContourMaker_module.cc
Go to the documentation of this file.
1 //CAF Contour Maker
2 //
3 //Prabhjot Singh May 30, 2017
4 //
5 //
6 #include <string>
7 #include <set>
8 #include <memory>
9 #include <sstream>
10 
11 #include "TH1.h"
12 #include "TH2.h"
13 #include "TGraph.h"
14 #include "TCanvas.h"
15 #include "TMath.h"
16 #include "TLegend.h"
17 #include "TMarker.h"
18 #include "TFile.h"
19 
20 // Framework includes
26 #include "fhiclcpp/ParameterSet.h"
29 
30 // NOvA includes
34 
35 namespace fnex {
36  class CAFContourMaker;
37 }
38 
40  public:
41  explicit CAFContourMaker(fhicl::ParameterSet const& pset);
42  virtual ~CAFContourMaker();
43 
44  void analyze(art::Event const & e) override;
45  void reconfigure (fhicl::ParameterSet const& p);
46 
47  private:
48 
49  void Make2DContours();
50 
51  void StoreContourGraphs(std::vector<double> const& percentiles,
52  std::vector< std::vector< std::unique_ptr<TGraph> > > & graphs,
53  TString analysis,
54  TH2F* hist) const;
55 
57 
58  std::string fInputRootFile; ///< The CAF .root file containing the histograms for making CAF contours
59  double fContourLevel1Sigma2D; ///< Level value for 2D upper contour
60  double fContourLevel2Sigma2D; ///< Level value for 2D upper contour
61  double fContourLevel3Sigma2D; ///< Level value for 3Sigma contour
62  std::string fContour1SigmaLabel; ///< Legend name for lower contour
63  std::string fContour2SigmaLabel; ///< Legend name for upper contour
64  std::string fContour3SigmaLabel; ///< Legend name for 3Sigma contour
65  std::string fHistogramName; ///< Histogram Name
66 
67  };
68 
69  //......................................................................
70 
72  :
73  EDAnalyzer(pset)
74  {
75  this->reconfigure(pset);
76 
77  return;
78  }
79 
80  //......................................................................
81 
83  {
84 
85  }
86 
87  //......................................................................
89  {
90 
91  fInputRootFile = pset.get<std::string >("InputRootFile", " " );
92  fContourLevel1Sigma2D = pset.get<double >("ContourLevel1Sigma2D" );
93  fContourLevel2Sigma2D = pset.get<double >("ContourLevel2Sigma2D" );
94  fContourLevel3Sigma2D = pset.get<double >("ContourLevel3Sigma2D" );
95  fContour1SigmaLabel = pset.get<std::string >("Contour1SigmaLabel", " " );
96  fContour2SigmaLabel = pset.get<std::string >("Contour2SigmaLabel", " " );
97  fContour3SigmaLabel = pset.get<std::string >("Contour3SigmaLabel", " " );
98  fHistogramName = pset.get<std::string >("HistogramName", " " );
99 
100  return;
101  }
102 
103  //......................................................................
104 
106  {
107  this->Make2DContours();
108  }
109  //......................................................................
110 
112 {
113 
115  TFile f(fInputRootFile.c_str(), "READ");
116 
118 
119  TH2F* hist_caf = (TH2F*)f.Get(histogram.c_str());
120  std::string title = ";"; title += hist_caf->GetXaxis()->GetTitle();
121  title += ";"; title += hist_caf->GetYaxis()->GetTitle();
122 
123  TH2F* hist_chisq = fTFS->make<TH2F>("hist_chisq",
124  title.c_str(),
125  hist_caf->GetXaxis()->GetNbins(),
126  hist_caf->GetXaxis()->GetXmin(),
127  hist_caf->GetXaxis()->GetXmax(),
128  hist_caf->GetYaxis()->GetNbins(),
129  hist_caf->GetYaxis()->GetXmin(),
130  hist_caf->GetYaxis()->GetXmax() );
131 
132  TH2F* hist = fTFS->make<TH2F>("hist_prob",
133  title.c_str(),
134  hist_caf->GetXaxis()->GetNbins(),
135  hist_caf->GetXaxis()->GetXmin(),
136  hist_caf->GetXaxis()->GetXmax(),
137  hist_caf->GetYaxis()->GetNbins(),
138  hist_caf->GetYaxis()->GetXmin(),
139  hist_caf->GetYaxis()->GetXmax() );
140 
141  for(int x = 0; x < hist->GetXaxis()->GetNbins() + 1; ++x){
142  for(int y = 0; y < hist->GetYaxis()->GetNbins()+ 1; ++y){
143  hist ->SetBinContent(x, y, 1.0-TMath::Prob(hist_caf->GetBinContent(x, y), 2) );
144  hist_chisq->SetBinContent(x, y, hist_caf->GetBinContent(x, y) );
145  }
146  }
147 
148  TString analysis = "CAF";
149  std::string histName = "CLContours";
150 
151  TH2F *hBackdrop = fTFS->make<TH2F>(histName.c_str(),
152  title.c_str(),
153  hist->GetXaxis()->GetNbins(),
154  hist->GetXaxis()->GetXmin(),
155  hist->GetXaxis()->GetXmax(),
156  hist->GetYaxis()->GetNbins(),
157  hist->GetYaxis()->GetXmin(),
158  hist->GetYaxis()->GetXmax());
159 
160  hBackdrop->GetXaxis()->CenterTitle();
161  hBackdrop->GetXaxis()->SetDecimals();
162 
163  hBackdrop->GetYaxis()->CenterTitle();
164  hBackdrop->GetYaxis()->SetDecimals();
165  hBackdrop->GetYaxis()->SetTitleOffset(1.0);
166 
167  // Grab the contours and color them appropriately
168  std::vector<double> cl_percentiles(3);
169  cl_percentiles[0] = fContourLevel1Sigma2D;
170  cl_percentiles[1] = fContourLevel2Sigma2D;
171  cl_percentiles[2] = fContourLevel3Sigma2D;
172 
173  std::vector< std::vector<std::unique_ptr<TGraph> > > cl_contours;
174  MakeGraphs(hist, cl_percentiles, cl_contours);
175 
176  this->StoreContourGraphs(cl_percentiles, cl_contours, analysis, hist);
177 
178  int num_to_draw_graphs = std::min(cl_contours.size(), cl_percentiles.size());
179 
180  for(int ialpha = 0; ialpha < num_to_draw_graphs; ++ialpha){
181  for(unsigned int igraph = 0; igraph < cl_contours[ialpha].size(); ++igraph){
182  if(ialpha == 0){
183  cl_contours[ialpha][igraph]->SetLineStyle(2);
184  cl_contours[ialpha][igraph]->SetLineColor(kBlue);
185  }
186  if(ialpha == 1){
187  cl_contours[ialpha][igraph]->SetLineStyle(1);
188  cl_contours[ialpha][igraph]->SetLineColor(kRed);
189  }
190  if(ialpha == 2){
191  cl_contours[ialpha][igraph]->SetLineStyle(1);
192  cl_contours[ialpha][igraph]->SetLineColor(6);
193  }
194  cl_contours[ialpha][igraph]->SetLineWidth(2);
195  }
196  }
197 
198  TString canName = "GaussianContoursCanv"+analysis; //canName += id->getUniqID();
199  TCanvas * can_gaus = fTFS->make<TCanvas>(canName.Data(), "Gaussian Contours", 1200, 800);
200  can_gaus->cd(1);
201  TString origSurfaceTitle = hBackdrop->GetTitle();
202  hBackdrop->SetTitle("Gaussian Probability Contours"+analysis);
203  hBackdrop->SetStats(false);
204  hBackdrop->Draw();
205  hBackdrop->SetTitle(origSurfaceTitle.Data());
206  for(int ialpha = 0; ialpha < num_to_draw_graphs; ++ialpha){
207  for(size_t igraph = 0; igraph < cl_contours[ialpha].size(); ++igraph){
208  cl_contours[ialpha][igraph]->Draw("C SAME");
209  }
210  }
211 //
212  TLegend legContours(0.76, 0.10, 0.90, 0.27);
213  if(cl_contours[0].size() > 0) cl_contours[0][0]->SetName("cl68"); // ROOT is inconsistent, and you must add graphs
214  if(cl_contours[1].size() > 0) cl_contours[1][0]->SetName("cl90"); // to a TLegend by name, rather than by pointer
215  if(cl_contours[2].size() > 0) cl_contours[2][0]->SetName("cl95");
216  for(int ialpha = 0; ialpha< num_to_draw_graphs; ++ialpha){
217  if(cl_contours[1].size() > 0 && ialpha == 1) legContours.AddEntry( "cl90", (fContour2SigmaLabel+ " CL") .c_str(), "l" );
218  if(cl_contours[0].size() > 0 && ialpha == 0) legContours.AddEntry( "cl68", (fContour1SigmaLabel+ " CL") .c_str(), "l" );
219  if(cl_contours[2].size() > 0 && ialpha == 2) legContours.AddEntry( "cl95", (fContour3SigmaLabel+ " CL") .c_str(), "l" );
220  }
221  legContours.Draw();
222 
223  can_gaus->Write();
224 
225  return;
226  }
227  //......................................................................
228  void fnex::CAFContourMaker::StoreContourGraphs(std::vector<double> const& percentiles,
229  std::vector< std::vector< std::unique_ptr<TGraph> > > & graphs,
230  TString analysis,
231  TH2F* hist) const
232  {
233  std::string XaxisTitle = hist->GetXaxis()->GetTitle();
234  std::string YaxisTitle = hist->GetYaxis()->GetTitle();
235 
236  std::string title = (";" + XaxisTitle +
237  ";" + YaxisTitle);
238 
239  for(size_t p = 0; p < graphs.size(); ++p){
240  for(size_t g = 0; g < graphs[p].size(); ++g){
241 
242  std::stringstream name;
243  name
244 // << std::setprecision(3)
245  << analysis
246  << "ContourGraph"
247  << "_CL"
248  << percentiles[p]
249  << "_Gr"
250  << g;
251  fTFS->makeAndRegister<TGraph>(name.str().c_str(),
252  title.c_str(),
253  graphs[p][g]->GetN(),
254  graphs[p][g]->GetX(),
255  graphs[p][g]->GetY());
256 
257  } // end loop over graphs for the percentile
258  } // end loop over different percentile graphs
259 
260  return;
261  }//end of StoreContourGraphs
262 
263 
264 
const XML_Char * name
Definition: expat.h:151
double fContourLevel2Sigma2D
Level value for 2D upper contour.
const char * p
Definition: xmltok.h:285
double fContourLevel1Sigma2D
Level value for 2D upper contour.
art::ServiceHandle< art::TFileService > fTFS
string filename
Definition: shutoffs.py:106
Create a list of fnex::Events to be used in fits.
void analyze(art::Event const &e) override
DEFINE_ART_MODULE(SpectrumTest)
void MakeGraphs(TH2F *in_hist, std::vector< double > vContours, std::vector< std::vector< std::unique_ptr< TGraph > > > &vGraphVectors)
Definition: PlotAssist.cxx:15
std::string fHistogramName
Histogram Name.
std::string fContour2SigmaLabel
Legend name for upper contour.
std::string fContour1SigmaLabel
Legend name for lower contour.
T get(std::string const &key) const
Definition: ParameterSet.h:231
CAFContourMaker(fhicl::ParameterSet const &pset)
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
T * makeAndRegister(char const *name, char const *title, ARGS...args) const
void reconfigure(fhicl::ParameterSet const &p)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
double fContourLevel3Sigma2D
Level value for 3Sigma contour.
std::string fContour3SigmaLabel
Legend name for 3Sigma contour.
T min(const caf::Proxy< T > &a, T b)
Float_t e
Definition: plot.C:35
void StoreContourGraphs(std::vector< double > const &percentiles, std::vector< std::vector< std::unique_ptr< TGraph > > > &graphs, TString analysis, TH2F *hist) const
std::string fInputRootFile
The CAF .root file containing the histograms for making CAF contours.