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(std::numeric_limits<double>::max(), std::numeric_limits<double>::lowest()))
39  , fYParamExtrema (std::make_pair(std::numeric_limits<double>::max(), 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" );
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  LOG_VERBATIM("PlotUtilities")
99  << "Plot Utilities: "
100  << "X Paramater: "
102  << " ("
103  << fXParamExtrema.first
104  << ", "
105  << fXParamExtrema.second
106  << ") numDivs: "
107  << fXParamDivs
108  << " Log axis: "
109  << fLogXAxis
110  << " Y Paramater: "
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
132  float logXParamExtremaLow = std::log10(fXParamExtrema.first);
133  float logXParamExtremaHigh = std::log10(fXParamExtrema.second);
134  paramDiv = (logXParamExtremaHigh - logXParamExtremaLow)/float(fXParamDivs);
135  paramStart = (logXParamExtremaLow - 0.5 * paramDiv);
136  paramStop = (logXParamExtremaHigh + 0.5 * paramDiv);
137  for(size_t i = 0; i < fXParamBinEdges.size() - 1; ++i){
138  fXParamBinEdges[i] = std::pow(10,(paramStart + (i * paramDiv)));
139  }
140  fXParamBinEdges.back() = std::pow(10, paramStop);
141  }
142  else {
143 
144  paramDiv = (fXParamExtrema.second - fXParamExtrema.first )/float(fXParamDivs);
145  paramStart = (fXParamExtrema.first - 0.5 * paramDiv) * fXParamScale.first;
146  paramStop = (fXParamExtrema.second + 0.5 * paramDiv) * fXParamScale.first;
147  for (size_t i = 0; i < fXParamBinEdges.size() - 1; ++i){
148  fXParamBinEdges[i] = paramStart + (i * paramDiv * fXParamScale.first);
149  }
150  fXParamBinEdges.back() = paramStop;
151  }
152 
153  // for(auto const& itr : fXParamBinEdges){
154  // LOG_VERBATIM("PlotUtilities")
155  // << "x bin edges: "
156  // << itr;
157  // }
158 
159  if(fLogYAxis)
160  {
161  // convert extrema to log values
162  float logYParamExtremaLow = std::log10(fYParamExtrema.first);
163  float logYParamExtremaHigh = std::log10(fYParamExtrema.second);
164  paramDiv = (logYParamExtremaHigh - logYParamExtremaLow)/float(fYParamDivs);
165  paramStart = (logYParamExtremaLow - 0.5 * paramDiv);
166  paramStop = (logYParamExtremaHigh + 0.5 * paramDiv);
167  for (size_t i = 0; i < fYParamBinEdges.size() - 1; ++i){
168  fYParamBinEdges[i] = std::pow(10,(paramStart + (i * paramDiv)));
169  }
170  fYParamBinEdges.back() = std::pow(10, paramStop);
171  }
172  else {
173  paramDiv = (fYParamExtrema.second - fYParamExtrema.first )/float(fYParamDivs);
174  paramStart = (fYParamExtrema.first - 0.5 * paramDiv) * fYParamScale.first;
175  paramStop = (fYParamExtrema.second + 0.5 * paramDiv) * fYParamScale.first;
176 
177  for (size_t i = 0; i < fYParamBinEdges.size() - 1; ++i){
178  fYParamBinEdges[i] = paramStart + (i * paramDiv * fYParamScale.first);
179  }
180  fYParamBinEdges.back() = paramStop;
181  }
182 
183  // for(auto const& itr : fYParamBinEdges){
184  // LOG_VERBATIM("PlotUtilities")
185  // << "y bin edges: "
186  // << itr;
187  // }
188 
189  }
190 
191  //......................................................................
193  std::string const& chiSqrName,
194  std::map<double, double> const& paramChiSqr,
195  bool xPar)
196  {
197  if(xPar) this->Make1DPlot(tfd,
198  chiSqrName,
199  paramChiSqr,
201  fXParamEnum,
202  fXParamScale);
203  else this->Make1DPlot(tfd,
204  chiSqrName,
205  paramChiSqr,
207  fYParamEnum,
208  fYParamScale);
209  }
210 
211  //......................................................................
213  std::string const& chiSqrName,
214  std::map<double, double> const& paramChiSqr,
215  std::vector<float> const& paramBins,
216  cmf::OscParm_t const& paramEnum,
217  std::pair<double, std::string> const& paramScale)
218  {
219 
220  std::string grName = chiSqrName + cmf::cOscParams_Strings[paramEnum] + "DeltaChiSqr";
221 
222  //make vectors holding the points for the graph
223  std::vector<double> param;
224  std::vector<double> deltaChiSqr;
225 
226  for(auto itr : paramChiSqr){
227  param .emplace_back(itr.first);
228  deltaChiSqr.emplace_back(itr.second);
229  }
230 
231  // make the graph
232  auto *gr = tfd.makeAndRegister<TGraph>(grName.c_str(),
233  grName.c_str(),
234  param.size(),
235  param.data(),
236  deltaChiSqr.data());
237 
238  // make a histogram to set the axis ranges, maximum DeltaChiSqr = 15
239  std::string histName = grName + "Hist";
241  if (paramEnum == cmf::kTh23 && fUseTrig) title = ";sin^{2}(#theta_{23})";
242  else if(paramEnum == cmf::kTh24 && fUseTrig) title = ";sin^{2}(#theta_{24})";
243  else if(paramEnum == cmf::kTh34 && fUseTrig) title = ";sin^{2}(#theta_{34})";
244  title += paramScale.second + ";#Delta#chi^{2}";
245 
246  TH1F *hist = tfd.make<TH1F>(histName.c_str(),
247  title.c_str(),
248  paramBins.size() - 1,
249  paramBins.data());
250  hist->SetMaximum(15.);
251  hist->GetXaxis()->CenterTitle();
252  hist->GetXaxis()->SetDecimals();
253 
254  hist->GetYaxis()->CenterTitle();
255  hist->GetYaxis()->SetDecimals();
256 
257  // Make TLines to show CLs
258  std::vector<Color_t> sigmaColors({kBlue, kRed, kMagenta});
259  std::vector<double> cl_percentiles({fContourLevel1Sigma1D,
262  std::vector<TLine*> sigmaLines(3);
263  for(size_t l = 0; l < sigmaLines.size(); ++l){
264  sigmaLines[l]= new TLine(paramBins.front(), cl_percentiles[l], paramBins.back(), cl_percentiles[l]);
265  sigmaLines[l]->SetLineStyle(2);
266  sigmaLines[l]->SetLineColor(sigmaColors[l]);
267  }
268 
269  LOG_DEBUG("PlotUtilities")
270  << " cl_percentile[0]: "
271  << cl_percentiles[0]
272  << " "
273  << fContourSigmaLabels[0]
274  << " cl_percentile[1]: "
275  << cl_percentiles[1]
276  << " "
277  << fContourSigmaLabels[1]
278  << " cl_percentile[2]: "
279  << cl_percentiles[2]
280  << " "
281  << fContourSigmaLabels[2];
282 
283  // Add a TLegend explain them
284  auto *legLines = new TLegend(0.8, 0.1, 0.9, 0.24);
285  for(size_t l = 0; l < sigmaLines.size(); ++l){
286  legLines->AddEntry(sigmaLines[l], (fContourSigmaLabels[l] + " CL").c_str(), "l");
287  }
288 
289  legLines->SetBorderSize(0);
290  legLines->SetFillStyle(0);
291 
292  // make a Canvas to hold it all
293  std::string canvName = grName + "Canv";
294  auto *canv = tfd.makeAndRegister<TCanvas>(canvName.c_str(),
295  canvName.c_str(),
296  2);
297  canv->cd();
298  hist->SetStats(false);
299  hist->Draw();
300  gr->Draw("lsame");
301  for(size_t l = 0; l < sigmaLines.size(); ++l) sigmaLines[l]->Draw();
302  legLines->Draw();
303 
304  // We're done here
305  canv->ForceUpdate();
306  }
307 
308  //......................................................................
310  std::string const& chiSqrName,
311  cmf::PointMap const& twoDChiSqr,
312  cmf::GridPoint const& bestFit,
313  bool useBestFit)
314  {
315  TH2F *hist = this->Make2DHist(tfd,
316  chiSqrName,
317  twoDChiSqr);
318 
319  std::vector<double> percentiles({fContourLevel1Sigma2D,
322  std::vector<std::string> contourNames({"1sigma",
323  "2sigma",
324  "3sigma"});
325  std::vector<std::string> contourLabels({fContourSigmaLabels[0],
327  fContourSigmaLabels[2]});
328 
329 
330  // check that the minimum for the histogram is less than the lowest contour level,
331  // if not, then get rid of that level
332  if(hist->GetMinimum() > fContourLevel1Sigma2D){
333  percentiles .erase(percentiles.begin() );
334  contourNames .erase(contourNames.begin() );
335  contourLabels.erase(contourLabels.begin());
336  }
337 
338 
339  // LOG_VERBATIM("PlotUtilities")
340  // << "draw contours for ";
341  // for(auto const& itr : percentiles)
342  // LOG_VERBATIM("PlotUtilities")
343  // << itr;
344 
345  hist->SetContour(percentiles.size(), percentiles.data());
346 
347  // This avoids disturbing canvases already
348  // Yes, we do need a canvas for this to work
349  // yes, we do need to call Update() on it
350  // No, ROOT is not very friendly.
351  TCanvas temp_can;
352  temp_can.cd();
353  hist->Draw("contlist");
354  temp_can.Update();
355 
356  auto * plah = (TObjArray *)(gROOT->GetListOfSpecials()->FindObject("contours"));
357  if(!plah) return;
358 
359  std::vector< std::vector<std::unique_ptr<TGraph> > > graphVec(percentiles.size());
360 
361  for(int i = 0; i < plah->GetSize(); ++i){
362  std::vector<std::unique_ptr<TGraph> > temp;
363  auto * list = (TList*) plah->At(i);
364  if(!list) break;
365  if(!(list->First())) break;
366  for(int igraph = 0; igraph < list->GetSize(); ++igraph){
367  if(list->At(igraph)) temp.emplace_back( std::make_unique<TGraph>( *((TGraph*) list->At(igraph))) );
368  }
369  graphVec[i] = std::move(temp);
370  }
371 
372  LOG_DEBUG("PlotAssist")
373  << "Found : "
374  << plah->GetSize()
375  << " sets of contours\n"
376  << "Paramater 1 and Enum 1: "
378  << " "
379  << "Paramater 2 and Enum 2: "
381  << " there are "
382  << graphVec.size()
383  << " vectors of graphs";
384 
385  std::string title = ";";
386  title += hist->GetXaxis()->GetTitle();
387  title += ";";
388  title += hist->GetYaxis()->GetTitle();
389 
390  int p = 0;
391  int g;
392  for(auto const& gritr : graphVec){
393 
394  LOG_DEBUG("PlotUtilities")
395  << "there are "
396  << gritr.size()
397  << " graphs for contour "
398  << p;
399 
400  g = 0;
401  for(auto const& itr : gritr){
402 
403  std::stringstream name;
404  name << std::setprecision(3)
405  << chiSqrName
406  << "ContourGraph"
409  << "_CL"
410  << contourLabels[p]
411  << "_Gr"
412  << g;
413 
414  LOG_DEBUG("PlotUtilities")
415  << "store graph with name "
416  << name.str();
417 
418  tfd.makeAndRegister<TGraph>(name.str().c_str(),
419  title.c_str(),
420  itr->GetN(),
421  itr->GetX(),
422  itr->GetY());
423 
424  ++g;
425  } // end loop over graphs for the percentile
426  ++p;
427  } // end loop over different percentile graphs
428 
429  //LOG_DEBUG << " Found " << graphVec.size() << " standard cl_contour sets\n";
430 
431  // Assign colors and line styles
432  // have to do some checks because we could have either 2 or 3 sets of graphs
433  //-------
434  g = 0;
435  for(auto & grv : graphVec){
436  for(auto & gr : grv){
437  if(percentiles.size() == 3 && g == 0){
438  gr->SetLineColor(kBlue);
439  }
440  if((percentiles.size() == 3 && g == 1) || (percentiles.size() == 2 && g == 0)){
441  gr->SetLineColor(kRed);
442  }
443  if((percentiles.size() == 3 && g == 2) || (percentiles.size() == 2 && g == 1)){
444  gr->SetLineColor(6);
445  }
446  gr->SetLineWidth(2);
447  }
448  ++g;
449  }
450 
451  // Make the 'backdrop' histogram against which the contours
452  // will be drawn
453 
454  LOG_DEBUG("PlotUtilities")
455  << "Make2DContours function: "
456  << "Paramater 1 and Enum 1: "
458  << " "
459  << "Paramater 2 and Enum 2: "
461 
462  std::string histName = (chiSqrName +
463  cmf::cOscParams_Strings[fXParamEnum] +
464  cmf::cOscParams_Strings[fYParamEnum] +
465  "CLContours");
466 
467  TH2F *hBackdrop = tfd.make<TH2F>(histName.c_str(),
468  title.c_str(),
469  hist->GetXaxis()->GetNbins(),
470  hist->GetXaxis()->GetXmin(),
471  hist->GetXaxis()->GetXmax(),
472  hist->GetYaxis()->GetNbins(),
473  hist->GetYaxis()->GetXmin(),
474  hist->GetYaxis()->GetXmax());
475 
476  hBackdrop->GetXaxis()->CenterTitle();
477  hBackdrop->GetXaxis()->SetDecimals();
478 
479  hBackdrop->GetYaxis()->CenterTitle();
480  hBackdrop->GetYaxis()->SetDecimals();
481  hBackdrop->GetYaxis()->SetTitleOffset(1.0);
482 
483  // Make a new canvas, plot the CL contours,
484  //-------
485  std::string canName(chiSqrName + "ContoursCanv");
486  auto *canv = tfd.make<TCanvas>(canName.c_str(),
487  (chiSqrName + " Contours").c_str(),
488  1200,
489  800);
490 
491  hBackdrop->Draw();
492  for(auto const& grv : graphVec){
493  for(auto const& itr : grv){
494  itr->Draw("C SAME");
495  }
496  }
497 
498  // Add a TLegend to describe the CL contours
499  // ROOT is inconsistent, and you must add graphs
500  // to a TLegend by name, rather than by pointer
501  //-------
502  TLegend legContours(0.8, 0.1, 0.9, 0.24);
503  legContours.SetBorderSize(0);
504  legContours.SetFillStyle(0);
505  g = 0;
506  for(auto & grv : graphVec){
507  if(!grv.empty()){
508  grv.front()->SetName(contourNames[g].c_str());
509  legContours.AddEntry(grv.front().get(),
510  (contourLabels[g] + " CL").c_str(),
511  "l");
512  }
513  ++g;
514  }
515  legContours.Draw();
516 
517  if(useBestFit){
518  // And a TLegend to describe the initial guess and
519  // best fit points
520  //-------
521  TLegend legPoints(0.12, 0.11, 0.27, 0.25);
522  legPoints.SetBorderSize(0);
523  legPoints.SetFillStyle(0);
524 
525  std::stringstream par1_legendentry;
526  std::stringstream par2_legendentry;
527 
529  if(fXParamEnum == cmf::kTh23 && fUseTrig) Parameter1Label = "sin^{2}(#theta_{23})";
530  if(fXParamEnum == cmf::kTh24 && fUseTrig) Parameter1Label = "sin^{2}(#theta_{24})";
531 
532  par1_legendentry
533  << Parameter1Label
534  << " = "
535  << std::setprecision(3)
536  << bestFit.X()
537  << " "
538  << fXParamScale.second;
539 
541  if(fYParamEnum == cmf::kTh23 && fUseTrig) Parameter2Label = "sin^{2}(#theta_{23})";
542  if(fYParamEnum == cmf::kTh24 && fUseTrig) Parameter2Label = "sin^{2}(#theta_{24})";
543 
544  par2_legendentry
545  << Parameter2Label
546  << " = "
547  << std::setprecision(3)
548  << bestFit.Y()
549  << " "
550  << fYParamScale.second;
551 
552  // add the initial and best fit point
553 
554  TMarker bestFitMarker(fXParamScale.first * bestFit.X(),
555  fYParamScale.first * bestFit.Y(),
556  8);
557  bestFitMarker.SetMarkerColor(1);
558  bestFitMarker.Draw();
559 
560  legPoints.AddEntry(&bestFitMarker, "Best Fit", "p");
561  legPoints.AddEntry(&bestFitMarker, par1_legendentry.str().data(), "");
562  legPoints.AddEntry(&bestFitMarker, par2_legendentry.str().data(), "");
563 
564  legPoints.Draw();
565  }
566 
567  canv->Update();
568  canv->Write();
569  }
570 
571  //......................................................................
573  std::string const& baseName,
574  std::vector<cmf::ChiSqrInfo> const& chiSqrInfoVec)
575  {
576 
577  // make a 2D histogram for each CL
578  // make a 2D histogram and fill it with the Delta chiSqr values
579  std::string name = (baseName +
582  "1sigmaHeatMap");
583 
584  std::vector<TH2F*> hmVec(3);
585 
586  for(size_t i = 0; i < hmVec.size(); ++i){
587 
588  if(i > 0)
589  name.replace(name.find(std::to_string(i) + "sigma"), 1, std::to_string(i+1));
590 
591  hmVec[i] = this->MakeGeneric2DHist(tfd,
592  name);
593  hmVec[i]->SetTitle((std::to_string(i + 1) + "#sigma").c_str());
594  }
595 
596  // loop over the chiSqrInfo structs to find the delta chiSqr at each point
597  LOG_DEBUG("PlotUtilities")
598  << "there are "
599  << chiSqrInfoVec.size()
600  << " chiSqrInfo objects to loop over";
601 
602  size_t uni = 0;
603  for(auto const& csItr : chiSqrInfoVec){
604  for(auto const& twoDItr : csItr.twoDChiSqr){
605 
606  LOG_DEBUG("PlotUtilities")
607  << "twoDChiSqr for universe "
608  << uni
609  << " "
610  << twoDItr.first
611  << " "
612  << twoDItr.second;
613 
614  if(1.0-TMath::Prob(twoDItr.second, 2) < fContourLevel1Sigma2D)
615  hmVec[0]->Fill(twoDItr.first.X() * fXParamScale.first,
616  twoDItr.first.Y() * fYParamScale.first);
617  if(1.0-TMath::Prob(twoDItr.second, 2) < fContourLevel2Sigma2D)
618  hmVec[1]->Fill(twoDItr.first.X() * fXParamScale.first,
619  twoDItr.first.Y() * fYParamScale.first);
620  if(1.0-TMath::Prob(twoDItr.second, 2) < fContourLevel3Sigma2D)
621  hmVec[2]->Fill(twoDItr.first.X() * fXParamScale.first,
622  twoDItr.first.Y() * fYParamScale.first);
623  }
624  ++uni;
625  } // end loop over ChiSqrInfos
626 
627  // make a canvas to hold these histograms
628  name = baseName +
631  "HeatMapCanv";
632 
633  auto *canv = tfd.make<TCanvas>(name.c_str(),
634  name.c_str(),
635  800,
636  800);
637 
638  canv->Divide(2,2);
639 
640  for(size_t i = 0; i < hmVec.size(); ++i){
641  canv->cd(i + 1);
642  hmVec[i]->Draw("colz");
643  }
644 
645  canv->Update();
646  canv->Write();
647 
648  return;
649  }
650 
651  //......................................................................
652  // Normalize the bin contents to be the equivalent to the desired width
654  double normToVal)
655  {
656  if(!fNormalizeBins) return;
657 
658  double norm;
659  for(int b = 0; b < hist->GetNbinsX(); ++b){
660  norm = normToVal / hist->GetBinWidth(b + 1);
661  hist->SetBinContent(b + 1, hist->GetBinContent(b + 1) * norm);
662  hist->SetBinError(b + 1, hist->GetBinError(b + 1) * norm);
663  } // end loop over bins
664  }
665 
666  //......................................................................
667  void PlotUtilities::MakeEnergySpectraFromBins(std::vector<float> const& binSpec,
668  std::vector<float> const& binCount,
669  std::map<long, TH1D*> & energySpecMap,
670  std::string const& name,
671  art::TFileDirectory & tfd,
672  std::string const& uniqueID)
673  {
675  std::string histTitle;
676  std::vector<double> bins;
677  TH1D *hist;
678 
679  // loop over the keys we are using to determine
680  // the offset for each selection
681  long offset;
683 
684  double binE;
685  double scale;
686 
687  for(auto const& dbsItr : cmf::SelectionUtility::Instance()->SelectionsToUse()){
688  md.detector = dbsItr.Detector();
689  md.period = (dbsItr.BeamType() == cmf::kFHC) ? 1 : 4;
690  scale = (dbsItr.Detector() == cmf::kNEARDET) ? 1.e-3 : 1.;
691 
692  for(auto const& itr : dbsItr.Keys()){
694 
695  // reverse engineer the metadata
697 
698  bins.clear();
700 
701  // check that we have a sane number of bins
702  if(offset + cmf::CovarianceBinUtility::Instance()->SelectionHighEdges(md).size() > binSpec.size())
703  throw cet::exception("CovarianceFitHelper")
704  << "Filling data energy spectra appears to want more bins"
705  << " than are in the data bin spectrum: "
707  << " vs "
708  << binSpec.size();
709 
710  histName = (name +
711  uniqueID +
713  md.DetectorString() +
715 
716  histTitle = (cmf::cBeamType_LatexStrings[md.BeamType()] + " " +
717  md.DetectorString() + " " +
719  ";Energy (GeV);");
720  histTitle += (scale < 1) ? "#times10^{3} Events" : "Events";
721 
722  if(md.isMC) histTitle.insert(0, "MC ");
723 
724  LOG_DEBUG("CovarianceFitHelper")
725  << "total Energy spectrum hist is "
726  << histName
727  << " "
728  << itr;
729 
730 // for(auto const& binitr : bins)
731 // LOG_VERBATIM("CovarianceFitHelper")
732 // << "bin edge: "
733 // << binitr;
734 
735  if(energySpecMap.count(itr) < 1)
736  hist = energySpecMap.emplace(itr, tfd.make<TH1D>(histName.c_str(),
737  histTitle.c_str(),
738  bins.size() - 1,
739  &bins[0])).first->second;
740  else
741  hist = energySpecMap.find(itr)->second;
742 
743  hist->GetXaxis()->CenterTitle();
744  hist->GetYaxis()->CenterTitle();
745  if(name.find("MC") != std::string::npos){
746  hist->SetMarkerColor(kRed);
747  hist->SetLineColor(kRed);
748  }
749 
750  for(size_t b = offset; b < offset + cmf::CovarianceBinUtility::Instance()->SelectionHighEdges(md).size(); ++b){
752  hist->Fill(binE, binSpec[b] * scale);
753 
754  // the uncertainty for each bin is the fractional uncertainty for the total count in the
755  // bin (1 / sqrt(binCount[b]) multiplied by the total weight in the bin (binSpec[b])
756  // only do it if the bin count > 0
757  if(binCount[b] > 0.)
758  hist->SetBinError(hist->FindBin(binE),
759  scale * (binSpec[b] / std::sqrt(binCount[b])));
760  }// end filling the spectrum
761 
762  LOG_VERBATIM("PlotUtilities")
763  << "histogram "
764  << hist->GetName()
765  << " has "
766  << hist->Integral()
767  << " selected events";
768 
769  // now normalize the bin contents if desired to the desired width
770  if(md.IsNuESelected())
771  this->NormalizeBinContents(hist, fNuENormVal);
772  else if(md.IsNuMuSelected())
773  this->NormalizeBinContents(hist, fNuMuNormVal);
774  else if (md.IsNCSelected())
775  this->NormalizeBinContents(hist, fNCNormVal);
776 
777  } // end loop over keys to use
778  } // end loop over selections being used in the analysis
779  }
780 
781  //......................................................................
783  std::string const& baseName,
784  cmf::Spectrum const& spectrum)
785  {
786  if(spectrum.size() < 1){
787  LOG_VERBATIM("PlotUtilities")
788  << "spectrum size for "
789  << baseName
790  << " is 0, are you sure you filled it?";
791  return;
792  }
793 
794  auto *spectrumHist = tfd.make<TH1F>(baseName.c_str(),
795  ";Bin;Events",
796  spectrum.size(),
797  0,
798  spectrum.size());
799 
800  for(size_t b = 0; b < spectrum.size(); ++b) spectrumHist->Fill(b, spectrum[b]);
801 
802  }
803 
804  //......................................................................
805  // The returned histogram has the axes and axis ranges defined by fXParam(2)
806  // only for now
808  std::string const& histName)
809  {
810  // Delta m^2 is always on the y axis for numu
811  // angle is always on y axis for nue contours
812 
815 
816  if(fXParamEnum == cmf::kTh23 && fUseTrig) XaxisTitle = "sin^{2}(#theta_{23})";
817  if(fXParamEnum == cmf::kTh24 && fUseTrig) XaxisTitle = "sin^{2}(#theta_{24})";
818  if(fYParamEnum == cmf::kTh23 && fUseTrig) YaxisTitle = "sin^{2}(#theta_{23})";
819  if(fYParamEnum == cmf::kTh24 && fUseTrig) YaxisTitle = "sin^{2}(#theta_{24})";
820 
821  std::string title = (";" + XaxisTitle +
822  ";" + YaxisTitle);
823 
824 
825  LOG_DEBUG("ContourMaker")
826  << "Creating 2D histogram with range "
827  << fXParamBinEdges.front() << " " << fXParamBinEdges.back() << " "
828  << fYParamBinEdges.front() << " " << fYParamBinEdges.back();
829 
830  auto* hist2D = tfd.make<TH2F>(histName.c_str(),
831  title.c_str(),
832  fXParamDivs+1,
833  fXParamBinEdges.data(),
834  fYParamDivs+1,
835  fYParamBinEdges.data());
836 
837  hist2D->GetXaxis()->CenterTitle();
838  hist2D->GetXaxis()->SetDecimals();
839  hist2D->GetYaxis()->CenterTitle();
840  hist2D->GetYaxis()->SetDecimals();
841 
842  return hist2D;
843  }
844 
845  //......................................................................
847  std::string const& chiSqrName,
848  cmf::PointMap const& twoDChiSqr)
849  {
850  // make a 2D histogram and fill it with the Delta chiSqr values
851  std::string histName = (chiSqrName +
854  "DeltaChiSqr");
855 
856  auto* hist_deltachisq = this->MakeGeneric2DHist(tfd,
857  histName);
858 
859  histName = (chiSqrName +
862  "DeltaChiSqr_Prob");
863 
864  auto* hist = this->MakeGeneric2DHist(tfd,
865  histName);
866 
867  double dcsq;
868  for(auto const& itr : twoDChiSqr){
869 
870  LOG_DEBUG("ContourMaker")
871  << itr.first.X()
872  << " "
873  << itr.first.Y()
874  << " "
875  << itr.second;
876 
877  dcsq = (itr.second < 0.1) ? 0.1 : itr.second;
878  hist_deltachisq->Fill(itr.first.X() * fXParamScale.first,
879  itr.first.Y() * fYParamScale.first,
880  dcsq);
881 
882  hist->Fill(itr.first.X() * fXParamScale.first,
883  itr.first.Y() * fYParamScale.first,
884  1.0-TMath::Prob(itr.second, 2));
885  }
886 
887 // for(int b = 0; b < hist_deltachisq->GetXaxis()->GetNbins(); ++ b)
888 // LOG_VERBATIM("CMFContourMaker")
889 // << "x axis "
890 // << b + 1
891 // << " "
892 // << hist_deltachisq->GetXaxis()->GetBinCenter(b + 1);
893 //
894 // for(int b = 0; b < hist_deltachisq->GetYaxis()->GetNbins(); ++ b)
895 // LOG_VERBATIM("CMFContourMaker")
896 // << "y axis "
897 // << b + 1
898 // << " "
899 // << hist_deltachisq->GetYaxis()->GetBinCenter(b + 1);
900 
901  return hist;
902  }
903 
904  //......................................................................
906  std::string const& baseName,
907  cmf::PointMap const& parVals)
908  {
909  if(parVals.size() < 1){
910  LOG_VERBATIM("PlotUtilities")
911  << "point map size for "
912  << baseName
913  << " is 0, are you sure you filled it?";
914  return;
915  }
916 
917  auto* parValHist = this->MakeGeneric2DHist(tfd,
918  baseName);
919 
920  for(auto const& itr : parVals){
921  parValHist->Fill(itr.first.X() * fXParamScale.first,
922  itr.first.Y() * fYParamScale.first,
923  itr.second);
924  }
925  }
926 
927 } // end cmf namespace
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
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
int KeyToOffset(long const &key, bool allSels=false)
void SelectionHistBinning(cmf::MetaData const &md, std::vector< double > &bins)
cmf::OscParm_t fFreeParamEnum
enumerated value of the oscillation parameter that floats freely
void MakeEnergySpectraFromBins(std::vector< float > const &binSpec, std::vector< float > const &binCount, std::map< long, TH1D * > &energySpecMap, std::string const &name, art::TFileDirectory &tfd, std::string const &uniqueID)
const std::string cSelectionType_Strings[12]
Definition: Constants.h:77
double fContourLevel2Sigma2D
Level value for 2 sigma 2D contour.
cmf::DetType_t detector
Definition: Structs.h:114
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
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
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()
void MakeSpectrumHistogram(art::TFileDirectory &tfd, std::string const &baseName, cmf::Spectrum const &spectrum)
TCanvas * canv
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
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
const std::string cOscParams_LatexScales[kNumOscParams]
Definition: Constants.h:233
const double cOscParams_Scales[kNumOscParams]
Definition: Constants.h:212
void Make1DPlot(art::TFileDirectory &tfd, std::string const &chiSqrName, std::map< double, double > const &paramChiSqr, bool xPar)
const std::string cOscParams_Strings[kNumOscParams]
Definition: Constants.h:253
void Make2DContours(art::TFileDirectory &tfd, std::string const &chiSqrName, cmf::PointMap const &twoDChiSqr, cmf::GridPoint const &bestFit, bool useBestFit=true)
std::string DetectorString() const
Definition: Structs.h:151
double BinToEnergy(int const &bin, bool allSels=false)
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)
cmf::SelectionType_t selectionType
Definition: Structs.h:116
double fNCNormVal
bin normalization in GeV
bool fLogXAxis
are bins on x axis logarithmic
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
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
bool IsNuMuSelected() const
Definition: Structs.cxx:219
std::pair< double, double > fXParamExtrema
first is the min, second is the max
const Binning bins
double fContourLevel1Sigma1D
Level value for 1 sigma 1D contour.
void Initialize(fhicl::ParameterSet const &plottingPars)
T * makeAndRegister(char const *name, char const *title, ARGS...args) const
std::vector< float > Spectrum
Definition: Constants.h:570
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 std::string cOscParams_Strings_Latex[kNumOscParams]
Definition: Constants.h:275
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
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
cmf::OscParm_t fXParamEnum
enumerated value of the parameter from CovarianceMatrixFit/dataProducts/Constants.h
T * make(ARGS...args) const
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
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:170
void MakeCLHeatMaps(art::TFileDirectory &tfd, std::string const &baseName, std::vector< cmf::ChiSqrInfo > const &chiSqrInfoVec)
void NormalizeBinContents(TH1D *hist, double normToVal)
long period
Definition: Structs.h:118
double fContourLevel3Sigma2D
Level value for 3 sigma 2D contour.
cmf::BeamType_t BeamType() const
Definition: Structs.cxx:178
#define LOG_VERBATIM(category)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const std::string cBeamType_LatexStrings[4]
Definition: Constants.h:38
Float_t e
Definition: plot.C:35
const std::string cBeamType_Strings[4]
Definition: Constants.h:33
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
const std::string cSelectionType_LatexStrings[12]
Definition: Constants.h:90
bool fUseTrig
Turn on plotting of trig functions for the angles.
cons_index_list< index_uni, nil_index_list > uni
void Make2DHiddenParHistogram(art::TFileDirectory &tfd, std::string const &baseName, cmf::PointMap const &parVals)
int fYParamDivs
number of divisions in the space for parameter 2
static CovarianceBinUtility * Instance()
bool IsNuESelected() const
Definition: Structs.cxx:225
bool IsNCSelected() const
Definition: Structs.cxx:231