MuonAlignmentAlgorithms/src/MuonCSCChamberResidual.cc

0001 /*
0002  * $Id: $
0003  */
0004
0005 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonCSCChamberResidual.h"
0006 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0007
0008 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0009 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0010
0011 MuonCSCChamberResidual::MuonCSCChamberResidual(edm::ESHandle<GlobalTrackingGeometry> globalGeometry,
0012                                                AlignableNavigator *navigator,
0013                                                DetId chamberId,
0014                                                AlignableDetOrUnitPtr chamberAlignable)
0015     : MuonHitsChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable) {
0016   m_type = MuonChamberResidual::kCSC;
0017   align::GlobalVector zDirection(0., 0., 1.);
0018   m_sign = m_globalGeometry->idToDet(m_chamberId)->toLocal(zDirection).z() > 0. ? 1. : -1.;
0019 }
0020
0021 void MuonCSCChamberResidual::addResidual(edm::ESHandle<Propagator> prop,
0022                                          const TrajectoryStateOnSurface *tsos,
0023                                          const TrackingRecHit *hit,
0024                                          double chamber_width,
0025                                          double chamber_length) {
0026   bool m_debug = false;
0027
0028   if (m_debug)
0029     std::cout << "MuonCSCChamberResidual::addResidual 1" << std::endl;
0030   DetId id = hit->geographicalId();
0031   if (m_debug)
0032     std::cout << "MuonCSCChamberResidual::addResidual 2" << std::endl;
0033   const CSCGeometry *cscGeometry = dynamic_cast<const CSCGeometry *>(m_globalGeometry->slaveGeometry(id));
0034   if (m_debug)
0035     std::cout << "MuonCSCChamberResidual::addResidual 3" << std::endl;
0036   assert(cscGeometry);
0037
0038   if (m_debug) {
0039     std::cout << " MuonCSCChamberResidual hit->localPosition() x: " << hit->localPosition().x()
0040               << " tsos->localPosition() x: " << tsos->localPosition().x() << std::endl;
0041     std::cout << "                        hit->localPosition() y: " << hit->localPosition().y()
0042               << " tsos->localPosition() y: " << tsos->localPosition().y() << std::endl;
0043     std::cout << "                        hit->localPosition() z: " << hit->localPosition().z()
0044               << " tsos->localPosition() z: " << tsos->localPosition().z() << std::endl;
0045   }
0046
0047   // hit->localPosition() is coordinate in local system of LAYER. Transfer it to coordiante in local system of chamber
0048   align::LocalPoint hitChamberPos =
0049       m_chamberAlignable->surface().toLocal(m_globalGeometry->idToDet(id)->toGlobal(hit->localPosition()));
0050   // TSOS->localPosition() is given in local system of CHAMBER (for segment-based reconstruction)
0051   //  align::LocalPoint tsosChamberPos = tsos->localPosition();
0052   align::LocalPoint tsosChamberPos =
0053       m_chamberAlignable->surface().toLocal(m_globalGeometry->idToDet(id)->toGlobal(tsos->localPosition()));
0054
0055   int strip = cscGeometry->layer(id)->geometry()->nearestStrip(hit->localPosition());
0056   double angle = cscGeometry->layer(id)->geometry()->stripAngle(strip) - M_PI / 2.;
0057   double sinAngle = sin(angle);
0058   double cosAngle = cos(angle);
0059
0060   double residual = cosAngle * (tsosChamberPos.x() - hitChamberPos.x()) +
0061                     sinAngle * (tsosChamberPos.y() - hitChamberPos.y());  // yes, that's +sin()
0062
0063   if (m_debug)
0064     std::cout << " MuonCSCChamberResidual residual: " << residual << std::endl;
0065
0066   double xx = hit->localPositionError().xx();
0067   double xy = hit->localPositionError().xy();
0068   double yy = hit->localPositionError().yy();
0069   double weight = 1. / (xx * cosAngle * cosAngle + 2. * xy * sinAngle * cosAngle + yy * sinAngle * sinAngle);
0070
0071   double layerPosition = tsosChamberPos.z();  // the layer's position in the chamber's coordinate system
0072   double layerHitPos = hitChamberPos.z();
0073
0074   m_numHits++;
0075
0076   // "x" is the layerPosition, "y" is the residual (this is a linear fit to residual versus layerPosition)
0077   m_residual_1 += weight;
0078   m_residual_x += weight * layerPosition;
0079   m_residual_y += weight * residual;
0080   m_residual_xx += weight * layerPosition * layerPosition;
0081   m_residual_xy += weight * layerPosition * residual;
0082
0083   // "x" is the layerPosition, "y" is chamberx (this is a linear fit to chamberx versus layerPosition)
0084   m_trackx_1 += weight;
0085   m_trackx_x += weight * layerPosition;
0086   m_trackx_y += weight * tsosChamberPos.x();
0087   m_trackx_xx += weight * layerPosition * layerPosition;
0088   m_trackx_xy += weight * layerPosition * tsosChamberPos.x();
0089
0090   // "x" is the layerPosition, "y" is chambery (this is a linear fit to chambery versus layerPosition)
0091   m_tracky_1 += weight;
0092   m_tracky_x += weight * layerPosition;
0093   m_tracky_y += weight * tsosChamberPos.y();
0094   m_tracky_xx += weight * layerPosition * layerPosition;
0095   m_tracky_xy += weight * layerPosition * tsosChamberPos.y();
0096
0097   m_hitx_1 += weight;
0098   m_hitx_x += weight * layerHitPos;
0099   m_hitx_y += weight * hitChamberPos.x();
0100   m_hitx_xx += weight * layerHitPos * layerHitPos;
0101   m_hitx_xy += weight * layerHitPos * hitChamberPos.x();
0102
0103   m_hity_1 += weight;
0104   m_hity_x += weight * layerHitPos;
0105   m_hity_y += weight * hitChamberPos.y();
0106   m_hity_xx += weight * layerHitPos * layerHitPos;
0107   m_hity_xy += weight * layerHitPos * hitChamberPos.y();
0108
0109   m_localIDs.push_back(id);
0110   m_localResids.push_back(residual);  // FIXME Check if this method is needed
0111   m_individual_x.push_back(layerPosition);
0112   m_individual_y.push_back(residual);
0113   m_individual_weight.push_back(weight);
0114
0115   if (m_numHits > 1)
0116     segment_fit();
0117 }