MuonAlignmentAlgorithms/src/MuonDT13ChamberResidual.cc

0001 /*
0002  * $Id: $
0003  */
0004
0005 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonDT13ChamberResidual.h"
0006
0007 MuonDT13ChamberResidual::MuonDT13ChamberResidual(edm::ESHandle<GlobalTrackingGeometry> globalGeometry,
0008                                                  AlignableNavigator *navigator,
0009                                                  DetId chamberId,
0010                                                  AlignableDetOrUnitPtr chamberAlignable)
0011     : MuonHitsChamberResidual(globalGeometry, navigator, chamberId, chamberAlignable) {
0012   m_type = MuonChamberResidual::kDT13;
0013   double rphiAngle = atan2(m_globalGeometry->idToDet(m_chamberId)->position().y(),
0014                            m_globalGeometry->idToDet(m_chamberId)->position().x()) +
0015                      M_PI / 2.;
0016   align::GlobalVector rphiDirection(cos(rphiAngle), sin(rphiAngle), 0.);
0017   m_sign = m_globalGeometry->idToDet(m_chamberId)->toLocal(rphiDirection).x() > 0. ? 1. : -1.;
0018 }
0019
0020 // void MuonDT13ChamberResidual::addResidual(const TrajectoryStateOnSurface *tsos, const TransientTrackingRecHit *hit)
0021
0022 void MuonDT13ChamberResidual::addResidual(edm::ESHandle<Propagator> prop,
0023                                           const TrajectoryStateOnSurface *tsos,
0024                                           const TrackingRecHit *hit,
0025                                           double chamber_width,
0026                                           double chamber_length) {
0027   bool m_debug = false;
0028
0029   m_chamber_width = chamber_width;
0030   m_chamber_length = chamber_length;
0031
0032   DetId id = hit->geographicalId();
0033
0034   align::LocalPoint hitChamberPos =
0035       m_chamberAlignable->surface().toLocal(m_globalGeometry->idToDet(id)->toGlobal(hit->localPosition()));
0036   align::LocalPoint tsosChamberPos =
0037       m_chamberAlignable->surface().toLocal(m_globalGeometry->idToDet(id)->toGlobal(tsos->localPosition()));
0038
0039   if (m_debug) {
0040     std::cout << " MuonDT13ChamberResidual hitChamberPos x: " << hitChamberPos.x()
0041               << " tsosChamberPos x: " << tsosChamberPos.x() << std::endl;
0042     std::cout << "                         hitChamberPos y: " << hitChamberPos.y()
0043               << " tsosChamberPos y: " << tsosChamberPos.y() << std::endl;
0044     std::cout << "                         hitChamberPos z: " << hitChamberPos.z()
0045               << " tsosChamberPos z: " << tsosChamberPos.z() << std::endl;
0046   }
0047
0048   double residual = tsosChamberPos.x() - hitChamberPos.x();  // residual is hit minus hit
0049   double weight = 1. / hit->localPositionError().xx();       // weight linear fit by hit-only local error
0050   double layerPosition = tsosChamberPos.z();                 // the layer's position in the chamber's coordinate system
0051   double layerHitPos = hitChamberPos.z();
0052
0053   m_numHits++;
0054
0055   // "x" is the layerPosition, "y" is the residual (this is a linear fit to residual versus layerPosition)
0056   m_residual_1 += weight;
0057   m_residual_x += weight * layerPosition;
0058   m_residual_y += weight * residual;
0059   m_residual_xx += weight * layerPosition * layerPosition;
0060   m_residual_xy += weight * layerPosition * residual;
0061
0062   // "x" is the layerPosition, "y" is chamberx (this is a linear fit to chamberx versus layerPosition)
0063   m_trackx_1 += weight;
0064   m_trackx_x += weight * layerPosition;
0065   m_trackx_y += weight * tsosChamberPos.x();
0066   m_trackx_xx += weight * layerPosition * layerPosition;
0067   m_trackx_xy += weight * layerPosition * tsosChamberPos.x();
0068
0069   // "x" is the layerPosition, "y" is chambery (this is a linear fit to chambery versus layerPosition)
0070   m_tracky_1 += weight;
0071   m_tracky_x += weight * layerPosition;
0072   m_tracky_y += weight * tsosChamberPos.y();
0073   m_tracky_xx += weight * layerPosition * layerPosition;
0074   m_tracky_xy += weight * layerPosition * tsosChamberPos.y();
0075
0076   m_hitx_1 += weight;
0077   m_hitx_x += weight * layerHitPos;
0078   m_hitx_y += weight * hitChamberPos.x();
0079   m_hitx_xx += weight * layerHitPos * layerHitPos;
0080   m_hitx_xy += weight * layerHitPos * hitChamberPos.x();
0081
0082   m_hity_1 += weight;
0083   m_hity_x += weight * layerHitPos;
0084   m_hity_y += weight * hitChamberPos.y();
0085   m_hity_xx += weight * layerHitPos * layerPosition;
0086   m_hity_xy += weight * layerHitPos * hitChamberPos.y();
0087
0088   m_localIDs.push_back(id);
0089   //  m_localResids.push_back(tsos->localPosition().x() - hit->localPosition().x()); //FIXME looks like this line is not used anywhere, moreover it is wrong for segment-based reconstruction, I changed it to the follwoing line
0090   m_localResids.push_back(residual);
0091   m_individual_x.push_back(layerPosition);
0092   m_individual_y.push_back(residual);
0093   m_individual_weight.push_back(weight);
0094
0095   if (m_numHits > 1)
0096     segment_fit();
0097 }