00001
00017 #include "RecoMuon/GlobalTrackingTools/interface/DynamicTruncation.h"
00018 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
00019 #include "TrackingTools/TrackFitters/interface/RecHitLessByDet.h"
00020 #include "RecoMuon/GlobalTrackingTools/interface/ChamberSegmentUtility.h"
00021 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
00022 #include "RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h"
00023 #include "RecoMuon/Records/interface/MuonRecoGeometryRecord.h"
00024 #include "RecoMuon/Navigation/interface/MuonNavigationSchool.h"
00025 #include "TrackingTools/DetLayers/interface/NavigationSetter.h"
00026 #include "RecoMuon/Navigation/interface/MuonNavigationPrinter.h"
00027 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00028 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
00029 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
00030 #include "DataFormats/TrackingRecHit/interface/AlignmentPositionError.h"
00031 #include "DataFormats/GeometryCommonDetAlgo/interface/ErrorFrameTransformer.h"
00032
00033 #include "CondFormats/Alignment/interface/Alignments.h"
00034 #include "CondFormats/Alignment/interface/AlignTransform.h"
00035 #include "CondFormats/AlignmentRecord/interface/DTAlignmentRcd.h"
00036 #include "CondFormats/AlignmentRecord/interface/CSCAlignmentRcd.h"
00037 #include "CondFormats/Alignment/interface/AlignmentErrors.h"
00038 #include "CondFormats/Alignment/interface/AlignTransformError.h"
00039 #include "CondFormats/AlignmentRecord/interface/DTAlignmentErrorRcd.h"
00040 #include "CondFormats/AlignmentRecord/interface/CSCAlignmentErrorRcd.h"
00041
00042
00043 #define MAX_THR 1e7
00044
00045 using namespace edm;
00046 using namespace std;
00047 using namespace reco;
00048
00049 DynamicTruncation::DynamicTruncation(const edm::Event& event, const MuonServiceProxy& theService):
00050 DTThr(0), CSCThr(0), useAPE(false)
00051 {
00052 theEvent = &event;
00053 theSetup = &theService.eventSetup();
00054 propagator = theService.propagator("SmartPropagator");
00055 propagatorCompatibleDet = theService.propagator("SteppingHelixPropagatorAny");
00056 theG = theService.trackingGeometry();
00057 theService.eventSetup().get<TransientRecHitRecord>().get("MuonRecHitBuilder",theMuonRecHitBuilder);
00058 theService.eventSetup().get<TrackingComponentsRecord>().get("KFUpdator",updatorHandle);
00059 theService.eventSetup().get<MuonGeometryRecord>().get(cscGeom);
00060 theService.eventSetup().get<MuonRecoGeometryRecord>().get(navMuon);
00061 theService.eventSetup().get<IdealMagneticFieldRecord>().get(magfield);
00062 navigation = new DirectMuonNavigation(theService.detLayerGeometry());
00063 }
00064
00065
00066
00067 DynamicTruncation::~DynamicTruncation() {
00068 if (navigation) delete navigation;
00069 }
00070
00071
00072
00073 TransientTrackingRecHit::ConstRecHitContainer DynamicTruncation::filter(const Trajectory& traj) {
00074 result.clear();
00075
00076
00077
00078 edm::ESHandle<AlignmentErrors> dtAlignmentErrors;
00079 theSetup->get<DTAlignmentErrorRcd>().get( dtAlignmentErrors );
00080 for ( std::vector<AlignTransformError>::const_iterator it = dtAlignmentErrors->m_alignError.begin();
00081 it != dtAlignmentErrors->m_alignError.end(); it++ ) {
00082 CLHEP::HepSymMatrix error = (*it).matrix();
00083 GlobalError glbErr(error[0][0], error[1][0], error[1][1], error[2][0], error[2][1], error[2][2]);
00084 DTChamberId DTid((*it).rawId());
00085 dtApeMap.insert( pair<DTChamberId, GlobalError> (DTid, glbErr) );
00086 }
00087
00088
00089
00090 edm::ESHandle<AlignmentErrors> cscAlignmentErrors;
00091 theSetup->get<CSCAlignmentErrorRcd>().get( cscAlignmentErrors );
00092 for ( std::vector<AlignTransformError>::const_iterator it = cscAlignmentErrors->m_alignError.begin();
00093 it != cscAlignmentErrors->m_alignError.end(); it++ ) {
00094 CLHEP::HepSymMatrix error = (*it).matrix();
00095 GlobalError glbErr(error[0][0], error[1][0], error[1][1], error[2][0], error[2][1], error[2][2]);
00096 CSCDetId CSCid((*it).rawId());
00097 cscApeMap.insert( pair<CSCDetId, GlobalError> (CSCid, glbErr) );
00098 }
00099
00100
00101 std::vector<TrajectoryMeasurement> muonMeasurements = traj.measurements();
00102 TrajectoryMeasurement lastTKm = muonMeasurements.front();
00103 for (std::vector<TrajectoryMeasurement>::const_iterator imT = muonMeasurements.begin(); imT != muonMeasurements.end(); imT++ ) {
00104 if ( !(*imT).recHit()->isValid() ) continue;
00105 const TransientTrackingRecHit* hit = &(*(*imT).recHit());
00106 if (hit->geographicalId().det() == DetId::Tracker) {
00107 result.push_back((*imT).recHit());
00108 if (!(*imT).forwardPredictedState().isValid()) continue;
00109 if ((*imT).forwardPredictedState().globalPosition().mag()
00110 > lastTKm.forwardPredictedState().globalPosition().mag()) lastTKm = *imT;
00111 }
00112 }
00113
00114
00115 currentState = lastTKm.forwardPredictedState();
00116
00117
00118 update(currentState, lastTKm.recHit());
00119
00120
00121 map<int, std::vector<DetId> > detMap;
00122 compatibleDets(currentState, detMap);
00123
00124
00125 filteringAlgo(detMap);
00126
00127 return result;
00128 }
00129
00130
00131
00132 void DynamicTruncation::setThr(int bThr, int eThr, int useAPE_) {
00133 if (useAPE_ == 1) useAPE = true;
00134 else useAPE = false;
00135 if (bThr <= MAX_THR && bThr >= 0) DTThr = bThr;
00136 else DTThr = MAX_THR;
00137 if (eThr <= MAX_THR && eThr >= 0) CSCThr = eThr;
00138 else CSCThr = MAX_THR;
00139 }
00140
00141
00142
00143 double DynamicTruncation::getBest(std::vector<CSCSegment>& segs, TrajectoryStateOnSurface& tsos, CSCSegment& bestCSCSeg) {
00144 unsigned int i = 0;
00145 double val = MAX_THR;
00146 std::vector<CSCSegment>::size_type sz = segs.size();
00147 for (i=0; i<sz; i++) {
00148 LocalError apeLoc;
00149 if (useAPE) apeLoc = ErrorFrameTransformer().transform(cscApeMap.find(segs[i].cscDetId())->second, theG->idToDet(segs[i].cscDetId())->surface());
00150 StateSegmentMatcher estim(&tsos, &segs[i], &apeLoc);
00151 double tmp = estim.value();
00152 if (tmp < val) {
00153 bestCSCSeg = segs[i];
00154 val = tmp;
00155 }
00156 }
00157 return val;
00158 }
00159
00160
00161
00162 double DynamicTruncation::getBest(std::vector<DTRecSegment4D>& segs, TrajectoryStateOnSurface& tsos, DTRecSegment4D& bestDTSeg) {
00163 unsigned int i = 0;
00164 double val = MAX_THR;
00165 std::vector<DTRecSegment4D>::size_type sz = segs.size();
00166 for (i=0; i<sz; i++) {
00167 LocalError apeLoc;
00168 if (useAPE) apeLoc = ErrorFrameTransformer().transform(dtApeMap.find(segs[i].chamberId())->second, theG->idToDet(segs[i].chamberId())->surface());
00169 StateSegmentMatcher estim(&tsos, &segs[i], &apeLoc);
00170 double tmp = estim.value();
00171 if (tmp < val) {
00172 bestDTSeg = segs[i];
00173 val = tmp;
00174 }
00175 }
00176 return val;
00177 }
00178
00179
00180
00181 void DynamicTruncation::compatibleDets(TrajectoryStateOnSurface& tsos, map<int, std::vector<DetId> >& detMap) {
00182
00183
00184 MeasurementEstimator *theEstimator = new Chi2MeasurementEstimator(100, 3);
00185 std::vector<const DetLayer *> navLayers;
00186 navLayers = navigation->compatibleLayers(*(currentState.freeState()), alongMomentum);
00187 for ( unsigned int ilayer=0; ilayer<navLayers.size(); ilayer++ ) {
00188
00189 if (navLayers[ilayer]->subDetector() == GeomDetEnumerators::RPCEndcap
00190 || navLayers[ilayer]->subDetector() == GeomDetEnumerators::RPCBarrel) continue;
00191 std::vector<DetLayer::DetWithState> comps = navLayers[ilayer]->compatibleDets(currentState,
00192 *propagatorCompatibleDet, *theEstimator);
00193
00194
00195 if (comps.size() > 0) {
00196 for ( unsigned int icomp=0; icomp<comps.size(); icomp++ ) {
00197 DetId id(comps[icomp].first->geographicalId().rawId());
00198 detMap[ilayer].push_back(id);
00199 }
00200 }
00201 }
00202 if (theEstimator) delete theEstimator;
00203 }
00204
00205
00206
00207 void DynamicTruncation::filteringAlgo(map<int, std::vector<DetId> >& detMap) {
00208 ChamberSegmentUtility getSegs(*theEvent, *theSetup);
00209 for (unsigned int iDet = 0; iDet < detMap.size(); ++iDet) {
00210 double bestLayerValue = MAX_THR;
00211 bool isDTorCSC = false;
00212 ConstRecHitContainer layerRH, layerSEG;
00213 std::vector<DetId> chamber = detMap[iDet];
00214 for (unsigned int j = 0; j < chamber.size(); ++j) {
00215 DetId id = chamber[j];
00216
00217 if (id.det() == DetId::Muon && id.subdetId() == MuonSubdetId::DT) {
00218 isDTorCSC = true;
00219
00220 DTChamberId DTid(id);
00221
00222 std::vector<DTRecSegment4D> allDTsegs;
00223 std::map<int, std::vector<DTRecSegment4D> >::const_iterator dtIter = getSegs.getDTlist().find(DTid.station());
00224 if (dtIter != getSegs.getDTlist().end()){
00225 allDTsegs = dtIter->second;
00226 }
00227 std::vector<DTRecSegment4D>::size_type sz = allDTsegs.size();
00228 for (unsigned int iSeg=0; iSeg<sz; ++iSeg) {
00229
00230
00231 TrajectoryStateOnSurface tsosdt = propagator->propagate(currentState,
00232 theG->idToDet(allDTsegs[iSeg].chamberId())->surface());
00233 if (!tsosdt.isValid()) continue;
00234
00235 std::vector<DTRecSegment4D> DTsegs;
00236 DTsegs.push_back(allDTsegs[iSeg]);
00237
00238 DTRecSegment4D bestDTSeg;
00239 double bestChamberValue = getBest(DTsegs, tsosdt, bestDTSeg);
00240 if (bestChamberValue < bestLayerValue) bestLayerValue = bestChamberValue;
00241
00242
00243 if (bestChamberValue >= DTThr || bestChamberValue > bestLayerValue) continue;
00244 layerRH.clear(); layerSEG.clear();
00245 layerSEG.push_back(theMuonRecHitBuilder->build(&bestDTSeg));
00246 std::vector<DTRecHit1D> DTrh = getSegs.getDTRHmap(bestDTSeg);
00247 for (std::vector<DTRecHit1D>::iterator it = DTrh.begin(); it != DTrh.end(); it++) {
00248 layerRH.push_back(theMuonRecHitBuilder->build(&*it));
00249 }
00250 }
00251 }
00252 if (id.det() == DetId::Muon && id.subdetId() == MuonSubdetId::CSC) {
00253 isDTorCSC = true;
00254 CSCDetId CSCid(id);
00255
00256 std::vector<CSCSegment> allCSCsegs;
00257 std::map<int, std::vector<CSCSegment> >::const_iterator cscIter = getSegs.getCSClist().find(CSCid.station());
00258 if (cscIter != getSegs.getCSClist().end()){
00259 allCSCsegs = cscIter->second;
00260 }
00261
00262 std::vector<CSCSegment>::size_type sz = allCSCsegs.size();
00263 for (unsigned int iSeg=0; iSeg<sz; ++iSeg) {
00264
00265
00266 TrajectoryStateOnSurface tsoscsc = propagator->propagate(currentState,
00267 theG->idToDet(allCSCsegs[iSeg].cscDetId())->surface());
00268 if (!tsoscsc.isValid()) continue;
00269
00270 std::vector<CSCSegment> CSCsegs;
00271 CSCsegs.push_back(allCSCsegs[iSeg]);
00272
00273 CSCSegment bestCSCSeg;
00274 double bestChamberValue = getBest(CSCsegs, tsoscsc, bestCSCSeg);
00275 if (bestChamberValue < bestLayerValue) bestLayerValue = bestChamberValue;
00276
00277
00278 if (bestChamberValue >= CSCThr || bestChamberValue > bestLayerValue) continue;
00279 layerRH.clear(); layerSEG.clear();
00280 layerSEG.push_back(theMuonRecHitBuilder->build(&bestCSCSeg));
00281 std::vector<CSCRecHit2D> CSCrh = getSegs.getCSCRHmap(bestCSCSeg);
00282 for (std::vector<CSCRecHit2D>::iterator it = CSCrh.begin(); it != CSCrh.end(); ++it) {
00283 layerRH.push_back(theMuonRecHitBuilder->build(&*it));
00284 }
00285 }
00286 }
00287 }
00288
00289 if (isDTorCSC) {
00290 if (bestLayerValue < MAX_THR) estimators.push_back(bestLayerValue);
00291 else estimators.push_back(-1);
00292 }
00293
00294
00295 if (layerSEG.size() > 0) {
00296 result.push_back(layerSEG.front());
00297 DetId id_ = layerSEG.front()->geographicalId();
00298 if (id_.subdetId() == MuonSubdetId::DT) {
00299 DTChamberId MuonChId(id_);
00300 TrajectoryStateOnSurface temp = propagator->propagate(currentState, theG->idToDet(MuonChId)->surface());
00301 if (temp.isValid())
00302 currentState = updatorHandle->update(temp, *layerSEG.front());
00303 }
00304 if (id_.subdetId() == MuonSubdetId::CSC) {
00305 CSCDetId MuonChId(id_);
00306 TrajectoryStateOnSurface temp = propagator->propagate(currentState, theG->idToDet(MuonChId)->surface());
00307 if (temp.isValid())
00308 currentState = updatorHandle->update(temp, *layerSEG.front());
00309 }
00310 }
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330 layerSEG.clear();
00331 layerRH.clear();
00332 }
00333 }
00334
00335
00336
00337 void DynamicTruncation::update(TrajectoryStateOnSurface& tsos, ConstRecHitPointer rechit) {
00338 TrajectoryStateOnSurface temp = updatorHandle->update(tsos, *rechit);
00339 if (temp.isValid()) currentState = updatorHandle->update(tsos, *rechit);
00340 }
00341
00342
00343
00344 void DynamicTruncation::updateWithDThits(ConstRecHitContainer& recHits) {
00345 for (ConstRecHitContainer::const_iterator it = recHits.begin(); it != recHits.end(); ++it) {
00346 DTLayerId layid((*it)->det()->geographicalId());
00347 TrajectoryStateOnSurface temp = propagator->propagate(currentState, theG->idToDet(layid)->surface());
00348 if (temp.isValid()) currentState = updatorHandle->update(temp, **it);
00349 }
00350 }
00351
00352
00353
00354 void DynamicTruncation::updateWithCSChits(ConstRecHitContainer& recHits) {
00355 for (ConstRecHitContainer::const_iterator it = recHits.begin(); it != recHits.end(); ++it) {
00356 const CSCLayer* cscChamber = cscGeom->layer((*it)->det()->geographicalId());
00357 CSCDetId layid = cscChamber->geographicalId();
00358 TrajectoryStateOnSurface temp = propagator->propagate(currentState, theG->idToDet(layid)->surface());
00359 if (temp.isValid()) currentState = updatorHandle->update(temp, **it);
00360 }
00361 }
00362
00363
00364
00365 TransientTrackingRecHit::ConstRecHitContainer DynamicTruncation::sort(ConstRecHitContainer& recHits) {
00366 unsigned int i=0;
00367 unsigned int j=0;
00368 ConstRecHitContainer::size_type n = recHits.size();
00369 for(i=1; i<n; ++i)
00370 for(j=n-1; j>=i; --j)
00371 {
00372 if(recHits[j-1]->globalPosition().mag() > recHits[j]->globalPosition().mag())
00373 {
00374 swap (recHits[j-1],recHits[j]);
00375 }
00376 }
00377 return recHits;
00378 }
00379