File indexing completed on 2024-04-06 12:23:38
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017 #include <vector>
0018
0019 #include "CondFormats/PhysicsToolsObjects/interface/Histogram.h"
0020 #include "PhysicsTools/MVAComputer/interface/VarProcessor.h"
0021 #include "PhysicsTools/MVAComputer/interface/Calibration.h"
0022 #include "PhysicsTools/MVAComputer/interface/Spline.h"
0023
0024 using namespace PhysicsTools;
0025
0026 namespace {
0027
0028 class ProcNormalize : public VarProcessor {
0029 public:
0030 typedef VarProcessor::Registry::Registry<ProcNormalize, Calibration::ProcNormalize> Registry;
0031
0032 ProcNormalize(const char *name, const Calibration::ProcNormalize *calib, const MVAComputer *computer);
0033 ~ProcNormalize() override {}
0034
0035 void configure(ConfIterator iter, unsigned int n) override;
0036 void eval(ValueIterator iter, unsigned int n) const override;
0037 std::vector<double> deriv(ValueIterator iter, unsigned int n) const override;
0038
0039 private:
0040 struct Map {
0041 Map(const Calibration::HistogramF &pdf) : min(pdf.range().min), width(pdf.range().width()) {
0042 std::vector<double> values(pdf.values().begin() + 1, pdf.values().end() - 1);
0043 spline.set(values.size(), &values.front());
0044 }
0045
0046 double min, width;
0047 Spline spline;
0048 };
0049
0050 void findMap(ValueIterator iter,
0051 unsigned int n,
0052 std::vector<Map>::const_iterator &begin,
0053 std::vector<Map>::const_iterator &end) const;
0054
0055 std::vector<Map> maps;
0056 int categoryIdx;
0057 unsigned int nCategories;
0058 };
0059
0060 ProcNormalize::Registry registry("ProcNormalize");
0061
0062 ProcNormalize::ProcNormalize(const char *name, const Calibration::ProcNormalize *calib, const MVAComputer *computer)
0063 : VarProcessor(name, calib, computer),
0064 maps(calib->distr.begin(), calib->distr.end()),
0065 categoryIdx(calib->categoryIdx),
0066 nCategories(1) {}
0067
0068 void ProcNormalize::configure(ConfIterator iter, unsigned int n) {
0069 if (categoryIdx >= 0) {
0070 if ((int)n < categoryIdx + 1)
0071 return;
0072 nCategories = maps.size() / (n - 1);
0073 if (nCategories * (n - 1) != maps.size())
0074 return;
0075 } else if (n != maps.size())
0076 return;
0077
0078 int i = 0;
0079 while (iter) {
0080 if (categoryIdx == i++)
0081 iter++(Variable::FLAG_NONE);
0082 else
0083 iter << iter++(Variable::FLAG_ALL);
0084 }
0085 }
0086
0087 void ProcNormalize::findMap(ValueIterator iter,
0088 unsigned int n,
0089 std::vector<Map>::const_iterator &begin,
0090 std::vector<Map>::const_iterator &end) const {
0091 if (categoryIdx >= 0) {
0092 ValueIterator iter2 = iter;
0093 for (int i = 0; i < categoryIdx; i++)
0094 ++iter2;
0095
0096 int cat = (int)*iter2;
0097 if (cat < 0 || (unsigned int)cat >= nCategories) {
0098 for (; iter; ++iter)
0099 iter();
0100 return;
0101 }
0102
0103 begin = maps.begin() + cat * (n - 1);
0104 end = begin + (n - 1);
0105 } else {
0106 begin = maps.begin();
0107 end = maps.end();
0108 }
0109 }
0110
0111 void ProcNormalize::eval(ValueIterator iter, unsigned int n) const {
0112 std::vector<Map>::const_iterator map, last;
0113 findMap(iter, n, map, last);
0114
0115 for (int i = 0; map != last; ++iter, i++) {
0116 if (i == categoryIdx)
0117 continue;
0118 for (double *value = iter.begin(); value < iter.end(); value++) {
0119 double val = *value;
0120 val = (val - map->min) / map->width;
0121 val = map->spline.integral(val);
0122 iter << val;
0123 }
0124 iter();
0125 ++map;
0126 }
0127 }
0128
0129 std::vector<double> ProcNormalize::deriv(ValueIterator iter, unsigned int n) const {
0130 std::vector<Map>::const_iterator map, last;
0131 findMap(iter, n, map, last);
0132
0133 unsigned int size = 0;
0134 for (ValueIterator iter2 = iter; iter2; ++iter2)
0135 size += iter2.size();
0136
0137 std::vector<double> result(size * size);
0138
0139 unsigned int j = 0;
0140 for (int i = 0; map != last; ++iter, i++) {
0141 if (i == categoryIdx) {
0142 j += iter.size();
0143 continue;
0144 }
0145
0146 for (double *value = iter.begin(); value < iter.end(); value++, j++) {
0147 double val = *value;
0148 val = (val - map->min) / map->width;
0149 val = map->spline.eval(val) * (map->spline.numberOfEntries() - 1) / (map->width * map->spline.getArea());
0150 result[j * size + j] = val;
0151 }
0152 ++map;
0153 }
0154
0155 return result;
0156 }
0157
0158 }