ContourMaker_plugin.cc
Go to the documentation of this file.
1 //
2 // ContourMaker_plugin.cpp
3 //
4 // Created by Brian Rebel on 3/23/16.
5 // Modified by Prabhjot Singh (prabhjot@fnal.gov) on Jan 25, 2017
6 //
7 
8 #include <string>
9 #include <set>
10 #include <memory>
11 #include <sstream>
12 
13 #include "TH1.h"
14 #include "TH2.h"
15 #include "TGraph.h"
16 #include "TCanvas.h"
17 #include "TMath.h"
18 #include "TLegend.h"
19 #include "TMarker.h"
20 
21 // Framework includes
29 #include "fhiclcpp/ParameterSet.h"
31 
32 // NOvA includes
38 
39 namespace fnex {
40 
41  typedef std::map<std::pair<double, double>, std::set<fnex::PointResult> > PointResMap;
42  typedef std::map<std::pair<double, double>, double > PointMap;
43 
45  public:
46  explicit ContourMaker(fhicl::ParameterSet const& pset);
47  virtual ~ContourMaker();
48 
49  void readResults (art::Results const& r);
50  void writeResults(art::Results & r);
51  void reconfigure (fhicl::ParameterSet const& p);
52  void clear();
53  void endJob();
54 
55  private:
56 
57  bool readSingleResults(art::Results const& r,
58  std::vector<fnex::PointResult> & prvec,
59  std::vector<fnex::InputPoint> & ipvec);
60 
61  bool readVectorResults(art::Results const& r,
62  std::vector<fnex::PointResult> & prvec,
63  std::vector<fnex::InputPoint> & ipvec);
64 
65  void AddPointToMap(fnex::InputPoint const& ip,
66  fnex::PointResult const& pr);
67 
68  void MakeDeltaChiSqrMaps(std::map<double, double> & param1ChiSqr,
69  std::map<double, double> & param2ChiSqr,
70  PointMap & twoDChiSqr);
71 
72  void Make1DPlot(std::map<double, double> const& paramChiSqr,
73  std::pair<double, double> const& paramExtrema,
74  fnex::OscParm_t const& paramEnum,
75  std::pair<double, std::string> const& paramScale) const;
76 
77  TH2F* Make2DHist(PointMap const& twoDChiSqr) const;
78  void Make2DContours(PointMap const& twoDChiSqr) const;
79  void StoreContourGraphs(std::vector<double> const& percentiles,
80  std::vector< std::vector< std::unique_ptr<TGraph> > > & graphs) const;
81  void Make2DSystematicHists(PointResult pr, TH2F* h) const;
82 
84  PointResMap fPointResMap; ///< Initial point passed to fitter, typedef above
85  std::string fPointLabel; ///< label of plugin storing the points
86  std::string fAltPointLabel; ///< reverse the order of the parameters in the label
87  fnex::OscParm_t fParam1Enum; ///< enumerated value of the parameter from FNEX/dataProducts/Constants.h
88  fnex::OscParm_t fParam2Enum; ///< enumerated value of the parameter from FNEX/dataProducts/Constants.h
89  std::pair<double, double> fParam1Extrema; ///< first is the min, second is the max
90  std::pair<double, double> fParam2Extrema; ///< first is the min, second is the max
91  std::pair<double, std::string> fParam1Scale; ///< scale the parameter if necessary ie is x10^{-3}, string is latex friendly
92  std::pair<double, std::string> fParam2Scale; ///< scale the parameter if necessary ie is x10^{-3}, string is latex friendly
93  fnex::PointResult fBestFit; ///< result with the minimum \chi^2 in the space
94  fnex::InputPoint fBestGuess; ///< initial guess that gave fBestFit
95  double fDeltaChiSqScale; ///< Rescale DeltaChiSq values (for sensitivity contours)
96  bool fUseTrig; ///< Turn on plotting of trig functions for the angles
97  int fNDivisions; ///< Force the number of divisions. If fNDivisions < 0, use the number of files
98  double fContourLevel1Sigma1D; ///< Level value for 1 sigma 1D contour
99  double fContourLevel2Sigma1D; ///< Level value for 2 sigma 1D contour
100  double fContourLevel3Sigma1D; ///< Level value for 3 sigma 1D contour
101  double fContourLevel1Sigma2D; ///< Level value for 1 sigma 2D contour
102  double fContourLevel2Sigma2D; ///< Level value for 2 sigma 2D contour
103  double fContourLevel3Sigma2D; ///< Level value for 3 sigma 2D contour
104  std::string fContour1SigmaLabel; ///< Legend name for 1 sigma contour
105  std::string fContour2SigmaLabel; ///< Legend name for 2 sigma contour
106  std::string fContour3SigmaLabel; ///< Legend name for 3 sigma contour
107 
108  };
109 
110  //......................................................................
114  , fParam1Extrema(std::make_pair(std::numeric_limits<double>::max(), std::numeric_limits<double>::lowest()))
115  , fParam2Extrema(std::make_pair(std::numeric_limits<double>::max(), std::numeric_limits<double>::lowest()))
116  , fParam1Scale(std::make_pair(1., ""))
117  , fParam2Scale(std::make_pair(1., ""))
118  {
119  this->reconfigure(pset);
120 
121  return;
122  }
123 
124  //......................................................................
126  {
127 
128  }
129 
130  //......................................................................
132  {
133 
134  fPointLabel = pset.get<std::string >("PointLabel" );
135  fDeltaChiSqScale = pset.get<double >("DeltaChiSqScale", 1.0 );
136  fUseTrig = pset.get<bool >("UseTrig" );
137  fNDivisions = pset.get<int >("NDivisions", -1 );
138  fContourLevel1Sigma1D = pset.get<double >("ContourLevel1Sigma1D", 1 );
139  fContourLevel2Sigma1D = pset.get<double >("ContourLevel2Sigma1D", 4 );
140  fContourLevel3Sigma1D = pset.get<double >("ContourLevel3Sigma1D", 9 );
141  fContourLevel1Sigma2D = pset.get<double >("ContourLevel1Sigma2D", 2.3 );
142  fContourLevel2Sigma2D = pset.get<double >("ContourLevel2Sigma2D", 6.18 );
143  fContourLevel3Sigma2D = pset.get<double >("ContourLevel3Sigma2D", 11.83 );
144  fContour1SigmaLabel = pset.get<std::string >("Contour1SigmaLabel", "1#sigma");
145  fContour2SigmaLabel = pset.get<std::string >("Contour2SigmaLabel", "2#sigma");
146  fContour3SigmaLabel = pset.get<std::string >("Contour3SigmaLabel", "3#sigma");
147 
148  auto hashPos = fPointLabel.find("#" );
149  auto Dmsq32Pos = fPointLabel.find("Dmsq32");
150  auto Th13Pos = fPointLabel.find("Th13" );
151  auto Th23Pos = fPointLabel.find("Th23" );
152  auto dCPPos = fPointLabel.find("dCP" );
153 
154  LOG_DEBUG("ContourMaker")
155  << "\nLabel : " << fPointLabel
156  << "\n\tHashPos : " << hashPos
157  << "\n\tDmsq32Pos : " << Dmsq32Pos
158  << "\n\tTh13Pos : " << Th13Pos
159  << "\n\tTh23Pos : " << Th23Pos
160  << "\n\tDCP : " << dCPPos;
161 
162  if(Th13Pos == Th23Pos && Dmsq32Pos == Th13Pos)
163  throw cet::exception("ContourMaker")
164  << "all possible parameters have the same position in the label"
165  << " - that cannot be good";
166 
167  std::string subPar1;
168  std::string subPar2;
169 
170  if(hashPos == std::string::npos)
171  throw cet::exception("ContourMaker")
172  << "cannot find the # in the input label, something is wrong";
173 
174  if(Dmsq32Pos != std::string::npos){
175  if(Dmsq32Pos > hashPos+2){
177  fParam2Scale = std::make_pair(1.e3, " (#times10^{-3} eV^{2})");
178  subPar2 = fPointLabel.substr(Dmsq32Pos, 6);
179  }
180  else{
182  fParam1Scale = std::make_pair(1.e3, " (#times10^{-3} eV^{2})");
183  subPar1 = fPointLabel.substr(Dmsq32Pos, 6);
184  }
185  }
186 
187  if(Th13Pos != std::string::npos){
188  if(Th13Pos > hashPos+2){
190  fParam2Scale = std::make_pair(1./3.14159, " (#pi rad)");
191  subPar2 = fPointLabel.substr(Th13Pos, 4);
192  }
193  else{
195  fParam1Scale = std::make_pair(1./3.14159, " (#pi rad)");
196  subPar1 = fPointLabel.substr(Th13Pos, 4);
197  }
198  }
199 
200  if(Th23Pos != std::string::npos){
201  if(Th23Pos > hashPos+2){
203  subPar2 = fPointLabel.substr(Th23Pos, 4);
204  if(!fUseTrig) fParam2Scale = std::make_pair(1./3.14159, " (#pi rad)");
205  }
206  else{
208  subPar1 = fPointLabel.substr(Th23Pos, 4);
209  if(!fUseTrig) fParam1Scale = std::make_pair(1./3.14159, " (#pi rad)");
210  }
211  }
212  if(dCPPos != std::string::npos){
213  if(dCPPos > hashPos+2){
215  fParam2Scale = std::make_pair(1./3.14159, " (#pi rad)");
216  subPar2 = fPointLabel.substr(dCPPos, 3);
217  }
218  else{
220  fParam1Scale = std::make_pair(1./3.14159, " (#pi rad)");
221  subPar1 = fPointLabel.substr(dCPPos, 3);
222  }
223  }
224  fAltPointLabel = fPointLabel.substr(0, hashPos+1) + subPar2 + subPar1;
225 
226  //for the numu analysis
227  // we always want Delta m^2 to be on the y axis, and the 2D drawing
228  // method always puts parameter 1 on the x axis...so make sure
229  // Delta m^2 is in the parameter 2 spot
230  //
231  // for the nue analysis
232  // we always want Th23 to be on the y axis, dCP on the x axis,
233  // and the 2D drawing method always puts parameter 1 on the x axis...
234  // so make sure Th23 is in the parameter 2 spot and dCP in
235  // the parameter 1 spot
236  //
237  //For any other analysis either use the paramter on the basis
238  //of their position in fPointLabel
239  //or create a special set for special parameters
240 
241  LOG_DEBUG("ContourMaker")
242  << "Reconfigure function: "
243  << "Before Swapping: "
244  << "Paramater 1 and Enum 1: "
246  << " "
247  << "Paramater 2 and Enum 2: "
249 
253  }
254 
258  }
259 
260  LOG_DEBUG("ContourMaker")
261  << "Reconfigure function: "
262  << "After Swapping: "
263  << "Paramater 1 and Enum 1: "
264  << fnex::cOscParams_Strings[fParam1Enum]
265  << " "
266  << "Paramater 2 and Enum 2: "
267  << fnex::cOscParams_Strings[fParam2Enum];
268 
269  LOG_DEBUG("ContourMaker")
270  << "Alternate label is "
271  << fAltPointLabel;
272 
273  //Over right and decide parameter1(2) and scale1(2) here
274  return;
275  }
276 
277  //......................................................................
278  // Method to clear out the collections of data products after the
279  // writeResults method is called.
281  {
282  return;
283  }
284 
285  //......................................................................
287  std::vector<fnex::PointResult> & prvec,
288  std::vector<fnex::InputPoint> & ipvec)
289  {
290  prvec.resize(1);
291  ipvec.resize(1);
292 
293  // get the FitPoints from the output files
295  r.getByLabel(fPointLabel, pr);
296 
297  // get the InputPoints from the output files - need these when
298  // attempting to make FC contours when there are many results for each
299  // InputPoint
301  r.getByLabel(fPointLabel, ip);
302 
303  // check that the handles are valid, add the points to the map
304  // and move along
305  if( !ip.isValid() || !pr.isValid() ){
306 
307  // see if the alternative label works
308  r.getByLabel(fAltPointLabel, pr);
309  r.getByLabel(fAltPointLabel, ip);
310 
311  if( !ip.isValid() || !pr.isValid() ) {
312  LOG_WARNING("ContourMaker")
313  << "unable to find any results in either\n"
314  << fPointLabel
315  << "\n or "
316  << fAltPointLabel;
317 
318  // no valid handles, clear the vectors.
319  prvec.clear();
320  ipvec.clear();
321 
322  return false;
323  }
324 
325  }
326 
327  // got here, so we have valid handles
328  prvec[0] = *pr;
329  ipvec[0] = *ip;
330 
331  return true;
332  }
333 
334  //......................................................................
336  std::vector <fnex::PointResult> & prvec,
337  std::vector <fnex::InputPoint> & ipvec)
338  {
339  prvec.clear();
340  ipvec.clear();
341 
342  // individual point handles are not valid, what about vectors of them?
344  r.getByLabel(fPointLabel, prhand);
345 
347  r.getByLabel(fPointLabel, iphand);
348 
349  if( !prhand.isValid() || !iphand.isValid() ){
350 
351  // see if the alternative label works
352  r.getByLabel(fAltPointLabel, prhand);
353  r.getByLabel(fAltPointLabel, iphand);
354 
355  if( !prhand.isValid() || !iphand.isValid() ){
356 
357  LOG_WARNING("ContourMaker")
358  << "PointResult handle isValid: "
359  << prhand.isValid()
360  << " InputPoint handle isValid: "
361  << iphand.isValid();
362 
363  return false;
364  } // end if the alternate label still does not make valid handles.
365  } // end if the handles are not valid
366 
367  // check that the vectors are the same size, we assume that the entries
368  // in each spot of the vectors go together
369  if( prhand->size() != iphand->size() )
370  throw cet::exception("ContourMaker")
371  << "PointResult and InputPoint vectors have different sizes, "
372  << "not sure we can do anything with these data so bail";
373 
374  // fill the vectors
375  prvec.resize(prhand->size());
376  ipvec.resize(iphand->size());
377 
378  for(size_t p = 0; p < prhand->size(); ++p){
379  prvec[p] = (*prhand)[p];
380  ipvec[p] = (*iphand)[p];
381  }
382 
383  return true;
384  }
385 
386  //......................................................................
388  {
389 
390  std::vector<fnex::PointResult> prvec;
391  std::vector<fnex::InputPoint> ipvec;
392 
393  // check to see if we have single points stored in the
394  // results, and if not look for vectors of points.
395  // if none of those either, give up
396  if( !this->readSingleResults(r, prvec, ipvec) ){
397 
398  if( !this->readVectorResults(r, prvec, ipvec) ){
399  LOG_WARNING("ContourMaker")
400  << "No valid handles to either single fnex::PointResult, fnex::InputPoint "
401  << "objects or vectors of them, skip this file";
402 
403  return;
404 
405  } // end if no vector point results
406 
407  } // end if no single point results
408 
409  // we found the results, so loop over the vectors and add them to the
410  // map
413  int fitStatus;
414  for(size_t p = 0; p < prvec.size(); ++p){
415 
416  pr = prvec[p];
417  ip = ipvec[p];
418 
419  fitStatus = pr.FitStatus();
420  LOG_DEBUG("ContourMaker")
421  << "Fit status of the fit is "
422  << fitStatus;
423 
424  if(fitStatus > -1 && fitStatus < 4){
425  this->AddPointToMap(ip, pr);
426  }
427  else{
428  // fit did not converge, so
429  // create a new PointResult with a very large chiSqr
430  // I am assuming that a value of -1 means the fit did not
431  // converge, so substitute it for a stupidly large chi^2
432  LOG_VERBATIM("ContourMaker")
433  << "switching chi^2 for "
434  << ip
435  << " "
436  << pr.ChiSqr()
437  << " and status "
438  << fitStatus
439  << " to "
441 
442  pr = fnex::PointResult(std::numeric_limits<double>::max(),
443  fitStatus,
445  pr.SystematicPulls());
446  this->AddPointToMap(ip, pr);
447  }
448 
449  } // end loop over points found in file
450 
451  return;
452  }
453 
454  //......................................................................
456  fnex::PointResult const& pr)
457  {
458  LOG_DEBUG("ContourMaker")
459  << " input point "
460  << ip
461  << " fit point "
462  << pr;
463 
464  // Add the point to the map, the first entry in the set is certain
465  // to be the lowest chi^2 value for the input point because that is
466  // how the set and < operator of PointResult work
467  // the key is the pair of the input parameter values
468  auto paramDetMap = pr.ParameterSpaceLocation();
469 
470  // for(auto itr : paramDetMap)
471  // LOG_VERBATIM("ContourMaker")
472  // << itr.first.first
473  // << " "
474  // << itr.first.second
475  // << " "
476  // << itr.second;
477 
478  LOG_DEBUG("ContourMaker")
479  << " have the paramMap with "
480  << paramDetMap.size()
481  << " entries";
482 
485 
486  if(paramDetMap.count(pd1) < 1 ||
487  paramDetMap.count(pd2) < 1)
488  throw cet::exception("ContourMaker")
489  << "Cannot find one of the parameters in the map: "
490  << "\n" << pd1.first
491  << " " << paramDetMap.count(pd1)
492  << "\n" << pd2.first
493  << " " << paramDetMap.count(pd2);
494 
495  std::pair<double, double> inputPoint = std::make_pair(fParam1Scale.first * paramDetMap.find(pd1)->second,
496  fParam2Scale.first * paramDetMap.find(pd2)->second);
497  LOG_VERBATIM("ContourMaker")
498  << "AddPointToMap function: "
499  << " Paramater 1 and Enum 1: "
501  << " "
502  << paramDetMap.find(pd1)->first.first
503  << " "
504  << paramDetMap.find(pd1)->first.second
505  << " "
506  << paramDetMap.find(pd1)->second
507  << " "
508  << paramDetMap.find(pd2)->first.first
509  << " "
510  << paramDetMap.find(pd2)->first.second
511  << " "
512  << paramDetMap.find(pd2)->second
513  << " "
514  << inputPoint.second;
515 
516  //If parameter 1 is Th23 then use sin(Th23)*sin(Th23) form
517  if(fUseTrig && fParam1Enum == fnex::kTh23){
518  inputPoint.first = std::sin(inputPoint.first) * std::sin(inputPoint.first);
519  }
520  //OR If parameter 2 is Th23 then use sin(Th23)*sin(Th23) form
521  if(fUseTrig && fParam2Enum == fnex::kTh23){
522  inputPoint.second = std::sin(inputPoint.second) * std::sin(inputPoint.second);
523  }
524 
525  fPointResMap[inputPoint].emplace(pr);
526 
527  // see if the chi^2 for this PointResult is lower than the minimum
528  if( pr.ChiSqr() < fBestFit.ChiSqr() ){
529  fBestFit = pr;
530  fBestGuess = ip;
531  }
532 
533  LOG_DEBUG("ContourMaker")
534  << "addding "
535  << inputPoint.first
536  << " "
537  << inputPoint.second
538  << " to map with chi^2 = "
539  << pr.ChiSqr()
540  << " size = "
541  << fPointResMap.size()
542  << " min chi^2 = "
543  << fBestFit.ChiSqr()
544  << " and fit status: "
545  << pr.FitStatus();
546 
547  return;
548  }
549 
550  //......................................................................
552  {
553  // not writing anything to the file for this plugin at this time
554 
555  return;
556  } // end writeResults
557 
558  //......................................................................
559  void ContourMaker::MakeDeltaChiSqrMaps(std::map<double, double> & param1ChiSqr,
560  std::map<double, double> & param2ChiSqr,
561  PointMap & twoDChiSqr)
562  {
563  // for each of the parameters from the grid need to project
564  // onto the 1D chi^2 space, will have to loop over all the points
565  // at each value for the parameter of interest and grab the lowest
566  // chi^2 for that value in the other parameter.
567 
568  // also make a delta chiSqr map for the whole space
569 
570  // Notes: itr.first is an input point, itr.second is a point result
571 
572  double param1Val = 0.;
573  double param2Val = 0.;
574  double minChiSqr = fBestFit.ChiSqr();
575 
576  LOG_DEBUG("ContourMaker")
577  << "There are "
578  << fPointResMap.size()
579  << " points in the map";
580 
581  std::pair<double, double> point;
582  for(auto itr : fPointResMap){
583 
584  param1Val = itr.first.first;
585  param2Val = itr.first.second;
586 
587  LOG_DEBUG("ContourMaker")
588  << param1Val
589  << " "
590  << param2Val
591  << " "
592  << itr.second.begin()->ChiSqr()
593  << " "
594  << itr.second.size();
595 
596  // find the edges of the parameter space represented by the input points
597  fParam1Extrema.first = std::min(fParam1Extrema.first, param1Val);
598  fParam1Extrema.second = std::max(fParam1Extrema.second, param1Val);
599  fParam2Extrema.first = std::min(fParam2Extrema.first, param2Val);
600  fParam2Extrema.second = std::max(fParam2Extrema.second, param2Val);
601 
602  // Have we seen this param1 yet? If not, set its param1ChiSqr = DeltaChiSq. wrt to the best fit chi2
603  // If so, see if this instance shows a smaller DeltaChiSq
604  // (would be at different param2 for this param1, hence the possibility of a different value)
605 
606  if(param1ChiSqr.count(param1Val) > 0)
607  param1ChiSqr[param1Val] = std::min(itr.second.begin()->ChiSqr() - minChiSqr,
608  param1ChiSqr.find(param1Val)->second);
609  else
610  param1ChiSqr[param1Val] = itr.second.begin()->ChiSqr() - minChiSqr;
611 
612  // Do the same for param2 (find minimimum DeltaChiSq for each value of Param2, by comparing
613  // to all InputPoint/PointResult pairs at Param2 that we see
614 
615  if(param2ChiSqr.count(param2Val) > 0)
616  param2ChiSqr[param2Val] = std::min(itr.second.begin()->ChiSqr() - minChiSqr,
617  param2ChiSqr.find(param2Val)->second);
618  else
619  param2ChiSqr[param2Val] = itr.second.begin()->ChiSqr() - minChiSqr;
620 
621  // At each point, update the 2D DeltaChiSq spot. It doesn't depend on having a particular
622  // value of param1 or param2, so we can just make a simple check down here at the bottom
623  twoDChiSqr[itr.first] = itr.second.begin()->ChiSqr() - minChiSqr;
624 
625  LOG_DEBUG("ContourMaker")
626  << "There are "
627  << twoDChiSqr.size()
628  << " points in the delta chisqr map";
629 
630  }
631 
632  //Rescale POT (for sensitivity contours)
633  LOG_DEBUG("ContourMaker")
634  << "fDeltaChiSqScale: "
635  << fDeltaChiSqScale;
636 
637  for(auto & itr: param1ChiSqr){ itr.second *= fDeltaChiSqScale; }
638  for(auto & itr: param2ChiSqr){ itr.second *= fDeltaChiSqScale; }
639  for(auto & itr: twoDChiSqr ){ itr.second *= fDeltaChiSqScale; }
640 
641  return;
642  }
643 
644  //......................................................................
645  void ContourMaker::Make1DPlot(std::map<double, double> const& paramChiSqr,
646  std::pair<double, double> const& paramExtrema,
647  fnex::OscParm_t const& paramEnum,
648  std::pair<double, std::string> const& paramScale) const
649  {
650 
651  std::string grName = fnex::cOscParams_Strings[paramEnum] + "DeltaChiSqr";
652 
653  //make vectors holding the points for the graph
654  std::vector<double> param;
655  std::vector<double> deltaChiSqr;
656 
657  for(auto itr : paramChiSqr){
658  param .push_back(itr.first);
659  deltaChiSqr.push_back(itr.second);
660  }
661 
662  // make the graph
663  TGraph *gr = fTFS->makeAndRegister<TGraph>(grName.c_str(),
664  grName.c_str(),
665  param.size(),
666  &param[0],
667  &deltaChiSqr[0]);
668 
669  // make a histogram to set the axis ranges, maximum DeltaChiSqr = 15
670  std::string histName = grName + "Hist";
672  if(paramEnum == fnex::kTh23 && fUseTrig) title = ";sin^{2}(#theta_{23})";
673  title += paramScale.second + ";#Delta#chi^{2}";
674 
675  TH1F *hist = fTFS->make<TH1F>(histName.c_str(),
676  title.c_str(),
677  2,
678  paramExtrema.first,
679  paramExtrema.second);
680  hist->SetMaximum(15.);
681  hist->GetXaxis()->CenterTitle();
682  hist->GetXaxis()->SetDecimals();
683 
684  hist->GetYaxis()->CenterTitle();
685  hist->GetYaxis()->SetDecimals();
686 
687  // Make TLines to show CLs
688  std::vector<double> cl_percentiles(3);
689  cl_percentiles[0] = fContourLevel1Sigma1D;
690  cl_percentiles[1] = fContourLevel2Sigma1D;
691  cl_percentiles[2] = fContourLevel3Sigma1D;
692  TLine * line_1sigma = new TLine(paramExtrema.first, cl_percentiles[0], paramExtrema.second, cl_percentiles[0]);
693  TLine * line_2sigma = new TLine(paramExtrema.first, cl_percentiles[1], paramExtrema.second, cl_percentiles[1]);
694  TLine * line_3sigma = new TLine(paramExtrema.first, cl_percentiles[2], paramExtrema.second, cl_percentiles[2]);
695 
696  LOG_DEBUG("ContourMaker")
697  << " cl_percentile[0]: "
698  << cl_percentiles[0]
699  << " "
701  << " cl_percentile[1]: "
702  << cl_percentiles[1]
703  << " "
705  << " cl_percentile[2]: "
706  << cl_percentiles[2]
707  << " "
709 
710 
711  line_1sigma->SetLineStyle(2);
712  line_1sigma->SetLineColor(kBlue);
713  line_2sigma->SetLineStyle(1);
714  line_2sigma->SetLineColor(kRed);
715  line_3sigma->SetLineStyle(1);
716  line_3sigma->SetLineColor(6);
717 
718  // Add a TLegend explain them
719  TLegend * legLines = new TLegend(0.8, 0.1, 0.9, 0.24);
720  legLines->AddEntry(line_1sigma, (fContour1SigmaLabel + " CL").c_str(), "l");
721  legLines->AddEntry(line_2sigma, (fContour2SigmaLabel + " CL").c_str(), "l");
722  legLines->AddEntry(line_3sigma, (fContour3SigmaLabel + " CL").c_str(), "l");
723 
724  legLines->SetBorderSize(0);
725  legLines->SetFillStyle(0);
726 
727  // make a Canvas to hold it all
728  std::string canvName = grName + "Canv";
729  TCanvas *canv = fTFS->makeAndRegister<TCanvas>(canvName.c_str(),
730  canvName.c_str(),
731  2);
732  canv->cd();
733  hist->SetStats(false);
734  hist->Draw();
735  gr->Draw("lsame");
736 
737  line_1sigma->Draw();
738  line_2sigma->Draw();
739  line_3sigma->Draw();
740 
741  legLines->Draw();
742 
743  // We're done here
744  canv->ForceUpdate();
745 
746  return;
747  }
748 
749  //......................................................................
750  TH2F* ContourMaker::Make2DHist(PointMap const& twoDChiSqr) const
751  {
752  // first make a histogram with the fit status code so we can tell
753  // if any fits failed
754 
757  "FitStatus");
758 
759  // Delta m^2 is always on the y axis for numu
760  // angle is always on y axis for nue contours
761 
764 
766  XaxisTitle = "sin^{2}(#theta_{23})";
768  YaxisTitle = "sin^{2}(#theta_{23})";
769 
770  std::string title = (";" + XaxisTitle +
771  ";" + YaxisTitle);
772 
773  // assume there are the same number of divisions for each parameter
774  // KM Aug 29 2016: changed ceil --> floor in following line so that adding a best
775  // fit point to the mix won't add an extra bin in each dimension
776  int numDivisions = (fNDivisions < 0) ? std::floor(std::sqrt(twoDChiSqr.size())) : fNDivisions;
777  LOG_DEBUG("ContourMaker") << "Using " << numDivisions << " divisions";
778 
779  double param1Div = ( fParam1Extrema.second - fParam1Extrema.first )/double(numDivisions);
780  double param2Div = ( fParam2Extrema.second - fParam2Extrema.first )/double(numDivisions);
781  double param1Start = fParam1Extrema.first - 0.5 * param1Div;
782  double param2Start = fParam2Extrema.first - 0.5 * param2Div;
783  double param1Stop = fParam1Extrema.second + 0.5 * param1Div;
784  double param2Stop = fParam2Extrema.second + 0.5 * param2Div;
785 
786 
787  LOG_DEBUG("ContourMaker")
788  << "param 1: "
789  << param1Start
790  << " "
791  << param1Stop
792  << " "
793  << param1Div
794  << " param 2: "
795  << param2Start
796  << " "
797  << param2Stop
798  << " "
799  << param2Div;
800 
801  TH2F *statusHist = fTFS->make<TH2F>(histName.c_str(),
802  title.c_str(),
803  numDivisions + 1,
804  param1Start,
805  param1Stop,
806  numDivisions + 1,
807  param2Start,
808  param2Stop);
809 
810  for(auto itr : fPointResMap)
811  statusHist->Fill(itr.first.first,
812  itr.first.second,
813  itr.second.begin()->FitStatus());
814 
815 
816  for(auto itr : fPointResMap){
817  this->Make2DSystematicHists(*(itr.second.begin()), statusHist);
818  break;
819  }
820 
821  // make a 2D histogram and fill it with the chiSqr values
824  "ChiSqr");
825  TH2F *hist_chisq = fTFS->make<TH2F>(histName.c_str(),
826  title.c_str(),
827  numDivisions + 1,
828  param1Start,
829  param1Stop,
830  numDivisions + 1,
831  param2Start,
832  param2Stop);
833 
834  for(auto itr : fPointResMap)
835  hist_chisq->Fill(itr.first.first,
836  itr.first.second,
837  itr.second.begin()->ChiSqr());
838 
839  // make a 2D histogram and fill it with the Delta chiSqr values
842  "DeltaChiSqr");
843 
844  LOG_VERBATIM("ContourMaker")
845  << "Creating 2D histogram with range "
846  << param1Start << " " << param1Stop << " "
847  << param2Start << " " << param2Stop << " "
848  << numDivisions;
849 
850  TH2F *hist_deltachisq = fTFS->make<TH2F>(histName.c_str(),
851  title.c_str(),
852  numDivisions + 1,
853  param1Start,
854  param1Stop,
855  numDivisions + 1,
856  param2Start,
857  param2Stop);
858 
859  for(auto itr : twoDChiSqr){
860 
861  LOG_DEBUG("ContourMaker")
862  << itr.first.first
863  << " "
864  << itr.first.second
865  << " "
866  << itr.second;
867 
868  hist_deltachisq->Fill(itr.first.first,
869  itr.first.second,
870  itr.second);
871 
872  }
873 
874  hist_deltachisq->GetXaxis()->CenterTitle();
875  hist_deltachisq->GetXaxis()->SetDecimals();
876 
877  hist_deltachisq->GetYaxis()->CenterTitle();
878  hist_deltachisq->GetYaxis()->SetDecimals();
879 
882  "DeltaChiSqr_Prob");
883 
884  TH2F *hist = fTFS->make<TH2F>(histName.c_str(),
885  title.c_str(),
886  numDivisions + 1,
887  param1Start,
888  param1Stop,
889  numDivisions + 1,
890  param2Start,
891  param2Stop);
892 
893  // If num divisions are not the same for each parameter, or if there are multiple entries per grid
894  // point, results shown will not really display what we want. Will have to change this behavior
895  // when we add a true '2D best fit' point to the lineup. For now will accomplish this by having two
896  // histograms (one counts entries, one counts sum of all entries, per grid point) and then dividing
897  // one by the other to get an average.
898  for(auto itr : twoDChiSqr){
899 
900  LOG_DEBUG("ContourMaker")
901  << itr.first.first
902  << " "
903  << itr.first.second
904  << " "
905  << itr.second;
906 
907  hist->Fill(itr.first.first,
908  itr.first.second,
909  1.0-TMath::Prob(itr.second, 2));
910  }
911 
912  hist->GetXaxis()->CenterTitle();
913  hist->GetXaxis()->SetDecimals();
914 
915  hist->GetYaxis()->CenterTitle();
916  hist->GetYaxis()->SetDecimals();
917 
918  return hist;
919  }
920 
921  //......................................................................
922  void ContourMaker::Make2DContours(PointMap const& twoDChiSqr) const
923  {
924  TH2F *hist = this->Make2DHist(twoDChiSqr);
925 
926  // Make the 'backdrop' histogram against which the contours
927  // will be drawn
928 
929  LOG_DEBUG("ContourMaker")
930  << "Make2DContours function: "
931  << "Paramater 1 and Enum 1: "
933  << " "
934  << "Paramater 2 and Enum 2: "
936 
937  std::string histName = (fnex::cOscParams_Strings[fParam1Enum] +
938  fnex::cOscParams_Strings[fParam2Enum] +
939  "CLContours");
940  std::string title = ";"; title += hist->GetXaxis()->GetTitle();
941  title += ";"; title += hist->GetYaxis()->GetTitle();
942 
943 
944  TH2F *hBackdrop = fTFS->make<TH2F>(histName.c_str(),
945  title.c_str(),
946  hist->GetXaxis()->GetNbins(),
947  hist->GetXaxis()->GetXmin(),
948  hist->GetXaxis()->GetXmax(),
949  hist->GetYaxis()->GetNbins(),
950  hist->GetYaxis()->GetXmin(),
951  hist->GetYaxis()->GetXmax());
952 
953  hBackdrop->GetXaxis()->CenterTitle();
954  hBackdrop->GetXaxis()->SetDecimals();
955 
956  hBackdrop->GetYaxis()->CenterTitle();
957  hBackdrop->GetYaxis()->SetDecimals();
958  hBackdrop->GetYaxis()->SetTitleOffset(1.0);
959 
960  // Grab the contours and color them appropriately
961  // There can be many TGraphs associated with a single cl_percentile value;
962  // this is the reason for the nested vectors structure (one vector of TGraphs
963  // for each cl_percentile value)
964  std::vector<double> cl_percentiles(3);
965  cl_percentiles[0] = fContourLevel1Sigma2D;
966  cl_percentiles[1] = fContourLevel2Sigma2D;
967  cl_percentiles[2] = fContourLevel3Sigma2D;
968 
969  LOG_DEBUG("ContourMaker")
970  << " 2D percentiles: "
971  << cl_percentiles[0]
972  << " "
973  << cl_percentiles[1]
974  << " "
975  << cl_percentiles[2];
976 
977  std::vector< std::vector<std::unique_ptr<TGraph> > > cl_contours;
978  fnex::MakeGraphs(hist, cl_percentiles, cl_contours);
979 
980  this->StoreContourGraphs(cl_percentiles, cl_contours);
981 
982  //LOG_DEBUG << " Found " << cl_contours.size() << " standard cl_contour sets\n";
983 
984  // Check to see how many contours we have to draw.
985  // (May 2 2016: I don't remember what condition this line is meant
986  // to correct. Maybe it can be removed?)
987  //-------
988  int num_to_draw_graphs = std::min(cl_contours.size(), cl_percentiles.size());
989 
990  // Assign their colors and line styles
991  //-------
992  for(int ialpha = 0; ialpha < num_to_draw_graphs; ++ialpha){
993  for(unsigned int igraph = 0; igraph < cl_contours[ialpha].size(); ++igraph){
994  if(ialpha == 0){
995  cl_contours[ialpha][igraph]->SetLineStyle(2);
996  cl_contours[ialpha][igraph]->SetLineColor(kBlue);
997  }
998  if(ialpha == 1){
999  cl_contours[ialpha][igraph]->SetLineStyle(1);
1000  cl_contours[ialpha][igraph]->SetLineColor(kRed);
1001  }
1002  if(ialpha == 2){
1003  cl_contours[ialpha][igraph]->SetLineStyle(1);
1004  cl_contours[ialpha][igraph]->SetLineColor(6);
1005  }
1006  cl_contours[ialpha][igraph]->SetLineWidth(2);
1007  }
1008  }
1009 
1010  // add the contour TGraphs to the TFileService outpu
1011 
1012  // Make a new canvas, plot the CL contours,
1013  //-------
1014  // FNEXUniqID * id = FNEXUniqID::getInstance();
1015  TString canName = "GaussianContoursCanv";
1016  TCanvas * can_gaus = fTFS->make<TCanvas>(canName.Data(), "Gaussian Contours", 1200, 800);
1017  can_gaus->cd(1);
1018  TString origSurfaceTitle = hBackdrop->GetTitle();
1019  hBackdrop->SetTitle("Gaussian Probability Contours");
1020  hBackdrop->SetStats(false);
1021  hBackdrop->Draw();
1022  hBackdrop->SetTitle(origSurfaceTitle.Data());
1023  for(int ialpha = 0; ialpha < num_to_draw_graphs; ++ialpha){
1024  for(size_t igraph = 0; igraph < cl_contours[ialpha].size(); ++igraph){
1025  cl_contours[ialpha][igraph]->Draw("C SAME");
1026  }
1027  }
1028 
1029 
1030  // add the initial and best fit point,
1031  //-------
1032  auto bestGuessParamDetMap = fBestGuess.ParameterSpaceLocation();
1033  auto bestFitParamDetMap = fBestFit.ParameterSpaceLocation();
1034  fnex::ParameterDet pd1(fnex::cOscParams_Strings[fParam1Enum], novadaq::cnv::kFARDET);
1035  fnex::ParameterDet pd2(fnex::cOscParams_Strings[fParam2Enum], novadaq::cnv::kFARDET);
1036 
1037  std::pair<double, double> bestGuessPoint
1038  = std::make_pair(fParam1Scale.first * bestGuessParamDetMap.find(pd1)->second,
1039  fParam2Scale.first * bestGuessParamDetMap.find(pd2)->second);
1040 
1041  std::pair<double, double> bestFitPoint
1042  = std::make_pair(fParam1Scale.first * bestFitParamDetMap.find(pd1)->second,
1043  fParam2Scale.first * bestFitParamDetMap.find(pd2)->second);
1044 
1045  if(fUseTrig && fParam1Enum == fnex::kTh23){
1046  bestFitPoint.first = std::sin(bestFitPoint.first) * std::sin(bestFitPoint.first);
1047  bestGuessPoint.first = std::sin(bestGuessPoint.first) * std::sin(bestGuessPoint.first);
1048  }
1049  if(fUseTrig && fParam2Enum == fnex::kTh23){
1050  bestFitPoint.second = std::sin(bestFitPoint.second) * std::sin(bestFitPoint.second);
1051  bestGuessPoint.second = std::sin(bestGuessPoint.second) * std::sin(bestGuessPoint.second);
1052  }
1053 
1054  TMarker bestFitMarker(bestFitPoint.first, bestFitPoint.second, 8);
1055  bestFitMarker.SetMarkerColor(1);
1056  bestFitMarker.Draw();
1057 
1058  // Add a TLegend to describe the CL contours
1059  //-------
1060  TLegend legContours(0.8, 0.1, 0.9, 0.24);
1061  if(cl_contours[0].size() > 0) cl_contours[0][0]->SetName("cl68"); // ROOT is inconsistent, and you must add graphs
1062  if(cl_contours[1].size() > 0) cl_contours[1][0]->SetName("cl90"); // to a TLegend by name, rather than by pointer
1063  if(cl_contours[2].size() > 0) cl_contours[2][0]->SetName("cl95");
1064  for(int ialpha = 0; ialpha< num_to_draw_graphs; ++ialpha){
1065  if(cl_contours[0].size() > 0 && ialpha == 0) legContours.AddEntry( "cl68", (fContour1SigmaLabel+ " CL") .c_str(), "l" );
1066  if(cl_contours[1].size() > 0 && ialpha == 1) legContours.AddEntry( "cl90", (fContour2SigmaLabel+ " CL") .c_str(), "l" );
1067  if(cl_contours[2].size() > 0 && ialpha == 2) legContours.AddEntry( "cl95", (fContour3SigmaLabel+ " CL") .c_str(), "l" );
1068  }
1069 
1070  legContours.SetBorderSize(0);
1071  legContours.SetFillStyle(0);
1072  legContours.Draw();
1073 
1074  // And a TLegend to describe the initial guess and
1075  // best fit points
1076  //-------
1077  TLegend legPoints(0.12, 0.11, 0.27, 0.25);
1078  legPoints.SetBorderSize(0);
1079  legPoints.SetFillStyle(0);
1080 
1081  std::stringstream par1_legendentry;
1082  std::stringstream par2_legendentry;
1083 
1084  std::string Parameter1Label(fnex::cOscParams_Strings_Latex[fParam1Enum]);
1085  if(fParam1Enum == fnex::kTh23 && fUseTrig)
1086  Parameter1Label = "sin^{2}(#theta_{23})";
1087 
1088  par1_legendentry
1089  << Parameter1Label
1090  << " = "
1091  << std::setprecision(3)
1092  << bestFitPoint.first
1093  << " "
1094  << fParam1Scale.second;
1095 
1096  std::string Parameter2Label(fnex::cOscParams_Strings_Latex[fParam2Enum]);
1097  if(fParam2Enum == fnex::kTh23 && fUseTrig)
1098  Parameter2Label = "sin^{2}(#theta_{23})";
1099 
1100  par2_legendentry
1101  << Parameter2Label
1102  << " = "
1103  << std::setprecision(3)
1104  << bestFitPoint.second
1105  << " "
1106  << fParam2Scale.second;
1107 
1108  legPoints.AddEntry(&bestFitMarker, "Best Fit", "p");
1109  legPoints.AddEntry(&bestFitMarker, par1_legendentry.str().data(), "" );
1110  legPoints.AddEntry(&bestFitMarker, par2_legendentry.str().data(), "" );
1111 
1112  legPoints.Draw();
1113 
1114  // Add canvas to the root output file
1115  //-------
1116  can_gaus->Write();
1117 
1118  // We're done here
1119  //-------
1120 
1121  // Then want to add CL lines to the 1D raster scan plots, as well
1122 
1123  return;
1124  }
1125 
1126 
1127  //......................................................................
1128  void ContourMaker::StoreContourGraphs(std::vector<double> const& percentiles,
1129  std::vector< std::vector< std::unique_ptr<TGraph> > > & graphs) const
1130  {
1131  LOG_DEBUG("ContourMaker")
1132  << "StoreContourGraphs function: "
1133  << "Paramater 1 and Enum 1: "
1135  << " "
1136  << "Paramater 2 and Enum 2: "
1138  << " there are "
1139  << graphs.size()
1140  << " vectors of graphs";
1141 
1142  //for the numu contours, Param1 is the angle
1143  //for the nue contours, Param2 is the angle
1146 
1147  if(fParam1Enum == fnex::kTh23 && fUseTrig)
1148  XaxisTitle = "sin^{2}(#theta_{23})";
1149 
1150  if(fParam2Enum == fnex::kTh23 && fUseTrig)
1151  YaxisTitle = "sin^{2}(#theta_{23})";
1152 
1153  std::string title = (";" + XaxisTitle +
1154  ";" + YaxisTitle);
1155 
1156  for(size_t p = 0; p < graphs.size(); ++p){
1157  LOG_DEBUG("ContourMaker")
1158  << "there are "
1159  << graphs[p].size()
1160  << " graphs for contour "
1161  << p;
1162 
1163  for(size_t g = 0; g < graphs[p].size(); ++g){
1164 
1165  std::stringstream name;
1166  name
1167  << std::setprecision(3)
1168  << "ContourGraph"
1171  << "_CL"
1172  << percentiles[p]
1173  << "_Gr"
1174  << g;
1175 
1176  LOG_DEBUG("ContourMaker")
1177  << "store graph with name "
1178  << name.str();
1179 
1180  fTFS->makeAndRegister<TGraph>(name.str().c_str(),
1181  title.c_str(),
1182  graphs[p][g]->GetN(),
1183  graphs[p][g]->GetX(),
1184  graphs[p][g]->GetY());
1185 
1186  } // end loop over graphs for the percentile
1187  } // end loop over different percentile graphs
1188 
1189  return;
1190  }
1191 
1192  //......................................................................
1194  {
1195  // get the maps for the 1 and 2D delta ChiSqr values
1196  std::map<double, double> param1ChiSqr;
1197  std::map<double, double> param2ChiSqr;
1198  PointMap twoDChiSqr;
1199 
1200  this->MakeDeltaChiSqrMaps(param1ChiSqr, param2ChiSqr, twoDChiSqr);
1201 
1202  // if there are no entries in the twoDChiSqr map, we are in trouble
1203  if(twoDChiSqr.size() < 1)
1204  throw cet::exception("ContourMaker")
1205  << "no entries in the map needed to make the contours";
1206 
1207  this->Make1DPlot(param1ChiSqr, fParam1Extrema, fParam1Enum, fParam1Scale);
1208  this->Make1DPlot(param2ChiSqr, fParam2Extrema, fParam2Enum, fParam2Scale);
1209  this->Make2DContours(twoDChiSqr);
1210 
1211  return;
1212  }
1213 
1214  //--------------------------------------------------------------------------------------
1215  //Make 2D histograms of the Systematic values in the oscillation space
1217 
1218  //subdirectory to store systematic histograms
1219  art::TFileDirectory itrDir = fTFS->mkdir("2DSystematicHists");
1220 
1221  //subdirectory to store 2D oscilaltion histograms
1222  art::TFileDirectory oscDir = fTFS->mkdir("2DOscillationHists");
1223 
1224  //get the list of systematic parameters
1225  auto sysParDets = pr.SystematicPulls(novadaq::cnv::kFARDET);
1226 
1227  //get the list of oscillation parameters
1228  auto oscParDets = pr.ParameterSpaceLocation(novadaq::cnv::kFARDET);
1229 
1230  //make a collection of histograms for systematics
1231  std::map < std::string, TH2F * > histMap;
1232 
1233  std::string XaxisTitle = h->GetXaxis()->GetTitle();
1234  std::string YaxisTitle = h->GetYaxis()->GetTitle();
1235  std::string title = (";" + XaxisTitle +
1236  ";" + YaxisTitle);
1237 
1238  for(auto itr : sysParDets){
1239 
1240  LOG_VERBATIM("ContourMaker")
1241  << " Make2DSystematicHists function: Systematic: "
1242  << itr.first;
1243  histMap[itr.first] = itrDir.make<TH2F>((itr.first).c_str(),
1244  title.c_str(),
1245  h->GetXaxis()->GetNbins(),
1246  h->GetXaxis()->GetXmin(),
1247  h->GetXaxis()->GetXmax(),
1248  h->GetYaxis()->GetNbins(),
1249  h->GetYaxis()->GetXmin(),
1250  h->GetYaxis()->GetXmax());
1251  histMap[itr.first]->SetTitle((itr.first).c_str());
1252  }//end of loop over systematics
1253 
1254  // Make the 2D histograms for the oscillation parameters
1255  for(auto itr : oscParDets){
1256 
1257  LOG_VERBATIM("ContourMaker")
1258  << " Make2DSystematicHists function: Oscillation Parameters: "
1259  << itr.first;
1260  histMap[itr.first] = oscDir.make<TH2F>((itr.first).c_str(),
1261  title.c_str(),
1262  h->GetXaxis()->GetNbins(),
1263  h->GetXaxis()->GetXmin(),
1264  h->GetXaxis()->GetXmax(),
1265  h->GetYaxis()->GetNbins(),
1266  h->GetYaxis()->GetXmin(),
1267  h->GetYaxis()->GetXmax());
1268  histMap[itr.first]->SetTitle((itr.first).c_str());
1269  }//end of loop over systematics
1270 
1271  //Now fill the histograms
1272  //loop over all points in the space
1273  for(auto itr_prm : fPointResMap){
1274 
1275  //loop over all systematic at a point and fill the histogram
1276  for(auto itr_pr: itr_prm.second.begin()->SystematicPulls(novadaq::cnv::kFARDET)){
1277 
1278  histMap[itr_pr.first]->Fill(itr_prm.first.first,
1279  itr_prm.first.second,
1280  itr_pr.second);
1281  }//end of loop over all systematic at a point and fill the histogram
1282  //loop over all systematic at a point and fill the histogram
1283  for(auto itr_osc_pr: itr_prm.second.begin()->ParameterSpaceLocation(novadaq::cnv::kFARDET)){
1284 
1285  histMap[itr_osc_pr.first]->Fill(itr_prm.first.first,
1286  itr_prm.first.second,
1287  itr_osc_pr.second);
1288  }//end of loop over all systematic at a point and fill the histogram
1289  }//loop over all points in the space
1290  }//end of Make2DSystematicHists function
1291 
1292 } // end fnex namespace
1293 
1294 namespace fnex{
1295 
1297 
1298 }
bool readVectorResults(art::Results const &r, std::vector< fnex::PointResult > &prvec, std::vector< fnex::InputPoint > &ipvec)
T max(const caf::Proxy< T > &a, T b)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
const XML_Char * name
Definition: expat.h:151
std::pair< double, double > fParam2Extrema
first is the min, second is the max
int fNDivisions
Force the number of divisions. If fNDivisions < 0, use the number of files.
art::ServiceHandle< art::TFileService > fTFS
TFileService.
int const & FitStatus() const
Definition: PointResult.h:71
std::map< std::pair< double, double >, std::set< fnex::PointResult > > PointResMap
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
TFileDirectory mkdir(std::string const &dir, std::string const &descr="")
double fContourLevel1Sigma1D
Level value for 1 sigma 1D contour.
TCanvas * canv
ContourMaker(fhicl::ParameterSet const &pset)
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void Make2DSystematicHists(PointResult pr, TH2F *h) const
void StoreContourGraphs(std::vector< double > const &percentiles, std::vector< std::vector< std::unique_ptr< TGraph > > > &graphs) const
TString ip
Definition: loadincs.C:5
ParameterDetMap const SystematicPulls() const
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
Create a list of fnex::Events to be used in fits.
double fDeltaChiSqScale
Rescale DeltaChiSq values (for sensitivity contours)
std::string fPointLabel
label of plugin storing the points
fnex::OscParm_t fParam2Enum
enumerated value of the parameter from FNEX/dataProducts/Constants.h
void MakeGraphs(TH2F *in_hist, std::vector< double > vContours, std::vector< std::vector< std::unique_ptr< TGraph > > > &vGraphVectors)
Definition: PlotAssist.cxx:15
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
bool isValid() const
Definition: Handle.h:189
TH2F * Make2DHist(PointMap const &twoDChiSqr) const
Far Detector at Ash River, MN.
std::string fContour3SigmaLabel
Legend name for 3 sigma contour.
PointResMap fPointResMap
Initial point passed to fitter, typedef above.
void MakeDeltaChiSqrMaps(std::map< double, double > &param1ChiSqr, std::map< double, double > &param2ChiSqr, PointMap &twoDChiSqr)
std::pair< double, std::string > fParam2Scale
scale the parameter if necessary ie is x10^{-3}, string is latex friendly
ParameterMap const ParameterSpaceLocation(novadaq::cnv::DetId const &det) const
Definition: InputPoint.cxx:330
#define DEFINE_ART_RESULTS_PLUGIN(klass)
T get(std::string const &key) const
Definition: ParameterSet.h:231
void AddPointToMap(fnex::InputPoint const &ip, fnex::PointResult const &pr)
void reconfigure(fhicl::ParameterSet const &p)
bool readSingleResults(art::Results const &r, std::vector< fnex::PointResult > &prvec, std::vector< fnex::InputPoint > &ipvec)
bool fUseTrig
Turn on plotting of trig functions for the angles.
double fContourLevel1Sigma2D
Level value for 1 sigma 2D contour.
double const & ChiSqr() const
Definition: PointResult.h:70
std::pair< double, std::string > fParam1Scale
scale the parameter if necessary ie is x10^{-3}, string is latex friendly
T * makeAndRegister(char const *name, char const *title, ARGS...args) const
ParameterDetMap const ParameterSpaceLocation() const
Definition: PointResult.cxx:82
double fContourLevel2Sigma2D
Level value for 2 sigma 2D contour.
fnex::InputPoint fBestGuess
initial guess that gave fBestFit
void Make1DPlot(std::map< double, double > const &paramChiSqr, std::pair< double, double > const &paramExtrema, fnex::OscParm_t const &paramEnum, std::pair< double, std::string > const &paramScale) const
#define LOG_WARNING(category)
enum fnex::osc_params OscParm_t
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
fnex::OscParm_t fParam1Enum
enumerated value of the parameter from FNEX/dataProducts/Constants.h
fnex::PointResult fBestFit
result with the minimum ^2 in the space
T * make(ARGS...args) const
T sin(T number)
Definition: d0nt_math.hpp:132
std::map< std::pair< double, double >, double > PointMap
std::string fContour1SigmaLabel
Legend name for 1 sigma contour.
void readResults(art::Results const &r)
std::string fAltPointLabel
reverse the order of the parameters in the label
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::pair< double, double > fParam1Extrema
first is the min, second is the max
void Make2DContours(PointMap const &twoDChiSqr) const
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
const std::string cOscParams_Strings[kNumOscParams]
Definition: Constants.h:196
std::pair< std::string, novadaq::cnv::DetId > ParameterDet
Definition: Constants.h:38
double fContourLevel3Sigma2D
Level value for 3 sigma 2D contour.
TRandom3 r(0)
std::string fContour2SigmaLabel
Legend name for 2 sigma contour.
const std::string cOscParams_Strings_Latex[kNumOscParams]
Definition: Constants.h:205
#define LOG_VERBATIM(category)
T min(const caf::Proxy< T > &a, T b)
double fContourLevel2Sigma1D
Level value for 2 sigma 1D contour.
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
void writeResults(art::Results &r)
double fContourLevel3Sigma1D
Level value for 3 sigma 1D contour.