File indexing completed on 2024-04-06 11:59:28
0001 #ifndef MinL3AlgoUniv_H
0002 #define MinL3AlgoUniv_H
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include <vector>
0015 #include <iostream>
0016 #include <map>
0017 #include <cmath>
0018
0019 template <class IDdet>
0020 class MinL3AlgoUniv {
0021 public:
0022 typedef std::map<IDdet, float> IDmap;
0023 typedef typename IDmap::value_type IDmapvalue;
0024 typedef typename IDmap::iterator iter_IDmap;
0025 typedef typename IDmap::const_iterator citer_IDmap;
0026
0027
0028
0029 MinL3AlgoUniv(float kweight_ = 0.);
0030
0031
0032 ~MinL3AlgoUniv();
0033
0034
0035
0036
0037
0038
0039 IDmap iterate(const std::vector<std::vector<float> >& eventMatrix,
0040 const std::vector<std::vector<IDdet> >& idMatrix,
0041 const std::vector<float>& energyVector,
0042 const int& nIter,
0043 const bool& normalizeFlag = false);
0044
0045
0046 void addEvent(const std::vector<float>& myCluster, const std::vector<IDdet>& idCluster, const float& energy);
0047
0048
0049 std::vector<float> recalibrateEvent(const std::vector<float>& myCluster,
0050 const std::vector<IDdet>& idCluster,
0051 const IDmap& newCalibration);
0052
0053
0054
0055 IDmap getSolution(const bool resetsolution = true);
0056
0057
0058 void resetSolution();
0059
0060 private:
0061 float kweight;
0062 int countEvents;
0063 IDmap wsum;
0064 IDmap Ewsum;
0065 };
0066
0067 template <class IDdet>
0068 MinL3AlgoUniv<IDdet>::MinL3AlgoUniv(float kweight_) : kweight(kweight_), countEvents(0) {
0069 resetSolution();
0070 }
0071
0072 template <class IDdet>
0073 MinL3AlgoUniv<IDdet>::~MinL3AlgoUniv() {}
0074
0075 template <class IDdet>
0076 typename MinL3AlgoUniv<IDdet>::IDmap MinL3AlgoUniv<IDdet>::iterate(const std::vector<std::vector<float> >& eventMatrix,
0077 const std::vector<std::vector<IDdet> >& idMatrix,
0078 const std::vector<float>& energyVector,
0079 const int& nIter,
0080 const bool& normalizeFlag) {
0081 int Nevents = eventMatrix.size();
0082
0083 IDmap totalSolution;
0084 IDmap iterSolution;
0085 std::vector<std::vector<float> > myEventMatrix(eventMatrix);
0086 std::vector<float> myEnergyVector(energyVector);
0087
0088 int i;
0089
0090
0091 for (int iter = 1; iter <= nIter; iter++) {
0092
0093 float sumOverEnergy;
0094 if (normalizeFlag) {
0095 float scale = 0.;
0096
0097 for (i = 0; i < Nevents; i++) {
0098 sumOverEnergy = 0.;
0099 for (unsigned j = 0; j < myEventMatrix[i].size(); j++) {
0100 sumOverEnergy += myEventMatrix[i][j];
0101 }
0102 sumOverEnergy /= myEnergyVector[i];
0103 scale += sumOverEnergy;
0104 }
0105 scale /= Nevents;
0106
0107 for (i = 0; i < Nevents; i++) {
0108 myEnergyVector[i] *= scale;
0109 }
0110 }
0111
0112
0113 for (int iEvt = 0; iEvt < Nevents; iEvt++) {
0114 addEvent(myEventMatrix[iEvt], idMatrix[iEvt], myEnergyVector[iEvt]);
0115 }
0116 iterSolution = getSolution();
0117 if (iterSolution.empty())
0118 return iterSolution;
0119
0120
0121 for (int ievent = 0; ievent < Nevents; ievent++) {
0122 myEventMatrix[ievent] = recalibrateEvent(myEventMatrix[ievent], idMatrix[ievent], iterSolution);
0123 }
0124
0125
0126 for (iter_IDmap i = iterSolution.begin(); i != iterSolution.end(); i++) {
0127 iter_IDmap itotal = totalSolution.find(i->first);
0128 if (itotal == totalSolution.end()) {
0129 totalSolution.insert(IDmapvalue(i->first, i->second));
0130 } else {
0131 itotal->second *= i->second;
0132 }
0133 }
0134
0135
0136 }
0137
0138 return totalSolution;
0139 }
0140
0141 template <class IDdet>
0142 void MinL3AlgoUniv<IDdet>::addEvent(const std::vector<float>& myCluster,
0143 const std::vector<IDdet>& idCluster,
0144 const float& energy) {
0145 countEvents++;
0146
0147 float w, invsumXmatrix;
0148 float eventw;
0149
0150 float sumXmatrix = 0.;
0151
0152 for (unsigned i = 0; i < myCluster.size(); i++) {
0153 sumXmatrix += myCluster[i];
0154 }
0155
0156
0157 eventw = 1 - fabs(1 - sumXmatrix / energy);
0158 eventw = pow(eventw, kweight);
0159
0160 if (sumXmatrix != 0.) {
0161 invsumXmatrix = 1 / sumXmatrix;
0162
0163 for (unsigned i = 0; i < myCluster.size(); i++) {
0164 w = myCluster[i] * invsumXmatrix;
0165
0166
0167 iter_IDmap iwsum = wsum.find(idCluster[i]);
0168 if (iwsum == wsum.end())
0169 wsum.insert(IDmapvalue(idCluster[i], w * eventw));
0170 else
0171 iwsum->second += w * eventw;
0172
0173 iter_IDmap iEwsum = Ewsum.find(idCluster[i]);
0174 if (iEwsum == Ewsum.end())
0175 Ewsum.insert(IDmapvalue(idCluster[i], (w * eventw * energy * invsumXmatrix)));
0176 else
0177 iEwsum->second += (w * eventw * energy * invsumXmatrix);
0178 }
0179 }
0180
0181 }
0182
0183 template <class IDdet>
0184 typename MinL3AlgoUniv<IDdet>::IDmap MinL3AlgoUniv<IDdet>::getSolution(const bool resetsolution) {
0185 IDmap solution;
0186
0187 for (iter_IDmap i = wsum.begin(); i != wsum.end(); i++) {
0188 iter_IDmap iEwsum = Ewsum.find(i->first);
0189 float myValue = 1;
0190 if (i->second != 0)
0191 myValue = iEwsum->second / i->second;
0192
0193 solution.insert(IDmapvalue(i->first, myValue));
0194 }
0195
0196 if (resetsolution)
0197 resetSolution();
0198
0199 return solution;
0200 }
0201
0202 template <class IDdet>
0203 void MinL3AlgoUniv<IDdet>::resetSolution() {
0204 wsum.clear();
0205 Ewsum.clear();
0206 }
0207
0208 template <class IDdet>
0209 std::vector<float> MinL3AlgoUniv<IDdet>::recalibrateEvent(const std::vector<float>& myCluster,
0210 const std::vector<IDdet>& idCluster,
0211 const IDmap& newCalibration) {
0212 std::vector<float> newCluster(myCluster);
0213
0214 for (unsigned i = 0; i < myCluster.size(); i++) {
0215 citer_IDmap icalib = newCalibration.find(idCluster[i]);
0216 if (icalib != newCalibration.end()) {
0217 newCluster[i] *= icalib->second;
0218 } else {
0219 std::cout << "No calibration available for this element." << std::endl;
0220 }
0221 }
0222
0223 return newCluster;
0224 }
0225
0226 #endif