CMFDecorrelator_plugin.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: CMFDecorrelator
3 // Plugin Type: resultsproducer (art v2_13_00)
4 // File: CMFDecorrelator_plugin.cc
5 //
6 // Generated at Thu Sep 17 14:31:24 2020 by Adam Lister using cetskelgen
7 // from cetlib version v3_06_01.
8 ////////////////////////////////////////////////////////////////////////
9 
16 #include "art_root_io/TFileService.h"
17 #include "art_root_io/TFileDirectory.h"
18 #include "fhiclcpp/ParameterSet.h"
20 
31 
32 #include "TH2.h"
33 
34 namespace cmf{
35  class CMFDecorrelator;
36 }
37 
38 
40 public:
41  explicit CMFDecorrelator(fhicl::ParameterSet const & p);
42  // The compiler-generated destructor is fine for non-base
43  // classes without bare pointers or other resource use.
44 
45  // Plugins should not be copied or assigned.
46  CMFDecorrelator(CMFDecorrelator const &) = delete;
47  CMFDecorrelator(CMFDecorrelator &&) = delete;
48  CMFDecorrelator & operator = (CMFDecorrelator const &) = delete;
50 
51  // Required functions.
52  void clear() override;
53  void writeResults(art::Results & res) override;
54 
55 private:
56 
57  int CantorPairingFunction(int k1, int k2);
58 
60 
61  void SetTemplateMatrix();
62 
63  void FillHistogramFromMatrix(Eigen::MatrixXd matrix, TH1D* hist);
64 
65  void FillHistogramFromVector(cmf::Spectrum const& spectrum, TH1D* hist);
66 
67  void FillFractionalHistograms(TH1D* ps, TH1D* init, TH1D* decor);
68 
69  void DecorrelateAgainstComponent(Eigen::MatrixXd matrix,
70  int nComponentBins,
71  int componentOffset);
72 
73  void AddEpsilonToDiagonal(Eigen::MatrixXd& matrix);
74 
75  Eigen::MatrixXd CycleMatrixComponents(Eigen::MatrixXd thisMat);
76 
78 
79  fhicl::ParameterSet fManipulatorParameters; ///< event list manipulator parameters
80  fhicl::ParameterSet fSAWParameters; ///< shifter and weighter parameters
81  bool fSystPlusStat; ///< include stats in decorrelation
82  cmf::Location loc; ///< parameter location
83 
84  TH2D* matrixHist = nullptr;
85  TH1D* initialError = nullptr;
86  TH1D* decorrelatedError = nullptr;
87  TH1D* predictedSpectra = nullptr;
88  TH1D* initialFractionalError = nullptr;
89  TH1D* decorrelatedFractionalError = nullptr;
90 
91  Eigen::MatrixXd covarianceMatrix;
92  Eigen::MatrixXd templateMatrix;
93  Eigen::MatrixXd decorrelatedMatrix;
94  Eigen::MatrixXd E_11;
95  Eigen::MatrixXd E_21;
96  Eigen::MatrixXd E_12;
97  Eigen::MatrixXd E_22;
98  Eigen::MatrixXd E_bar;
99 
100  std::map<int, int> componentInformation;
101  std::vector< cmf::SubMatrixData > thisSubMatrixDataVec;
102 
103 };
104 
105 //.............................................................................
107 // :
108 // Initialize member data here.
109 {
110 
111  fManipulatorParameters = p.get< fhicl::ParameterSet >("EventListManipulator" );
112  fSAWParameters = p.get< fhicl::ParameterSet >("ShifterAndWeighterConfig");
113  fSystPlusStat = p.get< bool >("SystPlusStat", false );
114 
118  cmf::ChiSqrCalculator::Instance() -> InitializeCovarianceMatrix(p.get<fhicl::ParameterSet>("ChiSqrCalculationParameters"));
119 
121 
123  fSAWParameters);
124 
125 
126  // this just pulls nbins for the configured selections
128 
129  matrixHist = fTFS->make<TH2D>("CovarianceMatrix",
130  ";Logical Bin;Logical Bin",
131  nbins, 0, nbins,
132  nbins, 0, nbins);
133 
134  predictedSpectra = fTFS->make<TH1D>("PredictedSpectra",
135  ";Logical Bin;Predicted Events",
136  nbins, 0, nbins);
137 
138  initialError = fTFS->make<TH1D>("InitialError",
139  ";Logical Bin;Total Uncertainty",
140  nbins, 0, nbins);
141 
142  decorrelatedError = fTFS->make<TH1D>("DecorrelatedError",
143  ";Logical Bin;Total Uncertainty",
144  nbins, 0, nbins);
145 
146  initialFractionalError = fTFS->make<TH1D>("InitialFractionalError",
147  ";Logical Bin; Fractional Uncertainty",
148  nbins, 0, nbins);
149 
150  decorrelatedFractionalError = fTFS->make<TH1D>("DecorrelatedFractionalError",
151  ";Logical Bin; Fractional Uncertainty",
152  nbins, 0, nbins);
153 
154 
155 }
156 
157 //.............................................................................
159 {
160  // Implementation of required member function here.
161 }
162 
163 //.............................................................................
165  return (((k1 + k2) * (k1 + k2 + 1))/2) + k2;
166 }
167 
168 
169 //.............................................................................
171  for (auto const& dbsItr : cmf::SelectionUtility::Instance()->SelectionsToUse()){
172  for (auto const& key : dbsItr.Keys()) componentInformation[key] = cmf::CovarianceBinUtility::Instance()->KeyToOffset(key);
173  }
174 }
175 
176 //.............................................................................
178 
179  int nComponents = componentInformation.size();
180  templateMatrix.resize(nComponents, nComponents);
181 
182  int i = 0;
183  for(const auto& it_i : componentInformation){
184 
185  int j = 0;
186  for(const auto& it_j : componentInformation){
187 
189 
190  MF_LOG_DEBUG("CMFDecorrelator")
191  << "\n cpn : " << templateMatrix(i,j)
192  << "\n key i : " << it_i.first
193  << "\n off i : " << it_i.second
194  << "\n pos i : " << i
195  << "\n nbin i: " << (int)cmf::CovarianceBinUtility::Instance()->SelectionHighEdges(it_i.first).size()
196  << "\n key j :" << it_j.first
197  << "\n off j :" << it_j.second
198  << "\n pos j :" << j
199  << "\n nbin j:" << (int)cmf::CovarianceBinUtility::Instance()->SelectionHighEdges(it_j.first).size();
200 
201  SubMatrixData thisSubMatrixData(templateMatrix(i,j),
202  it_i.first,
203  it_i.second,
204  i,
205  (int)cmf::CovarianceBinUtility::Instance()->SelectionHighEdges(it_i.first).size(),
206  it_j.first,
207  it_j.second,
208  j,
210 
211  thisSubMatrixDataVec.push_back(thisSubMatrixData);
212 
213  j++;
214  }
215  i++;
216  }
217 
218 }
219 
220 //.............................................................................
221 void cmf::CMFDecorrelator::AddEpsilonToDiagonal(Eigen::MatrixXd& thisMat){
222  for (int i = 0; i < thisMat.cols(); ++i){
223  thisMat(i,i) += 0.001;
224  }
225 }
226 
227 //.............................................................................
228 Eigen::MatrixXd cmf::CMFDecorrelator::CycleMatrixComponents(Eigen::MatrixXd thisMat){
229 
230  int matrixSize = thisMat.cols();
231 
232  Eigen::MatrixXd templateMatrixCycled(matrixSize,matrixSize);
233  cmf::SubMatrixData testSubMatrixData;
234 
235  // move c_last,last to c_0,0
236  templateMatrixCycled(0,0) = thisMat(matrixSize-1,matrixSize-1);
237 
238  // correct the 0th row
239  for (int i = 1; i < thisMat.cols(); i++){
240  for (size_t iv = 0; iv < thisSubMatrixDataVec.size(); iv++){
241  if (thisSubMatrixDataVec.at(iv).CantorPairingNumber == thisMat(matrixSize-1, i-1))
242  testSubMatrixData = thisSubMatrixDataVec.at(iv);
243  }
244  templateMatrixCycled(0,i) = thisMat(matrixSize-1,i-1);
245  }
246 
247  // correct the 1th column
248  for (int i = 1; i < thisMat.cols(); i++){
249  for (size_t iv = 0; iv < thisSubMatrixDataVec.size(); iv++){
250  if (thisSubMatrixDataVec.at(iv).CantorPairingNumber == thisMat(i-1, matrixSize-1))
251  testSubMatrixData = thisSubMatrixDataVec.at(iv);
252  }
253  templateMatrixCycled(i,0) = thisMat(i-1, matrixSize-1);
254  }
255 
256  // and now the rest
257  for (int i = 1; i < templateMatrixCycled.rows(); i++){
258  for (int j = 1; j < templateMatrixCycled.cols(); j++){
259 
260  for (size_t iv = 0; iv < thisSubMatrixDataVec.size(); iv++){
261  if (thisSubMatrixDataVec.at(iv).CantorPairingNumber == thisMat(i-1, j-1))
262  testSubMatrixData = thisSubMatrixDataVec.at(iv);
263  }
264 
265  templateMatrixCycled(i,j) = thisMat(i-1,j-1);
266  }
267 
268  }
269  return templateMatrixCycled;
270 }
271 
272 //.............................................................................
274  for (int i = 0; i < h->GetNbinsX(); ++i){
275  h->SetBinContent(i+1, std::sqrt(matrix(i,i)));
276  }
277 }
278 
279 //.............................................................................
281  for (size_t i = 0; i < s.size(); ++i){
282  h->SetBinContent(i+1, s[i]);
283  }
284 }
285 
286 //.............................................................................
288  TH1D* init,
289  TH1D* decor){
290 
291  // for whatever reason, art doesn't like it when you clone TH1s, so now do
292  // this dumb thing
293 
294  for (int i = 0; i < initialFractionalError->GetNbinsX(); ++i){
295  initialFractionalError->SetBinContent(i+1, init->GetBinContent(i+1)/ps->GetBinContent(i+1));
296  decorrelatedFractionalError->SetBinContent(i+1, decor->GetBinContent(i+1)/ps->GetBinContent(i+1));
297  }
298 }
299 
300 
302  int nComponentBins,
303  int componentOffset){
304 // now decorrelate each bin not in the 0,0 sample
305 
306  int binIterator = 1;
307  for (int i = 0; i < nComponentBins; i++){
308 
309  // firstly, size submatrices to be correct size
310  // and get the submatrices.
311 
312  int totCovMatNCols = matrix.cols()-1;
313 
314  E_11.resize(1 , 1 );
315  E_21.resize(totCovMatNCols , 1 );
316  E_12.resize(1 , totCovMatNCols );
317  E_22.resize(totCovMatNCols , totCovMatNCols );
318 
319  E_11 = matrix.topLeftCorner (1 , 1 ) ;
320  E_12 = matrix.topRightCorner (1 , totCovMatNCols ) ;
321  E_21 = matrix.bottomLeftCorner (totCovMatNCols , 1 ) ;
322  E_22 = matrix.bottomRightCorner(totCovMatNCols , totCovMatNCols ) ;
323 
324  E_bar = E_11 - (E_12 * (E_22.inverse() * E_21));
325 
326  MF_LOG_DEBUG("tmp")
327  << "E_11(" << E_11.rows() << "," << E_11.cols() << ") "
328  << "E_21(" << E_21.rows() << "," << E_21.cols() << ") "
329  << "E_12(" << E_12.rows() << "," << E_12.cols() << ") "
330  << "E_22(" << E_22.rows() << "," << E_22.cols() << ") "
331  << "E_bar(" << E_bar.rows() << "," << E_bar.cols() << ")";
332 
333  decorrelatedMatrix(i+componentOffset,i+componentOffset) = E_bar(0,0);
334 
335  matrix = E_22;
336 
337  binIterator++;
338 
339  }
340 }
341 
342 //.............................................................................
344 {
345 
346  // we're pulling in the elm so we can on-the-fly change the prediction
347  // by getting the correct "full" matrix for the fractional matrix
348 
350 
351  cmf::EventListColl eventLists;
352  manipulator.Deserialize(eventLists, cmf::kMC);
353 
354  cmf::ExposureMap exposureMap = manipulator.ExposureMap();
355 
356  cmf::Spectrum spectrum;
358 
359  cmf::FillSpectrum(eventLists,
360  exposureMap,
361  spectrum,
362  count);
363 
364  cmf::ChiSqrCalculator::Instance()->SetSpectrum(spectrum, false);
365 
366  // make a matrix for the calculation
367  covarianceMatrix = Eigen::MatrixXd::Zero(spectrum.size(), spectrum.size());
368  decorrelatedMatrix = Eigen::MatrixXd::Zero(spectrum.size(), spectrum.size());
369 
370  // get the covariance matrix for only the selections we care about,
371  // and using the predicted spectra from above
372 
375 
376  // fill the histogram to the file
377  for(int r = 0; r < matrixHist->GetNbinsX(); ++r){
378  for(int c = 0; c < matrixHist->GetNbinsX(); ++c){
379  matrixHist->SetBinContent(r+1, c+1, covarianceMatrix(r,c));
380  }
381  }
382 
383  // the decorrelation "rotates" the covariance matrix so that each iteration
384  // a different component takes the 0,0 place in order to decorrelate
385  // each component against each other component.
386  //
387  // the easiest way to do that is with a template matrix
388 
390  this->SetTemplateMatrix();
391 
392  for (size_t iComponent = 0; iComponent < componentInformation.size(); ++iComponent){
393 
394  Eigen::MatrixXd oldMat(spectrum.size(), spectrum.size());
395  Eigen::MatrixXd newMat(spectrum.size(), spectrum.size());
396 
397  // figure out the number of bins for this 0,0 component
398  int nComponentBins = 0;
399  int componentOffset = 0;
400  for ( size_t iv = 0; iv < thisSubMatrixDataVec.size(); ++iv ){
401  if (thisSubMatrixDataVec.at(iv).CantorPairingNumber == templateMatrix(0,0)){
402  nComponentBins = thisSubMatrixDataVec.at(iv).NumBinsI;
403  componentOffset = thisSubMatrixDataVec.at(iv).OffsetI;
404  }
405  }
406 
407  // construct oldMat using information from the templateMatrix
408  SubMatrixData thisSubMatrixData;
409  int offseti = 0;
410  for (int itmplt = 0; itmplt < templateMatrix.rows(); itmplt++){
411 
412  int offsetj = 0;
413  for (int jtmplt = 0; jtmplt < templateMatrix.cols(); jtmplt++){
414 
415  // find the corresponding SubMatrixData
416  for (size_t k = 0; k < thisSubMatrixDataVec.size(); k++){
417  if (thisSubMatrixDataVec.at(k).CantorPairingNumber == templateMatrix(itmplt,jtmplt))
418  thisSubMatrixData = thisSubMatrixDataVec.at(k);
419  }
420 
421  MF_LOG_VERBATIM("CovarianceMatrixDecorrelator")
422  << "Found component "
423  << thisSubMatrixData.ComponentKeyI
424  << ", "
425  << thisSubMatrixData.ComponentKeyJ
426  << " --- "
427  << nComponentBins
428  << " bins";
429 
430  // now loop over the relevant bins in the
431  // covariance matrix and fill in the oldMat with the correct bins
432  int iom = 0;
433  int finalBinI = thisSubMatrixData.OffsetI + thisSubMatrixData.NumBinsI;
434  int finalBinJ = thisSubMatrixData.OffsetJ + thisSubMatrixData.NumBinsJ;
435  for (int izrm = thisSubMatrixData.OffsetI; izrm < finalBinI; izrm++){
436  int jom = 0;
437  for (int jzrm = thisSubMatrixData.OffsetJ; jzrm < finalBinJ; jzrm++){
438  oldMat(offseti+iom,offsetj+jom) = covarianceMatrix(izrm, jzrm);
439  jom++;
440  }
441  iom++;
442  }
443  offsetj += thisSubMatrixData.NumBinsJ;
444  }
445  offseti+= thisSubMatrixData.NumBinsI;
446  }
447 
448  newMat = oldMat;
449 
450  this->DecorrelateAgainstComponent(newMat, nComponentBins, componentOffset);
452  }
453 
454  // fill histograms
455  this->FillHistogramFromVector(spectrum, predictedSpectra);
456  this->FillHistogramFromMatrix(decorrelatedMatrix, decorrelatedError);
459 }
460 
std::map< cmf::MetaData, cmf::SpillSummary > ExposureMap
Definition: Structs.h:215
Eigen::MatrixXd decorrelatedMatrix
void FillSpectrum(cmf::EventListColl const &eventLists, cmf::ExposureMap const &exposureMap, cmf::Spectrum &spectrum, cmf::Spectrum &count, cmf::InteractionType_t intType, bool applyExposureNorm)
int KeyToOffset(long const &key, bool allSels=false)
void AddEpsilonToDiagonal(Eigen::MatrixXd &matrix)
void writeResults(art::Results &res) override
fhicl::ParameterSet fManipulatorParameters
event list manipulator parameters
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
static SelectionUtility * Instance()
std::vector< double > Spectrum
Definition: Constants.h:746
cmf::Location loc
parameter location
void FillHistogramFromMatrix(Eigen::MatrixXd matrix, TH1D *hist)
void Initialize(Bool_t useTMVAStyle=kTRUE)
Definition: tmvaglob.C:184
void DecorrelateAgainstComponent(Eigen::MatrixXd matrix, int nComponentBins, int componentOffset)
art::ServiceHandle< art::TFileService > fTFS
TFileService.
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}))
Eigen::MatrixXd CycleMatrixComponents(Eigen::MatrixXd thisMat)
const XML_Char * s
Definition: expat.h:262
static ParameterUtility * Instance()
const int nbins
Definition: cellShifts.C:15
std::vector< cmf::SubMatrixData > thisSubMatrixDataVec
void FillHistogramFromVector(cmf::Spectrum const &spectrum, TH1D *hist)
#define DEFINE_ART_RESULTS_PLUGIN(klass)
void FillFractionalHistograms(TH1D *ps, TH1D *init, TH1D *decor)
T get(std::string const &key) const
Definition: ParameterSet.h:231
fhicl::ParameterSet fSAWParameters
shifter and weighter parameters
cmf::ExposureMap const & ExposureMap() const
const double j
Definition: BetheBloch.cxx:29
bool fSystPlusStat
include stats in decorrelation
static ShifterAndWeighter * Instance()
CMFDecorrelator(fhicl::ParameterSet const &p)
static cmf::ChiSqrCalculator * Instance()
size_t TotalBins(bool allSels=false)
int CantorPairingFunction(int k1, int k2)
Module to combine a set of results into a single file currently only does one data product type at a ...
Definition: Event.cxx:24
std::vector< cmf::EventList > EventListColl
Definition: Event.h:150
void InitShiftsAndWeightsToUse(cmf::Location const &loc, fhicl::ParameterSet const &parset)
std::map< double, int > const & SelectionHighEdges(cmf::MetaData const &md)
void Zero()
#define MF_LOG_VERBATIM(category)
Eigen::MatrixXd covarianceMatrix
#define MF_LOG_DEBUG(id)
cmf::Location ParametersAsLocation()
TRandom3 r(0)
std::map< int, int > componentInformation
CMFDecorrelator & operator=(CMFDecorrelator const &)=delete
void FillCovarianceMatrix(Eigen::MatrixXd &matrix, bool isSystOnly=false)
void SetSpectrum(cmf::Spectrum const &spectrum, bool isData)
static CovarianceBinUtility * Instance()