CMS 3D CMS Logo

LinearCalibrator.cc
Go to the documentation of this file.
2 #include <cassert>
3 #include <cmath>
4 #include "TVector.h"
5 #include "TDecompLU.h"
6 #include "TDecompSVD.h"
7 
8 #include "TDecompBK.h"
9 
10 #include <iostream>
11 
12 using namespace pftools;
13 
14 /* These utility functions allow you to print matrices and vectors to a stream,
15  * and copy and paste the output directly into Octave or Matlab
16  */
17 void printMat(std::ostream& s, const TMatrixD& input) {
18  s << "\t[";
19  for (int i(0); i < input.GetNrows(); i++) {
20  for (int j(0); j < input.GetNcols(); j++) {
21  s << input[i][j]<< ", \t";
22  }
23  if (i != (input.GetNrows() - 1)) {
24  s << ";\n";
25  s << "\t";
26  }
27  }
28  s << "\t]\n";
29 }
30 
31 void printVec(std::ostream& s, const TVectorD& input) {
32  s << "\t[";
33  for (int i(0); i < input.GetNrows(); i++) {
34  s << input[i];
35  if (i != (input.GetNrows() - 1)) {
36  s << ",\n\t";
37  }
38  }
39  s << "\t]\n";
40 }
41 
43 
44 }
45 
47 }
49  return new LinearCalibrator(*this);
50 }
51 
55 }
56 
58  return new LinearCalibrator();
59 }
60 
61 std::map<DetectorElementPtr, double> LinearCalibrator::getCalibrationCoefficientsCore()
62  noexcept(false) {
63 // std::cout << __PRETTY_FUNCTION__
64 // << ": determining linear calibration coefficients...\n";
65  if (!hasParticles()) {
66  //I have no particles to calibrate to - throw exception.
67  PFToolsException me("Calibrator has no particles for calibration!");
68  throw me;
69  }
70  //std::cout << "\tGetting eij matrix...\n";
71  TMatrixD eij;
72  TVectorD truthE;
73  initEijMatrix(eij, truthE);
74 
75 // std::cout << "\tEij matrix:\n";
76 // printMat(std::cout, eij);
77 
78  //std::cout << "\tGetting projections...\n";
79  TVectorD proj;
80  TMatrixD hess;
81 
82  getProjections(eij, proj, truthE);
83  //std::cout << "\tProjections:\n";
84  //printVec(std::cout, proj);
85  getHessian(eij, hess, truthE);
86 
87  TDecompLU lu;
88  lu.SetMatrix(hess);
89  //std::cout << "\tHessian:\n";
90  //printMat(std::cout, hess);
91 
92  lu.SetTol(1e-25);
93  TMatrixD hessInv = lu.Invert();
94  //std::cout <<"\tInverse Hessian:\n";
95  //printMat(std::cout, hessInv);
96  TVectorD calibsSolved(eij.GetNcols());
97 
98  bool ok(true);
99  calibsSolved = lu.Solve(proj, ok);
100  if (ok) {
101  //std::cout << "\tLU reports ok.\n";
102  }
103  else {
104  std::cout << "\tWARNING: LU reports NOT ok.\n";
105  //This usually happens when you've asked the calibrator to solve for the 'a' term, without including
106  //a dummy 'a' term in each particle deposit.
107  //Make sure you do
108  //Deposition dOffset(offset, eta, phi, 1.0);
109  //particle->addRecDeposition(dOffset);
110  //particle->addTruthDeposition(dOffset);
111  std::cout
112  << "\tDid you forget to add a dummy offset deposition to each particle?\n";
113  std::cout << "\tThrowing an exception!"<< std::endl;
115  me("TDecompLU did not converge successfully when finding calibrations. Did you forget to add a dummy offset deposition to each particle?");
116  throw me;
117  }
118 
119  //std::cout << "\tCalibrations: \n";
120  //printVec(std::cout, calibsSolved);
121 
122  std::map<DetectorElementPtr, double> answers;
123  for (std::map<DetectorElementPtr, unsigned>::iterator it =
124  myDetElIndex.begin(); it != myDetElIndex.end(); ++it) {
125  DetectorElementPtr de = (*it).first;
126  answers[de] = calibsSolved[(*it).second];
127  }
128  return answers;
129 }
130 
131 void LinearCalibrator::initEijMatrix(TMatrixD& eij, TVectorD& truthE) {
132  //std::cout << __PRETTY_FUNCTION__ << "\n";
133  //std::cout << "\tGetting detector element indices...\n";
135  eij.Clear();
136  eij.Zero();
137 
138  truthE.ResizeTo(myParticleDeposits.size());
139  truthE.Zero();
140 
141  //First determine number of calibration constants.
142 
143  eij.ResizeTo(myParticleDeposits.size(), myDetElIndex.size());
144 
145  //loop over all particle deposits
146  unsigned index(0);
147  for (std::vector<ParticleDepositPtr>::const_iterator cit =
148  myParticleDeposits.begin(); cit != myParticleDeposits.end(); ++cit) {
149  ParticleDepositPtr p = *cit;
150  //get each of the relevant detector elements
151 
152  for (std::vector<DetectorElementPtr>::const_iterator detElIt =
153  myDetectorElements.begin(); detElIt != myDetectorElements.end(); ++detElIt) {
154  DetectorElementPtr de = *detElIt;
155  eij[index][myDetElIndex[de]] = p->getRecEnergy(de);
156  //truthE[p->getId()] += p->getTruthEnergy(de);
157  }
158  truthE[index] = p->getTruthEnergy();
159  ++index;
160  }
161 
162 }
163 
164 TVectorD& LinearCalibrator::getProjections(const TMatrixD& eij, TVectorD& proj,
165  const TVectorD& truthE) const {
166  //std::cout << __PRETTY_FUNCTION__ << "\n";
167  proj.ResizeTo(eij.GetNcols());
168  proj.Zero();
169 
170  for (int j(0); j < eij.GetNcols(); ++j) {
171  for (int i(0); i < eij.GetNrows(); ++i) {
172  proj[j] += eij[i][j] / truthE[i];
173  }
174  }
175 
176  return proj;
177 }
178 
179 TMatrixD& LinearCalibrator::getHessian(const TMatrixD& eij, TMatrixD& hess,
180  const TVectorD& truthE) const {
181  //std::cout << __PRETTY_FUNCTION__ << "\n";
182  unsigned nCalibs(eij.GetNcols());
183  hess.ResizeTo(nCalibs, nCalibs);
184  hess.Zero();
185 
186  for (unsigned i(0); i < nCalibs; ++i) {
187  for (unsigned j(0); j < nCalibs; ++j) {
188  for (int n(0); n < eij.GetNrows(); ++n) {
189  hess[i][j] += eij[n][i] * eij[n][j]/ pow(truthE[n], 2.0);
190  }
191  }
192  }
193 
194  return hess;
195 }
196 
198  //reserve index = 0 for the constant term, if we're told to compute it
199  unsigned index(0);
200 
201  myDetElIndex.clear();
202  //loop over known detector elements, and assign a unique row/column index to each
203  for (std::vector<DetectorElementPtr>::const_iterator cit =
204  myDetectorElements.begin(); cit != myDetectorElements.end(); ++cit) {
205  DetectorElementPtr de = *cit;
206  //std::cout << "\tGot element: "<< *de;
207  //check we don't have duplicate detector elements
208  if (myDetElIndex.count(de) == 0) {
209  myDetElIndex[de] = index;
210  ++index;
211  }
212  }
213 
214 }
215 
std::vector< ParticleDepositPtr > myParticleDeposits
Definition: Calibrator.h:61
LinearCalibrator * create() const override
virtual TVectorD & getProjections(const TMatrixD &eij, TVectorD &proj, const TVectorD &truthE) const
#define noexcept
static std::string const input
Definition: EdmProvDump.cc:44
void printVec(std::ostream &s, const TVectorD &input)
LinearCalibrator * clone() const override
boost::shared_ptr< ParticleDeposit > ParticleDepositPtr
int hasParticles() const
Definition: Calibrator.h:46
boost::shared_ptr< DetectorElement > DetectorElementPtr
std::map< DetectorElementPtr, double > getCalibrationCoefficientsCore() noexcept(false) override
General option file parser.
Definition: Calibratable.h:15
General purpose exception class for use by classes in the pftools namespace.
std::vector< DetectorElementPtr > myDetectorElements
Definition: Calibrator.h:60
virtual TMatrixD & getHessian(const TMatrixD &eij, TMatrixD &hess, const TVectorD &truthE) const
virtual void initEijMatrix(TMatrixD &eij, TVectorD &truthE)
std::map< DetectorElementPtr, unsigned > myDetElIndex
void printMat(std::ostream &s, const TMatrixD &input)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40