File indexing completed on 2024-04-06 11:59:29
0001
0002
0003
0004
0005
0006
0007
0008 #include "Calibration/Tools/interface/GenericMinL3Algorithm.h"
0009
0010 GenericMinL3Algorithm::GenericMinL3Algorithm(bool normalise)
0011 : normaliseFlag(normalise)
0012
0013 {}
0014
0015 GenericMinL3Algorithm::~GenericMinL3Algorithm() {}
0016
0017 std::vector<float> GenericMinL3Algorithm::iterate(const std::vector<std::vector<float> >& eventMatrix,
0018 const std::vector<float>& energyVector,
0019 const int nIter) {
0020 std::vector<float> solution;
0021 std::vector<std::vector<float> > myEventMatrix(eventMatrix);
0022 int Nevents = eventMatrix.size();
0023 int Nchannels = eventMatrix[0].size();
0024 std::vector<float> theCalibVector(Nchannels, 1.);
0025
0026
0027 for (int iter = 1; iter <= nIter; iter++) {
0028
0029 solution = iterate(myEventMatrix, energyVector);
0030
0031 if (solution.empty())
0032 return solution;
0033
0034
0035
0036 for (int i = 0; i < Nchannels; i++) {
0037 for (int ievent = 0; ievent < Nevents; ievent++) {
0038 myEventMatrix[ievent][i] *= solution[i];
0039 }
0040
0041 theCalibVector[i] *= solution[i];
0042 }
0043
0044 }
0045
0046 return theCalibVector;
0047 }
0048
0049 std::vector<float> GenericMinL3Algorithm::iterate(const std::vector<std::vector<float> >& eventMatrix,
0050 const std::vector<float>& energyVector) {
0051 std::vector<float> solution;
0052 std::vector<float> myEnergyVector(energyVector);
0053
0054 int Nevents = eventMatrix.size();
0055 int Nchannels = eventMatrix[0].size();
0056
0057
0058 if (Nevents != int(myEnergyVector.size())) {
0059 std::cout << "GenericMinL3Algorithm::iterate(): Error: bad matrix dimensions. Dropping out." << std::endl;
0060 return solution;
0061 }
0062
0063
0064 solution.assign(Nchannels, 1.);
0065
0066 int ievent, i, j;
0067
0068
0069 float sumOverEnergy;
0070 if (normaliseFlag) {
0071 float scale = 0.;
0072
0073 std::cout << "GenericMinL3Algorithm::iterate(): Normalising event data" << std::endl;
0074
0075 for (i = 0; i < Nevents; i++) {
0076 sumOverEnergy = 0.;
0077 for (j = 0; j < Nchannels; j++) {
0078 sumOverEnergy += eventMatrix[i][j];
0079 }
0080 sumOverEnergy /= myEnergyVector[i];
0081 scale += sumOverEnergy;
0082 }
0083 scale /= Nevents;
0084 std::cout << " Normalisation = " << scale << std::endl;
0085
0086 for (i = 0; i < Nevents; i++) {
0087 myEnergyVector[i] *= scale;
0088 }
0089 }
0090
0091
0092 float sum25, invsum25;
0093 float w;
0094 std::vector<float> wsum(Nchannels, 0.);
0095 std::vector<float> Ewsum(Nchannels, 0.);
0096
0097
0098 for (ievent = 0; ievent < Nevents; ievent++) {
0099
0100 sum25 = 0.;
0101
0102 for (i = 0; i < Nchannels; i++) {
0103 sum25 += eventMatrix[ievent][i];
0104 }
0105
0106 if (sum25 != 0.) {
0107 invsum25 = 1 / sum25;
0108
0109 for (i = 0; i < Nchannels; i++) {
0110 w = eventMatrix[ievent][i] * invsum25;
0111 wsum[i] += w;
0112 Ewsum[i] += (w * myEnergyVector[ievent] * invsum25);
0113 }
0114 } else {
0115 std::cout << " Debug: dropping null event: " << ievent << std::endl;
0116 }
0117
0118 }
0119
0120
0121 for (i = 0; i < Nchannels; i++) {
0122 if (wsum[i] != 0.) {
0123 solution[i] *= Ewsum[i] / wsum[i];
0124 } else {
0125 std::cout << "warning - no event data for crystal index " << i << std::endl;
0126 }
0127 }
0128
0129 return solution;
0130 }