CMFSpectraFromPredictions_plugin.cc
Go to the documentation of this file.
1 //
2 // Plugin to create histograms of the library prediction spectra.
3 // Pass the art job a list of input library prediction files for the space
4 //
5 // Created by Brian Rebel on 3/16/21.
6 //
7 
8 #include <string>
9 #include <set>
10 #include <map>
11 #include <memory>
12 #include <sstream>
13 
14 //#include "TStopwatch.h"
15 #include "TMath.h"
16 #include "TGraph2D.h"
17 
18 // Framework includes
24 #include "art_root_io/TFileService.h"
25 #include "art_root_io/TFileDirectory.h"
26 #include "fhiclcpp/ParameterSet.h"
28 #include "cetlib_except/exception.h"
29 
30 // NOvA includes
39 
40 namespace cmf{
41 
43  public:
44  explicit SpectraFromPredictions(fhicl::ParameterSet const& pset);
45 
46  ~SpectraFromPredictions() override = default;
47 
48  void readResults(art::Results const& r) override;
49 
50  void writeResults(art::Results& r) override;
51 
52  void reconfigure(fhicl::ParameterSet const& p);
53 
54  void clear() override;
55  void endJob() override;
56 
57  private:
58 
59  float fNDFHCExposureScale; ///< how much the ND FHC bins in the prediction libraries should be scaled by
60  float fNDRHCExposureScale; ///< how much the ND RHC bins in the prediction libraries should be scaled by
61  float fFDFHCExposureScale; ///< how much the FD FHC bins in the prediction libraries should be scaled by
62  float fFDRHCExposureScale; ///< how much the FD RHC bins in the prediction libraries should be scaled by
63  std::string fPredictionModule; ///< label for module containing the prediction point library
64  cmf::OscParm_t fXParamEnum; ///< enumerated value of the parameter from CovarianceMatrixFit/dataProducts/Constants.h
65  cmf::OscParm_t fYParamEnum; ///< enumerated value of the parameter from CovarianceMatrixFit/dataProducts/Constants.h
66  fhicl::ParameterSet fPlotConfig; ///< configuration for the plot utilities
67  bool fLogXAxis; ///< are parameters on a log scale on the X axis
68  bool fLogYAxis; ///< are parameters on a log scale on the Y axis
69  bool fScaleExposure; ///< should the exposure of the prediction libraries be scaled?
70  std::vector<cmf::PointSpectra> fPointSpectraVec; ///< the point spectra for all the oscillation points
71 
72  std::map<cmf::OscParm_t, std::vector<float> > fOscParValuesToKeep; ///< which values of osc parameters do we want to keep?
73  std::vector<cmf::DetBeamBins> fDBBinsVec; ///< which bins are for what detector/beam type combination
74  std::map<std::pair<cmf::DetType_t, cmf::BeamType_t>, float> fDBScaleMap; ///< how do we want to scale our exposure?
75 
76  };
77 
78  //......................................................................
79  // the X and Y parameter enumerations are only there for
80  // ordering the predictions so it does not matter if they are 3F or 4F
81  // pick Th23 and Dmsq32 as those are always present for 3F, 4F and NSI
85  , fPlotConfig(pset.get<fhicl::ParameterSet>("PlotUtilities"))
86  {
87  this->reconfigure(pset);
88  }
89 
90  //......................................................................
92  {
95  cmf::ParameterUtility::Instance() ->Initialize(pset.get< fhicl::ParameterSet >("ParametersToUse" ));
97 
98  fNDFHCExposureScale = pset.get< float >("NDFHCExposureScale", 1.0 );
99  fFDFHCExposureScale = pset.get< float >("FDFHCExposureScale", 1.0 );
100  fNDRHCExposureScale = pset.get< float >("NDRHCExposureScale", 1.0 );
101  fFDRHCExposureScale = pset.get< float >("FDRHCExposureScale", 1.0 );
102  fScaleExposure = pset.get< bool >("ScaleExposure", false);
103  fPredictionModule = pset.get<std::string >("PredictionLibraryModule" );
104 
106 
108  fDBScaleMap.emplace(std::make_pair(cmf::kFARDET, cmf::kFHC), fFDFHCExposureScale);
109  fDBScaleMap.emplace(std::make_pair(cmf::kNEARDET, cmf::kRHC), fNDRHCExposureScale);
110  fDBScaleMap.emplace(std::make_pair(cmf::kFARDET, cmf::kRHC), fFDRHCExposureScale);
111 
112  // figure out the points in space we want to use
113  auto pointConfig = pset.get<fhicl::ParameterSet>("SpectraPoints", fhicl::ParameterSet());
114 
115  for(auto const& itr : pointConfig.get_names()){
116 
117  auto parConfig = pointConfig.get<fhicl::ParameterSet>(itr);
118 
119  MF_LOG_VERBATIM("SpectraFromPredictions")
120  << parConfig.to_indented_string();
121 
122  fOscParValuesToKeep.emplace(cmf::StringToOscParmType(parConfig.get<std::string>("OscParName")),
123  parConfig.get<std::vector<float>>("ParVals"));
124  }
125 
126  MF_LOG_VERBATIM("SFP")
127  << "Write spectra for the following points to file:";
128  for(auto const& itr : fOscParValuesToKeep){
130  for(auto const& vItr : itr.second){
131  values += " " + std::to_string(vItr);
132  }
133  MF_LOG_VERBATIM("SFP")
134  << cmf::cOscParams_Strings[itr.first]
135  << values;
136  }
137  }
138 
139  //......................................................................
140  // Method to clear out the collections of data products after the
141  // writeResults method is called.
143  {}
144 
145  //......................................................................
146  // We need to read in several types of files - the best fit,
147  // prediction library and random universe files. The latter is what we
148  // want if we are making median sensitivities and also has the Asimov
149  // predictions for the fake data.
151  {
152 
154  if( r.getByLabel(fPredictionModule, predLib) &&
155  predLib.isValid() ){
156  MF_LOG_VERBATIM("CovarianceMatrixFitter")
157  << "Valid handle to std::vector<cmf::OscParamPoint> "
158  << "objects found in "
160  << " "
161  << " vector size "
162  << predLib->size();
163 
164  auto useTrig = fPlotConfig.get<bool>("UseTrig");
165  auto psvItr = fPointSpectraVec.begin();
166  float xVal;
167  float yVal;
168 
169  // maybe we should first figure out how many grid points we need in the space
170  // and how many entries we need for the hidden parameter combinations at each
171  // point and reserve that space
172 
173  for(auto const& itrUnscaled : *predLib){
174 
175  cmf::OscParamPoint itr = (fScaleExposure) ? itrUnscaled.ScaleSpectra(fDBBinsVec, fDBScaleMap) : itrUnscaled;
176 
177  MF_LOG_DEBUG("CMFSpectraFromPredictions")
178  << itr;
179 
180  auto tmpMap = itr.OscPoint();
181 
182  // set collect to decide if we want this point or not.
183  bool keepPoint = true;
184 
185  // keep all points where the desired parameters have the desired values
186  for(auto const& keepParItr : fOscParValuesToKeep){
187  auto tmpItr = tmpMap.find(keepParItr.first);
188  if(tmpItr != tmpMap.end()){
189  auto const& parVals = keepParItr.second;
190 
191  auto opvItr = std::find(parVals.begin(), parVals.end(), tmpItr->second);
192  if(opvItr == parVals.end()){
193  keepPoint = false;
194  }
195  } // end if the parameter we care about is in the space for these predictions
196  } // end loop over parameters to check
197 
198  // if keepPoint is false, we don't want this prediction
199  if(!keepPoint) continue;
200 
201  xVal = (useTrig && cmf::IsAngleParameter(fXParamEnum)) ? std::pow(std::sin(itr.OscPoint().find(fXParamEnum)->second), 2.) : itr.OscPoint().find(fXParamEnum)->second;
202  yVal = (useTrig && cmf::IsAngleParameter(fYParamEnum)) ? std::pow(std::sin(itr.OscPoint().find(fYParamEnum)->second), 2.) : itr.OscPoint().find(fYParamEnum)->second;
203 
204  auto gridPoint = cmf::GridPoint(xVal,
205  yVal,
206  fXParamEnum,
207  fYParamEnum);
208  psvItr = std::find(fPointSpectraVec.begin(),
209  fPointSpectraVec.end(),
210  PointSpectra(gridPoint));
211  if(psvItr == fPointSpectraVec.end()){
212  fPointSpectraVec.emplace_back(gridPoint);
213  fPointSpectraVec.back().oscPoints.emplace_back(itr);
214  }
215  else
216  psvItr->oscPoints.emplace_back(itr);
217  }
218 
219  // release space needed to read in the prediction library
220  r.removeCachedProduct(predLib);
221  }
222 
223  // for(auto const& itr : fPointSpectraVec)
224  // MF_LOG_VERBATIM("SpectraFromPredictions")
225  // << itr.point
226  // << " has "
227  // << itr.oscPoints.size()
228  // << " spectra associated with it";
229 
230  }
231 
232  //----------------------------------------------------------------------------
234  {
235  // make sure we aren't taking up too much space
236  fPointSpectraVec.shrink_to_fit();
237 
238  // do the same for the osc points in the vector
239  for(auto & itr : fPointSpectraVec) itr.oscPoints.shrink_to_fit();
240 
241  MF_LOG_VERBATIM("SpectraFromPredictions")
242  << " There are "
243  << fPointSpectraVec.size()
244  << " library points in the space with "
245  << fPointSpectraVec.front().oscPoints.size()
246  << " predictions at each point ";
247 
248  // now sort the universe point results in the library points vector
249  std::sort(fPointSpectraVec.begin(), fPointSpectraVec.end());
250 
251  for(auto & itr : fPointSpectraVec){
252  // sort the vector of OscParamPoints
253  std::sort(itr.oscPoints.begin(), itr.oscPoints.end());
254 
255  MF_LOG_DEBUG("SpectraFromPredictions")
256  << itr;
257  }
258 
259  // for(auto const& oitr : fPointSpectraVec.front().oscPoints)
260  // MF_LOG_VERBATIM("SpectraFromPredictions")
261  // << "check osc point ordering: "
262  // << oitr;
263 
265 
266  // make a TFileDirectory to store the histograms for this calculation
267  art::TFileDirectory curDir = tfs->mkdir("predictedSpectra");
268 
269  std::string baseName;
270  for(auto & itr : fPointSpectraVec){
271  for(auto & spItr : itr.oscPoints){
272  baseName.clear();
273  for(auto & oscParItr : spItr.OscPointAsVec()){
274  baseName += cmf::cOscParams_Strings[oscParItr.Param()] + "_" + std::to_string(oscParItr.Value()) + "_";
275  }
277  baseName + "Spectrum",
278  ";Bin;Events",
279  spItr.PredictedSpectrum());
280  }
281  }
282 
283  }
284 
285  //----------------------------------------------------------------------------
287  {
288  }
289 
291 
292 } // end namespace
cmf::OscillationParameterMap OscPoint() const
enum cmf::osc_params OscParm_t
void Initialize(fhicl::ParameterSet const &pset)
std::string fPredictionModule
label for module containing the prediction point library
const char * p
Definition: xmltok.h:285
constexpr T pow(T x)
Definition: pow.h:72
fhicl::ParameterSet fPlotConfig
configuration for the plot utilities
static SelectionUtility * Instance()
float fFDFHCExposureScale
how much the FD FHC bins in the prediction libraries should be scaled by
std::vector< cmf::DetBeamBins > const & DetectorBeamBins()
TH1F * MakeSpectrumHistogram(art::TFileDirectory &tfd, std::string const &baseName, std::string const &baseTitle, cmf::Spectrum const &spectrum)
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
std::map< std::pair< cmf::DetType_t, cmf::BeamType_t >, float > fDBScaleMap
how do we want to scale our exposure?
static bool IsAngleParameter(cmf::OscParm_t const &par)
Definition: StaticFuncs.h:354
cmf::OscParamPoint ScaleSpectra(std::vector< cmf::DetBeamBins > const &dbBinsVec, std::map< cmf::DetBeamPair, float > const &dBScaleMap) const
bool isValid() const
Definition: Handle.h:183
float fFDRHCExposureScale
how much the FD RHC bins in the prediction libraries should be scaled by
static ParameterUtility * Instance()
void readResults(art::Results const &r) override
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
bool fLogXAxis
are parameters on a log scale on the X axis
std::map< cmf::OscParm_t, std::vector< float > > fOscParValuesToKeep
which values of osc parameters do we want to keep?
#define DEFINE_ART_RESULTS_PLUGIN(klass)
T get(std::string const &key) const
Definition: ParameterSet.h:231
void Initialize(fhicl::ParameterSet const &pset)
bool fScaleExposure
should the exposure of the prediction libraries be scaled?
void Initialize(fhicl::ParameterSet const &plottingPars)
bool removeCachedProduct(Handle< PROD > &) const
Definition: DataViewImpl.h:971
float fNDRHCExposureScale
how much the ND RHC bins in the prediction libraries should be scaled by
Module to combine a set of results into a single file currently only does one data product type at a ...
Definition: Event.cxx:24
float fNDFHCExposureScale
how much the ND FHC bins in the prediction libraries should be scaled by
T sin(T number)
Definition: d0nt_math.hpp:132
#define MF_LOG_VERBATIM(category)
#define MF_LOG_DEBUG(id)
void reconfigure(fhicl::ParameterSet const &p)
static PlotUtilities * Instance()
TRandom3 r(0)
cmf::OscParm_t fYParamEnum
enumerated value of the parameter from CovarianceMatrixFit/dataProducts/Constants.h
cmf::OscParm_t fXParamEnum
enumerated value of the parameter from CovarianceMatrixFit/dataProducts/Constants.h
~SpectraFromPredictions() override=default
void Initialize(fhicl::ParameterSet const &pset)
void writeResults(art::Results &r) override
bool fLogYAxis
are parameters on a log scale on the Y axis
std::string to_string(ModuleType const mt)
Definition: ModuleType.h:34
std::vector< cmf::DetBeamBins > fDBBinsVec
which bins are for what detector/beam type combination
std::vector< cmf::PointSpectra > fPointSpectraVec
the point spectra for all the oscillation points
static cmf::OscParm_t StringToOscParmType(std::string const &str)
SpectraFromPredictions(fhicl::ParameterSet const &pset)
static CovarianceBinUtility * Instance()
const std::vector< std::string > cOscParams_Strings({"L","Rho","Dmsq21","Dmsq32","Th12","Th13","Th23","dCP","Th14","Th24","Th34","Dmsq41","d24","Eps_ee","Eps_emu","Eps_etau","Eps_mumu","Eps_mutau","Eps_tautau","Delta_emu","Delta_etau","Delta_mutau"})
enum BeamMode string