GridPointResult.cxx
Go to the documentation of this file.
1 //
2 // Created by Brian Rebel on 12/6/19.
3 //
4 
6 #include "cetlib_except/exception.h"
7 
11 
12 namespace cmf{
13 
14  //----------------------------------------------------------------------------
16  {
17  if(this->X() < other.X()) return true;
18  else if(this->X() == other.X()){
19  if(this->Y() < other.Y()) return true;
20  }
21 
22  return false;
23  }
24 
25  //----------------------------------------------------------------------------
27  {
28  if(this->X() == other.X() && this->Y() == other.Y()) return true;
29 
30  return false;
31  }
32 
33  //----------------------------------------------------------------------------
34  std::ostream & operator<<(std::ostream & o, cmf::GridPoint const& gp)
35  {
36 
37  if(cmf::IsAngleParameter(gp.XPar()))
38  o << "sin^{2}("
40  << ") : ";
41  else
42  o << cmf::cOscParams_Strings[gp.XPar()]
43  << " : ";
44 
45  o << gp.X()
46  << " ";
47 
48  if(cmf::IsAngleParameter(gp.YPar()))
49  o << "sin^{2}("
50  << cmf::cOscParams_Strings[gp.YPar()]
51  << ") : ";
52  else
53  o << cmf::cOscParams_Strings[gp.YPar()]
54  << " : ";
55 
56  o << gp.Y();
57 
58  return o;
59  }
60 
61  //----------------------------------------------------------------------------
62  OscParamPoint::OscParamPoint(std::vector<cmf::OscParm_t> const& oscParms,
63  std::vector<float> const& point,
64  cmf::Spectrum const& mcSpectrum)
65  : fSpectrum(mcSpectrum)
66  {
67  // We want to order the oscillation parameters for future use according to
68  // which analysis the prediction is for:
69  // For three flavor analyses, order by Th23, Dmsq32, Th13, dCP
70  // For four flavor analyses, order by Th24, Dmsq41, Th34, Th23, Dmsq32
71 
72  std::map<cmf::OscParm_t, size_t> parToIdx;
73  for(size_t p = 0; p < oscParms.size(); ++p) parToIdx.emplace(oscParms[p], p);
74 
75  std::vector<cmf::OscParm_t> parOrder;
76 
77  if(parToIdx.find(cmf::kTh24) != parToIdx.end() ){
78  parOrder = std::vector<cmf::OscParm_t>({cmf::kTh24, cmf::kDmsq41, cmf::kTh34, cmf::kTh23, cmf::kDmsq32, cmf::kd24});
79  }
80  // look for an NSI parameter - this one needs to be fixed up
81  else if(parToIdx.find(cmf::kEps_mumu) != parToIdx.end()){
82  parOrder = std::vector<cmf::OscParm_t>({cmf::kEps_mumu, cmf::kEps_ee, cmf::kTh23, cmf::kDmsq32});
83  }
84  else{
85  parOrder = std::vector<cmf::OscParm_t>({cmf::kTh23, cmf::kDmsq32, cmf::kTh13, cmf::kdCP});
86  }
87 
88  for(auto const& itr : parOrder){
89  if(parToIdx.count(itr) > 0){
90  size_t idx = parToIdx.find(itr)->second;
91  fOscParVec.emplace_back(itr, point[idx]);
92  MF_LOG_DEBUG("GridPointResult")
93  << "adding "
94  << fOscParVec.back();
95  }
96  else
97  MF_LOG_WARNING("GridPointResult")
98  << "expected parameter "
100  << " not in parToIdx map";
101  }
102 
103  }
104 
105  //----------------------------------------------------------------------------
107  cmf::Spectrum const& mcSpectrum)
108  : fSpectrum(mcSpectrum)
109  {
110  // We want to order the oscillation parameters for future use according to
111  // which analysis the prediction is for:
112  // For three flavor analyses, order by Th23, Dmsq32, Th13, dCP
113  // For four flavor analyses, order by Th24, Dmsq41, Th34, Th23, Dmsq32
114 
115  std::vector<cmf::OscParm_t> parOrder;
116  if(point.find(cmf::kTh24) != point.end() ){
117  parOrder = std::vector<cmf::OscParm_t>({cmf::kTh24, cmf::kDmsq41, cmf::kTh34, cmf::kTh23, cmf::kDmsq32, cmf::kd24});
118  }
119  // look for an NSI parameter - this one needs to be fixed up
120  else if(point.find(cmf::kEps_mumu) != point.end()){
121  parOrder = std::vector<cmf::OscParm_t>({cmf::kEps_mumu, cmf::kEps_ee, cmf::kTh23, cmf::kDmsq32});
122  }
123  else{
124  parOrder = std::vector<cmf::OscParm_t>({cmf::kTh23, cmf::kDmsq32, cmf::kTh13, cmf::kdCP});
125  }
126 
127  for(auto const& itr : parOrder){
128  if(point.count(itr) > 0){
129  fOscParVec.emplace_back(point.find(itr)->first, point.find(itr)->second);
130  MF_LOG_DEBUG("GridPointResult")
131  << "adding "
132  << fOscParVec.back();
133  }
134  else
135  MF_LOG_WARNING("GridPointResult")
136  << "expected parameter "
138  << " not in parToIdx map";
139  }
140 
141  }
142 
143  //----------------------------------------------------------------------------
145  {
147  for(auto const& itr : fOscParVec)
148  pointMap.emplace(itr.Param(), itr.Value());
149 
150  return std::move(pointMap);
151  }
152 
153  //----------------------------------------------------------------------------
154  cmf::OscParamPoint OscParamPoint::ScaleSpectra(std::vector<cmf::DetBeamBins> const& dbBinsVec,
155  std::map<cmf::DetBeamPair, float> const& dbScaleMap) const
156  {
157  cmf::Spectrum spect = fSpectrum;
158  for(auto const& dbBins : dbBinsVec)
159  {
160  cmf::DetBeamPair binDetBeam;
161  binDetBeam.first = dbBins.fDet;
162  binDetBeam.second = dbBins.fBeam;
163 
164  auto findScale = dbScaleMap.find(binDetBeam);
165  if (findScale != dbScaleMap.end())
166  {
167  MF_LOG_VERBATIM("OscParamPoint")
168  << "attempting to scale spectra for "
169  << cmf::cDetType_Strings[binDetBeam.first]
170  << " "
171  << cmf::cBeamType_Strings[binDetBeam.second]
172  << " by "
173  << findScale->second;
174 
175  if (dbBins.fHighBin > fSpectrum.size() - 1)
176  throw cet::exception("OscParamPoint")
177  << "attempting to scale spectrum out of range. There is no bin "
178  << dbBins.fHighBin
179  << " in spectrum of size "
180  << fSpectrum.size();
181  for (size_t bin = dbBins.fLowBin; bin < dbBins.fHighBin + 1; ++bin)
182  {
183  spect[bin] *= findScale->second;
184  }
185  }
186  }
187  return OscParamPoint(this->OscPoint(), spect);
188  }
189 
190  //----------------------------------------------------------------------------
192  {
193  // look for a four flavor parameter to see if we want that ordering
194  std::vector<float> parVals;
195  std::vector<float> otherVals;
196 
197  for(auto const& itr : fOscParVec) parVals .emplace_back(itr.Value());
198  for(auto const& itr : other.OscPointAsVec()) otherVals.emplace_back(itr.Value());
199 
200  if(parVals[0] < otherVals[0]) return true;
201  else if(parVals[0] == otherVals[0]){
202  if(parVals[1] < otherVals[1]) return true;
203  else if(parVals[1] == otherVals[1]){
204  if(parVals[2] < otherVals[2]) return true;
205  else if(parVals[2] == otherVals[2]){
206  if(parVals[3] < otherVals[3]) return true;
207  if(parVals.size() == 5){
208  if(parVals[3] == otherVals[3]){
209  if(parVals[4] < otherVals[4])
210  return true;
211  }
212  } // end if we are in the four flavor analysis
213  }
214  }
215  }
216 
217  return false;
218  }
219 
220  //----------------------------------------------------------------------------
222  {
223  // if we can't find one of the parameters in other or any of the parameters
224  // have different values these OscParamPoints are different
225  for(auto const& itr : fOscParVec){
226  auto const otherItr = std::find(other.OscPointAsVec().begin(), other.OscPointAsVec().end(), itr);
227  if(otherItr == other.OscPointAsVec().end()) return false;
228  if(otherItr->Value() != itr.Value()) return false;
229  }
230 
231  return true;
232  }
233 
234  //----------------------------------------------------------------------------
236  {
237  if(fSpectrum.size() != other.PredictedSpectrum().size()){
238  MF_LOG_WARNING("OscParamPoint")
239  << "We have spectra of different sizes: "
240  << fSpectrum.size()
241  << " vs other "
242  << other.PredictedSpectrum().size()
243  << " so bail";
244  return;
245  }
246 
247  // we got here, so we can add the spectra
248  for(size_t b = 0; b < fSpectrum.size(); ++b){
249  fSpectrum[b] += other.PredictedSpectrum()[b];
250  }
251  }
252 
253  //----------------------------------------------------------------------------
254  std::ostream& operator<< (std::ostream & o, cmf::OscParamPoint const& hpp)
255  {
256  for(auto const& itr : hpp.OscPointAsVec())
257  o << cmf::cOscParams_Strings[itr.Param()]
258  << " "
259  << itr.Value()
260  << " ";
261 
262  o << std::endl;
263 
264  return o;
265  }
266 
267  //----------------------------------------------------------------------------
268  void UniversePointBestFits::AddPointsToUniverse(std::vector<cmf::PointBestFit> const& randomColl,
269  std::vector<cmf::PointBestFit> const& asimovColl)
270  {
271 
272  for(auto const& itr : randomColl) fRandomColl.emplace_back(itr);
273  for(auto const& itr : asimovColl) fAsimovColl.emplace_back(itr);
274 
275  std::sort(fRandomColl.begin(), fRandomColl.end());
276  std::sort(fAsimovColl.begin(), fAsimovColl.end());
277  }
278 
279  //----------------------------------------------------------------------------
281  {
282  std::sort(fRandomColl.begin(), fRandomColl.end());
283  std::sort(fAsimovColl.begin(), fAsimovColl.end());
284  }
285 }
cmf::OscillationParameterMap OscPoint() const
cmf::Spectrum fSpectrum
spectrum for this combination of parameters
const char * p
Definition: xmltok.h:285
const std::vector< std::string > cDetType_Strings({"UnknownDet","NearDet","FarDet","MINOSNear","MINOSFar","AllDetectors"})
std::vector< cmf::OscPar > fOscParVec
the point in oscillation space
bool operator==(cmf::OscParamPoint const &other) const
std::vector< double > Spectrum
Definition: Constants.h:743
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
friend std::ostream & operator<<(std::ostream &o, cmf::OscParamPoint const &hpp)
bool operator<(cmf::OscParamPoint const &other) const
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
void operator+=(cmf::OscParamPoint const &other)
std::map< cmf::OscParm_t, float > OscillationParameterMap
Definition: Constants.h:745
cmf::Spectrum const & PredictedSpectrum() const
bool operator==(cmf::GridPoint const &other) const
const std::vector< std::string > cBeamType_Strings({"FHC","RHC","0HC","UnknownBeam"})
float bin[41]
Definition: plottest35.C:14
float const & X() const
Module to combine a set of results into a single file currently only does one data product type at a ...
Definition: Event.cxx:24
bool operator<(cmf::GridPoint const &other) const
std::vector< cmf::OscPar > const & OscPointAsVec() const
cmf::OscParm_t const & XPar() const
void AddPointsToUniverse(std::vector< cmf::PointBestFit > const &randomColl, std::vector< cmf::PointBestFit > const &asimovColl)
#define MF_LOG_VERBATIM(category)
#define MF_LOG_DEBUG(id)
float const & Y() const
OscParamPoint()=default
const hit & b
Definition: hits.cxx:21
friend std::ostream & operator<<(std::ostream &o, cmf::GridPoint const &gp)
#define MF_LOG_WARNING(category)
cmf::OscParm_t const & YPar() const
std::pair< cmf::DetType_t, cmf::BeamType_t > DetBeamPair
Definition: Constants.h:738
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"})