|
|
|
#include <Alignment/MuonHIPAlignmentRefitter/src/MuonHIPAlignmentRefitter.cc>
Public Member Functions | |
| MuonHIPAlignmentRefitter (const edm::ParameterSet &) | |
| ~MuonHIPAlignmentRefitter () | |
Private Member Functions | |
| virtual void | beginJob (const edm::EventSetup &) |
| virtual void | endJob () |
| virtual void | produce (edm::Event &, const edm::EventSetup &) |
Private Attributes | |
| edm::InputTag | m_muonSource |
| std::string | m_propagatorSource |
| TrackTransformer * | m_trackTransformer |
Implementation: <Notes on="" implementation>="">
Definition at line 57 of file MuonHIPAlignmentRefitter.cc.
| MuonHIPAlignmentRefitter::MuonHIPAlignmentRefitter | ( | const edm::ParameterSet & | iConfig | ) | [explicit] |
Definition at line 86 of file MuonHIPAlignmentRefitter.cc.
References edm::ParameterSet::getParameter(), m_muonSource, m_propagatorSource, and m_trackTransformer.
00087 { 00088 m_muonSource = iConfig.getParameter<edm::InputTag>("MuonSource"); 00089 m_propagatorSource = iConfig.getParameter<std::string>("MuonPropagator"); 00090 00091 m_trackTransformer = new TrackTransformer(iConfig.getParameter<edm::ParameterSet>("TrackerTrackTransformer")); 00092 00093 produces<std::vector<Trajectory> >(); 00094 produces<TrajTrackAssociationCollection>(); 00095 }
| MuonHIPAlignmentRefitter::~MuonHIPAlignmentRefitter | ( | ) |
Definition at line 98 of file MuonHIPAlignmentRefitter.cc.
00099 { 00100 00101 // do anything here that needs to be done at desctruction time 00102 // (e.g. close files, deallocate resources etc.) 00103 00104 }
| void MuonHIPAlignmentRefitter::beginJob | ( | const edm::EventSetup & | ) | [private, virtual] |
| void MuonHIPAlignmentRefitter::produce | ( | edm::Event & | iEvent, | |
| const edm::EventSetup & | iSetup | |||
| ) | [private, virtual] |
Implements edm::EDProducer.
Definition at line 113 of file MuonHIPAlignmentRefitter.cc.
References CSCDetId::chamber(), TrajectoryStateOnSurface::charge(), MuonSubdetId::CSC, TrajectoryStateOnSurface::curvilinearError(), MuonSubdetId::DT, CSCDetId::endcap(), error, TrajectoryMeasurement::forwardPredictedState(), edm::EventSetup::get(), edm::Event::getByLabel(), TrajectoryStateOnSurface::globalDirection(), TrajectoryStateOnSurface::globalPosition(), i, TrajectoryStateOnSurface::isValid(), Trajectory::lastMeasurement(), TrajectoryStateOnSurface::localPosition(), m_muonSource, m_propagatorSource, m_trackTransformer, PV3DBase< T, PVType, FrameType >::mag(), CurvilinearTrajectoryError::matrix(), metsig::muon, DetId::Muon, muons_cfi::muons, funct::pow(), edm::Event::put(), DetId::rawId(), CSCDetId::ring(), MuonSubdetId::RPC, TrackTransformer::setServices(), TrajectoryStateOnSurface::signedInverseMomentum(), funct::sqrt(), CSCDetId::station(), TrajectoryStateOnSurface::surface(), Surface::toGlobal(), GloballyPositioned< T >::toLocal(), TrackTransformer::transform(), PV3DBase< T, PVType, FrameType >::x(), x, TkRotation< T >::xx(), TkRotation< T >::xy(), TkRotation< T >::xz(), PV3DBase< T, PVType, FrameType >::y(), y, TkRotation< T >::yx(), TkRotation< T >::yy(), TkRotation< T >::yz(), PV3DBase< T, PVType, FrameType >::z(), and z.
00114 { 00115 m_trackTransformer->setServices(iSetup); 00116 00117 edm::Handle<reco::MuonCollection> muons; 00118 iEvent.getByLabel(m_muonSource, muons); 00119 00120 edm::ESHandle<Propagator> propagator; 00121 iSetup.get<TrackingComponentsRecord>().get(m_propagatorSource, propagator); 00122 00123 edm::ESHandle<DTGeometry> dtGeometry; 00124 iSetup.get<MuonGeometryRecord>().get(dtGeometry); 00125 00126 edm::ESHandle<CSCGeometry> cscGeometry; 00127 iSetup.get<MuonGeometryRecord>().get(cscGeometry); 00128 00129 edm::ESHandle<MagneticField> magneticField; 00130 iSetup.get<IdealMagneticFieldRecord>().get(magneticField); 00131 00132 for (reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) { 00133 std::vector<Trajectory> trackerTrajectories = m_trackTransformer->transform(*muon->track()); 00134 00135 if (trackerTrajectories.size() == 1) { 00136 const Trajectory trackerTrajectory = *(trackerTrajectories.begin()); 00137 TrajectoryStateOnSurface tracker_tsos = trackerTrajectory.lastMeasurement().forwardPredictedState(); 00138 TrajectoryStateOnSurface last_tsos = trackerTrajectory.lastMeasurement().forwardPredictedState(); 00139 00140 int last_chamber = 0; 00141 std::vector<double> last_x, last_xerr2, last_zforx; 00142 std::vector<double> last_y, last_yerr2, last_zfory; 00143 std::vector<const TrackingRecHit*> last_hits; 00144 00145 for (trackingRecHit_iterator hit = muon->standAloneMuon()->recHitsBegin(); hit != muon->standAloneMuon()->recHitsEnd(); ++hit) { 00146 DetId id = (*hit)->geographicalId(); 00147 00148 if (id.det() == DetId::Muon && id.subdetId() != MuonSubdetId::RPC) { 00149 int chamberId = -1; 00150 if (id.subdetId() == MuonSubdetId::DT) { 00151 DTChamberId dtChamberId(id.rawId()); 00152 chamberId = dtChamberId.rawId(); 00153 } 00154 else if (id.subdetId() == MuonSubdetId::CSC) { 00155 CSCDetId cscChamberId(id.rawId()); 00156 cscChamberId = CSCDetId(cscChamberId.endcap(), cscChamberId.station(), cscChamberId.ring(), cscChamberId.chamber(), 0); 00157 chamberId = cscChamberId.rawId(); 00158 } 00159 else assert(false); 00160 00161 // First chamber: set "last_chamber" variable 00162 if (last_chamber == 0) { 00163 last_chamber = chamberId; 00164 } 00165 // New chamber: re-orient propagator based on linear fit to chamber's 1-12 dof 00166 else if (last_chamber != chamberId) { 00167 GlobalPoint globalPosition = last_tsos.globalPosition(); 00168 GlobalVector globalDirection = last_tsos.globalDirection(); 00169 00170 const Surface* surface; 00171 if (DetId(last_chamber).subdetId() == MuonSubdetId::DT) { 00172 surface = &(dtGeometry->idToDet(DetId(last_chamber))->surface()); 00173 } 00174 else { 00175 surface = &(cscGeometry->idToDet(DetId(last_chamber))->surface()); 00176 } 00177 00178 LocalPoint localPosition = surface->toLocal(globalPosition); 00179 LocalVector localDirection = surface->toLocal(globalDirection); 00180 // we want a normalization in which dz/dz == 1 00181 localDirection /= localDirection.z(); 00182 00183 if (last_zforx.size() > 3) { 00184 // weighted linear fit to x(z) 00185 double x_S = 0; 00186 double x_Sz = 0; 00187 double x_Sx = 0; 00188 double x_Szz = 0; 00189 double x_Szx = 0; 00190 00191 for (unsigned int i = 0; i < last_zforx.size(); i++) { 00192 double x = last_x[i]; 00193 double z = last_zforx[i]; 00194 double x_weight = 1./last_xerr2[i]; 00195 00196 x_S += x_weight; 00197 x_Sz += z * x_weight; 00198 x_Sx += x * x_weight; 00199 x_Szz += z*z * x_weight; 00200 x_Szx += z*x * x_weight; 00201 } 00202 00203 double x_delta = x_S * x_Szz - x_Sz*x_Sz; 00204 double x_intercept = (x_Szz * x_Sx - x_Sz * x_Szx) / x_delta; 00205 // double x_intercept_err2 = x_Szz / x_delta; 00206 double x_slope = (x_S * x_Szx - x_Sz * x_Sx) / x_delta; 00207 // double x_slope_err2 = x_S / x_delta; 00208 // double x_covariance = -x_Sz / x_delta; 00209 00210 double x_on_surface = x_intercept + x_slope * localPosition.z(); 00211 // double x_on_surface_err2 = x_intercept_err2 + localPosition.z()*localPosition.z() * x_slope_err2 + localPosition.z() * x_covariance; 00212 00213 localPosition = LocalPoint(x_on_surface, localPosition.y(), localPosition.z()); 00214 localDirection = LocalVector(x_slope, localDirection.y(), 1); 00215 } 00216 00217 if (last_zfory.size() > 3) { 00218 // weighted linear fit to y(z) 00219 double y_S = 0; 00220 double y_Sz = 0; 00221 double y_Sy = 0; 00222 double y_Szz = 0; 00223 double y_Szy = 0; 00224 for (unsigned int i = 0; i < last_zfory.size(); i++) { 00225 double y = last_y[i]; 00226 double z = last_zfory[i]; 00227 double y_weight = 1./last_yerr2[i]; 00228 00229 y_S += y_weight; 00230 y_Sz += z * y_weight; 00231 y_Sy += y * y_weight; 00232 y_Szz += z*z * y_weight; 00233 y_Szy += z*y * y_weight; 00234 } 00235 00236 double y_delta = y_S * y_Szz - y_Sz*y_Sz; 00237 double y_intercept = (y_Szz * y_Sy - y_Sz * y_Szy) / y_delta; 00238 // double y_intercept_err2 = y_Szz / y_delta; 00239 double y_slope = (y_S * y_Szy - y_Sz * y_Sy) / y_delta; 00240 // double y_slope_err2 = y_S / y_delta; 00241 // double y_covariance = -y_Sz / y_delta; 00242 00243 double y_on_surface = y_intercept + y_slope * localPosition.z(); 00244 // double y_on_surface_err2 = y_intercept_err2 + localPosition.z()*localPosition.z() * y_slope_err2 + localPosition.z() * y_covariance; 00245 00246 localPosition = LocalPoint(localPosition.x(), y_on_surface, localPosition.z()); 00247 localDirection = LocalVector(localDirection.x(), y_slope, 1); 00248 } 00249 00250 if (last_zforx.size() > 3 || last_zfory.size() > 3) { 00251 // make the direction normalized again 00252 localDirection /= sqrt(localDirection.mag()); 00253 00254 globalPosition = surface->toGlobal(localPosition); 00255 globalDirection = surface->toGlobal(localDirection); 00256 GlobalVector globalMomentum = globalDirection * fabs(1./last_tsos.signedInverseMomentum()); 00257 00258 AlgebraicSymMatrix55 error(last_tsos.curvilinearError().matrix()); 00259 GlobalTrajectoryParameters globalTrajectoryParameters(globalPosition, globalMomentum, last_tsos.charge(), &*magneticField); 00260 FreeTrajectoryState freeTrajectoryState(globalTrajectoryParameters, CurvilinearTrajectoryError(error)); 00261 last_tsos = TrajectoryStateOnSurface(freeTrajectoryState, last_tsos.surface()); 00262 } 00263 00264 last_x.clear(); 00265 last_y.clear(); 00266 last_zforx.clear(); 00267 last_zfory.clear(); 00268 last_xerr2.clear(); 00269 last_yerr2.clear(); 00270 last_chamber = chamberId; 00271 00272 last_hits.clear(); 00273 } 00274 00275 // Collect hits and propagate to surface 00276 LocalPoint chamberPoint; 00277 align::RotationType layerRot, chamberRot; 00278 00279 last_hits.push_back(&**hit); 00280 GlobalPoint global_hit_position; 00281 TrajectoryStateOnSurface extrapolation, fromtracker; 00282 if (id.subdetId() == MuonSubdetId::DT) { 00283 extrapolation = propagator->propagate(last_tsos, dtGeometry->idToDet(id)->surface()); 00284 fromtracker = propagator->propagate(tracker_tsos, dtGeometry->idToDet(id)->surface()); 00285 LocalPoint local = (*hit)->localPosition(); 00286 if (extrapolation.isValid()) { 00287 // The extrapolated track is used to set a y position for the hit; 00288 // y positions are always small corrections to chamber coordinates, 00289 // applicable only when layers are rotated inside the chamber 00290 local = LocalPoint(local.x(), extrapolation.localPosition().y(), 0); 00291 } 00292 00293 GlobalPoint globalPoint = dtGeometry->idToDet(id)->surface().toGlobal(local); 00294 chamberPoint = dtGeometry->idToDet(chamberId)->surface().toLocal(globalPoint); 00295 layerRot = dtGeometry->idToDet(id)->surface().rotation(); 00296 chamberRot = dtGeometry->idToDet(chamberId)->surface().rotation(); 00297 global_hit_position = globalPoint; 00298 } 00299 else if (id.subdetId() == MuonSubdetId::CSC) { 00300 extrapolation = propagator->propagate(last_tsos, cscGeometry->idToDet(id)->surface()); 00301 fromtracker = propagator->propagate(tracker_tsos, cscGeometry->idToDet(id)->surface()); 00302 GlobalPoint globalPoint = cscGeometry->idToDet(id)->surface().toGlobal((*hit)->localPosition()); 00303 chamberPoint = cscGeometry->idToDet(chamberId)->surface().toLocal(globalPoint); 00304 layerRot = cscGeometry->idToDet(id)->surface().rotation(); 00305 chamberRot = cscGeometry->idToDet(chamberId)->surface().rotation(); 00306 global_hit_position = globalPoint; 00307 } 00308 else assert(false); 00309 00310 double lRxx = layerRot.xx(); 00311 double lRxy = layerRot.xy(); 00312 double lRxz = layerRot.xz(); 00313 double lRyx = layerRot.yx(); 00314 double lRyy = layerRot.yy(); 00315 double lRyz = layerRot.yz(); 00316 // double lRzx = layerRot.zx(); 00317 // double lRzy = layerRot.zy(); 00318 // double lRzz = layerRot.zz(); 00319 double cRxx = chamberRot.xx(); 00320 double cRxy = chamberRot.xy(); 00321 double cRxz = chamberRot.xz(); 00322 double cRyx = chamberRot.yx(); 00323 double cRyy = chamberRot.yy(); 00324 double cRyz = chamberRot.yz(); 00325 // double cRzx = chamberRot.zx(); 00326 // double cRzy = chamberRot.zy(); 00327 // double cRzz = chamberRot.zz(); 00328 00329 if (id.subdetId() == MuonSubdetId::DT) { 00330 if (fabs(cRxx*lRxx + cRxy*lRxy + cRxz*lRxz) > fabs(cRxx*lRyx + cRxy*lRyy + cRxz*lRyz)) { 00331 last_x.push_back(chamberPoint.x()); 00332 last_xerr2.push_back(pow((cRxx*lRxx + cRxy*lRxy + cRxz*lRxz), 2) * (*hit)->localPositionError().xx() + 00333 pow((cRxx*lRyx + cRxy*lRyy + cRxz*lRyz), 2) * (*hit)->localPositionError().yy()); 00334 last_zforx.push_back(chamberPoint.z()); 00335 } 00336 else { 00337 last_y.push_back(chamberPoint.y()); 00338 last_yerr2.push_back(pow((cRyx*lRxx + cRyy*lRxy + cRyz*lRxz), 2) * (*hit)->localPositionError().xx() + 00339 pow((cRyx*lRyx + cRyy*lRyy + cRyz*lRyz), 2) * (*hit)->localPositionError().yy()); 00340 last_zfory.push_back(chamberPoint.z()); 00341 } 00342 } 00343 else if (id.subdetId() == MuonSubdetId::CSC) { 00344 last_x.push_back(chamberPoint.x()); 00345 last_xerr2.push_back(pow((cRxx*lRxx + cRxy*lRxy + cRxz*lRxz), 2) * (*hit)->localPositionError().xx() + 00346 pow((cRxx*lRyx + cRxy*lRyy + cRxz*lRyz), 2) * (*hit)->localPositionError().yy()); 00347 last_zforx.push_back(chamberPoint.z()); 00348 00349 last_y.push_back(chamberPoint.y()); 00350 last_yerr2.push_back(pow((cRyx*lRxx + cRyy*lRxy + cRyz*lRxz), 2) * (*hit)->localPositionError().xx() + 00351 pow((cRyx*lRyx + cRyy*lRyy + cRyz*lRyz), 2) * (*hit)->localPositionError().yy()); 00352 last_zfory.push_back(chamberPoint.z()); 00353 } 00354 else assert(false); 00355 00356 if (extrapolation.isValid()) { 00357 last_tsos = extrapolation; 00358 } 00359 } // end if Muon and not RPC 00360 } // end loop over standAlone hits 00361 } // end if successfully refit tracker track 00362 } // end loop over muons 00363 00364 std::auto_ptr<std::vector<Trajectory> > trajectoryCollection(new std::vector<Trajectory>); 00365 std::auto_ptr<TrajTrackAssociationCollection> trajTrackMap(new TrajTrackAssociationCollection); 00366 iEvent.put(trajectoryCollection); 00367 iEvent.put(trajTrackMap); 00368 }
Definition at line 69 of file MuonHIPAlignmentRefitter.cc.
Referenced by MuonHIPAlignmentRefitter(), and produce().
std::string MuonHIPAlignmentRefitter::m_propagatorSource [private] |
Definition at line 70 of file MuonHIPAlignmentRefitter.cc.
Referenced by MuonHIPAlignmentRefitter(), and produce().
Definition at line 72 of file MuonHIPAlignmentRefitter.cc.
Referenced by MuonHIPAlignmentRefitter(), and produce().
1.5.4