PlotUtilities.cxx
Go to the documentation of this file.
1 //
2 // Created by Brian Rebel on 11/19/19.
3 //
4 #include <sstream>
5 #include <iomanip>
6 
7 #include "TCanvas.h"
8 #include "TMath.h"
9 #include "TLegend.h"
10 #include "TMarker.h"
11 #include "TList.h"
12 #include "TLine.h"
13 #include "TROOT.h"
14 #include "TKey.h"
15 #include "TObjArray.h"
16 
20 
21 namespace cmf{
22 
23  static PlotUtilities* gPlotUtil = nullptr;
24 
25  //-----------------------------------------------------------------------------
27  {
28  if(gPlotUtil == nullptr) gPlotUtil = new PlotUtilities();
29 
30  return gPlotUtil;
31  }
32 
33  //-----------------------------------------------------------------------------
38  , fXParamExtrema (std::make_pair(cmf::kGarbageDouble, std::numeric_limits<double>::lowest()))
39  , fYParamExtrema (std::make_pair(cmf::kGarbageDouble, std::numeric_limits<double>::lowest()))
40  , fXParamScale (std::make_pair(1., ""))
41  , fYParamScale (std::make_pair(1., ""))
42  , fContourSigmaLabels(std::vector<std::string>(3))
43  {
44  }
45 
46  //-----------------------------------------------------------------------------
48  {
49  fUseTrig = plottingPars.get<bool >("UseTrig", false );
50  fContourLevel1Sigma1D = plottingPars.get<double >("ContourLevel1Sigma1D", 1 );
51  fContourLevel2Sigma1D = plottingPars.get<double >("ContourLevel2Sigma1D", 4 );
52  fContourLevel3Sigma1D = plottingPars.get<double >("ContourLevel3Sigma1D", 9 );
53  fContourLevel1Sigma2D = plottingPars.get<double >("ContourLevel1Sigma2D", 2.3 );
54  fContourLevel2Sigma2D = plottingPars.get<double >("ContourLevel2Sigma2D", 6.18 );
55  fContourLevel3Sigma2D = plottingPars.get<double >("ContourLevel3Sigma2D", 11.83 );
56  fContourSigmaLabels[0] = plottingPars.get<std::string>("Contour1SigmaLabel", "1#sigma");
57  fContourSigmaLabels[1] = plottingPars.get<std::string>("Contour2SigmaLabel", "2#sigma");
58  fContourSigmaLabels[2] = plottingPars.get<std::string>("Contour3SigmaLabel", "3#sigma");
59  fNormalizeBins = plottingPars.get<bool >("NormalizeBinWidths", false );
60  fNuMuNormVal = plottingPars.get<double >("NuMuNormVal" , 0.1 );
61  fNuENormVal = plottingPars.get<double >("NuENormVal" , 0.5 );
62  fNCNormVal = plottingPars.get<double >("NCNormVal" , 1.0 );
63  }
64 
65  //-----------------------------------------------------------------------------
67  cmf::OscParm_t const& xParEnum,
68  cmf::OscParm_t const& yParEnum,
69  size_t const& numXParPoints,
70  size_t const& numYParPoints,
71  std::pair<double, double> const& parXExtrema,
72  std::pair<double, double> const& parYExtrema,
73  bool const& logXAxis,
74  bool const& logYAxis)
75  {
76  this->Initialize(plottingPars);
77 
78  fXParamExtrema = parXExtrema;
79  fYParamExtrema = parYExtrema;
80  fXParamDivs = numXParPoints;
81  fYParamDivs = numYParPoints;
82  fLogXAxis = logXAxis;
83  fLogYAxis = logYAxis;
84  fXParamEnum = xParEnum;
85  fYParamEnum = yParEnum;
86 
88  throw cet::exception("ContourMaker")
89  << "x and y parameters are the same "
91  << " - that cannot be good, bail";
92 
94  cmf::cOscParams_LatexScales[fXParamEnum]);
96  cmf::cOscParams_LatexScales[fYParamEnum]);
97 
98  MF_LOG_VERBATIM("PlotUtilities")
99  << "Plot Utilities: "
100  << "X Parameter: "
102  << " ("
103  << fXParamExtrema.first
104  << ", "
105  << fXParamExtrema.second
106  << ") numDivs: "
107  << fXParamDivs
108  << " Log axis: "
109  << fLogXAxis
110  << " Y Parameter: "
112  << " ("
113  << fYParamExtrema.first
114  << ", "
115  << fYParamExtrema.second
116  << ") numDivs: "
117  << fYParamDivs
118  << " Log axis: "
119  << fLogYAxis;
120 
121  float paramDiv;
122  float paramStart;
123  float paramStop;
124 
125  // set the size of the vectors for the parameter bins to be the number of
126  // divisions + 2 to keep track of the end points
127  fXParamBinEdges = std::vector<float>(fXParamDivs + 2);
128  fYParamBinEdges = std::vector<float>(fYParamDivs + 2);
129 
130  if(fLogXAxis){
131  // convert extrema to log values - if we are logging a parameter we
132  // have a scale of 1, so don't bother using the f*ParamScale
133  float logXParamExtremaLow = std::log10(fXParamExtrema.first);
134  float logXParamExtremaHigh = std::log10(fXParamExtrema.second);
135  paramDiv = (logXParamExtremaHigh - logXParamExtremaLow)/float(fXParamDivs);
136  paramStart = (logXParamExtremaLow - 0.5 * paramDiv);
137  paramStop = (logXParamExtremaHigh + 0.5 * paramDiv);
138  for(size_t i = 0; i < fXParamBinEdges.size() - 1; ++i){
139  fXParamBinEdges[i] = std::pow(10,(paramStart + (i * paramDiv)));
140  }
141  fXParamBinEdges.back() = std::pow(10, paramStop);
142  }
143  else {
144 
145  paramDiv = (fXParamExtrema.second - fXParamExtrema.first )/float(fXParamDivs);
146  paramStart = (fXParamExtrema.first - 0.5 * paramDiv) * fXParamScale.first;
147  paramStop = (fXParamExtrema.second + 0.5 * paramDiv) * fXParamScale.first;
148  for (size_t i = 0; i < fXParamBinEdges.size() - 1; ++i){
149  fXParamBinEdges[i] = paramStart + (i * paramDiv * fXParamScale.first);
150  }
151  fXParamBinEdges.back() = paramStop;
152  }
153 
154  // for(auto const& itr : fXParamBinEdges){
155  // MF_LOG_VERBATIM("PlotUtilities")
156  // << "x bin edges: "
157  // << itr;
158  // }
159 
160  if(fLogYAxis)
161  {
162  // convert extrema to log values - if we are logging a parameter we
163  // have a scale of 1, so don't bother using the f*ParamScale
164  float logYParamExtremaLow = std::log10(fYParamExtrema.first);
165  float logYParamExtremaHigh = std::log10(fYParamExtrema.second);
166  paramDiv = (logYParamExtremaHigh - logYParamExtremaLow)/float(fYParamDivs);
167  paramStart = (logYParamExtremaLow - 0.5 * paramDiv);
168  paramStop = (logYParamExtremaHigh + 0.5 * paramDiv);
169  for (size_t i = 0; i < fYParamBinEdges.size() - 1; ++i){
170  fYParamBinEdges[i] = std::pow(10,(paramStart + (i * paramDiv)));
171  }
172  fYParamBinEdges.back() = std::pow(10, paramStop);
173  }
174  else {
175  paramDiv = (fYParamExtrema.second - fYParamExtrema.first )/float(fYParamDivs);
176  paramStart = (fYParamExtrema.first - 0.5 * paramDiv) * fYParamScale.first;
177  paramStop = (fYParamExtrema.second + 0.5 * paramDiv) * fYParamScale.first;
178 
179  for (size_t i = 0; i < fYParamBinEdges.size() - 1; ++i){
180  fYParamBinEdges[i] = paramStart + (i * paramDiv * fYParamScale.first);
181  }
182  fYParamBinEdges.back() = paramStop;
183  }
184 
185  // for(auto const& itr : fYParamBinEdges){
186  // MF_LOG_VERBATIM("PlotUtilities")
187  // << "y bin edges: "
188  // << itr;
189  // }
190 
191  }
192 
193  //......................................................................
194  void PlotUtilities::Make1DPlot(art::TFileDirectory & tfd,
195  std::string const& chiSqrName,
196  std::map<double, double> const& paramChiSqr,
197  bool xPar)
198  {
199 
200  std::vector<float> paramBins = fYParamBinEdges;
201  cmf::OscParm_t paramEnum = fYParamEnum;
202  std::pair<double, std::string> paramScale = fYParamScale;
203  if(xPar){
204  paramBins = fXParamBinEdges;
205  paramEnum = fXParamEnum;
206  paramScale = fXParamScale;
207  }
208 
209  std::string grName = chiSqrName + cmf::cOscParams_Strings[paramEnum] + "DeltaChiSqr";
210 
211  //make vectors holding the points for the graph
212  std::vector<double> param;
213  std::vector<double> deltaChiSqr;
214 
215  for(auto itr : paramChiSqr){
216  param .emplace_back(itr.first * paramScale.first);
217  deltaChiSqr.emplace_back(itr.second);
218 
219  MF_LOG_DEBUG("PlotUtilities")
220  << grName
221  << " "
222  << param.back()
223  << " "
224  << deltaChiSqr.back();
225  }
226 
227  // make the graph
228  auto *gr = tfd.makeAndRegister<TGraph>(grName.c_str(),
229  grName.c_str(),
230  param.size(),
231  param.data(),
232  deltaChiSqr.data());
233 
234  // make a histogram to set the axis ranges, maximum DeltaChiSqr = 15
235  std::string histName = grName + "Hist";
237  if (paramEnum == cmf::kTh23 && fUseTrig) title = ";sin^{2}(#theta_{23})";
238  else if(paramEnum == cmf::kTh24 && fUseTrig) title = ";sin^{2}(#theta_{24})";
239  else if(paramEnum == cmf::kTh34 && fUseTrig) title = ";sin^{2}(#theta_{34})";
240  title += paramScale.second + ";#Delta#chi^{2}";
241 
242  TH1F *hist = tfd.make<TH1F>(histName.c_str(),
243  title.c_str(),
244  paramBins.size() - 1,
245  paramBins.data());
246  hist->SetMaximum(15.);
247  hist->GetXaxis()->CenterTitle();
248  hist->GetXaxis()->SetDecimals();
249 
250  hist->GetYaxis()->CenterTitle();
251  hist->GetYaxis()->SetDecimals();
252 
253  // Make TLines to show CLs
254  std::vector<Color_t> sigmaColors({kBlue, kRed, kMagenta});
255  std::vector<double> cl_percentiles({fContourLevel1Sigma1D,
258  std::vector<TLine*> sigmaLines(3);
259  for(size_t l = 0; l < sigmaLines.size(); ++l){
260  sigmaLines[l]= new TLine(paramBins.front(), cl_percentiles[l], paramBins.back(), cl_percentiles[l]);
261  sigmaLines[l]->SetLineStyle(2);
262  sigmaLines[l]->SetLineColor(sigmaColors[l]);
263  }
264 
265  MF_LOG_DEBUG("PlotUtilities")
266  << " cl_percentile[0]: "
267  << cl_percentiles[0]
268  << " "
269  << fContourSigmaLabels[0]
270  << " cl_percentile[1]: "
271  << cl_percentiles[1]
272  << " "
273  << fContourSigmaLabels[1]
274  << " cl_percentile[2]: "
275  << cl_percentiles[2]
276  << " "
277  << fContourSigmaLabels[2];
278 
279  // Add a TLegend explain them
280  auto *legLines = new TLegend(0.8, 0.1, 0.9, 0.24);
281  for(size_t l = 0; l < sigmaLines.size(); ++l){
282  legLines->AddEntry(sigmaLines[l], (fContourSigmaLabels[l] + " CL").c_str(), "l");
283  }
284 
285  legLines->SetBorderSize(0);
286  legLines->SetFillStyle(0);
287 
288  // make a Canvas to hold it all
289  std::string canvName = grName + "Canv";
290  auto *canv = tfd.makeAndRegister<TCanvas>(canvName.c_str(),
291  canvName.c_str(),
292  2);
293 
294  // Log the x axis if this parameter is to be logged - this is a 1D plot
295  // so even for the Y parameter we log the x axis
296  if ( xPar && fLogXAxis) canv->SetLogx();
297  else if(!xPar && fLogYAxis) canv->SetLogx();
298 
299  canv->cd();
300  hist->SetStats(false);
301  hist->Draw();
302  gr->Draw("lsame");
303  for(size_t l = 0; l < sigmaLines.size(); ++l) sigmaLines[l]->Draw();
304  legLines->Draw();
305 
306  // We're done here
307  canv->ForceUpdate();
308  }
309 
310  //......................................................................
311  void PlotUtilities::Make2DContours(art::TFileDirectory & tfd,
312  std::string const& chiSqrName,
313  cmf::PointMap const& twoDChiSqr,
314  cmf::GridPoint const& bestFit,
315  bool useBestFit)
316  {
317  TH2F *hist = this->Make2DHist(tfd,
318  chiSqrName,
319  twoDChiSqr);
320 
321  std::vector<double> percentiles({fContourLevel1Sigma2D,
324  std::vector<std::string> contourNames({"1sigma",
325  "2sigma",
326  "3sigma"});
327  std::vector<std::string> contourLabels({fContourSigmaLabels[0],
329  fContourSigmaLabels[2]});
330 
331 
332  // check that the minimum for the histogram is less than the lowest contour level,
333  // if not, then get rid of that level
334  if(hist->GetMinimum() > fContourLevel1Sigma2D){
335  percentiles .erase(percentiles.begin() );
336  contourNames .erase(contourNames.begin() );
337  contourLabels.erase(contourLabels.begin());
338  }
339 
340 
341  // MF_LOG_VERBATIM("PlotUtilities")
342  // << "draw contours for ";
343  // for(auto const& itr : percentiles)
344  // MF_LOG_VERBATIM("PlotUtilities")
345  // << itr;
346 
347  hist->SetContour(percentiles.size(), percentiles.data());
348 
349  // This avoids disturbing canvases already
350  // Yes, we do need a canvas for this to work
351  // yes, we do need to call Update() on it
352  // No, ROOT is not very friendly.
353  TCanvas temp_can;
354  temp_can.cd();
355  hist->Draw("contlist");
356  temp_can.Update();
357 
358  auto * plah = (TObjArray *)(gROOT->GetListOfSpecials()->FindObject("contours"));
359  if(!plah) return;
360 
361  std::vector< std::vector<std::unique_ptr<TGraph> > > graphVec(percentiles.size());
362 
363  for(int i = 0; i < plah->GetSize(); ++i){
364  std::vector<std::unique_ptr<TGraph> > temp;
365  auto * list = (TList*) plah->At(i);
366  if(!list) break;
367  if(!(list->First())) break;
368  for(int igraph = 0; igraph < list->GetSize(); ++igraph){
369  if(list->At(igraph)) temp.emplace_back( std::make_unique<TGraph>( *((TGraph*) list->At(igraph))) );
370  }
371  graphVec[i] = std::move(temp);
372  }
373 
374  MF_LOG_DEBUG("PlotAssist")
375  << "Found : "
376  << plah->GetSize()
377  << " sets of contours\n"
378  << "Parameter 1 and Enum 1: "
380  << " "
381  << "Parameter 2 and Enum 2: "
383  << " there are "
384  << graphVec.size()
385  << " vectors of graphs";
386 
387  std::string title = ";";
388  title += hist->GetXaxis()->GetTitle();
389  title += ";";
390  title += hist->GetYaxis()->GetTitle();
391 
392  int p = 0;
393  int g;
394  for(auto const& gritr : graphVec){
395 
396  MF_LOG_DEBUG("PlotUtilities")
397  << "there are "
398  << gritr.size()
399  << " graphs for contour "
400  << p;
401 
402  g = 0;
403  for(auto const& itr : gritr){
404 
405  std::stringstream name;
406  name << std::setprecision(3)
407  << chiSqrName
408  << "ContourGraph"
411  << "_CL"
412  << contourLabels[p]
413  << "_Gr"
414  << g;
415 
416  MF_LOG_DEBUG("PlotUtilities")
417  << "store graph with name "
418  << name.str();
419 
420  tfd.makeAndRegister<TGraph>(name.str().c_str(),
421  title.c_str(),
422  itr->GetN(),
423  itr->GetX(),
424  itr->GetY());
425 
426  ++g;
427  } // end loop over graphs for the percentile
428  ++p;
429  } // end loop over different percentile graphs
430 
431  //MF_LOG_DEBUG << " Found " << graphVec.size() << " standard cl_contour sets\n";
432 
433  // Assign colors and line styles
434  // have to do some checks because we could have either 2 or 3 sets of graphs
435  //-------
436  g = 0;
437  for(auto & grv : graphVec){
438  for(auto & gr : grv){
439  if(percentiles.size() == 3 && g == 0){
440  gr->SetLineColor(kBlue);
441  }
442  if((percentiles.size() == 3 && g == 1) || (percentiles.size() == 2 && g == 0)){
443  gr->SetLineColor(kRed);
444  }
445  if((percentiles.size() == 3 && g == 2) || (percentiles.size() == 2 && g == 1)){
446  gr->SetLineColor(6);
447  }
448  gr->SetLineWidth(2);
449  }
450  ++g;
451  }
452 
453  // Make the 'backdrop' histogram against which the contours
454  // will be drawn
455 
456  MF_LOG_DEBUG("PlotUtilities")
457  << "Make2DContours function: "
458  << "Parameter 1 and Enum 1: "
460  << " "
461  << "Parameter 2 and Enum 2: "
463 
464  std::string histName = (chiSqrName +
465  cmf::cOscParams_Strings[fXParamEnum] +
466  cmf::cOscParams_Strings[fYParamEnum] +
467  "CLContours");
468 
469  TH2F *hBackdrop = tfd.make<TH2F>(histName.c_str(),
470  title.c_str(),
471  hist->GetXaxis()->GetNbins(),
472  hist->GetXaxis()->GetXmin(),
473  hist->GetXaxis()->GetXmax(),
474  hist->GetYaxis()->GetNbins(),
475  hist->GetYaxis()->GetXmin(),
476  hist->GetYaxis()->GetXmax());
477 
478  hBackdrop->GetXaxis()->CenterTitle();
479  hBackdrop->GetXaxis()->SetDecimals();
480 
481  hBackdrop->GetYaxis()->CenterTitle();
482  hBackdrop->GetYaxis()->SetDecimals();
483  hBackdrop->GetYaxis()->SetTitleOffset(1.0);
484 
485  // Make a new canvas, plot the CL contours,
486  //-------
487  std::string canName(chiSqrName + "ContoursCanv");
488  auto *canv = tfd.make<TCanvas>(canName.c_str(),
489  (chiSqrName + " Contours").c_str(),
490  1200,
491  800);
492 
493  hBackdrop->Draw();
494  for(auto const& grv : graphVec){
495  for(auto const& itr : grv){
496  itr->Draw("C SAME");
497  }
498  }
499 
500  // Add a TLegend to describe the CL contours
501  // ROOT is inconsistent, and you must add graphs
502  // to a TLegend by name, rather than by pointer
503  //-------
504  TLegend legContours(0.8, 0.1, 0.9, 0.24);
505  legContours.SetBorderSize(0);
506  legContours.SetFillStyle(0);
507  g = 0;
508  for(auto & grv : graphVec){
509  if(!grv.empty()){
510  grv.front()->SetName(contourNames[g].c_str());
511  legContours.AddEntry(grv.front().get(),
512  (contourLabels[g] + " CL").c_str(),
513  "l");
514  }
515  ++g;
516  }
517  legContours.Draw();
518 
519  if(useBestFit){
520  // And a TLegend to describe the initial guess and
521  // best fit points
522  //-------
523  TLegend legPoints(0.12, 0.11, 0.27, 0.25);
524  legPoints.SetBorderSize(0);
525  legPoints.SetFillStyle(0);
526 
527  std::stringstream par1_legendentry;
528  std::stringstream par2_legendentry;
529 
531  if(fXParamEnum == cmf::kTh23 && fUseTrig) Parameter1Label = "sin^{2}(#theta_{23})";
532  if(fXParamEnum == cmf::kTh24 && fUseTrig) Parameter1Label = "sin^{2}(#theta_{24})";
533 
534  par1_legendentry
535  << Parameter1Label
536  << " = "
537  << std::setprecision(3)
538  << bestFit.X()
539  << " "
540  << fXParamScale.second;
541 
543  if(fYParamEnum == cmf::kTh23 && fUseTrig) Parameter2Label = "sin^{2}(#theta_{23})";
544  if(fYParamEnum == cmf::kTh24 && fUseTrig) Parameter2Label = "sin^{2}(#theta_{24})";
545 
546  par2_legendentry
547  << Parameter2Label
548  << " = "
549  << std::setprecision(3)
550  << bestFit.Y()
551  << " "
552  << fYParamScale.second;
553 
554  // add the initial and best fit point
555 
556  TMarker bestFitMarker(fXParamScale.first * bestFit.X(),
557  fYParamScale.first * bestFit.Y(),
558  8);
559  bestFitMarker.SetMarkerColor(1);
560  bestFitMarker.Draw();
561 
562  legPoints.AddEntry(&bestFitMarker, "Best Fit", "p");
563  legPoints.AddEntry(&bestFitMarker, par1_legendentry.str().data(), "");
564  legPoints.AddEntry(&bestFitMarker, par2_legendentry.str().data(), "");
565 
566  legPoints.Draw();
567  }
568 
569  if(fLogXAxis) canv->SetLogx();
570  if(fLogYAxis) canv->SetLogy();
571 
572  canv->Update();
573  canv->Write();
574  }
575 
576  //......................................................................
577  void PlotUtilities::MakeCLHeatMaps(art::TFileDirectory & tfd,
578  std::string const& baseName,
579  std::vector<cmf::ChiSqrInfo> const& chiSqrInfoVec)
580  {
581 
582  // make a 2D histogram for each CL
583  // make a 2D histogram and fill it with the Delta chiSqr values
584  std::string name = (baseName +
587  "1sigmaHeatMap");
588 
589  std::vector<TH2F*> hmVec(3);
590 
591  for(size_t i = 0; i < hmVec.size(); ++i){
592 
593  if(i > 0)
594  name.replace(name.find(std::to_string(i) + "sigma"), 1, std::to_string(i+1));
595 
596  hmVec[i] = this->MakeGeneric2DHist(tfd,
597  name);
598  hmVec[i]->SetTitle((std::to_string(i + 1) + "#sigma").c_str());
599  }
600 
601  // loop over the chiSqrInfo structs to find the delta chiSqr at each point
602  MF_LOG_DEBUG("PlotUtilities")
603  << "there are "
604  << chiSqrInfoVec.size()
605  << " chiSqrInfo objects to loop over";
606 
607  size_t uni = 0;
608  for(auto const& csItr : chiSqrInfoVec){
609  for(auto const& twoDItr : csItr.twoDChiSqr){
610 
611  MF_LOG_DEBUG("PlotUtilities")
612  << "twoDChiSqr for universe "
613  << uni
614  << " "
615  << twoDItr.first
616  << " "
617  << twoDItr.second;
618 
619  if(1.0-TMath::Prob(twoDItr.second, 2) < fContourLevel1Sigma2D)
620  hmVec[0]->Fill(twoDItr.first.X() * fXParamScale.first,
621  twoDItr.first.Y() * fYParamScale.first);
622  if(1.0-TMath::Prob(twoDItr.second, 2) < fContourLevel2Sigma2D)
623  hmVec[1]->Fill(twoDItr.first.X() * fXParamScale.first,
624  twoDItr.first.Y() * fYParamScale.first);
625  if(1.0-TMath::Prob(twoDItr.second, 2) < fContourLevel3Sigma2D)
626  hmVec[2]->Fill(twoDItr.first.X() * fXParamScale.first,
627  twoDItr.first.Y() * fYParamScale.first);
628  }
629  ++uni;
630  } // end loop over ChiSqrInfos
631 
632  // make a canvas to hold these histograms
633  name = baseName +
636  "HeatMapCanv";
637 
638  auto *canv = tfd.make<TCanvas>(name.c_str(),
639  name.c_str(),
640  800,
641  800);
642 
643  if(fLogXAxis) canv->SetLogx();
644  if(fLogYAxis) canv->SetLogy();
645  canv->Divide(2,2);
646 
647  for(size_t i = 0; i < hmVec.size(); ++i){
648  canv->cd(i + 1);
649  hmVec[i]->Draw("colz");
650  }
651 
652  canv->Update();
653  canv->Write();
654 
655  return;
656  }
657 
658  //......................................................................
659  // Normalize the bin contents to be the equivalent to the desired width
661  double normToVal)
662  {
663  if(!fNormalizeBins) return;
664 
665  double norm;
666  for(int b = 0; b < hist->GetNbinsX(); ++b){
667  norm = normToVal / hist->GetBinWidth(b + 1);
668  hist->SetBinContent(b + 1, hist->GetBinContent(b + 1) * norm);
669  hist->SetBinError(b + 1, hist->GetBinError(b + 1) * norm);
670  } // end loop over bins
671  }
672 
673  //......................................................................
675  cmf::Spectrum const& binCount,
676  std::map<long, TH1D*> & energySpecMap,
677  std::string const& name,
678  art::TFileDirectory & tfd,
679  std::string const& uniqueID)
680  {
682  std::string histTitle;
684  TH1D *hist;
685 
686  // loop over the keys we are using to determine
687  // the offset for each selection
688  long offset;
689 
690  double binE;
691  double scale;
692 
696 
697  for(auto const& itr : cmf::SelectionUtility::Instance()->SelectionsToUseAsKeys()){
699  sel = cmf::KeyToSelectionType(itr);
700  det = cmf::KeyToDetectorType(itr);
701  beam = cmf::KeyToBeamType(itr);
702  scale = (det == cmf::kNEARDET) ? 1.e-3 : 1.;
703 
704  bins.clear();
706  det,
707  bins);
708 
709  // check that we have a sane number of bins
710  if(offset + cmf::CovarianceBinUtility::Instance()->SelectionHighEdges(itr).size() > binSpec.size())
711  throw cet::exception("CovarianceFitHelper")
712  << "Filling data energy spectra appears to want more bins"
713  << " than are in the data bin spectrum: "
715  << " vs "
716  << binSpec.size();
717 
718  histName = (name +
719  uniqueID +
723 
724  histTitle = (cmf::cBeamType_LatexStrings[beam] + " " +
725  cmf::cDetType_Strings[det] + " " +
726  cmf::cSelectionType_LatexStrings[sel] + ";Energy (GeV);");
727 
728  histTitle += (scale < 1) ? "#times10^{3} Events" : "Events";
729 
730  if(name.find("MC") != std::string::npos) histTitle.insert(0, "MC ");
731 
732  MF_LOG_DEBUG("CovarianceFitHelper")
733  << "total Energy spectrum hist is "
734  << histName
735  << " "
736  << itr;
737 
738 // for(auto const& binitr : bins)
739 // LOG_VERBATIM("CovarianceFitHelper")
740 // << "bin edge: "
741 // << binitr;
742 
743  if(energySpecMap.count(itr) < 1)
744  hist = energySpecMap.emplace(itr, tfd.make<TH1D>(histName.c_str(),
745  histTitle.c_str(),
746  bins.size() - 1,
747  bins.data())).first->second;
748  else
749  hist = energySpecMap.find(itr)->second;
750 
751  hist->GetXaxis()->CenterTitle();
752  hist->GetYaxis()->CenterTitle();
753  if(name.find("MC") != std::string::npos){
754  hist->SetMarkerColor(kRed);
755  hist->SetLineColor(kRed);
756  }
757 
758  for(size_t b = offset; b < offset + cmf::CovarianceBinUtility::Instance()->SelectionHighEdges(itr).size(); ++b){
760  hist->Fill(binE, binSpec[b] * scale);
761 
762  // the uncertainty for each bin is the fractional uncertainty for the total count in the
763  // bin (1 / sqrt(binCount[b]) multiplied by the total weight in the bin (binSpec[b])
764  // only do it if the bin count > 0
765  if(binCount[b] > 0.)
766  hist->SetBinError(hist->FindBin(binE),
767  scale * (binSpec[b] / std::sqrt(binCount[b])));
768  }// end filling the spectrum
769 
770  MF_LOG_DEBUG("PlotUtilities")
771  << "histogram "
772  << hist->GetName()
773  << " has "
774  << hist->Integral()
775  << " selected events";
776 
777  // now normalize the bin contents if desired to the desired width
778  if(cmf::IsNuESelected(sel))
779  this->NormalizeBinContents(hist, fNuENormVal);
780  else if(cmf::IsNuMuSelected(sel))
781  this->NormalizeBinContents(hist, fNuMuNormVal);
782  else if (cmf::IsNCSelected(sel))
783  this->NormalizeBinContents(hist, fNCNormVal);
784 
785  } // end loop over selections being used in the analysis
786  }
787 
788  //......................................................................
789  TH1F* PlotUtilities::MakeSpectrumHistogram(art::TFileDirectory & tfd,
790  std::string const& baseName,
791  std::string const& baseTitle,
792  cmf::Spectrum const& spectrum)
793  {
794  if(spectrum.size() < 1){
795  MF_LOG_VERBATIM("PlotUtilities")
796  << "spectrum size for "
797  << baseName
798  << " is 0, are you sure you filled it?";
799  return nullptr;
800  }
801 
802  auto spectrumHist = tfd.make<TH1F>(baseName.c_str(),
803  baseTitle.c_str(),
804  spectrum.size(),
805  0,
806  spectrum.size());
807 
808  spectrumHist->Sumw2();
809 
810  for(size_t b = 0; b < spectrum.size(); ++b) spectrumHist->Fill(b, spectrum[b]);
811 
812  return spectrumHist;
813  }
814 
815  //......................................................................
816  // The returned histogram has the axes and axis ranges defined by fXParam(2)
817  // only for now
818  TH2F* PlotUtilities::MakeGeneric2DHist(art::TFileDirectory & tfd,
819  std::string const& histName)
820  {
821  // Delta m^2 is always on the y axis for numu
822  // angle is always on y axis for nue contours
823 
826 
829 
830  std::string title = (";" + XaxisTitle +
831  ";" + YaxisTitle);
832 
833 
834  MF_LOG_DEBUG("ContourMaker")
835  << "Creating 2D histogram with range "
836  << fXParamBinEdges.front() << " " << fXParamBinEdges.back() << " "
837  << fYParamBinEdges.front() << " " << fYParamBinEdges.back();
838 
839  auto* hist2D = tfd.make<TH2F>(histName.c_str(),
840  title.c_str(),
841  fXParamDivs+1,
842  fXParamBinEdges.data(),
843  fYParamDivs+1,
844  fYParamBinEdges.data());
845 
846  hist2D->GetXaxis()->CenterTitle();
847  hist2D->GetXaxis()->SetDecimals();
848  hist2D->GetYaxis()->CenterTitle();
849  hist2D->GetYaxis()->SetDecimals();
850 
851  return hist2D;
852  }
853 
854  //......................................................................
855  TH2F* PlotUtilities::Make2DHist(art::TFileDirectory & tfd,
856  std::string const& chiSqrName,
857  cmf::PointMap const& twoDChiSqr)
858  {
859  // make a 2D histogram and fill it with the Delta chiSqr values
860  std::string histName = (chiSqrName +
863  "DeltaChiSqr");
864 
865  auto* hist_deltachisq = this->MakeGeneric2DHist(tfd,
866  histName);
867 
868  histName = (chiSqrName +
871  "DeltaChiSqr_Prob");
872 
873  auto* hist = this->MakeGeneric2DHist(tfd,
874  histName);
875 
876  double dcsq;
877  for(auto const& itr : twoDChiSqr){
878 
879  MF_LOG_DEBUG("ContourMaker")
880  << itr.first.X()
881  << " "
882  << itr.first.Y()
883  << " "
884  << itr.second;
885 
886  dcsq = (itr.second < 0.1) ? 0.1 : itr.second;
887  hist_deltachisq->Fill(itr.first.X() * fXParamScale.first,
888  itr.first.Y() * fYParamScale.first,
889  dcsq);
890 
891  hist->Fill(itr.first.X() * fXParamScale.first,
892  itr.first.Y() * fYParamScale.first,
893  1.0-TMath::Prob(itr.second, 2));
894  }
895 
896 // for(int b = 0; b < hist_deltachisq->GetXaxis()->GetNbins(); ++ b)
897 // MF_LOG_VERBATIM("CMFContourMaker")
898 // << "x axis "
899 // << b + 1
900 // << " "
901 // << hist_deltachisq->GetXaxis()->GetBinCenter(b + 1);
902 //
903 // for(int b = 0; b < hist_deltachisq->GetYaxis()->GetNbins(); ++ b)
904 // MF_LOG_VERBATIM("CMFContourMaker")
905 // << "y axis "
906 // << b + 1
907 // << " "
908 // << hist_deltachisq->GetYaxis()->GetBinCenter(b + 1);
909 
910  return hist;
911  }
912 
913  //......................................................................
914  void PlotUtilities::Make2DHiddenParHistogram(art::TFileDirectory & tfd,
915  std::string const& baseName,
916  cmf::PointMap const& parVals)
917  {
918  if(parVals.size() < 1){
919  MF_LOG_VERBATIM("PlotUtilities")
920  << "point map size for "
921  << baseName
922  << " is 0, are you sure you filled it?";
923  return;
924  }
925 
926  auto* parValHist = this->MakeGeneric2DHist(tfd,
927  baseName);
928 
929  for(auto const& itr : parVals){
930  parValHist->Fill(itr.first.X() * fXParamScale.first,
931  itr.first.Y() * fYParamScale.first,
932  itr.second);
933  }
934  }
935 
936 } // end cmf namespace
const XML_Char * name
Definition: expat.h:151
tree Draw("slc.nhit")
enum cmf::osc_params OscParm_t
double fContourLevel3Sigma1D
Level value for 3 sigma 1D contour.
std::pair< double, std::string > fXParamScale
scale the parameter if necessary ie is x10^{-3}, string is latex friendly
enum BeamMode kRed
int KeyToOffset(long const &key, bool allSels=false)
static bool IsNCSelected(cmf::SelectionType_t const &sel)
Definition: StaticFuncs.h:387
cmf::OscParm_t fFreeParamEnum
enumerated value of the oscillation parameter that floats freely
double fContourLevel2Sigma2D
Level value for 2 sigma 2D contour.
std::pair< double, std::string > fYParamScale
scale the parameter if necessary ie is x10^{-3}, string is latex friendly
std::vector< std::string > fContourSigmaLabels
labels for contour sigma levels
static PlotUtilities * gPlotUtil
const char * p
Definition: xmltok.h:285
const std::vector< std::string > cDetType_Strings({"UnknownDet","NearDet","FarDet","MINOSNear","MINOSFar","AllDetectors"})
static const double kGarbageDouble
Definition: Constants.h:18
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:72
std::pair< double, double > fYParamExtrema
first is the min, second is the max
bool fNormalizeBins
whether to normalize bins with different sizes
static SelectionUtility * Instance()
std::vector< double > Spectrum
Definition: Constants.h:746
TCanvas * canv
enum cmf::det_type DetType_t
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void SelectionHistBinning(cmf::MetaData const &md, cmf::Spectrum &bins)
TH1F * MakeSpectrumHistogram(art::TFileDirectory &tfd, std::string const &baseName, std::string const &baseTitle, cmf::Spectrum const &spectrum)
const std::vector< std::string > cOscParams_LatexScales({" (m)"," (g/cc)"," (#times10^{-5} eV^{2})"," (#times10^{-3} eV^{2})","","",""," (#pi rad)","","","","(eV^{2})"," (#pi rad)","","","","","","","","",""})
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
int fXParamDivs
number of divisions in the space for parameter 1
static cmf::BeamType_t KeyToBeamType(long const &key)
Definition: StaticFuncs.h:135
enum cmf::sel_type SelectionType_t
static bool IsAngleParameter(cmf::OscParm_t const &par)
Definition: StaticFuncs.h:354
static bool IsNuESelected(cmf::SelectionType_t const &sel)
Definition: StaticFuncs.h:380
void Make1DPlot(art::TFileDirectory &tfd, std::string const &chiSqrName, std::map< double, double > const &paramChiSqr, bool xPar)
void Make2DContours(art::TFileDirectory &tfd, std::string const &chiSqrName, cmf::PointMap const &twoDChiSqr, cmf::GridPoint const &bestFit, bool useBestFit=true)
double BinToEnergy(int const &bin, bool allSels=false)
static bool IsNuMuSelected(cmf::SelectionType_t const &sel)
Definition: StaticFuncs.h:366
std::vector< float > fYParamBinEdges
binning for the y parameter
Double_t scale
Definition: plot.C:25
TH2F * MakeGeneric2DHist(art::TFileDirectory &tfd, std::string const &histName)
const std::vector< std::string > cOscParams_Strings_Latex({"L","#rho","#Delta m^{2}_{21}","#Delta m^{2}_{32}","#theta_{12}","#theta_{13}","#theta_{23}","#delta_{CP}","#theta_{14}","#theta_{24}","#theta_{34}","#Delta m^{2}_{41}","#delta_{24}","#epsilon_{ee}","#epsilon_{e#mu}","#epsilon_{e#tau}","#epsilon_{#mu#mu}","#epsilon_{#mu#tau}","#epsilon_{#tau#tau}","#delta_{e#mu}","#delta_{e#tau}","#delta_{#mu#tau}"})
double fNCNormVal
bin normalization in GeV
bool fLogXAxis
are bins on x axis logarithmic
void MakeEnergySpectraFromBins(cmf::Spectrum const &binSpec, cmf::Spectrum const &binCount, std::map< long, TH1D * > &energySpecMap, std::string const &name, art::TFileDirectory &tfd, std::string const &uniqueID)
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
static cmf::DetType_t KeyToDetectorType(long const &key)
Definition: StaticFuncs.h:69
T get(std::string const &key) const
Definition: ParameterSet.h:231
double fContourLevel1Sigma2D
Level value for 1 sigma 2D contour.
cmf::OscParm_t fYParamEnum
enumerated value of the parameter from CovarianceMatrixFit/dataProducts/Constants.h
std::pair< double, double > fXParamExtrema
first is the min, second is the max
const std::vector< std::string > cSelectionType_Strings({"NuESel_AllPID","NuESel_LowPID","NuESel_MidPID","NuESel_HighPID","NuESel_Peripheral","NuMuSel","NuMuSelQ1","NuMuSelQ2","NuMuSelQ3","NuMuSelQ4","NCSel","UnknownSel"})
const std::vector< std::string > cBeamType_Strings({"FHC","RHC","0HC","UnknownBeam"})
const Binning bins
double fContourLevel1Sigma1D
Level value for 1 sigma 1D contour.
void Initialize(fhicl::ParameterSet const &plottingPars)
bool fLogYAxis
are bins on y axis logarithmic
float const & X() const
std::vector< float > fXParamBinEdges
binning for the x parameter
TH2F * Make2DHist(art::TFileDirectory &tfd, std::string const &chiSqrName, cmf::PointMap const &twoDChiSqr)
const Binning binE
Definition: SINCpi0_Bins.h:18
double fNuMuNormVal
bin normalization in GeV
Module to combine a set of results into a single file currently only does one data product type at a ...
Definition: Event.cxx:24
double fNuENormVal
bin normalization in GeV
cmf::OscParm_t fXParamEnum
enumerated value of the parameter from CovarianceMatrixFit/dataProducts/Constants.h
double fContourLevel2Sigma1D
Level value for 2 sigma 1D contour.
std::map< cmf::GridPoint, double > PointMap
Float_t norm
std::map< double, int > const & SelectionHighEdges(cmf::MetaData const &md)
T log10(T number)
Definition: d0nt_math.hpp:120
#define MF_LOG_VERBATIM(category)
#define MF_LOG_DEBUG(id)
enum cmf::beam_type BeamType_t
float const & Y() const
const hit & b
Definition: hits.cxx:21
static PlotUtilities * Instance()
static cmf::SelectionType_t KeyToSelectionType(long const &key)
Definition: StaticFuncs.h:169
const std::vector< double > cOscParams_Scales({1., 1., 1.e5, 1.e3, 1., 1., 1., 1./3.14159, 1., 1., 1., 1., 1./3.14159, 1., 1., 1., 1., 1., 1., 1., 1., 1.})
void MakeCLHeatMaps(art::TFileDirectory &tfd, std::string const &baseName, std::vector< cmf::ChiSqrInfo > const &chiSqrInfoVec)
const std::vector< std::string > cSelectionType_LatexStrings({"#nu_{e} Selection","#nu_{e} Selection - Low PID","#nu_{e} Selection - Mid PID","#nu_{e} Selection - High PID","#nu_{e} Selection - Peripheral","#nu_{#mu} Selection","#nu_{#mu} Selection - Quantile 1","#nu_{#mu} Selection - Quantile 2","#nu_{#mu} Selection - Quantile 3","#nu_{#mu} Selection - Quantile 4","NC Selection","Unknown Selection"})
double fContourLevel3Sigma2D
Level value for 3 sigma 2D contour.
Float_t e
Definition: plot.C:35
enum BeamMode kBlue
void NormalizeBinContents(TH1 *hist, double normToVal)
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
bool fUseTrig
Turn on plotting of trig functions for the angles.
void Make2DHiddenParHistogram(art::TFileDirectory &tfd, std::string const &baseName, cmf::PointMap const &parVals)
int fYParamDivs
number of divisions in the space for parameter 2
const std::vector< std::string > cBeamType_LatexStrings({"FHC","RHC","0HC""Unknown Beam"})
static CovarianceBinUtility * Instance()
const std::vector< std::string > cOscParams_Strings({"L","Rho","Dmsq21","Dmsq32","Th12","Th13","Th23","dCP","Th14","Th24","Th34","Dmsq41","d24","Eps_ee","Eps_emu","Eps_etau","Eps_mumu","Eps_mutau","Eps_tautau","Delta_emu","Delta_etau","Delta_mutau"})
enum BeamMode string