Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GblData.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GblData.cc
1 /*
2  * GblData.cpp
3  *
4  * Created on: Aug 18, 2011
5  * Author: kleinwrt
6  */
7 
30 #include "GblData.h"
31 
33 namespace gbl {
34 
36 
41 GblData::GblData(unsigned int aLabel, double aValue, double aPrec) :
42  theLabel(aLabel), theValue(aValue), thePrecision(aPrec), theDownWeight(
43  1.), thePrediction(0.), theParameters(), theDerivatives(), globalLabels(), globalDerivatives() {
44 
45 }
46 
48 }
49 
51 
63 void GblData::addDerivatives(unsigned int iRow,
64  const std::vector<unsigned int> &labDer, const SMatrix55 &matDer,
65  unsigned int iOff, const TMatrixD &derLocal,
66  const std::vector<int> &labGlobal, const TMatrixD &derGlobal,
67  unsigned int extOff, const TMatrixD &extDer) {
68 
69  unsigned int nParMax = 5 + derLocal.GetNcols() + extDer.GetNcols();
70  theParameters.reserve(nParMax); // have to be sorted
71  theDerivatives.reserve(nParMax);
72 
73  for (int i = 0; i < derLocal.GetNcols(); ++i) // local derivatives
74  {
75  if (derLocal(iRow - iOff, i)) {
76  theParameters.push_back(i + 1);
77  theDerivatives.push_back(derLocal(iRow - iOff, i));
78  }
79  }
80 
81  for (int i = 0; i < extDer.GetNcols(); ++i) // external derivatives
82  {
83  if (extDer(iRow - iOff, i)) {
84  theParameters.push_back(extOff + i + 1);
85  theDerivatives.push_back(extDer(iRow - iOff, i));
86  }
87  }
88 
89  for (unsigned int i = 0; i < 5; ++i) // curvature, offset derivatives
90  {
91  if (labDer[i] and matDer(iRow, i)) {
92  theParameters.push_back(labDer[i]);
93  theDerivatives.push_back(matDer(iRow, i));
94  }
95  }
96 
97  globalLabels = labGlobal;
98  for (int i = 0; i < derGlobal.GetNcols(); ++i) // global derivatives
99  globalDerivatives.push_back(derGlobal(iRow - iOff, i));
100 }
101 
103 
111 void GblData::addDerivatives(unsigned int iRow,
112  const std::vector<unsigned int> &labDer, const SMatrix27 &matDer,
113  unsigned int extOff, const TMatrixD &extDer) {
114 
115  unsigned int nParMax = 7 + extDer.GetNcols();
116  theParameters.reserve(nParMax); // have to be sorted
117  theDerivatives.reserve(nParMax);
118 
119  for (int i = 0; i < extDer.GetNcols(); ++i) // external derivatives
120  {
121  if (extDer(iRow, i)) {
122  theParameters.push_back(extOff + i + 1);
123  theDerivatives.push_back(extDer(iRow, i));
124  }
125  }
126 
127  for (unsigned int i = 0; i < 7; ++i) // curvature, offset derivatives
128  {
129  if (labDer[i] and matDer(iRow, i)) {
130  theParameters.push_back(labDer[i]);
131  theDerivatives.push_back(matDer(iRow, i));
132  }
133  }
134 }
135 
137 
142 void GblData::addDerivatives(const std::vector<unsigned int> &index,
143  const std::vector<double> &derivatives) {
144  for (unsigned int i = 0; i < derivatives.size(); ++i) // any derivatives
145  {
146  if (derivatives[i]) {
147  theParameters.push_back(index[i]);
148  theDerivatives.push_back(derivatives[i]);
149  }
150  }
151 }
152 
154 void GblData::setPrediction(const VVector &aVector) {
155 
156  thePrediction = 0.;
157  for (unsigned int i = 0; i < theDerivatives.size(); ++i) {
158  thePrediction += theDerivatives[i] * aVector(theParameters[i] - 1);
159  }
160 }
161 
163 
166 double GblData::setDownWeighting(unsigned int aMethod) {
167 
168  double aWeight = 1.;
169  double scaledResidual = fabs(theValue - thePrediction) * sqrt(thePrecision);
170  if (aMethod == 1) // Tukey
171  {
172  if (scaledResidual < 4.6851) {
173  aWeight = (1.0 - 0.045558 * scaledResidual * scaledResidual);
174  aWeight *= aWeight;
175  } else {
176  aWeight = 0.;
177  }
178  } else if (aMethod == 2) //Huber
179  {
180  if (scaledResidual >= 1.345) {
181  aWeight = 1.345 / scaledResidual;
182  }
183  } else if (aMethod == 3) //Cauchy
184  {
185  aWeight = 1.0 / (1.0 + (scaledResidual * scaledResidual / 5.6877));
186  }
187  theDownWeight = aWeight;
188  return aWeight;
189 }
190 
192 
195 double GblData::getChi2() const {
196  double aDiff = theValue - thePrediction;
197  return aDiff * aDiff * thePrecision * theDownWeight;
198 }
199 
201 void GblData::printData() const {
202 
203  std::cout << " measurement at label " << theLabel << ": " << theValue
204  << ", " << thePrecision << std::endl;
205  std::cout << " param " << theParameters.size() << ":";
206  for (unsigned int i = 0; i < theParameters.size(); ++i) {
207  std::cout << " " << theParameters[i];
208  }
209  std::cout << std::endl;
210  std::cout << " deriv " << theDerivatives.size() << ":";
211  for (unsigned int i = 0; i < theDerivatives.size(); ++i) {
212  std::cout << " " << theDerivatives[i];
213  }
214  std::cout << std::endl;
215 }
216 
218 
224 void GblData::getLocalData(double &aValue, double &aWeight,
225  std::vector<unsigned int>* &indLocal, std::vector<double>* &derLocal) {
226 
227  aValue = theValue;
228  aWeight = thePrecision * theDownWeight;
229  indLocal = &theParameters;
230  derLocal = &theDerivatives;
231 }
232 
234 
242 void GblData::getAllData(double &aValue, double &aErr,
243  std::vector<unsigned int>* &indLocal, std::vector<double>* &derLocal,
244  std::vector<int>* &labGlobal, std::vector<double>* &derGlobal) {
245  aValue = theValue;
246  aErr = 1.0 / sqrt(thePrecision);
247  indLocal = &theParameters;
248  derLocal = &theDerivatives;
249  labGlobal = &globalLabels;
250  derGlobal = &globalDerivatives;
251 }
252 
254 
261 void GblData::getResidual(double &aResidual, double &aVariance,
262  double &aDownWeight, std::vector<unsigned int>* &indLocal,
263  std::vector<double>* &derLocal) {
264  aResidual = theValue - thePrediction;
265  aVariance = 1.0 / thePrecision;
266  aDownWeight = theDownWeight;
267  indLocal = &theParameters;
268  derLocal = &theDerivatives;
269 }
270 }