00001 #include <math.h>
00002 #include "DataFormats/TrackReco/interface/HitPattern.h"
00003 #include "DataFormats/TrackReco/interface/TrackResiduals.h"
00004 using namespace reco;
00005
00006 TrackResiduals::TrackResiduals () : residualType(X_Y_RESIDUALS)
00007 {
00008 memset(residuals_, 0, sizeof(residuals_));
00009 }
00010
00011 TrackResiduals::TrackResiduals (enum ResidualType type) : residualType(type)
00012 {
00013 memset(residuals_, 0, sizeof(residuals_));
00014 }
00015
00016 void TrackResiduals::setResidualType (enum ResidualType type)
00017 {
00018 residualType = type;
00019 }
00020
00021 void TrackResiduals::setResidualXY (int idx, double residualX, double residualY)
00022 {
00023 assert(residualType == X_Y_RESIDUALS);
00024 residuals_[idx] = (pack_residual(residualX) << 4) | pack_residual(residualY);
00025 }
00026
00027 double TrackResiduals::residualX (int i) const
00028 {
00029 switch (residualType) {
00030 case X_Y_RESIDUALS:
00031 return unpack_residual(residuals_[i] >> 4);
00032 case X_Y_PULLS:
00033 return unpack_pull(residuals_[i] >> 4);
00034 default:
00035 assert(0);
00036 }
00037 return 0;
00038 }
00039
00040 double TrackResiduals::residualY (int i) const
00041 {
00042 switch (residualType) {
00043 case X_Y_RESIDUALS:
00044 return unpack_residual(residuals_[i] & 0x0f);
00045 case X_Y_PULLS:
00046 return unpack_pull(residuals_[i] & 0x0f);
00047 default:
00048 assert(0);
00049 }
00050 return 0;
00051 }
00052
00053 static int index_to_hitpattern (int i_hitpattern, const HitPattern &h)
00054 {
00055 int i_residuals = 0;
00056 assert(i_hitpattern < h.numberOfHits());
00057 if (!h.validHitFilter(h.getHitPattern(i_hitpattern)))
00058
00059 return -999;
00060 for (int i = 0; i <= i_hitpattern; i++) {
00061 if (h.validHitFilter(h.getHitPattern(i)))
00062 i_residuals++;
00063 }
00064 assert(i_residuals > 0);
00065 return i_residuals - 1;
00066 }
00067
00068 double TrackResiduals::residualX (int i, const HitPattern &h) const
00069 {
00070 int idx = index_to_hitpattern(i, h);
00071 if (idx == -999)
00072 return -999;
00073 return residualX(idx);
00074 }
00075
00076 double TrackResiduals::residualY (int i, const HitPattern &h) const
00077 {
00078 int idx = index_to_hitpattern(i, h);
00079 if (idx == -999)
00080 return -999;
00081 return residualY(idx);
00082 }
00083
00084 void TrackResiduals::setPullXY (int idx, double pullX, double pullY)
00085 {
00086 assert(residualType == X_Y_PULLS);
00087 residuals_[idx] = (pack_pull(pullX) << 4) | pack_pull(pullY);
00088 }
00089
00090 static const double pull_char_to_double[8][2] = {
00091 { 0, 0.5 },
00092 { 0.5, 1 },
00093 { 1, 1.5 },
00094 { 1.5, 2 },
00095 { 2, 2.5 },
00096 { 2.5, 3.5 },
00097 { 3.5, 4.5 },
00098 { 4.5, 5.5 },
00099 };
00100
00101 double TrackResiduals::unpack_pull (unsigned char pull)
00102 {
00103 int sgn = 1 - 2 * ((pull & 0x08) >> 3);
00104 unsigned char mag = pull & 0x07;
00105 return sgn *
00106 (pull_char_to_double[mag][0] + pull_char_to_double[mag][1]) / 2;
00107 }
00108
00109 unsigned char TrackResiduals::pack_pull (double pull)
00110 {
00111 unsigned char sgn = (pull < 0) * 0x08;
00112 int mag = -1;
00113 while (++mag < 8 && pull_char_to_double[mag][1] < fabs(pull));
00114 return sgn + mag;
00115 }
00116
00117 static const double residual_char_to_double[8][2] = {
00118 { 0, 0.5 },
00119 { 0.5, 1 },
00120 { 1, 1.5 },
00121 { 1.5, 2 },
00122 { 2, 2.5 },
00123 { 2.5, 3.5 },
00124 { 3.5, 4.5 },
00125 { 4.5, 5.5 },
00126 };
00127
00128 double TrackResiduals::unpack_residual (unsigned char pull)
00129 {
00130 int sgn = 1 - 2 * ((pull & 0x08) >> 3);
00131 unsigned char mag = pull & 0x07;
00132 return sgn *
00133 (pull_char_to_double[mag][0] + pull_char_to_double[mag][1]) / 2;
00134 }
00135
00136 unsigned char TrackResiduals::pack_residual (double pull)
00137 {
00138 unsigned char sgn = (pull < 0) * 0x08;
00139 int mag = -1;
00140 while (++mag < 8 && pull_char_to_double[mag][1] < fabs(pull));
00141 return sgn + mag;
00142 }
00143
00144 void TrackResiduals::print (std::ostream &stream) const
00145 {
00146 stream << "TrackResiduals" << std::endl;
00147 std::ios_base::fmtflags flags = stream.flags();
00148 stream.setf ( std::ios_base::hex, std::ios_base::basefield );
00149 stream.setf ( std::ios_base::showbase );
00150 for (int i = 0; i < numResiduals; i++) {
00151 unsigned char residual = residuals_[i];
00152 printf("0x%x\n", residual);
00153
00154 }
00155 stream.flags(flags);
00156 }
00157
00158 void TrackResiduals::print (const HitPattern &h, std::ostream &stream) const
00159 {
00160 stream << "TrackResiduals" << std::endl;
00161 for (int i = 0; i < h.numberOfHits(); i++) {
00162 stream << (h.validHitFilter(h.getHitPattern(i)) ?
00163 "valid hit: " : "invalid hit: ") <<
00164 "( " << residualX(i, h) << " , " << residualY(i, h) << " )" <<
00165 std::endl;
00166 }
00167 }