CovarianceMatrixFitter_plugin.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 /// \brief Use covariance matrices in fits
3 /// \author brebel@fnal.gov
4 ///////////////////////////////////////////////////////////////////////
5 #include <cmath>
6 #include <string>
7 #include <map>
8 
9 // Framework includes
13 #include "art_root_io/TFileDirectory.h"
14 #include "fhiclcpp/ParameterSet.h"
16 
17 // NOvA includes
29 
30 // ROOT includes
31 #include "Math/Functor.h"
32 
33 namespace cmf {
34 
36  public:
37 
38  explicit CovarianceMatrixFitter(fhicl::ParameterSet const& pset);
39  virtual ~CovarianceMatrixFitter();
40 
41  // Plugins should not be copied or assigned.
46 
47  void readResults(art::Results const& r) override;
48  void writeResults(art::Results & r) override;
49  void clear() override;
50  void beginJob() override;
51 
52  void reconfigure(fhicl::ParameterSet const& pset);
53 
54  private:
55 
56  void FillDataSpectrum(cmf::Spectrum & spec);
57 
58  fhicl::ParameterSet fManipulatorPars; ///< event list manipulator parameters
59  fhicl::ParameterSet fFakeParameters; ///< the fake point parameters
60  fhicl::ParameterSet fChiSqrCalcPars; ///< configuration for the ChiSqrCalculator singleton
61  fhicl::ParameterSet fHelperPars; ///< event list manipulator parameters
62  fhicl::ParameterSet fSAWParameters; ///< ShifterAndWeighter parameters
63  fhicl::ParameterSet fMinimizerPars; ///< ROOT minimizer configuration
64  std::vector<std::string> fChiSqrToUse; ///< vector of chi^2 calculation methods to use for fitting
65  std::vector<cmf::OscParamPoint> fPredictionLib; ///< library of predictions
66  std::vector<cmf::FakeUniverse> fFakeUniverses; ///< library of predictions
67  std::string fPredictionModule; ///< label for module containing the prediction point library
68  std::string fInputUniverseLabel; ///< label for module containing the fake universe we want to fit
69  size_t fUniverseToFit; ///< which universe to fit, if std::numeric_limits<unsigned int>::max(), then Asimov
70  float fFourFlavord24Val; ///< value of delta_24 to use
71  bool fPerformMinimization; ///< tell fitter whether to perform minimisation or not
72  };
73 
74  //......................................................................
76  : fManipulatorPars(pset.get<fhicl::ParameterSet>("EventListManipulator" ))
77  , fFakeParameters (pset.get<fhicl::ParameterSet>("FakeParameters", fhicl::ParameterSet()))
78  , fChiSqrCalcPars (pset.get<fhicl::ParameterSet>("ChiSqrCalculationParameters" ))
79  , fHelperPars (pset.get<fhicl::ParameterSet>("FitHelperParameters" ))
80  , fSAWParameters (pset.get<fhicl::ParameterSet>("ShifterAndWeighterParameters" ))
81  , fMinimizerPars (pset.get<fhicl::ParameterSet>("ROOTMinimizerParameters" ))
82  {
83  this->reconfigure(pset);
84 
85  for(auto const& itr : fChiSqrToUse)
86  produces<cmf::PointResult>(itr);
87 
88  if( !fFakeParameters.is_empty() )
89  produces<cmf::InputPoint>("FakeInput");
90  }
91 
92  //......................................................................
94  {
95  }
96 
97  //......................................................................
99  {
102  cmf::ParameterUtility::Instance() ->Initialize(pset.get< fhicl::ParameterSet >("ParametersToUse" ));
103  cmf::PlotUtilities::Instance() ->Initialize(pset.get< fhicl::ParameterSet >("PlotUtilities" ));
104 
105  fChiSqrToUse = pset.get<std::vector<std::string>>("ChiSqrCalcsToUse");
106 
107  fFourFlavord24Val = pset.get<float >("FourFlavord24Value", 0.0 );
108  fPerformMinimization = pset.get<bool >("PerformMinimization" );
109  fPredictionModule = pset.get<std::string >("PredictionLibraryModule", std::string() );
110  fInputUniverseLabel = pset.get<std::string >("InputUniverseModule", std::string() );
111  fUniverseToFit = pset.get<unsigned int>("UniverseToFit", std::numeric_limits<unsigned int>::max());
112 
113  MF_LOG_DEBUG("CovarianceMatrixFitter")
114  << pset.to_indented_string();
115 
116  }
117 
118  //......................................................................
119  // Method to clear out the collections of data products after the
120  // writeResults method is called.
122  {
123  }
124 
125  //......................................................................
127  {
128  }
129 
130  //......................................................................
132  {
134  if( r.getByLabel(fInputUniverseLabel, universes) &&
135  universes.isValid() ){
136 
137  fFakeUniverses.reserve(universes->size() + fFakeUniverses.size());
138  for(auto const& itr : *universes) fFakeUniverses.emplace_back(itr);
139  MF_LOG_VERBATIM("FitFeldmanCousinsPoint")
140  << "Valid handle to std::vector<cmf::FakeUniverse> "
141  << "objects found in "
143  << " vector size "
144  << fFakeUniverses.size();
145  r.removeCachedProduct(universes);
146  }
147 
148 
149  if(fPredictionModule.empty()){
150  MF_LOG_VERBATIM("CovarianceMatrixFitter")
151  << "no prediction module label specified, will not attempt to read in a library";
152  return;
153  }
154 
156  if( r.getByLabel(fPredictionModule, predLib) &&
157  predLib.isValid() ){
158  MF_LOG_VERBATIM("CovarianceMatrixFitter")
159  << "Valid handle to std::vector<cmf::OscParamPoint> "
160  << "objects found in "
162  << " "
163  << predLib->size();
164 
165  for (auto const &itr : *predLib){
166 
167  MF_LOG_DEBUG("CMFContourFromLibrary")
168  << itr;
169 
170  auto tmpMap = itr.OscPoint();
171  if(tmpMap.count(cmf::kd24) > 0){
172  if(tmpMap.find(cmf::kd24)->second != fFourFlavord24Val){
173  MF_LOG_DEBUG("CMFContourFromLibrary")
174  << "don't use this point";
175 
176  continue;
177  }
178  }
179 
180  // add the predictions to the full list
181  fPredictionLib.emplace_back(itr);
182  } // end loop over predictions
183 
184  // release space needed to read in the prediction library
185  r.removeCachedProduct(predLib);
186  }
187  }
188 
189  //......................................................................
191  {
192 
193  // if we were passed fake universes, we don't need to do anything with the
194  // event list manipulator
195  if(!fFakeUniverses.empty()){
197  spec = fFakeUniverses.front().AsimovSpectrum();
198  else if(fUniverseToFit < fFakeUniverses.size())
199  spec = fFakeUniverses[fUniverseToFit].PoissonSpectrum();
200  else
201  throw cet::exception("CovarianceMatrixFitter")
202  << "Requested to fit fake universe "
203  << fUniverseToFit
204  << " but vector only contains "
205  << fFakeUniverses.size()
206  << " universes";
207 
208  return;
209  }
210 
212  cmf::EventListColl evLists;
213 
214  if(!fFakeParameters.is_empty()){
217  cmf::Location fakeLoc;
218  // loop over the systematic and oscillation parameters
219  auto const& partypes = fFakeParameters.get_names();
220 
221  MF_LOG_DEBUG("CovarianceMatrixFitter")
223 
224  for(auto const& itr : partypes){
225  MF_LOG_DEBUG("CovarianceMatrixFitter")
226  << itr;
227 
228  for(auto const& detItr : cmf::SelectionUtility::Instance()->DetectorsToUse()){
230  detItr,
231  fakeLoc);
232  }
233  }
234 
235  MF_LOG_VERBATIM("CovarianceMatrixFitter")
236  << "We have fake parameters, so making fake data generated at\n"
237  << fakeLoc;
238 
240 
241  elm.Deserialize(evLists, cmf::kMC);
242 
243  } // end if we want fake data
244  else{
245  elm.Deserialize(evLists, cmf::kData);
246  }
247 
249 
250  cmf::FillSpectrum(evLists,
251  elm.ExposureMap(),
252  spec,
253  count);
254 
255 
256  // for(size_t b = 0; b < spec.size(); ++b)
257  // MF_LOG_VERBATIM("CovarianceFitHelper")
258  // << "Data "
259  // << b
260  // << " "
261  // << spec[b];
262 
263  }
264 
265  //......................................................................
266  // Should be able to both do the fit and make the non-FC corrected contours
267  // in this method. The contours have to be super easy as we just need to
268  // go on the grid and find the chi^2 at each point - there is no fitting
269  // of systematic uncertainties to slow us down and force a separate set
270  // of jobs for that.
272  {
273  // get the data spectrum and exposure map to pass to the fit helper
274  cmf::Spectrum dataSpectrum;
275  this->FillDataSpectrum(dataSpectrum);
276 
282  dataSpectrum);
283 
284  if(fPredictionLib.size() > 0 && fPerformMinimization)
286 
287  for(auto const& chiSqrCalc : fChiSqrToUse){
288 
291 
292  // create plots showing the minimizer progress
293  cmf::CovarianceFitHelper::Instance()->MakeResultPlots("FitResults", chiSqrCalc);
294 
295  // create the point result to store
296  std::unique_ptr<cmf::PointResult> pr = cmf::CovarianceFitHelper::Instance()->Result();
297 
298  r.put(std::move(pr), chiSqrCalc);
299  } // end loop over chi^2 calculation methods
300 
301  if( !fFakeParameters.is_empty() ){
302 
304  // loop over the systematic and oscillation parameters
305  auto const& partypes = fFakeParameters.get_names();
306  for(auto const& itr : partypes){
307  for(auto const& detItr : cmf::SelectionUtility::Instance()->DetectorsToUse()){
309  detItr,
310  loc);
311  }
312  }
313 
314  std::unique_ptr<cmf::InputPoint> fakeInputPoint = std::make_unique<InputPoint>(loc);
315  MF_LOG_VERBATIM("FitPoint")
316  << "Fake Input Point\n"
317  << *fakeInputPoint;
318 
319  r.put(std::move(fakeInputPoint), "FakeInput");
320  } // end if we are making the fake input point to store
321  }
322 
324 
325 } // end cmf namespace
bool fPerformMinimization
tell fitter whether to perform minimisation or not
void FillSpectrum(cmf::EventListColl const &eventLists, cmf::ExposureMap const &exposureMap, cmf::Spectrum &spectrum, cmf::Spectrum &count, cmf::InteractionType_t intType, bool applyExposureNorm)
CovarianceMatrixFitter(fhicl::ParameterSet const &pset)
fhicl::ParameterSet fChiSqrCalcPars
configuration for the ChiSqrCalculator singleton
void Initialize(fhicl::ParameterSet const &pset)
bool is_empty() const
void reconfigure(fhicl::ParameterSet const &pset)
static SelectionUtility * Instance()
std::vector< double > Spectrum
Definition: Constants.h:746
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::vector< std::string > fChiSqrToUse
vector of chi^2 calculation methods to use for fitting
std::string fInputUniverseLabel
label for module containing the fake universe we want to fit
void Deserialize(cmf::EventListColl &eventLists, cmf::DataMC_t dataMC=cmf::kBoth, std::set< cmf::DetType_t > const &detectors=std::set< cmf::DetType_t >({cmf::kNEARDET, cmf::kFARDET}))
bool isValid() const
Definition: Handle.h:183
static ParameterUtility * Instance()
void Minimize(cmf::ChiSqrCalc_t const &chiSqrCalc=cmf::kCovMat)
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
void readResults(art::Results const &r) override
fhicl::ParameterSet fMinimizerPars
ROOT minimizer configuration.
void writeResults(art::Results &r) override
#define DEFINE_ART_RESULTS_PLUGIN(klass)
T get(std::string const &key) const
Definition: ParameterSet.h:231
CovarianceMatrixFitter & operator=(CovarianceMatrixFitter const &)=delete
void FindInitialGuess(std::vector< cmf::OscParamPoint > const &library)
cmf::ExposureMap const & ExposureMap() const
void Initialize(fhicl::ParameterSet const &pset)
static ShifterAndWeighter * Instance()
std::vector< std::string > get_names() const
void Initialize(fhicl::ParameterSet const &plottingPars)
float fFourFlavord24Val
value of delta_24 to use
bool removeCachedProduct(Handle< PROD > &) const
Definition: DataViewImpl.h:971
static CovarianceFitHelper * Instance()
std::string fPredictionModule
label for module containing the prediction point library
fhicl::ParameterSet fHelperPars
event list manipulator parameters
static cmf::ChiSqrCalc_t StringToChiSqrType(std::string const &str)
Module to combine a set of results into a single file currently only does one data product type at a ...
Definition: Event.cxx:24
void ParameterSetToLocation(fhicl::ParameterSet const &pset, cmf::DetType_t const &det, cmf::Location &loc)
void MakeResultPlots(std::string const &dirBaseName, std::string const &uniqueID=std::string(), bool fillFromLists=true)
std::vector< cmf::EventList > EventListColl
Definition: Event.h:150
std::string to_indented_string() const
void InitShiftsAndWeightsToUse(cmf::Location const &loc, fhicl::ParameterSet const &parset)
std::vector< cmf::OscParamPoint > fPredictionLib
library of predictions
#define MF_LOG_VERBATIM(category)
fhicl::ParameterSet fFakeParameters
the fake point parameters
void Initialize(fhicl::ParameterSet const &helperPars, fhicl::ParameterSet const &chiSqrCalcPars, fhicl::ParameterSet const &manipulatorPars, fhicl::ParameterSet const &minimizerPars, fhicl::ParameterSet const &sawPars)
std::vector< cmf::FakeUniverse > fFakeUniverses
library of predictions
#define MF_LOG_DEBUG(id)
static PlotUtilities * Instance()
TRandom3 r(0)
size_t fUniverseToFit
which universe to fit, if std::numeric_limits<unsigned int>::max(), then Asimov
void FillDataSpectrum(cmf::Spectrum &spec)
void Initialize(fhicl::ParameterSet const &pset)
fhicl::ParameterSet fManipulatorPars
event list manipulator parameters
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:66
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
void SetCurrentVals(cmf::Location const &loc)
std::unique_ptr< cmf::PointResult > Result(bool loadFromEventList=true)
fhicl::ParameterSet fSAWParameters
ShifterAndWeighter parameters.
static CovarianceBinUtility * Instance()
std::string to_string() const
Definition: ParameterSet.h:137
enum BeamMode string