CovarianceFitHelper.cxx
Go to the documentation of this file.
1 //
2 // Created by Brian Rebel on 12/15/18.
3 //
4 
5 #include <sstream>
6 
10 
11 
12 // ROOT includes
13 #include "TCanvas.h"
14 #include "TFile.h"
15 #include "TGraph.h"
16 #include "TGraphErrors.h"
17 #include "TH1.h"
18 #include "TH2.h"
19 #include "TLegend.h"
20 
21 namespace fnex{
22 
23  //----------------------------------------------------------------------------
25  : fBinUtil(std::make_unique<fnex::CovarianceBinUtility>())
26  , fCovarianceMatrixStat(nullptr)
27  , fCovarianceMatrixSyst(nullptr)
28  , fCovarianceMatrixHist(nullptr)
29  , fIterations(std::make_unique<std::vector<fnex::CovFitIteration>>() )
30  , fPrevDiffs (std::make_unique<std::vector<double> >(378, 0.))
31  , fKeysToUse (std::make_unique<std::set<int> >() )
32 
33  {
34  }
35 
36  //----------------------------------------------------------------------------
38  {
39  fInvCovarianceMatrix.reset();
40  fRowDiffMatrix .reset();
41  fColDiffMatrix .reset();
42  fInputPoint .reset();
43  fDataSpectrum .reset();
44  fMCSpectrum .reset();
45  fIterations .reset();
46 
47  }
48 
49  //----------------------------------------------------------------------------
51  {
52  fEventLists.clear();
53  }
54 
55  //----------------------------------------------------------------------------
57  {
58  // create the manipulator
59  auto manipulatorParameters = pset.get<fhicl::ParameterSet>("EventListManipulator" );
60  auto inputPointParameters = pset.get<fhicl::ParameterSet>("InitialGuessInputPoint");
61 
62  // set up the initial guess for the fit
63  // the input point should have no systematic shifts except for the
64  // ones that are always on and fixed
65  fInputPoint = std::make_unique<fnex::InputPoint>(inputPointParameters);
66 
67  // make the EventListManipulator
68  fManipulator = std::make_unique<fnex::EventListManipulator>(manipulatorParameters,
69  *fInputPoint);
70 
71  // get data samples to use from the manipulatorParameters
72  int key = 0;
73  int keyBase = 0;
74  std::string bmStr;
75  std::string detStr;
76 
77  for(auto const& itr : manipulatorParameters.get<std::vector<fhicl::ParameterSet>>("SelectionsToUse")){
78  keyBase = 0;
79 
80  bmStr = itr.get<std::string>("Beam");
81  detStr = itr.get<std::string>("Detector");
82 
83  keyBase += (bmStr.find("FHC") != std::string::npos) ? 1000 * (1 + fnex::kFHC) : 1000 * (1 + fnex::kRHC);
84  keyBase += (detStr.find("Far") != std::string::npos) ? 100 * novadaq::cnv::kFARDET : 100 * novadaq::cnv::kNEARDET;
85 
86  auto selStrs = itr.get<std::vector<std::string>>("Selections");
87  for(auto& sel: selStrs){
88  LOG_DEBUG("CovarianceFitHelper")
89  << bmStr
90  << " "
91  << detStr
92  << " "
93  << sel
94  << " "
95  << keyBase;
96 
97  if (sel.find("NuESel_Peripheral") != std::string::npos) key = *(fKeysToUse->insert(keyBase + fnex::kNuESelectionPeripheral).first);
98  else if(sel.find("NuESel_LowPID") != std::string::npos) key = *(fKeysToUse->insert(keyBase + fnex::kNuESelectionLowPID ).first);
99  else if(sel.find("NuESel_HighPID") != std::string::npos) key = *(fKeysToUse->insert(keyBase + fnex::kNuESelectionHighPID ).first);
100  else if(sel.find("NuMuSelQ1") != std::string::npos) key = *(fKeysToUse->insert(keyBase + fnex::kNuMuSelectionQ1 ).first);
101  else if(sel.find("NuMuSelQ2") != std::string::npos) key = *(fKeysToUse->insert(keyBase + fnex::kNuMuSelectionQ2 ).first);
102  else if(sel.find("NuMuSelQ3") != std::string::npos) key = *(fKeysToUse->insert(keyBase + fnex::kNuMuSelectionQ3 ).first);
103  else if(sel.find("NuMuSelQ4") != std::string::npos) key = *(fKeysToUse->insert(keyBase + fnex::kNuMuSelectionQ4 ).first);
104 
105  LOG_DEBUG("CovarianceFitHelper")
106  << "use key "
107  << key
108  << " in fit";
109 
110  } // end loop over selection strings
111 
112  } // end loop over samples
113 
114  fStatOnlyFit = pset.get<bool>("StatOnlyFit", false);
115 
116  // get the information for where the covariance matrix is stored
117  // assume that the matrices for individual systematics have already
118  // been summed by this point
119  fCovarianceHistName = pset.get<std::string>("CovarianceMatrixHistName", "totalCovariance");
120 
121  // xrootd path to file in /pnfs/nova/scratch/
122  fCovarianceHistFileName = pset.get<std::string>("CovarianceMatrixHistFile");
123 
124  // the Nuisance Parameters were put into a singleton when the InputPoint was made
125 
126  return;
127 
128  }
129 
130  //----------------------------------------------------------------------------
132  {
133  // get the covariance matrix from the input file - it is assumed to already be summed
134  // for the individual systematic uncertainties.
135 
136  // Load the covariance matrix from the input histogram file.
137  // Assume that the matrix has been properly normalized for the total
138  // number of shifts tested before we get it
139  TFile *covHistFile = TFile::Open(fCovarianceHistFileName.c_str());
140  TH2D* inputCovarianceSyst = dynamic_cast<TH2D*>(covHistFile->Get(fCovarianceHistName.c_str()));
141 
142  // here we need to check which keys we have in the map of s
143 
144  int numBins = inputCovarianceSyst->GetNbinsX();
145 
146  LOG_DEBUG("CovarianceFitHelper")
147  << "number of bins: "
148  << numBins;
149 
150  fCovarianceMatrixSyst = fTFS->make<TH2D>("CovarianceMatrixSyst",
151  ";Energy Bin;Energy Bin",
152  numBins,
153  1,
154  numBins + 1,
155  numBins,
156  1,
157  numBins + 1);
158 
159  fCovarianceMatrixSyst->Add(inputCovarianceSyst);
160 
161  covHistFile->Close();
162 
163  fCovarianceMatrixSystScaled = fTFS->make<TH2D>("CovarianceMatrixSystScaled",
164  ";Energy Bin;Energy Bin",
165  numBins,
166  1,
167  numBins + 1,
168  numBins,
169  1,
170  numBins + 1);
171 
172  // make a covariance matrix for the statistical covariance
173  // and one for the combined syst + stat
174  fCovarianceMatrixStat = fTFS->make<TH2D>("CovarianceMatrixStat",
175  ";Energy Bin;Energy Bin",
176  numBins,
177  1,
178  numBins + 1,
179  numBins,
180  1,
181  numBins + 1);
182 
183  fCovarianceMatrixHist = fTFS->make<TH2D>("CovarianceMatrixTot",
184  ";Energy Bin;Energy Bin",
185  numBins,
186  1,
187  numBins + 1,
188  numBins,
189  1,
190  numBins + 1);
191 
192  fMCSpectrum = std::make_unique<std::vector<double>>(numBins, 0.);
193  fDataSpectrum = std::make_unique<std::vector<double>>(numBins, 0.);
194 
195  LOG_DEBUG("CovarianceFitHelper")
196  << "size of spectra: "
197  << fMCSpectrum->size()
198  << " "
199  << fDataSpectrum->size();
200 
201 
202  // make histograms to store the MC and data spectrum in energy for each beam/detector/selection/file type
203  // first make vectors of the bin boundaries
204  std::vector<double> numuBins(1, 0.);
205  std::vector<double> nueBins (1, 0.);
206  auto numuHighEdge = fBinUtil->NuMuHighEdges();
207  auto nueHighEdge = fBinUtil->NuEHighEdges();
208 
209  for(auto const& itr : numuHighEdge) numuBins.push_back(itr.first);
210  for(auto const& itr : nueHighEdge ) nueBins.push_back(itr.first);
211 
212  // for(auto const & itr : numuBins)
213  // LOG_VERBATIM("CovarianceFitHelper")
214  // << "numu bins: "
215  // << itr;
216 
217  // for(auto const & itr : nueBins)
218  // LOG_VERBATIM("CovarianceFitHelper")
219  // << "nue bins: "
220  // << itr;
221 
222  fDataSpectraE = std::make_unique<std::map<int, TH1D*>>();
223  fMCSpectraE = std::make_unique<std::map<int, TH1D*>>();
224 
225  std::string beamDetFileName;
227  int key = 0;
228  std::set <fnex::FileType_t> fileTypes({fnex::kDataFile,
229  fnex::kBeam,
230  fnex::kSwap,
233  for(int b = 0; b < 2; ++b){
234  for(int d = 1; d < 3; ++d){
235 
236  // loop over file types
237  for(auto fitr : fileTypes){
238 
239  // in the ND, we only care about the nonswap file type
240  if(d == 0 &&
241  (fitr != fnex::kBeam && fitr != fnex::kDataFile)) continue;
242 
243  beamDetFileName = (fnex::cBeamType_Strings[b] +
245  fnex::kFileTypeStrings[fitr]);
246 
248 
249  // MidPID selection is not used after 2017
250  // We only have peripheral samples in the FD
252  (e == fnex::kNuESelectionPeripheral && d == 0))
253  continue;
254 
255  key = (b + 1) * 1000 + d * 100 + e + (fitr + 1) * 10 ;
256 
257  if(fitr != fnex::kDataFile){
258 
259  histName = "MCSpectraEnergy" + beamDetFileName + fnex::cSelectionType_Strings[e];
260 
261  fMCSpectraE->emplace(key, fTFS->make<TH1D>(histName.c_str(),
262  ";Energy (GeV);Events",
263  nueBins.size() - 1,
264  &nueBins[0]));
265  } // end if MC file
266  else{
267 
268  histName = "DataSpectraEnergy" + beamDetFileName + fnex::cSelectionType_Strings[e];
269 
270  fDataSpectraE->emplace(key, fTFS->make<TH1D>(histName.c_str(),
271  ";Energy (GeV);Events",
272  nueBins.size() - 1,
273  &nueBins[0]));
274  } // end if data file
275 
276  } // end loop over nue selections
277 
278  // now the numu bins
279  for(size_t m = fnex::kNuMuSelectionQ1; m < fnex::kNuMuSelectionQ4 + 1; ++m){
280 
281  key = (b + 1) * 1000 + d * 100 + m + (fitr + 1) * 10;
282 
283  if(fitr != fnex::kDataFile){
284 
285  histName = "MCSpectraEnergy" + beamDetFileName + fnex::cSelectionType_Strings[m];
286 
287  fMCSpectraE->emplace(key, fTFS->make<TH1D>(histName.c_str(),
288  ";Energy (GeV);Events",
289  numuBins.size() - 1,
290  &numuBins[0]));
291  } // end if MC file
292  else{
293 
294  histName = "DataSpectraEnergy" + beamDetFileName + fnex::cSelectionType_Strings[m];
295 
296  fDataSpectraE->emplace(key, fTFS->make<TH1D>(histName.c_str(),
297  ";Energy (GeV);Events",
298  numuBins.size() - 1,
299  &numuBins[0]));
300  } // end if data file
301 
302  } // end loop over numu selections
303 
304  } // end loop over file types
305  } // end loop over detectors
306  } // end loop over beam types
307 
308  // for(auto const& itr : fDataSpectraE)
309  // LOG_VERBATIM("CovarianceMatrixFitHelper")
310  // << "Data spectra map key: "
311  // << itr.first
312  // << " "
313  // << itr.second;
314 
315  // for(auto const& itr : fMCSpectraE)
316  // LOG_VERBATIM("CovarianceMatrixFitHelper")
317  // << "MC spectra map key: "
318  // << itr.first
319  // << " "
320  // << itr.second;
321 
322  LOG_DEBUG("CovarianceMatrixFitHelper")
323  << "histograms made";
324 
325  // The shifter and weighter stuff was initialized in the EventListManipulator::Deserialize method
326 
327  return;
328  }
329 
330  // Function to fill the spectra for each meta data type
331  // only fill the data once, don't bother looking at data
332  // event lists if the fillData argument is false
333  // Also, use fillData to determine if we need to normalize the
334  // the event lists to the data POT or not
335  //......................................................................
336  void CovarianceFitHelper::FillSpectrum(std::vector<double> & spectrum,
337  bool data)
338  {
339  if(data) fDataPOT.clear();
340 
341  // first reset the spectrum
342  spectrum.assign(spectrum.size(), 0.);
343 
344  // declare some variables we will be using a lot in the for loops below
345  double weight = 1.;
346  double energy = 0.;
347  double norm = 1.;
348  int bin = 0;
349  int offset = 0;
350  int key = 0;
351 
352  // loop over each list and event in the list
353  for(auto const& itr : fEventLists){
354  auto const md = itr.first;
355 
356  // check that the metadata and the requested
357  // spectrum to fill are consistent
358  if ( data && md.isMC) continue;
359  else if(!data && !md.isMC) continue;
360 
361  // check if we should be using this metadata
362  key = fBinUtil->MetaDataToKey(md);
363  if(fKeysToUse->find(key) == fKeysToUse->end()) continue;
364 
365  // we want to use these events, so determine the offset.
366  offset = fBinUtil->KeyToOffset(key);
367 
368  if(!md.isMC){
369  // need to include an epoch component to the
370  // key for the data because the MC is done by
371  // epoch. MetaData::Epoch is the period * 1000,
372  // so another factor of 10 is enough to get us
373  // unique keys
374  fDataPOT[key + md.epoch * 10] = itr.second.ListSpillSummary().goodPOT;
375  norm = 1.;
376  }
377  else{
378  norm = fDataPOT[key + md.epoch * 10]/itr.second.ListSpillSummary().goodPOT;
379  }
380 
381  LOG_DEBUG("CovarianceMatrixFitter")
382  << md.ToString()
383  << " key: "
384  << key
385  << " bin offset: "
386  << offset
387  << " metadata fileType "
388  << md.fileType
389  << " isMC "
390  << md.isMC;
391 
392  for(auto const& evt : itr.second){
393 
394  energy = evt->dataVals->val_at(fnex::kNu_RecoE, md);
395  bin = fBinUtil->EnergyToBin(energy, md);
397 
398  // now fill the spectrum with the proper weight and
399  // normalization
400  spectrum[bin + offset - 1] += weight * norm;
401 
402  } // end loop over events
403  } // end loop over lists of events
404 
405  return;
406  }
407 
408  //......................................................................
410  {
411 
412  // first get all the data events
413  fEventLists = fManipulator->DeserializeData();
414 
415  this->FillSpectrum(*fDataSpectrum, true);
416 
417  // the data spectrum is filled, so drop all the data events
418  std::set<fnex::MetaData> dataMD;
419  for(auto const& itr : fEventLists) if(!itr.first.isMC) dataMD.insert(itr.first);
420 
421  // fill the energy spectra for the data
422  this->FillEnergySpectra(dataMD);
423 
424  for(auto const& itr : dataMD) fEventLists.erase(itr);
425 
426  return;
427  }
428 
429  //......................................................................
431  {
432  // first check if we have any event lists, deserialize them
433  // if we don't
434  if(fEventLists.size() < 1){
435  fManipulator->DeserializeMC();
436  }
437 
438  // check that the nominal spectrum has been filled
439  if(fMCSpectrumNominal.get() == nullptr){
440  fMCSpectrumNominal = std::make_unique<std::vector<double>>(fMCSpectrum->size(), 0.);
442 
443  // the imported systematic histogram is a relative histogram and needs to have
444  // the entries scaled by the number of events in the i and j bins
446 
447  double count = 0.;
448 
449  for(size_t i = 0; i < fMCSpectrumNominal->size(); ++i){
450  for(size_t j = 0; j < fMCSpectrumNominal->size(); ++j){
451 
452  count = (*fMCSpectrumNominal)[i] * (*fMCSpectrumNominal)[j] * fCovarianceMatrixSyst->GetBinContent(i + 1, j + 1);
453 
454  fCovarianceMatrixSystScaled->SetBinContent(i + 1, j + 1, count);
455 
456  } // end loop over j bins
457  } // end loop over i bins
458  }// end if we need to make the nominal spectrum and fill it
459 
460  // now fill the MC spectrum with the current point
461  this->FillSpectrum(*fMCSpectrum);
462 
463  return;
464  }
465 
466  // Function to make TH1D from the data and MC 1D spectra vectors
467  // Ideally this function is called after the fit is finished and
468  // the fit values have been set in the InputPoint
469  //......................................................................
471  {
472  size_t numBins = fDataSpectrum->size();
473 
474  auto data1DHist = fTFS->make<TH1D>("dataSpectrumBins",
475  ";Logical Bin;Events",
476  numBins,
477  0,
478  numBins);
479 
480  auto mc1DHist = fTFS->make<TH1D>("mcSpectrumBins",
481  ";Logical Bin;Events",
482  numBins,
483  0,
484  numBins);
485 
486  for(size_t b = 0; b < numBins; ++b){
487  data1DHist->SetBinContent(b+1, (*fDataSpectrum)[b]);
488  mc1DHist ->SetBinContent(b+1, (*fMCSpectrum)[b]);
489  }
490 
491  return;
492  }
493 
494  // This function extracts the energy spectrum for each MetaData type and
495  // puts it into a histogram
496  //......................................................................
497  void CovarianceFitHelper::FillEnergySpectra(std::set<fnex::MetaData> const& mdToFill)
498  {
499  int offset = 0;
500  int key = 0;
501  size_t binsForMD = 0;
502  double energy = 0.;
503  double events = 0.;
504 
505  // check the first entry in the set to determine if we are filling
506  // MC or data spectra
507  bool isMC = mdToFill.begin()->isMC;
508 
509  auto & spectMap = (isMC) ? *fMCSpectraE : *fDataSpectraE;
510  auto const& spectrum = (isMC) ? *fMCSpectrum : *fDataSpectrum;
511 
512  // first clear out the spectra histograms
513  for(auto const& itr : spectMap ) itr.second->Reset();
514 
515  // the event lists only have the MC lists at this point because
516  // we dropped the data ones to conserve memory early on
517  for(auto const& md : mdToFill){
518 
519  key = fBinUtil->MetaDataToKey(md);
520  offset = fBinUtil->KeyToOffset(key);
521 
522  key += 10 * (md.fileType + 1);
523 
524  binsForMD = (md.IsNuESelected()) ? fBinUtil->NuEHighEdges().size() : fBinUtil->NuMuHighEdges().size();
525  if(md.selectionType == fnex::kNuESelectionPeripheral) binsForMD = 1;
526 
527  // we only want to fill the spectrum for each key once - given the different interaction
528  // types we could end up filling the same information multiple times because we are
529  // filling from the total spectrum
530  // check to see if the histogram for this key has a non-zero integral and if so,
531  // skip it
532  if(spectMap[key]->Integral() > 0.) continue;
533 
534  // loop over the spectrum to get the bins
535  for(size_t b = (size_t)offset; b < offset + binsForMD; ++b){
536  energy = fBinUtil->BinToEnergy(b, md);
537  events = spectrum[b];
538 
539  spectMap[key]->Fill(energy, events);
540 
541  LOG_DEBUG("CovarianceFitHelper")
542  << md.ToString()
543  << " filling spectrum for key "
544  << key
545  << " energy "
546  << energy
547  << " events "
548  << events
549  << " total events "
550  << spectMap[key]->GetBinContent(spectMap[key]->FindBin(energy));
551  }
552  } // end loop over event lists
553 
554  return;
555  }
556 
557  //......................................................................
559  {
560  // fill the statistical covariance matrix only if it does not already
561  // have entries.
562  if(fCovarianceMatrixStat->GetEntries() < 1){
563  LOG_DEBUG("CovarianceFitHelper")
564  << "data spectrum size: "
565  << fDataSpectrum->size();
566 
567  double count = 0.;
568 
569  // for(auto const& itr : *fDataSpectrum)
570  // LOG_VERBATIM("CovarianceFitHelper")
571  // << itr;
572 
573  // fill the statistical covariance matrix, then combine it with the systematic one
574  // if we have 0 events in a bin, then the uncertainty is given by the 90% CL
575  // from the Feldman & Cousins paper, ie 1
576  for(size_t b = 0; b < fDataSpectrum->size(); ++b){
577  count = ((*fDataSpectrum)[b] != 0) ? (*fDataSpectrum)[b] : 1.;
578  fCovarianceMatrixStat->SetBinContent(b, b, count);
579  }
580  } // end if we need to fill the statistical covariance matrix
581 
582  // next iteration of the fit, have to reset the histogram.
583  fCovarianceMatrixHist->Reset();
584 
585  // the stat only fit is easy
587  else{
589  }
590 
591  // now take the covariance histogram and turn it into a matrix
592  TMatrixD covmat = TMatrixD(fCovarianceMatrixHist->GetNbinsX() + 2,
593  fCovarianceMatrixHist->GetNbinsY() + 2,
594  fCovarianceMatrixHist->GetArray(),
595  "F");
596 
597  covmat.ResizeTo(0, 377, 0, 377);
598 
599  // LOG_VERBATIM("CovarianceFitHelper")
600  // << "Print the matrix from the histogram: ";
601  // covmat.Print();
602 
603  // invert the matrix because that is really what we want
604  fInvCovarianceMatrix.reset();
605  fInvCovarianceMatrix = std::make_unique<TMatrixD>(covmat.Invert());
606 
607  // LOG_VERBATIM("CovarianceFitHelper")
608  // << "Print the inverse covariance matrix: ";
609  // fInvCovarianceMatrix->Print();
610 
611  return;
612  }
613 
614  //......................................................................
616  {
617  // get the difference of the data and MC stored in a vector
618  std::vector<double> diffs(fDataSpectrum->size(), 0.);
619 
620  for(size_t i = 0; i < diffs.size(); ++i){
621  diffs[i] = (*fDataSpectrum)[i] - (*fMCSpectrum)[i];
622  LOG_DEBUG("CovarianceFitHelper")
623  << i
624  << " data "
625  << (*fDataSpectrum)[i]
626  << " mc "
627  << (*fMCSpectrum)[i]
628  << " delta diffs: "
629  << diffs[i]
630  << " "
631  << (*fPrevDiffs)[i]
632  << " "
633  << diffs[i] - (*fPrevDiffs)[i];
634  (*fPrevDiffs)[i] = diffs[i];
635  }
636 
637  // now make a row matrix out of the diffs and a column matrix
638  // TMatrixD rowDiffMatrix(1, diffs.size(), &diffs[0]);
639  // TMatrixD colDiffMatrix(diffs.size(), 1, &diffs[0]);
640 
641  if(!fRowDiffMatrix){
642  fRowDiffMatrix = std::make_unique<TMatrixD>(1, diffs.size(), &diffs[0]);
643  fColDiffMatrix = std::make_unique<TMatrixD>(diffs.size(), 1, &diffs[0]);
644 
645  // fRowDiffMatrix->Print();
646  // fColDiffMatrix->Print();
647  }
648  else{
649  fRowDiffMatrix->SetMatrixArray(&diffs[0]);
650  fColDiffMatrix->SetMatrixArray(&diffs[0]);
651  }
652 
653  // apparently in place multiplication only works when the matrix is square, so
654  // fColDiffMatrix is not square and we have to do it this way
655  auto intermed = (*fInvCovarianceMatrix) * (*fColDiffMatrix);
656 
657  //(*fInvCovarianceMatrix).GetSub(98, 183, 98, 183).Print();
658  // Print out the matrix for the numu FD neutrino section
659  // as a csv
660  // for(int r = 98; r < 184; ++r){
661  // std::stringstream covmatss;
662  // for(int c = 98; c < 184; ++c){
663  // covmatss << (*fInvCovarianceMatrix)(r, c);
664  // if(c < 183) covmatss << ", ";
665  // }
666  // LOG_VERBATIM("CovarianceFitHelper")
667  // << covmatss.str();
668  // }
669 
670 
671  // LOG_VERBATIM("CovarianceFitHelper")
672  // << "number of intermediate matrix elements is "
673  // << intermed.GetNoElements();
674  // intermed.Print();
675 
676  auto chiSqrMat = (*fRowDiffMatrix) * intermed;
677 
678  // LOG_VERBATIM("CovarianceFitHelper")
679  // << "number of chi^2 matrix elements is "
680  // << chiSqrMat.GetNoElements();
681 
682  // chiSqrMat.Print();
683 
684  // we want to get any chi^2 contribution from the nuisance parameters now
685  double chiSqrNP = 0.;
686  double curNPVal = 0.;
687 
688  LOG_DEBUG("CovarianceFitHelper")
689  << "there are "
691  << " nuisance parameters";
692 
693  for(auto const& itr : fnex::NuisanceParameters::Instance()->NuisanceParameterInfo()){
694 
695  curNPVal = fInputPoint->ParameterValue(std::make_pair(itr.first, novadaq::cnv::kFARDET));
696  chiSqrNP += std::pow((curNPVal - itr.second.fCentralValue) / itr.second.fSigma, 2.);
697 
698  LOG_DEBUG("CovarianceFitHelper")
699  << "nuisance parameter "
700  << itr.first
701  << " adds "
702  << std::pow((curNPVal - itr.second.fCentralValue) / itr.second.fSigma, 2.)
703  << " to the chi^2, current value "
704  << curNPVal
705  << " central value "
706  << itr.second.fCentralValue
707  << " sigma "
708  << itr.second.fSigma;
709 
710  }
711 
712  // the chiSqrMat should be a 1 x 1 matrix, so the sum is just that single entry
713  return chiSqrMat.Sum() + chiSqrNP;
714  }
715 
716  //......................................................................
718  {
719  *fInputPoint = point;
720 
721  return;
722  }
723 
724  //......................................................................
726  {
727  // create histograms for the fit parameters and the chi^2 to see
728  // how they changed during the minimization
729  auto oscParDets = (*fIterations)[0].point.FitParameters();
730 
731  // make a collection of TGraphs for the oscillation parameters
733  std::map<ParameterDet, TGraph*> grMap;
734  for(auto itr : oscParDets){
735 
736  // do nothing if this parameter is not in the FD
737  if( itr.second != fd ) continue;
738 
739  grMap[itr] = fTFS->make<TGraph>(fIterations->size());
740  grMap[itr]->SetName((itr.first + "IterGraph").c_str());
741  grMap[itr]->SetTitle(itr.first.c_str());
742  grMap[itr]->SetMarkerStyle(20);
743  }
744 
745  // Make the TGraph for the chi^2 values
746  TGraph* chisq_graph = fTFS->make<TGraph>(fIterations->size());
747  chisq_graph->SetName("IterGraphChiSqr");
748  chisq_graph->SetTitle("#chi^{2}");
749  chisq_graph->SetMarkerStyle(20);
750 
751  //loop over the iterations to fill the graphs
752  // the last value in fIterations is not the minimum, it is
753  // the value that Minuit used to confirm it had found a minimum,
754  // so don't use it in the graphs
755  double parVal = 0.;
756  for(size_t i = 0; i < fIterations->size() - 1; ++i){
757 
758  // get the value for each parameter we care about and put it into its
759  // graph
760  for(auto mapItr : grMap){
761  parVal = (*fIterations)[i].point.ParameterValue(mapItr.first);
762 
763  mapItr.second->SetPoint(i, i * 1., parVal);
764  }
765 
766  // fill the chi^2 graph
767  chisq_graph->SetPoint(i, i * 1., (*fIterations)[i].chiSqr);
768  } // end loop over iterations
769 
770  // now loop over the graphs and create a canvas for each
771  TCanvas* can = nullptr;
772 
773  std::string canvasname;
774  for(auto mapItr : grMap){
775  canvasname = "ParameterIterationCanv_" + mapItr.first.first;
776  can = fTFS->make<TCanvas>(canvasname.data(), canvasname.data(), 800, 800);
777 
778  mapItr.second->Draw("ACP");
779 
780  // Save the canvas to file
781  can->Update();
782  can->Write();
783 
784  }//end of loop over canvases
785 
786  // now for the chi^2 iterations
787  can = fTFS->make<TCanvas>("chisq_can", "chisq", 800, 800);
788  chisq_graph->Draw("ACP");
789 
790  can->Update();
791 
792  // Save the canvas to file
793  can->Write();
794 
795  return;
796  }
797 
798  //......................................................................
800  {
801  // normalize the data and MC histogram bin contents to be equivalent to 1 GeV wide bins
802  for(auto itr : *fDataSpectraE){
803  for(int b = 0; b < itr.second->GetNbinsX(); ++b){
804  itr.second->SetBinContent(b+1, itr.second->GetBinContent(b+1)/itr.second->GetBinWidth(b+1));
805  }
806  }
807 
808  for(auto itr : *fMCSpectraE){
809  for(int b = 0; b < itr.second->GetNbinsX(); ++b){
810  itr.second->SetBinContent(b+1, itr.second->GetBinContent(b+1)/itr.second->GetBinWidth(b+1));
811  }
812  }
813 
814  return;
815  }
816 
817  //......................................................................
819  {
820  this->NormalizeHistograms();
821 
822  // make 1 canvas each for the numu and nue data/MC comparison
823  // The spectra maps have what we want and we can use the keys
824  // to determine how to make the canvases
825 
826  auto tfd = fTFS->mkdir("results");
827 
828  // for now we need 12 pads for the nue in the ND and FD, FHC and RHC
829  // arrange them ND above FD, FHC grouped about RHC
830 
831  std::vector<fnex::SelectionType_t> numuSels = {fnex::kNuMuSelectionQ1, fnex::kNuMuSelectionQ2, fnex::kNuMuSelectionQ3, fnex::kNuMuSelectionQ4};
832  std::vector<fnex::SelectionType_t> nueSels = {fnex::kNuESelectionLowPID, fnex::kNuESelectionHighPID, fnex::kNuESelectionPeripheral};
833 
834  std::set <fnex::FileType_t> fileTypes({fnex::kBeam,
835  fnex::kSwap,
838 
839  // make vectors of the bin boundaries for the new histogram
840  std::vector<double> numuBins(1, 0.);
841  for(auto const& itr : fBinUtil->NuMuHighEdges()){
842  numuBins.push_back(itr.first);
843  LOG_DEBUG("CovarianceFitHelper")
844  << "bins for total MC and ratios "
845  << numuBins.back()
846  << " "
847  << numuBins.size();
848  }
849 
850  // make vectors of the bin boundaries for the new histogram
851  std::vector<double> nueBins(1, 0.);
852  for(auto const& itr : fBinUtil->NuEHighEdges()){
853  nueBins.push_back(itr.first);
854  LOG_DEBUG("CovarianceFitHelper")
855  << "bins for total MC and ratios "
856  << nueBins.back()
857  << " "
858  << nueBins.size();
859  }
860 
861  int canvDiv = 0;
862  std::string nameBase;
863  for(int bt = 0; bt < 2; ++bt){
864  for(int d = 1; d < 3; ++d){
865 
866  for(int s = 0; s < 4; ++s){
867 
869 
870  // for now we need 16 pads for the numu in the ND and FD, FHC and RHC
871  // arrange them ND above FD, FHC grouped about RHC
872  TCanvas *numuCanv = tfd.make<TCanvas>(("Numu" + nameBase + "FitCanvas").c_str(), "#nu_{#mu} Canvas", 600, 600);
873  //numuCanv->Divide(4, 4);
874 
875  canvDiv = (bt * 8) + (d - 1) * 4 + s + 1;
876 
877  int keyBase = (bt + 1) * 1000 + d * 100 + numuSels[s];
878 
879  auto dataHist = (*fDataSpectraE)[keyBase + (fnex::kDataFile + 1) * 10];
880 
881 
882  // make the histogram
883  TH1D *totMC = tfd.make<TH1D>(("MC" + nameBase + "Summed").c_str(),
884  ";Energy (Gev);Events",
885  numuBins.size() - 1,
886  &numuBins[0]);
887 
888  TH1D *ratio = tfd.make<TH1D>(("DataMCRatio" + nameBase).c_str(),
889  ";Energy (Gev);Ratio",
890  numuBins.size() - 1,
891  &numuBins[0]);
892 
893  totMC->GetXaxis()->CenterTitle();
894  totMC->GetYaxis()->CenterTitle();
895  totMC->SetLineColor(kRed);
896  totMC->SetMarkerColor(kRed);
897 
898  ratio->GetXaxis()->CenterTitle();
899  ratio->GetXaxis()->SetDecimals();
900  ratio->GetYaxis()->CenterTitle();
901  ratio->GetYaxis()->SetDecimals();
902  ratio->GetXaxis()->SetLabelSize(0.1);
903  ratio->GetYaxis()->SetLabelSize(0.1);
904 
905  // fill it
906  for(auto const& itr : fileTypes){
907  totMC->Add((*fMCSpectraE)[keyBase + (itr + 1) * 10]);
908  LOG_DEBUG("CovarianceFitHelper")
909  << (*fMCSpectraE)[keyBase + (itr + 1) * 10]->GetNbinsX()
910  << " "
911  << keyBase + (itr + 1) * 10;
912  }
913 
914  ratio->Divide(dataHist, totMC);
915 
916  double max = std::max(dataHist->GetMaximum(), totMC->GetMaximum());
917  totMC->SetMaximum(1.2 * max);
918 
919  //numuCanv->cd(canvDiv);
920 
921  LOG_DEBUG("CovarianceFitHelper")
922  << "current subcanvas is "
923  << canvDiv;
924 
925  // make a couple of pads in the canvas
926  TPad *padUpper = tfd.make<TPad>(("Upper" + nameBase).c_str(), ("Upper" + nameBase).c_str(), 0, 0.3, 1, 1.0);
927  padUpper->SetBottomMargin(0.016);
928  padUpper->SetBorderSize(0);
929 
930  TPad *padLower = tfd.make<TPad>(("Lower" + nameBase).c_str(), ("Lower" + nameBase).c_str(), 0, 0.05, 1, 0.3);
931  padLower->SetTopMargin(0);
932  padLower->SetBottomMargin(0.25);
933  padLower->SetBorderSize(0);
934 
935  padUpper->Draw();
936  padUpper->cd();
937 
938  totMC ->Draw();
939  dataHist->Draw("psame");
940 
941  TLegend leg(0.55, 0.3, 0.9, 0.9);
942  leg.SetBorderSize(0);
943  leg.SetFillStyle(0);
944  leg.AddEntry(dataHist, "Data", "lep");
945  leg.AddEntry(totMC, "MC", "l");
946  leg.Draw();
947 
948  //numuCanv->cd(canvDiv);
949  numuCanv->cd();
950  padLower->Draw();
951  padLower->cd();
952 
953  ratio->SetMaximum(1.15);
954  ratio->SetMinimum(0.85);
955  ratio->Draw("e1");
956 
957  numuCanv->Update();
958  numuCanv->Write();
959  } // end loop over numu selections
960 
961 
962  for(int s = 0; s < 3; ++s){
963 
964  if(s == 2 && d == 1) continue;
965 
967 
968  TCanvas *nueCanv = tfd.make<TCanvas>(("Nue" + nameBase + "FitCanvas").c_str(), "#nu_{e} Canvas", 600, 800);
969  //nueCanv->Divide(4, 3);
970 
971  int keyBase = (bt + 1) * 1000 + d * 100 + nueSels[s];
972 
973  auto dataHist = (*fDataSpectraE)[keyBase + (fnex::kDataFile + 1) * 10];
974 
975  // make the histogram
976  TH1D *totMC = tfd.make<TH1D>(("MC" + nameBase + "Summed").c_str(),
977  ";Energy (Gev);Events",
978  nueBins.size() - 1,
979  &nueBins[0]);
980 
981  TH1D *ratio = tfd.make<TH1D>(("DataMCRatio" + nameBase).c_str(),
982  ";Energy (Gev);Ratio",
983  nueBins.size() - 1,
984  &nueBins[0]);
985 
986  totMC->GetXaxis()->CenterTitle();
987  totMC->GetYaxis()->CenterTitle();
988  totMC->SetLineColor(kRed);
989  totMC->SetMarkerColor(kRed);
990 
991  ratio->GetXaxis()->CenterTitle();
992  ratio->GetYaxis()->CenterTitle();
993  ratio->GetYaxis()->CenterTitle();
994  ratio->GetYaxis()->SetDecimals();
995  ratio->GetXaxis()->SetLabelSize(0.1);
996  ratio->GetYaxis()->SetLabelSize(0.1);
997 
998  // fill it
999  for(auto const& itr : fileTypes){
1000  totMC->Add((*fMCSpectraE)[keyBase + (itr + 1) * 10]);
1001  LOG_DEBUG("CovarianceFitHelper")
1002  << (*fMCSpectraE)[keyBase + (itr + 1) * 10]->GetNbinsX()
1003  << " "
1004  << keyBase + (itr + 1) * 10;
1005  }
1006 
1007  ratio->Divide(dataHist, totMC);
1008 
1009  double max = std::max(dataHist->GetMaximum(), totMC->GetMaximum());
1010  totMC->SetMaximum(1.2 * max);
1011 
1012  canvDiv = (bt * 8) + (d - 1) * 3 + s + 1;
1013  //nueCanv->cd(canvDiv);
1014 
1015  // make a couple of pads in the canvas
1016  TPad *padUpper = tfd.make<TPad>(("Upper" + nameBase).c_str(), ("Upper" + nameBase).c_str(), 0, 0.3, 1, 1.0);
1017  padUpper->SetBottomMargin(0.016);
1018  padUpper->SetBorderSize(0);
1019 
1020  TPad *padLower = tfd.make<TPad>(("Lower" + nameBase).c_str(), ("Lower" + nameBase).c_str(), 0, 0.05, 1, 0.3);
1021  padLower->SetTopMargin(0);
1022  padLower->SetBottomMargin(0.25);
1023  padLower->SetBorderSize(0);
1024 
1025  padUpper->Draw();
1026  padUpper->cd();
1027 
1028  totMC ->Draw();
1029  dataHist->Draw("psame");
1030 
1031  TLegend leg(0.55, 0.3, 0.9, 0.9);
1032  leg.SetBorderSize(0);
1033  leg.SetFillStyle(0);
1034  leg.AddEntry(dataHist, "Data", "lep");
1035  leg.AddEntry(totMC, "MC", "l");
1036  leg.Draw();
1037 
1038  nueCanv->cd();
1039  padLower->Draw();
1040  padLower->cd();
1041 
1042  ratio->SetMaximum(1.15);
1043  ratio->SetMinimum(0.85);
1044  ratio->Draw("e1");
1045 
1046  nueCanv->Update();
1047  nueCanv->Write();
1048 
1049  } // end loop over nue selections
1050  } // end loop over detectors
1051  } // end loop over beams
1052 
1053  return;
1054  }
1055 
1056  //......................................................................
1058  TCanvas *canv,
1059  int div,
1060  fnex::BeamType_t const& bt,
1061  fnex::SelectionType_t const& sel,
1062  novadaq::cnv::DetId const& det)
1063  {
1064  canv->cd(div);
1065 
1066  // grab all of the MC histograms for this selection and detector combination and
1067  // create a summed histogram
1068  std::set <fnex::FileType_t> fileTypes({fnex::kBeam,
1069  fnex::kSwap,
1072 
1073 
1074  return;
1075  }
1076 
1077 } // end namespace
static std::string GetName(int id)
std::map< int, double > fDataPOT
Data POT values for each detector/beam/selection.
T max(const caf::Proxy< T > &a, T b)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
std::unique_ptr< TMatrixD > fInvCovarianceMatrix
inverted covariance matrix
static const unsigned char kNu_RecoE
Definition: VarVals.h:36
std::unique_ptr< fnex::EventListManipulator > fManipulator
event list manipulator
std::unique_ptr< TMatrixD > fColDiffMatrix
column-wise matrix for the difference between data and MC
const Var weight
bool fStatOnlyFit
Statistics only fit.
double Integral(const Spectrum &s, const double pot, int cvnregion)
TH1 * ratio(TH1 *h1, TH1 *h2)
TFileDirectory mkdir(std::string const &dir, std::string const &descr="")
constexpr T pow(T x)
Definition: pow.h:75
std::unique_ptr< std::vector< fnex::CovFitIteration > > fIterations
keep track of the minimizer iterations
TCanvas * canv
void FillSpectrum(std::vector< double > &spectrum, bool data=false)
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.
void reconfigure(fhicl::ParameterSet const &pset)
std::string fCovarianceHistFileName
name of the covariance histogram file
static ShifterAndWeighter * Instance()
std::unique_ptr< std::map< int, TH1D * > > fDataSpectraE
Data spectrum for each detector/file/beam/selection in energy.
const XML_Char const XML_Char * data
Definition: expat.h:268
const XML_Char * s
Definition: expat.h:262
void UpdateInputPoint(fnex::InputPoint const &point)
Far Detector at Ash River, MN.
std::unique_ptr< fnex::InputPoint > fInputPoint
input point used to set the ShifterAndWeighter before each iteration
enum fnex::sel_type SelectionType_t
double Weight(fnex::EventPtr const &curEvent, fnex::MetaData const &md, fnex::WeightType_t const wgtType=kAllWgt)
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
std::unique_ptr< std::vector< double > > fMCSpectrumNominal
Nominal MC spectrum in bins.
const std::string kFileTypeStrings[8]
Definition: Constants.h:386
std::unique_ptr< std::vector< double > > fDataSpectrum
data spectrum in bins
double energy
Definition: plottest35.C:25
std::string fCovarianceHistName
full path to covariance histogram in file
Near Detector in the NuMI cavern.
Float_t d
Definition: plot.C:236
const double j
Definition: BetheBloch.cxx:29
std::unique_ptr< std::set< int > > fKeysToUse
keys to include in the fit
fnex::EventListMap fEventLists
the lists of events we are going to use, data and MC
TMatrixT< double > TMatrixD
Definition: Utilities.h:18
float bin[41]
Definition: plottest35.C:14
std::unique_ptr< TMatrixD > fRowDiffMatrix
row-wise matrix for the difference between data and MC
std::unique_ptr< std::map< int, TH1D * > > fMCSpectraE
MC spectrum for each detector/file/beam/selection in energy.
const std::string cBeamType_Strings[4]
Definition: Constants.h:76
void FillEnergySpectra(std::set< fnex::MetaData > const &mdToFill)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
TH2D * fCovarianceMatrixSyst
Covariance matrix.
void FillSelectionPad(art::TFileDirectory &tfd, TCanvas *canv, int div, fnex::BeamType_t const &bt, fnex::SelectionType_t const &sel, novadaq::cnv::DetId const &det)
Float_t norm
TH2D * fCovarianceMatrixStat
Covariance matrix.
const hit & b
Definition: hits.cxx:21
enum fnex::beam_type BeamType_t
std::map< std::string, NPInfo > const & NuisanceParameterInfo() const
const std::string cSelectionType_Strings[11]
Definition: Constants.h:101
Float_t e
Definition: plot.C:35
std::unique_ptr< std::vector< double > > fMCSpectrum
MC spectrum in bins.
TH2D * fCovarianceMatrixSystScaled
Covariance matrix.
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
art::ServiceHandle< art::TFileService > fTFS
TFileService.
TH2D * fCovarianceMatrixHist
Covariance matrix.
void events(int which)
Definition: Cana.C:52
std::unique_ptr< fnex::CovarianceBinUtility > fBinUtil
useful functions for mapping bins to energy and back
static NuisanceParameters * Instance()