StarPlot_plugin.cc
Go to the documentation of this file.
1 //
2 // StarPlot_plugin.cpp
3 //
4 // Created by Prabhjot Singh (prabhjot) on 05/04/2017
5 //
6 // Purpose: Make star plot and overlay star plot on the contours.
7 //
8 
9 #include <string>
10 #include <set>
11 #include <memory>
12 #include <sstream>
13 
14 #include "TH1.h"
15 #include "TH2.h"
16 #include "TGraph.h"
17 #include "TCanvas.h"
18 #include "TMath.h"
19 #include "TLegend.h"
20 #include "TMarker.h"
21 #include "TPaveText.h"
22 #include "TFile.h"
23 
24 // Framework includes
32 #include "fhiclcpp/ParameterSet.h"
34 
35 // NOvA includes
41 
42 namespace fnex {
43 
44  typedef std::map<std::pair<double,double>, std::set<fnex::PointResult> > PointResMap;
45  typedef std::map<std::pair<double,double>, double > PointMap;
46 
47  class StarPlot : public art::ResultsProducer {
48  public:
49  explicit StarPlot(fhicl::ParameterSet const& pset);
50  virtual ~StarPlot();
51 
52  void readResults (art::Results const& r);
53  void writeResults(art::Results & r);
54  void reconfigure (fhicl::ParameterSet const& p);
55  void clear();
56  void endJob();
57 
58  private:
59 
60  bool readSingleResults(art::Results const& r,
61  std::vector<fnex::PointResult> & prvec,
62  std::vector<fnex::InputPoint> & ipvec,
63  std::vector<fnex::InputPoint> & ipvec_fake);
64 
65  void AddPointToMap(fnex::InputPoint const& ip,
66  fnex::PointResult const& pr,
67  fnex::InputPoint const& ip_fake);
68 
69  void MakeStarPlot();
70 
71  void MakeTable();
72 
73  float ParameterDiff(float, float);
74 
75  void Legends(TLine* l);
76 
78  bool fUseTrig; ///< Turn on plotting of trig functions for the angles
79  PointResMap fPointResMap; ///< Initial point passed to fitter, typedef above
80  std::string fPointLabel; ///< label of plugin storing the points
81  std::string fAltPointLabel; ///< reverse the order of the parameters in the label
82  fnex::OscParm_t fParam1Enum; ///< enumerated value of the parameter from FNEX/dataProducts/Constants.h
83  fnex::OscParm_t fParam2Enum; ///< enumerated value of the parameter from FNEX/dataProducts/Constants.h
84  std::pair<double, double> fParam1Extrema; ///< first is the min, second is the max
85  std::pair<double, double> fParam2Extrema; ///< first is the min, second is the max
86  std::pair<double, std::string> fParam1Scale; ///< scale the parameter if necessary ie is x10^{-3}, string is latex friendly
87  std::pair<double, std::string> fParam2Scale; ///< scale the parameter if necessary ie is x10^{-3}, string is latex friendly
88  float fMinParX; ///< Minimum X range in Star Plot
89  float fMaxParX; ///< Maximum X range in Star Plot
90  float fMinParY; ///< Minimum Y range in Star Plot
91  float fMaxParY; ///< Maximum Y range in Star Plot
92  float fDisplayMinX; ///< Display minimum X range in Star Plot
93  float fDisplayMaxX; ///< Display minimum X range in Star Plot
94  float fDisplayMinY; ///< Display minimum X range in Star Plot
95  float fDisplayMaxY; ///< Display minimum X range in Star Plot
96  std::vector<float> fparX_vec; ///< A vector that stores values to be plotted on X axis;
97  std::vector<float> fparY_vec; ///< A vector that stores values to be plotted on Y axis.
98  std::vector<float> fplusminus_shift; ///< A vector that stores the plus or minus values of the systematic shift in the fake input point
99  std::vector<TString> fsyst_vec; ///< A vector that stores systematics names
100  // std::vector<std::string> fsyst_vec;
101  std::map<TString, TString> fmap_syst_name; ///< map to store human readable systematic names.
102  std::map<TString, int> fmap_syst_color; ///< map to store line colors for same systematics.
103  bool fOneSystOnly; ///< only looks for one shifter or weighter per input file
104  bool fOverlayContour; ///< true if we want to overlay the results on a contour
105  std::string fContourFile; ///< filename of contour to overlay plots on
106 
107  };
111  , fParam1Extrema(std::make_pair(std::numeric_limits<double>::max(), std::numeric_limits<double>::min()))
112  , fParam2Extrema(std::make_pair(std::numeric_limits<double>::max(), std::numeric_limits<double>::min()))
113  , fParam1Scale(std::make_pair(1., ""))
114  , fParam2Scale(std::make_pair(1., ""))
115  {
116  this->reconfigure(pset);
117 
118  return;
119  }
120 
121  //......................................................................
123  {
124 
125  }
126 
127  //......................................................................
129  {
130 
131  fPointLabel = pset.get<std::string >("PointLabel" );
132  fUseTrig = pset.get<bool >("UseTrig" );
133  fMinParX = pset.get<float >("MinParX" );
134  fMaxParX = pset.get<float >("MaxParX" );
135  fMinParY = pset.get<float >("MinParY" );
136  fMaxParY = pset.get<float >("MaxParY" );
137  fDisplayMinX = pset.get<float >("DisplayMinX" );
138  fDisplayMaxX = pset.get<float >("DisplayMaxX" );
139  fDisplayMinY = pset.get<float >("DisplayMinY" );
140  fDisplayMaxY = pset.get<float >("DisplayMaxY" );
141  fOneSystOnly = pset.get<bool >("OneSystOnly", true );
142  fOverlayContour = pset.get<bool >("OverlayContour", false );
143  fContourFile = pset.get<std::string >("ContourFile" );
144 
145  auto hashPos = fPointLabel.find("#" );
146  auto Dmsq32Pos = fPointLabel.find("Dmsq32");
147  auto Th13Pos = fPointLabel.find("Th13" );
148  auto Th23Pos = fPointLabel.find("Th23" );
149  auto dCPPos = fPointLabel.find("dCP" );
150 
151  // these lines are for backwards compatibility and can be removed after
152  // July 2016
153  if(Dmsq32Pos == std::string::npos) Dmsq32Pos = fPointLabel.find("deltaMSqr32");
154  if(Th23Pos == std::string::npos) Th23Pos = fPointLabel.find("theta23" );
155  if(dCPPos == std::string::npos) dCPPos = fPointLabel.find("deltaCP" );
156 
157  if(Th13Pos == Th23Pos && Dmsq32Pos == Th13Pos)
158  throw cet::exception("StarPlot")
159  << "all possible parameters have the same position in the label"
160  << " - that cannot be good";
161 
162  LOG_DEBUG("StarPlot")
163  << "\nLabel : " << fPointLabel
164  << "\n\tHashPos : " << hashPos
165  << "\n\tDmsq32Pos : " << Dmsq32Pos
166  << "\n\tTh13Pos : " << Th13Pos
167  << "\n\tTh23Pos : " << Th23Pos
168  << "\n\tDCP : " << dCPPos;
169 
170  std::string subPar1;
171  std::string subPar2;
172 
173  if(hashPos == std::string::npos)
174  throw cet::exception("StarPlot")
175  << "cannot find the # in the input label, something is wrong";
176 
177  if(Dmsq32Pos != std::string::npos){
178  if(Dmsq32Pos > hashPos+2){
180  fParam2Scale = std::make_pair(1.e3, " (#times10^{-3} eV^{2})");
181  subPar2 = fPointLabel.substr(Dmsq32Pos, 6);
182  }
183  else{
185  fParam1Scale = std::make_pair(1.e3, " (#times10^{-3} eV^{2})");
186  subPar1 = fPointLabel.substr(Dmsq32Pos, 6);
187  }
188  }
189 
190  if(Th13Pos != std::string::npos){
191  if(Th13Pos > hashPos+2){
193  fParam2Scale = std::make_pair(1./3.14159, " (#pi rad)");
194  subPar2 = fPointLabel.substr(Th13Pos, 4);
195  }
196  else{
198  fParam1Scale = std::make_pair(1./3.14159, " (#pi rad)");
199  subPar1 = fPointLabel.substr(Th13Pos, 4);
200  }
201  }
202 
203  if(Th23Pos != std::string::npos){
204  if(Th23Pos > hashPos+2){
206  subPar2 = fPointLabel.substr(Th23Pos, 4);
207  if(!fUseTrig) fParam2Scale = std::make_pair(1./3.14159, " (#pi rad)");
208  }
209  else{
211  subPar1 = fPointLabel.substr(Th23Pos, 4);
212  if(!fUseTrig) fParam1Scale = std::make_pair(1./3.14159, " (#pi rad)");
213  }
214  }
215  if(dCPPos != std::string::npos){
216  if(dCPPos > hashPos+2){
218  fParam2Scale = std::make_pair(1./3.14159, " (#pi rad)");
219  subPar2 = fPointLabel.substr(dCPPos, 3);
220  }
221  else{
223  fParam1Scale = std::make_pair(1./3.14159, " (#pi rad)");
224  subPar1 = fPointLabel.substr(dCPPos, 3);
225  }
226  }
227  fAltPointLabel = fPointLabel.substr(0, hashPos+1) + subPar2 + subPar1;
228 
229  //for the numu analysis
230  // we always want Delta m^2 to be on the y axis, and the 2D drawing
231  // method always puts parameter 1 on the x axis...so make sure
232  // Delta m^2 is in the parameter 2 spot
233  //
234  // for the nue analysis
235  // we always want Th23 to be on the y axis, dCP on the x axis,
236  // and the 2D drawing method always puts parameter 1 on the x axis...
237  // so make sure Th23 is in the parameter 2 spot and dCP in
238  // the parameter 1 spot
239  //
240  //For any other analysis either use the paramter on the basis
241  //of their position in fPointLabel
242  //or create a special set for special parameters
243 
244  LOG_DEBUG("StarPlot")
245  << "Reconfigure function: "
246  << "Before Swapping: "
247  << "Paramater 1 and Enum 1: "
249  << " "
250  << "Paramater 2 and Enum 2: "
252 
256  }
257 
261  }
262 
263  LOG_DEBUG("StarPlot")
264  << "Reconfigure function: "
265  << "After Swapping: "
266  << "Paramater 1 and Enum 1: "
267  << fnex::cOscParams_Strings[fParam1Enum]
268  << " "
269  << "Paramater 2 and Enum 2: "
270  << fnex::cOscParams_Strings[fParam2Enum];
271 
272  LOG_DEBUG("StarPlot")
273  << "Alternate label is "
274  << fAltPointLabel;
275 
276  //store systematic names in a map into a human readable form
277  fmap_syst_name["AbsNormWgt"] = "Abs. Norm.";
278  fmap_syst_name["RelNormWgt"] = "Rel. Norm.";
279  fmap_syst_name["TauNormWgt"] = "Tau Norm.";
280  fmap_syst_name["CosmicMuNormWgt"] = "Cosmic #mu Norm.";
281  fmap_syst_name["RockNormOnFidWgt"] = "Rock norm.";
282  fmap_syst_name["XsecCVWgt"] = "XSecCV Weights";
283  fmap_syst_name["RPACCQEWgt"] = "RPA CC QE Wgt";
284  fmap_syst_name["MECNormWgt"] = "MEC";
285  fmap_syst_name["NearBeamSystWgt"] = "Near Beam";
286  fmap_syst_name["FarBeamSystWgt"] = "Beam Syst.";
287  fmap_syst_name["NCNormWgt"] = "NC Norm.";
288  fmap_syst_name["RecoHadEShifter"] = "Reco. E_{Had} Shift";
289  fmap_syst_name["RecoLepEShifter"] = "Reco. E_{#mu} Shift";
290  fmap_syst_name["RelRecoHadEShifter"] = "Rel. Reco. E_{Had} Shift";
291  fmap_syst_name["RelRecoLepEShifter"] = "Rel. Reco. E_{#mu} Shift";
292  fmap_syst_name["RecoNuEShifter"] = "Reco. E_{#nu} Shift";
293  fmap_syst_name["RecoNuEBirksFarSystWgt"] = "Birks Syst.";
294  fmap_syst_name["RecoNuEBirksFarShifter"] = "Birks Syst.";
295  fmap_syst_name["RecoNuECalibFarSystWgt"] = "Calib. Syst.";
296  fmap_syst_name["RecoNuECalibFarShifter"] = "Calib. Syst.";
297  fmap_syst_name["RecoNuENoiseFarSystWgt"] = "Noise Syst.";
298  fmap_syst_name["RecoNuENoiseFarShifter"] = "Noise Syst.";
299  fmap_syst_name["MaCCQE_shapeWgt"] = "GENIE Ma CC QE Shape";
300  fmap_syst_name["NormCCQEWgt"] = "GENIE CC QE Norm.";
301  fmap_syst_name["MaCCRESWgt"] = "GENIE Ma CC Resonance";
302  fmap_syst_name["MvCCRESWgt"] = "GENIE Mv CC Resonance";
303  fmap_syst_name["MaNCRESWgt"] = "GENIE Ma NC Resonance";
304  fmap_syst_name["MvNCRESWgt"] = "GENIE Mv NC Resonance";
305  fmap_syst_name["MaNCELWgt"] = "GENIE Ma NC EL";
306  fmap_syst_name["MFP_NWgt"] = "GENIE MFP";
307  fmap_syst_name["ReweightCCQEPauliSupViaKFWgt"] = "GENIE CC QE Pauli Sup Via KF";
308  fmap_syst_name["OtherGENIEWgt"] = "Other GENIE";
309  fmap_syst_name["NueBirksCSystCVNWgt"] = "Birks C";
310  fmap_syst_name["NueShiftXSystCVNWgt"] = "Calib. X Shift";
311  fmap_syst_name["NueShiftYSystCVNWgt"] = "Calib. Y Shift";
312  fmap_syst_name["NueXYAbsSystCVNWgt"] = "Abs. Calib.";
313  fmap_syst_name["NueXYRelSystCVNWgt"] = "Rel. Calib.";
314  fmap_syst_name["NueAbsNormWgt"] = "Abs. Norm.";
315  fmap_syst_name["NueRelNormWgt"] = "Rel. Norm.";
316  fmap_syst_name["RadCorrNueWgt"] = "Rad. Corr. #nu_{e}";
317  fmap_syst_name["RadCorrNueBarWgt"] = "Rad. Corr. #bar{#nu}_{e}";
318  fmap_syst_name["SecondClassCurrWgt"] = "Second Class Curr.";
319  fmap_syst_name["NueExtrapBkg2017Wgt"] = "#nu_{e} Bkg. Extrap.";
320  fmap_syst_name["NueExtrapSig2017Wgt"] = "#nu_{e} Signal Extrap.";
321  fmap_syst_name["NueNormBothWgt"] = "#nu_{e} Norm.";
322  fmap_syst_name["MECEnuShapeNuWgt"] = "MEC E_{#nu} Shape: #nu";
323  fmap_syst_name["MECEnuShapeAntiNuWgt"] = "MEC E_{#nu} Shape: anti #nu";
324  fmap_syst_name["MECShapeNuWgt"] = "MEC Shape: #nu";
325  fmap_syst_name["MECShapeAntiNuWgt"] = "MEC Shape: anti #nu";
326  fmap_syst_name["MECInitNPFracNuWgt"] = "MEC Init. NP Frac.: #nu";
327  fmap_syst_name["MECInitNPFracAntiNuWgt"] = "MEC Init. NP Frac.: anti #nu";
328  fmap_syst_name["RPARESWgt"] = "RPA Resonance";
329  fmap_syst_name["RPACCQEshapeEnhWgt"] = "RPA CCQE Shape Enh.";
330 
331  //map to store line colors for same systematics.
332  fmap_syst_color["AbsNormWgt"] = 1;
333  fmap_syst_color["CosmicMuNormWgt"] = 2;
334  fmap_syst_color["RockNormOnFidWgt"] = 205;
335  fmap_syst_color["NCNormWgt"] = 225;
336  fmap_syst_color["RelNormWgt"] = 7;
337  fmap_syst_color["TauNormWgt"] = 51;
338  fmap_syst_color["XsecCVWgt"] = 6;
339  fmap_syst_color["RPACCQEWgt"] = 96;
340  fmap_syst_color["MECNormWgt"] = 29;
341  fmap_syst_color["FarBeamSystWgt"] = 209;
342  fmap_syst_color["MaCCQE_shapeWgt"] = 15;
343  fmap_syst_color["NormCCQEWgt"] = 205;
344  fmap_syst_color["MaCCRESWgt"] = 3;
345  fmap_syst_color["MvCCRESWgt"] = 4;
346  fmap_syst_color["MaNCRESWgt"] = 30;
347  fmap_syst_color["MvNCRESWgt"] = 28;
348  fmap_syst_color["MaNCELWgt"] = 90;
349  fmap_syst_color["MFP_NWgt"] = 42;
350  fmap_syst_color["ReweightCCQEPauliSupViaKFWgt"] = 43;
351  fmap_syst_color["OtherGENIEWgt"] = 44;
352  fmap_syst_color["RecoHadEShifter"] = 12;
353  fmap_syst_color["RecoLepEShifter"] = 208;
354  fmap_syst_color["RelRecoHadEShifter"] = 8;
355  fmap_syst_color["RelRecoLepEShifter"] = 4;
356  fmap_syst_color["RecoNuEShifter"] = 6;
357  fmap_syst_color["RecoNuEBirksFarShifter"] = 219;
358  fmap_syst_color["RecoNuECalibFarShifter"] = 94;
359  fmap_syst_color["RecoNuENoiseFarShifter"] = 213;
360  fmap_syst_color["NueBirksCSystCVNWgt"] = 5;
361  fmap_syst_color["NueShiftXSystCVNWgt"] = 6;
362  fmap_syst_color["NueShiftYSystCVNWgt"] = 7;
363  fmap_syst_color["NueXYAbsSystCVNWgt"] = 8;
364  fmap_syst_color["NueXYRelSystCVNWgt"] = 9;
365  fmap_syst_color["NueAbsNormWgt"] = 28;
366  fmap_syst_color["NueRelNormWgt"] = 39;
367  fmap_syst_color["RadCorrNueWgt"] = 1;
368  fmap_syst_color["RadCorrNueBarWgt"] = 2;
369  fmap_syst_color["SecondClassCurrWgt"] = 3;
370  fmap_syst_color["NueExtrapBkg2017Wgt"] = 4;
371  fmap_syst_color["NueExtrapSig2017Wgt"] = 5;
372  fmap_syst_color["NueNormBothWgt"] = 6;
373  fmap_syst_color["MECEnuShapeNuWgt"] = 99;
374  fmap_syst_color["MECEnuShapeAntiNuWgt"] = 6;
375  fmap_syst_color["MECShapeNuWgt"] = 65;
376  fmap_syst_color["MECShapeAntiNuWgt"] = 1;
377  fmap_syst_color["MECInitNPFracNuWgt"] = 51;
378  fmap_syst_color["MECInitNPFracAntiNuWgt"] = kGray+2;
379  fmap_syst_color["RPARESWgt"] = 5;
380  fmap_syst_color["RPACCQEshapeEnhWgt"] = kMagenta+1;
381 
382  return;
383  }
384  //......................................................................
385  // Method to clear out the collections of data products after the
386  // writeResults method is called.
388  {
389  return;
390  }
391 
392  //......................................................................
394  std::vector<fnex::PointResult> & prvec,
395  std::vector<fnex::InputPoint> & ipvec,
396  std::vector<fnex::InputPoint> & ipvec_fake)
397  {
398  prvec .resize(1);
399  ipvec .resize(1);
400  ipvec_fake .resize(1);
401 
402  // get the FitPoints from the output files
404  r.getByLabel(fPointLabel, pr);
405 
406  // get the InputPoints from the output files - need these when
407  // attempting to make FC contours when there are many results for each
408  // InputPoint
410  r.getByLabel(fPointLabel, ip);
411 
412  // get the Fake InputPoints from the output files
414  r.getByLabel(fPointLabel, "FakeInput", ip_fake);
415 
416  // check that the handles are valid, add the points to the map
417  // and move along
418  if( !ip.isValid() || !pr.isValid() ){
419 
420  // see if the alternative label works
421  r.getByLabel(fAltPointLabel, pr);
422  r.getByLabel(fAltPointLabel, ip);
423  r.getByLabel(fAltPointLabel, "FakeInput", ip_fake);
424 
425  if( !ip.isValid() || !pr.isValid() ) {
426  LOG_WARNING("StarPlot")
427  << "unable to find any results in either\n"
428  << fPointLabel
429  << "\n or "
430  << fAltPointLabel;
431 
432  // no valid handles, clear the vectors.
433  prvec.clear();
434  ipvec.clear();
435 
436  return false;
437  }
438 
439  if( !ip_fake.isValid()) {
440  LOG_VERBATIM("StarPlot")
441  << "Fake Input Point not available: No worries";
442  }
443 
444  }
445 
446  // got here, so we have valid handles
447  prvec[0] = *pr;
448  ipvec[0] = *ip;
449  ipvec_fake[0] = *ip_fake;
450 
451  return true;
452  }
453 
454  //......................................................................
456  {
457 
458  std::vector<fnex::PointResult> prvec;
459  std::vector<fnex::InputPoint> ipvec;
460  std::vector<fnex::InputPoint> ipvec_fake;
461 
462  // check to see if we have single points stored in the
463  // results, and if not look for vectors of points.
464  // if none of those either, give up
465  if( !this->readSingleResults(r, prvec, ipvec, ipvec_fake) ){
466 
467  LOG_WARNING("StarPlot")
468  << "No valid handles to either single fnex::PointResult, fnex::InputPoint "
469  << "objects of them, skip this file";
470 
471  return;
472 
473  } // end if no single point results
474 
475  // we found the results, so loop over the vectors and add them to the
476  // map
479  fnex::InputPoint ip_fake;
480  int fitStatus;
481  for(size_t p = 0; p < prvec.size(); ++p){
482 
483  pr = prvec[p];
484  ip = ipvec[p];
485  ip_fake = ipvec_fake[p];
486 
487  fitStatus = pr.FitStatus();
488 
489  if(fitStatus > -1 && fitStatus < 4){
490  this->AddPointToMap(ip, pr, ip_fake);
491  }
492  else{
493  // fit did not converge, so
494  // create a new PointResult with a very large chiSqr
495  // I am assuming that a value of -1 means the fit did not
496  // converge, so substitute it for a stupidly large chi^2
497  LOG_VERBATIM("StarPlot")
498  << "switching chi^2 for "
499  << ip
500  << " "
501  << pr.ChiSqr()
502  << " and status "
503  << fitStatus
504  << " to "
506 
507  pr = fnex::PointResult(std::numeric_limits<double>::max(),
508  fitStatus,
510  pr.SystematicPulls());
511  this->AddPointToMap(ip, pr, ip_fake);
512  }
513 
514  } // end loop over points found in file
515 
516  return;
517  }
518  //......................................................................
520  fnex::PointResult const& pr,
521  fnex::InputPoint const& ip_fake)
522  {
523  LOG_DEBUG("StarPlot")
524  << ip
525  << pr;
526 
527  // Add the point to the map, the first entry in the set is certain
528  // to be the lowest chi^2 value for the input point because that is
529  // how the set and < operator of PointResult work
530  // the key is the pair of the input parameter values
531  auto paramDetMap = pr.ParameterSpaceLocation();
532 
533  for(auto itr : paramDetMap)
534  LOG_DEBUG("StarPlot")
535  << itr.first.first
536  << " "
537  << itr.first.second
538  << " "
539  << itr.second;
540 
541  LOG_DEBUG("StarPlot")
542  << " have the paramMap with " << paramDetMap.size() << " entries";
543 
546 
547  if(paramDetMap.count(pd1) < 1 ||
548  paramDetMap.count(pd2) < 1)
549  throw cet::exception("StarPlot")
550  << "Cannot find one of the parameters in the map: "
551  << "\n" << pd1.first
552  << " " << paramDetMap.count(pd1)
553  << "\n" << pd2.first
554  << " " << paramDetMap.count(pd2);
555 
556  std::pair<double, double> inputPoint = std::make_pair(fParam1Scale.first * paramDetMap.find(pd1)->second,
557  fParam2Scale.first * paramDetMap.find(pd2)->second);
558  LOG_DEBUG("StarPlot")
559  << "AddPointToMap function: "
560  << "Parameter 1 and Enum 1: "
562  << " "
563  << "Parameter 2 and Enum 2: "
565 
566  //If parameter 1 is Th23 then use sin(Th23)*sin(Th23) form
567  if(fUseTrig && fParam1Enum == fnex::kTh23){
568  inputPoint.first = std::sin(inputPoint.first) * std::sin(inputPoint.first);
569  }
570  //OR If parameter 2 is Th23 then use sin(Th23)*sin(Th23) form
571  if(fUseTrig && fParam2Enum == fnex::kTh23){
572  inputPoint.second = std::sin(inputPoint.second) * std::sin(inputPoint.second);
573  }
574 
575  fPointResMap[inputPoint].emplace(pr);
576 
577  if(inputPoint.first < 0) inputPoint.first = inputPoint.first + 2.;
578  if(inputPoint.second < 0) inputPoint.second = inputPoint.second + 2.;
579 
580  LOG_VERBATIM("StarPlot")
581  << "Adding "
582  << inputPoint.first
583  << " "
584  << inputPoint.second
585  << " to map with chi^2 = "
586  << pr.ChiSqr()
587  << " size = "
588  << fPointResMap.size()
589  << " chi^2 = "
590  << pr.ChiSqr();
591 
592  LOG_DEBUG("StarPlot")
593  << "Point Result: "
594  << pr;
595 
596  //Store the values of the x and y parameters
597  fparX_vec.push_back(inputPoint.first);
598  fparY_vec.push_back(inputPoint.second);
599 
600  LOG_VERBATIM("StarPlot") << "Getting systematic names";
601 
602  bool gotshifter = false;
603  //Store the names of the systematic names
604  auto sysParDets = ip_fake.SystematicPulls(novadaq::cnv::kFARDET);
605  for(auto itr : sysParDets){
606  if(itr.first=="XsecCVWgt" || (gotshifter && fOneSystOnly)) continue;
607  if(itr.first=="CalibFarSystWgt" ||
608  itr.first=="NoiseFarSystWgt" ||
609  itr.first=="BirksFarSystWgt") continue;
610  LOG_VERBATIM("StarPlot") << "adding systematic " << itr.first;
611  fsyst_vec.push_back(itr.first);
612  gotshifter = true;
613  }
614 
615  if(!gotshifter){
616  sysParDets = ip_fake.SystematicPulls(novadaq::cnv::kNEARDET);
617  for(auto itr : sysParDets){
618  if(itr.first=="XsecCVWgt" || (gotshifter && fOneSystOnly)) continue;
619  LOG_VERBATIM("StarPlot") << "adding systematic " << itr.first;
620  fsyst_vec.push_back(itr.first);
621  gotshifter = true;
622  }
623  }
624 
625  //Store the values of the systematic shift in the fake input point
626  gotshifter = false;
627  auto sysParDets_fake = ip_fake.SystematicPulls(novadaq::cnv::kFARDET);
628  for(auto itr : sysParDets_fake){
629  if(itr.first=="XsecCVWgt" || (gotshifter && fOneSystOnly)) continue;
630  if(itr.first=="CalibFarSystWgt" ||
631  itr.first=="NoiseFarSystWgt" ||
632  itr.first=="BirksFarSystWgt") continue;
633  fplusminus_shift.push_back(itr.second);
634  gotshifter = true;
635  }
636  if (!gotshifter){
637  sysParDets_fake = ip_fake.SystematicPulls(novadaq::cnv::kNEARDET);
638  for(auto itr : sysParDets_fake){
639  if(itr.first=="XsecCVWgt" || (gotshifter && fOneSystOnly)) continue;
640  fplusminus_shift.push_back(itr.second);
641  gotshifter = true;
642  }
643  }
644 
645  return;
646  }
647  //......................................................................
649  {
650  // not writing anything to the file for this plugin at this time
651 
652  return;
653  } // end writeResults
654 
655  //......................................................................
657  {
658 
659  LOG_VERBATIM("StarPlot") << "Making star plot";
660  LOG_VERBATIM("StarPlot") << "Making histogram";
661 
664  "_StarPlot");
665 
668 
669  LOG_VERBATIM("StarPlot") << "Do we need to sin th23?";
670 
672  XaxisTitle = "sin^{2}(#theta_{23})";
674  YaxisTitle = "sin^{2}(#theta_{23})";
675 
676  std::string title = (";" + XaxisTitle +
677  ";" + YaxisTitle);
678  LOG_VERBATIM("StarPlot") << "histogram title: " << title;
679 
680  TH2F* h_starplot = fTFS->make<TH2F>(histName.c_str(),
681  title.c_str(),
682  500, fMinParX, fMaxParX,
683  500, fMinParY, fMaxParY );
684  h_starplot->SetStats(false);
685  h_starplot->GetXaxis()->SetRangeUser(fDisplayMinX, fDisplayMaxX);
686  h_starplot->GetYaxis()->SetRangeUser(fDisplayMinY, fDisplayMaxY);
687 
688  LOG_VERBATIM("StarPlot") << "Getting stat-only values: ";
689 
690  //Stat only value
691  float stat_par1 = fPointResMap.begin()->first.first;
692  float stat_par2 = fPointResMap.begin()->first.second;
693 
694  LOG_VERBATIM("StarPlot") << "\t" << stat_par1 << "\t" << stat_par2;
695 
696  //loop over one less than number of input files
697  //to make the TLine objects and draw them
698  //assuming the first file is always Stat only
699 
700  //Canvas to draw the Star Plot
701  TCanvas* fcan = fTFS->make<TCanvas>("StarPlot",
702  "Star Plot",
703  1200 ,
704  800 );
705  fcan ->Divide(2, 2);
706 
707  TCanvas* fcan2 = fTFS->make<TCanvas>("StarPlot2",
708  "Star Plot",
709  1600 ,
710  800 );
711  fcan2->Divide(2);
712 
713  fcan->cd(1);
714  h_starplot->Draw();
715  fcan2->cd(1);
716  h_starplot->Draw();
717 
718  LOG_VERBATIM("StarPlot") << "Vector X has " << fparX_vec.size() << " entries, Y has " << fparY_vec.size();
719 
720  //Marker for presenting the Stat only point
721  TMarker* marker_stat = new TMarker(fparX_vec[0], fparY_vec[0], 20);
722  marker_stat ->SetMarkerColor(46);
723 
724  //Legends
725  TLegend *legend = new TLegend(.22,.01,.77,.90);
726  legend->SetBorderSize(0);
727  legend->SetFillStyle(0);
728  legend->SetTextSize(0.04);
729  legend->AddEntry((TObject*)0, "Systematics Effects", " " );
730  legend->AddEntry(marker_stat, "Input Value", "p" );
731 
732  for(unsigned int i = 1; i < fparX_vec.size(); ++i)
733  {
734  LOG_VERBATIM("StarPlot") << "entry " << i << " of " << fparX_vec.size();
735 
736  TLine* l = new TLine(fparX_vec[0], fparY_vec[0], fparX_vec[i], fparY_vec[i]);
737  l->SetLineWidth(5);
738  l->SetLineColor(fmap_syst_color[fsyst_vec[i-1]]);
739  // if(i < 10) l->SetLineColor(i);
740  // else l->SetLineColor(i+30);
741  if(fplusminus_shift[i-1]==1) l->SetLineStyle(1);
742  else if(fplusminus_shift[i-1]==-1) l->SetLineStyle(2);
743  else l->SetLineStyle(10);
744  fcan->cd(1);
745  l->Draw("same");
746  fcan2->cd(1);
747  l->Draw("same");
748  //Markers for each of the Systematics
749  TMarker* marker = new TMarker(fparX_vec[i], fparY_vec[i], 20);
750  marker->SetMarkerColor(l->GetLineColor());
751 
752  fcan ->cd(1);
753  marker->Draw("same");
754  fcan2 ->cd(1);
755  marker->Draw("same");
756 
757  legend->AddEntry(l, (fmap_syst_name[fsyst_vec[i-1]]+": "+TString::Format("%g #sigma", fplusminus_shift[i-1])).Data() , "l");
758 
759  LOG_VERBATIM("StarPlot")
760  << "Systematic Names: "
761  << fmap_syst_name[ fsyst_vec[i-1] ]
762  << " "
763  << "Fake Input Value: "
764  << fplusminus_shift[i-1];
765 
766  }//end of looping over number of lines drawn
767 
768  fcan ->cd(1);
769  marker_stat->Draw("same");
770  fcan2->cd(1);
771  marker_stat->Draw("same");
772 
773  //Legends
774  fcan ->cd(2);
775  legend ->Draw();
776  fcan2 ->cd(2);
777  legend ->Draw();
778 
779  //Display Plus Systematic values
780  fcan ->cd(3);
781 
782  TPaveText *pt_plus = new TPaveText(.05,.1,.95,.8);
783  pt_plus ->SetFillColor(kWhite);
784  pt_plus ->SetBorderSize(0);
785 
786  std::stringstream stat_only;
787  std::stringstream syst_heading;
788  std::stringstream syst_all;
789  syst_heading
790  << std::setw(25) << std::right
791  << "Case: "
792  << std::setw(10) << std::right
793  << h_starplot->GetXaxis()->GetTitle()
794  << std::setw(25) << std::right
795  << h_starplot->GetYaxis()->GetTitle();
796 
797  stat_only
798  << std::setw(25) << std::right
799  << "Input Point: "
800  << std::setw(25) << std::right
801  << std::setprecision(4)
802  << marker_stat->GetX()
803  << std::setw(25) << std::right
804  << std::setprecision(3)
805  << marker_stat->GetY();
806 
807  pt_plus->AddText(syst_heading .str().data());
808  pt_plus->AddText(stat_only .str().data());
809 
810  for(unsigned int i = 1; i < fparX_vec.size(); ++i)
811  {
812  if(fplusminus_shift[i-1]==-1) continue;
813  syst_all.str("");
814 
815  syst_all
816  << std::setw(25) << std::right
817  << fmap_syst_name[ fsyst_vec[i-1] ]
818  << " "
819  << TString::Format("%g #sigma", fplusminus_shift[i-1])
820  << std::setw(10) << std::right
821  << std::setprecision(4)
822  << fparX_vec[i]
823  << std::setw(10) << std::right
824  << " ( "
825  << std::setprecision(3)
826  << ParameterDiff( marker_stat->GetX(), fparX_vec[i])
827  << " %) "
828  << std::setw(12) << std::right
829  << std::setprecision(3)
830  << fparY_vec[i]
831  << std::setw(10) << std::right
832  << " ( "
833  << std::setprecision(3)
834  << ParameterDiff(marker_stat->GetY(), fparY_vec[i])
835  << " %) ";
836 
837  pt_plus->AddText(syst_all .str().data());
838  pt_plus->SetTextSize(0.035);
839 
840  }
841  pt_plus->Draw();
842 
843  //Display Minus Systematic values
844  fcan ->cd(4);
845 
846  TPaveText *pt_minus = new TPaveText(.05,.1,.95,.8);
847  pt_minus->SetFillColor(kWhite);
848  pt_minus->SetBorderSize(0);
849  pt_minus->AddText(syst_heading .str().data());
850  pt_minus->AddText(stat_only .str().data());
851 
852  for(unsigned int i = 1; i < fparX_vec.size(); ++i)
853  {
854  if(fplusminus_shift[i-1]==1) continue;
855  syst_all.str("");
856 
857  syst_all
858  << std::setw(25) << std::right
859  << fmap_syst_name[ fsyst_vec[i-1] ]
860  << " "
861  << TString::Format("%g #sigma", fplusminus_shift[i-1])
862  << std::setw(10) << std::right
863  << std::setprecision(4)
864  << fparX_vec[i]
865  << std::setw(10) << std::right
866  << " ( "
867  << std::setprecision(3)
868  << ParameterDiff( marker_stat->GetX(), fparX_vec[i])
869  << " %) "
870  << std::setw(12) << std::right
871  << std::setprecision(3)
872  << fparY_vec[i]
873  << std::setw(10) << std::right
874  << " ( "
875  << std::setprecision(3)
876  << ParameterDiff(marker_stat->GetY(), fparY_vec[i])
877  << " %) ";
878 
879  pt_minus->AddText(syst_all .str().data());
880  pt_minus->SetTextSize(0.035);
881 
882  }
883  pt_minus->Draw();
884 
885  fcan ->Write();
886  fcan2->Write();
887 
888  //If we want to overlay the star over a contour, do so
889  if(fOverlayContour)
890  {
891  LOG_VERBATIM("StarPlot") << "\nStarting to try to overlay with a contour. Oh joy...";
892  TFile * contourfile = TFile::Open(fContourFile.data());
893  LOG_VERBATIM("StarPlot") << "Got file";
894  TCanvas * contourcanv = fTFS->make<TCanvas>("contourcanv", "contourcanv", 1200, 800);
895  contourcanv = (TCanvas*)contourfile->Get("contours/GaussianContoursCanv")->Clone();
896  LOG_VERBATIM("StarPlot") << "Got canvas";
897 
898  for(unsigned int i = 1; i < fparX_vec.size(); ++i)
899  {
900  LOG_VERBATIM("StarPlot") << "entry " << i << " of " << fparX_vec.size();
901 
902  TLine* l = new TLine(fparX_vec[0], fparY_vec[0], fparX_vec[i], fparY_vec[i]);
903  l->SetLineWidth(5);
904  l->SetLineColor(fmap_syst_color[fsyst_vec[i-1]]);
905  // if(i < 10) l->SetLineColor(i);
906  // else l->SetLineColor(i+30);
907  if(fplusminus_shift[i-1]==1) l->SetLineStyle(1);
908  else if(fplusminus_shift[i-1]==-1) l->SetLineStyle(2);
909  else l->SetLineStyle(10);
910  l->Draw("same");
911  //Markers for each of the Systematics
912  TMarker* marker = new TMarker(fparX_vec[i], fparY_vec[i], 20);
913  marker->SetMarkerColor(l->GetLineColor());
914 
915  marker->Draw("same");
916 
917  LOG_VERBATIM("StarPlot")
918  << "Systematic Names: "
919  << fmap_syst_name[ fsyst_vec[i-1] ]
920  << " "
921  << "Fake Input Value: "
922  << fplusminus_shift[i-1];
923 
924  }//end of looping over number of lines drawn
925  contourcanv->Write();
926 
927  TCanvas * legcanv = fTFS->make<TCanvas>("legcanv", "legend", 400, 800);
928  legcanv->cd();
929  legend->Draw();
930  legcanv->Write();
931  }//end of if(fOverlayContour)
932  }
933  //......................................................................
934 
936  {
937  LOG_VERBATIM("StarPlot") << "Making table??";
938 
939  }
940 
941  //......................................................................
942 
943  float StarPlot::ParameterDiff(float x1, float x2)
944  {
945  return (x2-x1)/x1 * 100.;
946  }
947 
948  //......................................................................
950  {
951 
952  this->MakeStarPlot();
953  this->MakeTable();
954 
955  return;
956  }
957 
958 
959 } // end fnex namespace
960 
961 namespace fnex{
962 
964 
965 }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
void writeResults(art::Results &r)
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:112
std::pair< double, double > fParam1Extrema
first is the min, second is the max
float fDisplayMinY
Display minimum X range in Star Plot.
int const & FitStatus() const
Definition: PointResult.h:71
std::map< TString, int > fmap_syst_color
map to store line colors for same systematics.
fnex::OscParm_t fParam1Enum
enumerated value of the parameter from FNEX/dataProducts/Constants.h
std::map< std::pair< double, double >, std::set< fnex::PointResult > > PointResMap
void reconfigure(fhicl::ParameterSet const &p)
Float_t x1[n_points_granero]
Definition: compare.C:5
const char * p
Definition: xmltok.h:285
void readResults(art::Results const &r)
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
TString ip
Definition: loadincs.C:5
ParameterDetMap const SystematicPulls() const
art::ServiceHandle< art::TFileService > fTFS
TFileService.
float fMinParX
Minimum X range in Star Plot.
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.
float fDisplayMaxX
Display minimum X range in Star Plot.
bool fOneSystOnly
only looks for one shifter or weighter per input file
const XML_Char const XML_Char * data
Definition: expat.h:268
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
bool isValid() const
Definition: Handle.h:189
float fMaxParY
Maximum Y range in Star Plot.
Far Detector at Ash River, MN.
float ParameterDiff(float, float)
std::vector< TString > fsyst_vec
A vector that stores systematics names.
#define DEFINE_ART_RESULTS_PLUGIN(klass)
PointResMap fPointResMap
Initial point passed to fitter, typedef above.
T get(std::string const &key) const
Definition: ParameterSet.h:231
std::vector< float > fparX_vec
A vector that stores values to be plotted on X axis;.
float fDisplayMaxY
Display minimum X range in Star Plot.
Near Detector in the NuMI cavern.
std::string fAltPointLabel
reverse the order of the parameters in the label
bool fUseTrig
Turn on plotting of trig functions for the angles.
std::string fContourFile
filename of contour to overlay plots on
double const & ChiSqr() const
Definition: PointResult.h:70
std::vector< float > fparY_vec
A vector that stores values to be plotted on Y axis.
float fDisplayMinX
Display minimum X range in Star Plot.
ParameterDetMap const ParameterSpaceLocation() const
Definition: PointResult.cxx:82
std::pair< double, std::string > fParam2Scale
scale the parameter if necessary ie is x10^{-3}, string is latex friendly
#define LOG_WARNING(category)
enum fnex::osc_params OscParm_t
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
float fMaxParX
Maximum X range in Star Plot.
std::vector< float > fplusminus_shift
A vector that stores the plus or minus values of the systematic shift in the fake input point...
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
bool fOverlayContour
true if we want to overlay the results on a contour
std::pair< double, std::string > fParam1Scale
scale the parameter if necessary ie is x10^{-3}, string is latex friendly
T * make(ARGS...args) const
T sin(T number)
Definition: d0nt_math.hpp:132
std::map< std::pair< double, double >, double > PointMap
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void Legends(TLine *l)
std::map< TString, TString > fmap_syst_name
map to store human readable systematic names.
StarPlot(fhicl::ParameterSet const &pset)
const std::string cOscParams_Strings[kNumOscParams]
Definition: Constants.h:196
std::pair< std::string, novadaq::cnv::DetId > ParameterDet
Definition: Constants.h:38
TRandom3 r(0)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const std::string cOscParams_Strings_Latex[kNumOscParams]
Definition: Constants.h:205
std::string fPointLabel
label of plugin storing the points
#define LOG_VERBATIM(category)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
std::pair< double, double > fParam2Extrema
first is the min, second is the max
float fMinParY
Minimum Y range in Star Plot.
bool readSingleResults(art::Results const &r, std::vector< fnex::PointResult > &prvec, std::vector< fnex::InputPoint > &ipvec, std::vector< fnex::InputPoint > &ipvec_fake)
fnex::OscParm_t fParam2Enum
enumerated value of the parameter from FNEX/dataProducts/Constants.h
void AddPointToMap(fnex::InputPoint const &ip, fnex::PointResult const &pr, fnex::InputPoint const &ip_fake)
ParameterMap const SystematicPulls(novadaq::cnv::DetId const &det) const
Definition: InputPoint.cxx:364