MuonAlignmentAlgorithms/src/MuonDT2ChamberResidual.cc

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