GeometryCommonDetAlgo/interface/GlobalErrorBase.h

0001 #ifndef GlobalErrorType_H
0002 #define GlobalErrorType_H
0003
0004 #include "DataFormats/GeometryCommonDetAlgo/interface/DeepCopyPointer.h"
0005 #include "DataFormats/Math/interface/AlgebraicROOTObjects.h"
0006 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0007 //
0008 // Exceptions
0009 //
0010 #include "FWCore/Utilities/interface/Exception.h"
0011
0012 /**
0013    * Templated class representing a symmetric 3*3 matrix describing,  according
0014    * to the ErrorWeightType tag, a (cartesian) covariance matrix or the weight
0015    * matrix (the inverse of the covariance matrix).
0016    * \li To have a covariance matrix, the ErrorMatrixTag has to be used, and a
0017    * typedef is available as GlobalError
0018    * \li To have a weight matrix, the WeightMatrixTag has to be used, and a
0019    * typedef is available as Globalweight
0020    *
0021    * The typedefs should be used in the code.
0022    */
0023
0024 template <class T, class ErrorWeightType>
0025 class GlobalErrorBase {
0026 public:
0027   /// Tag to request a null error matrix
0028   class NullMatrix {};
0029
0030   /**
0031    * Default constructor, creating a null 3*3 matrix (all values are 0)
0032    */
0033   GlobalErrorBase() {}
0034
0035   /**
0036    * Obsolete  Constructor that allocates a null GlobalErrorBase (it does not create the error matrix at all)
0037    */
0038   GlobalErrorBase(const NullMatrix&) {}
0039
0040   /**
0041    * Constructor.
0042    * The symmetric matrix stored as a lower triangular matrix
0043    */
0044   GlobalErrorBase(T c11, T c21, T c22, T c31, T c32, T c33) {
0045     theCartesianError(0, 0) = c11;
0046     theCartesianError(1, 0) = c21;
0047     theCartesianError(1, 1) = c22;
0048     theCartesianError(2, 0) = c31;
0049     theCartesianError(2, 1) = c32;
0050     theCartesianError(2, 2) = c33;
0051     theCartesianError(3, 0) = 0.;
0052     theCartesianError(3, 1) = 0.;
0053     theCartesianError(3, 2) = 0.;
0054     theCartesianError(3, 3) = 0.;
0055   }
0056
0057   /**
0058    * Constructor.
0059    * The symmetric matrix stored as a lower triangular matrix (4D)
0060    */
0061   GlobalErrorBase(T c11, T c21, T c22, T c31, T c32, T c33, T c41, T c42, T c43, T c44) {
0062     theCartesianError(0, 0) = c11;
0063     theCartesianError(1, 0) = c21;
0064     theCartesianError(1, 1) = c22;
0065     theCartesianError(2, 0) = c31;
0066     theCartesianError(2, 1) = c32;
0067     theCartesianError(2, 2) = c33;
0068     theCartesianError(3, 0) = c41;
0069     theCartesianError(3, 1) = c42;
0070     theCartesianError(3, 2) = c43;
0071     theCartesianError(3, 3) = c44;
0072   }
0073
0074   /**
0075    * Constructor from SymMatrix. The original matrix has to be a 3*3 matrix.
0076    */
0077   GlobalErrorBase(const AlgebraicSymMatrix33& err) {
0078     theCartesianError(0, 0) = err(0, 0);
0079     theCartesianError(1, 0) = err(1, 0);
0080     theCartesianError(1, 1) = err(1, 1);
0081     theCartesianError(2, 0) = err(2, 0);
0082     theCartesianError(2, 1) = err(2, 1);
0083     theCartesianError(2, 2) = err(2, 2);
0084     theCartesianError(3, 0) = 0.;
0085     theCartesianError(3, 1) = 0.;
0086     theCartesianError(3, 2) = 0.;
0087     theCartesianError(3, 3) = 0.;
0088   }
0089
0090   /**
0091    * Constructor from SymMatrix. The original matrix has to be a 4*4 matrix.
0092    */
0093   GlobalErrorBase(const AlgebraicSymMatrix44& err) : theCartesianError(err) {}
0094
0095   ~GlobalErrorBase() {}
0096
0097   T cxx() const { return theCartesianError(0, 0); }
0098
0099   T cyx() const { return theCartesianError(1, 0); }
0100
0101   T cyy() const { return theCartesianError(1, 1); }
0102
0103   T czx() const { return theCartesianError(2, 0); }
0104
0105   T czy() const { return theCartesianError(2, 1); }
0106
0107   T czz() const { return theCartesianError(2, 2); }
0108
0109   T ctx() const { return theCartesianError(3, 0); }
0110
0111   T cty() const { return theCartesianError(3, 1); }
0112
0113   T ctz() const { return theCartesianError(3, 2); }
0114
0115   T ctt() const { return theCartesianError(3, 3); }
0116
0117   /**
0118    * Access method to the matrix,
0119    * /return The SymMatrix
0120    */
0121   const AlgebraicSymMatrix33 matrix() const {
0122     AlgebraicSymMatrix33 out;
0123     out(0, 0) = theCartesianError(0, 0);
0124     out(1, 0) = theCartesianError(1, 0);
0125     out(1, 1) = theCartesianError(1, 1);
0126     out(2, 0) = theCartesianError(2, 0);
0127     out(2, 1) = theCartesianError(2, 1);
0128     out(2, 2) = theCartesianError(2, 2);
0129     return out;
0130   }
0131
0132   /**
0133    * Access method to the matrix,
0134    * /return The SymMatrix 4x4
0135    */
0136   const AlgebraicSymMatrix44& matrix4D() const { return theCartesianError; }
0137
0138   T rerr(const GlobalPoint& aPoint) const {
0139     T r2 = aPoint.perp2();
0140     T x2 = aPoint.x() * aPoint.x();
0141     T y2 = aPoint.y() * aPoint.y();
0142     T xy = aPoint.x() * aPoint.y();
0143     if (r2 != 0)
0144       return std::max<T>(0, (1. / r2) * (x2 * cxx() + 2. * xy * cyx() + y2 * cyy()));
0145     else
0146       return 0.5 * (cxx() + cyy());
0147   }
0148
0149   T phierr(const GlobalPoint& aPoint) const {
0150     T r2 = aPoint.perp2();
0151     T x2 = aPoint.x() * aPoint.x();
0152     T y2 = aPoint.y() * aPoint.y();
0153     T xy = aPoint.x() * aPoint.y();
0154     if (r2 != 0)
0155       return std::max<T>(0, (1. / (r2 * r2)) * (y2 * cxx() - 2. * xy * cyx() + x2 * cyy()));
0156     else
0157       return 0;
0158   }
0159
0160   GlobalErrorBase operator+(const GlobalErrorBase& err) const {
0161     return GlobalErrorBase(theCartesianError + err.theCartesianError);
0162   }
0163   GlobalErrorBase operator-(const GlobalErrorBase& err) const {
0164     return GlobalErrorBase(theCartesianError - err.theCartesianError);
0165   }
0166
0167 private:
0168   AlgebraicSymMatrix44 theCartesianError;
0169 };
0170
0171 #endif