CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoMuon/GlobalTrackingTools/src/DynamicTruncation.cc

Go to the documentation of this file.
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   // Fill the map with the APE chamber by chamber (DT)
00077   // IT ALSO WORKS IF APE IS 0
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   // Fill the map with the APE chamber by chamber (CSC)
00089   // IT ALSO WORKS IF APE IS 0
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   // Put the tracker hits in the final vector and get the last tracker valid measure
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   // get the last (forward) predicted state for the tracker
00115   currentState = lastTKm.forwardPredictedState();
00116   
00117   // update the state with the last tracker measure
00118   update(currentState, lastTKm.recHit());
00119 
00120   // use the navigation to get the list of compatible dets
00121   map<int, std::vector<DetId> > detMap;
00122   compatibleDets(currentState, detMap);
00123   
00124   // selects the muon hits for the final refit
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; // DT thr 
00136   else DTThr = MAX_THR;
00137   if (eThr <= MAX_THR && eThr >= 0) CSCThr = eThr; // CSC thr
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; //default constructor is all zeroes, OK
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; //default constructor is all zeroes, OK
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   //  SteppingHelixPropagator prop(magfield.product(), anyDirection);
00183   //  MuonPatternRecoDumper dumper;
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     // Skip RPC layers
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     //    cout << comps.size() << " compatible Dets with " << navLayers[ilayer]->subDetector() << " Layer " << ilayer << " "
00194     //     << dumper.dumpLayer(navLayers[ilayer]) << " " << endl;
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           // Propagate the state to the current chamber
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           // Check if the best estimator value is below the THR and then get the RH componing the segment
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           // Propagate the state to the current chamber
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           // Check if the best estimator value is below the THR and then get the RH componing the segment
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     // For segment based fits
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       // For hit based fits 
00317     if (layerRH.size() > 0) {                                                                                                                                          
00318       for (ConstRecHitContainer::iterator it = layerRH.begin(); it != layerRH.end(); ++it) result.push_back((*it)); 
00319 
00320       // sort the vector                                                                                                                                             
00321       layerRH = sort(layerRH);                                                                                                                                         
00322       
00323       // update the currentState using all rh                                                                                                             
00324       DetId id = layerRH.front()->geographicalId();                                                                                                      
00325       if (id.subdetId() == MuonSubdetId::DT) updateWithDThits(layerRH);                                                                                 
00326       else updateWithCSChits(layerRH);                                                                                                                 
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