00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include <memory>
00015
00016
00017 #include "FWCore/Framework/interface/Frameworkfwd.h"
00018 #include "FWCore/Framework/interface/EDProducer.h"
00019
00020 #include "FWCore/Framework/interface/Event.h"
00021 #include "FWCore/Framework/interface/EventSetup.h"
00022
00023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00024
00025 #include "DataFormats/Common/interface/Handle.h"
00026 #include "DataFormats/TrackReco/interface/Track.h"
00027 #include "DataFormats/MuonReco/interface/Muon.h"
00028 #include "DataFormats/MuonReco/interface/CaloMuon.h"
00029 #include "DataFormats/MuonReco/interface/MuonTime.h"
00030 #include "DataFormats/MuonReco/interface/MuonTimeExtra.h"
00031 #include "DataFormats/MuonReco/interface/MuonTimeExtraMap.h"
00032 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
00033 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
00034
00035 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
00036 #include "Utilities/Timing/interface/TimerStack.h"
00037
00038 #include <boost/regex.hpp>
00039 #include "RecoMuon/MuonIdentification/plugins/MuonIdProducer.h"
00040 #include "RecoMuon/MuonIdentification/interface/MuonIdTruthInfo.h"
00041 #include "RecoMuon/MuonIdentification/interface/MuonArbitrationMethods.h"
00042 #include "RecoMuon/MuonIdentification/interface/MuonMesh.h"
00043
00044 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositExtractorFactory.h"
00045 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00046
00047 #include <algorithm>
00048
00049 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00050 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
00051 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
00052 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
00053
00054 #include "DataFormats/Common/interface/ValueMap.h"
00055
00056 #include "RecoMuon/MuonIdentification/interface/MuonKinkFinder.h"
00057
00058 MuonIdProducer::MuonIdProducer(const edm::ParameterSet& iConfig):
00059 muIsoExtractorCalo_(0),muIsoExtractorTrack_(0),muIsoExtractorJet_(0)
00060 {
00061 produces<reco::MuonCollection>();
00062 produces<reco::CaloMuonCollection>();
00063 produces<reco::MuonTimeExtraMap>("combined");
00064 produces<reco::MuonTimeExtraMap>("dt");
00065 produces<reco::MuonTimeExtraMap>("csc");
00066
00067 minPt_ = iConfig.getParameter<double>("minPt");
00068 minP_ = iConfig.getParameter<double>("minP");
00069 minPCaloMuon_ = iConfig.getParameter<double>("minPCaloMuon");
00070 minNumberOfMatches_ = iConfig.getParameter<int>("minNumberOfMatches");
00071 addExtraSoftMuons_ = iConfig.getParameter<bool>("addExtraSoftMuons");
00072 maxAbsEta_ = iConfig.getParameter<double>("maxAbsEta");
00073 maxAbsDx_ = iConfig.getParameter<double>("maxAbsDx");
00074 maxAbsPullX_ = iConfig.getParameter<double>("maxAbsPullX");
00075 maxAbsDy_ = iConfig.getParameter<double>("maxAbsDy");
00076 maxAbsPullY_ = iConfig.getParameter<double>("maxAbsPullY");
00077 fillCaloCompatibility_ = iConfig.getParameter<bool>("fillCaloCompatibility");
00078 fillEnergy_ = iConfig.getParameter<bool>("fillEnergy");
00079 fillMatching_ = iConfig.getParameter<bool>("fillMatching");
00080 fillIsolation_ = iConfig.getParameter<bool>("fillIsolation");
00081 writeIsoDeposits_ = iConfig.getParameter<bool>("writeIsoDeposits");
00082 fillGlobalTrackQuality_ = iConfig.getParameter<bool>("fillGlobalTrackQuality");
00083 fillGlobalTrackRefits_ = iConfig.getParameter<bool>("fillGlobalTrackRefits");
00084
00085 fillTrackerKink_ = false;
00086 if (fillGlobalTrackQuality_) fillTrackerKink_ = iConfig.getParameter<bool>("fillTrackerKink");
00087
00088 ptThresholdToFillCandidateP4WithGlobalFit_ = iConfig.getParameter<double>("ptThresholdToFillCandidateP4WithGlobalFit");
00089 sigmaThresholdToFillCandidateP4WithGlobalFit_ = iConfig.getParameter<double>("sigmaThresholdToFillCandidateP4WithGlobalFit");
00090 caloCut_ = iConfig.getParameter<double>("minCaloCompatibility");
00091 arbClean_ = iConfig.getParameter<bool>("runArbitrationCleaner");
00092
00093
00094 edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
00095 parameters_.loadParameters( parameters );
00096
00097
00098 edm::ParameterSet timingParameters = iConfig.getParameter<edm::ParameterSet>("TimingFillerParameters");
00099 theTimingFiller_ = new MuonTimingFiller(timingParameters);
00100
00101 if (fillCaloCompatibility_){
00102
00103 parameters = iConfig.getParameter<edm::ParameterSet>("MuonCaloCompatibility");
00104 muonCaloCompatibility_.configure( parameters );
00105 }
00106
00107 if (fillIsolation_){
00108
00109 edm::ParameterSet caloExtractorPSet = iConfig.getParameter<edm::ParameterSet>("CaloExtractorPSet");
00110 std::string caloExtractorName = caloExtractorPSet.getParameter<std::string>("ComponentName");
00111 muIsoExtractorCalo_ = IsoDepositExtractorFactory::get()->create( caloExtractorName, caloExtractorPSet);
00112
00113 edm::ParameterSet trackExtractorPSet = iConfig.getParameter<edm::ParameterSet>("TrackExtractorPSet");
00114 std::string trackExtractorName = trackExtractorPSet.getParameter<std::string>("ComponentName");
00115 muIsoExtractorTrack_ = IsoDepositExtractorFactory::get()->create( trackExtractorName, trackExtractorPSet);
00116
00117 edm::ParameterSet jetExtractorPSet = iConfig.getParameter<edm::ParameterSet>("JetExtractorPSet");
00118 std::string jetExtractorName = jetExtractorPSet.getParameter<std::string>("ComponentName");
00119 muIsoExtractorJet_ = IsoDepositExtractorFactory::get()->create( jetExtractorName, jetExtractorPSet);
00120 }
00121 if (fillIsolation_ && writeIsoDeposits_){
00122 trackDepositName_ = iConfig.getParameter<std::string>("trackDepositName");
00123 produces<reco::IsoDepositMap>(trackDepositName_);
00124 ecalDepositName_ = iConfig.getParameter<std::string>("ecalDepositName");
00125 produces<reco::IsoDepositMap>(ecalDepositName_);
00126 hcalDepositName_ = iConfig.getParameter<std::string>("hcalDepositName");
00127 produces<reco::IsoDepositMap>(hcalDepositName_);
00128 hoDepositName_ = iConfig.getParameter<std::string>("hoDepositName");
00129 produces<reco::IsoDepositMap>(hoDepositName_);
00130 jetDepositName_ = iConfig.getParameter<std::string>("jetDepositName");
00131 produces<reco::IsoDepositMap>(jetDepositName_);
00132 }
00133
00134 inputCollectionLabels_ = iConfig.getParameter<std::vector<edm::InputTag> >("inputCollectionLabels");
00135 inputCollectionTypes_ = iConfig.getParameter<std::vector<std::string> >("inputCollectionTypes");
00136 if (inputCollectionLabels_.size() != inputCollectionTypes_.size())
00137 throw cms::Exception("ConfigurationError") << "Number of input collection labels is different from number of types. " <<
00138 "For each collection label there should be exactly one collection type specified.";
00139 if (inputCollectionLabels_.size()>7 ||inputCollectionLabels_.empty())
00140 throw cms::Exception("ConfigurationError") << "Number of input collections should be from 1 to 7.";
00141
00142 debugWithTruthMatching_ = iConfig.getParameter<bool>("debugWithTruthMatching");
00143 if (debugWithTruthMatching_) edm::LogWarning("MuonIdentification")
00144 << "========================================================================\n"
00145 << "Debugging mode with truth matching is turned on!!! Make sure you understand what you are doing!\n"
00146 << "========================================================================\n";
00147 if (fillGlobalTrackQuality_){
00148 globalTrackQualityInputTag_ = iConfig.getParameter<edm::InputTag>("globalTrackQualityInputTag");
00149 }
00150
00151 if (fillTrackerKink_) {
00152 trackerKinkFinder_.reset(new MuonKinkFinder(iConfig.getParameter<edm::ParameterSet>("TrackerKinkFinderParameters")));
00153 }
00154
00155
00156 meshAlgo_ = new MuonMesh(iConfig.getParameter<edm::ParameterSet>("arbitrationCleanerOptions"));
00157 }
00158
00159
00160 MuonIdProducer::~MuonIdProducer()
00161 {
00162 if (muIsoExtractorCalo_) delete muIsoExtractorCalo_;
00163 if (muIsoExtractorTrack_) delete muIsoExtractorTrack_;
00164 if (muIsoExtractorJet_) delete muIsoExtractorJet_;
00165 if (theTimingFiller_) delete theTimingFiller_;
00166 if (meshAlgo_) delete meshAlgo_;
00167
00168 }
00169
00170 void MuonIdProducer::init(edm::Event& iEvent, const edm::EventSetup& iSetup)
00171 {
00172
00173
00174
00175 innerTrackCollectionHandle_.clear();
00176 outerTrackCollectionHandle_.clear();
00177 linkCollectionHandle_.clear();
00178 muonCollectionHandle_.clear();
00179
00180 tpfmsCollectionHandle_.clear();
00181 pickyCollectionHandle_.clear();
00182 dytCollectionHandle_.clear();
00183
00184
00185
00186 edm::ESHandle<Propagator> propagator;
00187 iSetup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagator);
00188 trackAssociator_.setPropagator(propagator.product());
00189
00190 if (fillTrackerKink_) trackerKinkFinder_->init(iSetup);
00191
00192
00193 for ( unsigned int i = 0; i < inputCollectionLabels_.size(); ++i ) {
00194 if ( inputCollectionTypes_[i] == "inner tracks" ) {
00195 iEvent.getByLabel(inputCollectionLabels_[i], innerTrackCollectionHandle_);
00196 if (! innerTrackCollectionHandle_.isValid())
00197 throw cms::Exception("FatalError") << "Failed to get input track collection with label: " << inputCollectionLabels_[i];
00198 LogTrace("MuonIdentification") << "Number of input inner tracks: " << innerTrackCollectionHandle_->size();
00199 continue;
00200 }
00201 if ( inputCollectionTypes_[i] == "outer tracks" ) {
00202 iEvent.getByLabel(inputCollectionLabels_[i], outerTrackCollectionHandle_);
00203 if (! outerTrackCollectionHandle_.isValid())
00204 throw cms::Exception("FatalError") << "Failed to get input track collection with label: " << inputCollectionLabels_[i];
00205 LogTrace("MuonIdentification") << "Number of input outer tracks: " << outerTrackCollectionHandle_->size();
00206 continue;
00207 }
00208 if ( inputCollectionTypes_[i] == "links" ) {
00209 iEvent.getByLabel(inputCollectionLabels_[i], linkCollectionHandle_);
00210 if (! linkCollectionHandle_.isValid())
00211 throw cms::Exception("FatalError") << "Failed to get input link collection with label: " << inputCollectionLabels_[i];
00212 LogTrace("MuonIdentification") << "Number of input links: " << linkCollectionHandle_->size();
00213 continue;
00214 }
00215 if ( inputCollectionTypes_[i] == "muons" ) {
00216 iEvent.getByLabel(inputCollectionLabels_[i], muonCollectionHandle_);
00217 if (! muonCollectionHandle_.isValid())
00218 throw cms::Exception("FatalError") << "Failed to get input muon collection with label: " << inputCollectionLabels_[i];
00219 LogTrace("MuonIdentification") << "Number of input muons: " << muonCollectionHandle_->size();
00220 continue;
00221 }
00222 if ( fillGlobalTrackRefits_ && inputCollectionTypes_[i] == "tev firstHit" ) {
00223 iEvent.getByLabel(inputCollectionLabels_[i], tpfmsCollectionHandle_);
00224 if (! tpfmsCollectionHandle_.isValid())
00225 throw cms::Exception("FatalError") << "Failed to get input muon collection with label: " << inputCollectionLabels_[i];
00226 LogTrace("MuonIdentification") << "Number of input muons: " << tpfmsCollectionHandle_->size();
00227 continue;
00228 }
00229
00230 if ( fillGlobalTrackRefits_ && inputCollectionTypes_[i] == "tev picky" ) {
00231 iEvent.getByLabel(inputCollectionLabels_[i], pickyCollectionHandle_);
00232 if (! pickyCollectionHandle_.isValid())
00233 throw cms::Exception("FatalError") << "Failed to get input muon collection with label: " << inputCollectionLabels_[i];
00234 LogTrace("MuonIdentification") << "Number of input muons: " << pickyCollectionHandle_->size();
00235 continue;
00236 }
00237
00238 if ( fillGlobalTrackRefits_ && inputCollectionTypes_[i] == "tev dyt" ) {
00239 iEvent.getByLabel(inputCollectionLabels_[i], dytCollectionHandle_);
00240 if (! dytCollectionHandle_.isValid())
00241 throw cms::Exception("FatalError") << "Failed to get input muon collection with label: " << inputCollectionLabels_[i];
00242 LogTrace("MuonIdentification") << "Number of input muons: " << dytCollectionHandle_->size();
00243 continue;
00244 }
00245 throw cms::Exception("FatalError") << "Unknown input collection type: " << inputCollectionTypes_[i];
00246 }
00247 }
00248
00249 reco::Muon MuonIdProducer::makeMuon(edm::Event& iEvent, const edm::EventSetup& iSetup,
00250 const reco::TrackRef& track, MuonIdProducer::TrackType type)
00251 {
00252 LogTrace("MuonIdentification") << "Creating a muon from a track " << track.get()->pt() <<
00253 " Pt (GeV), eta: " << track.get()->eta();
00254 reco::Muon aMuon( makeMuon( *(track.get()) ) );
00255
00256 aMuon.setMuonTrack(type,track);
00257 return aMuon;
00258 }
00259
00260
00261 reco::CaloMuon MuonIdProducer::makeCaloMuon( const reco::Muon& muon )
00262 {
00263 reco::CaloMuon aMuon;
00264 aMuon.setInnerTrack( muon.innerTrack() );
00265
00266 if (muon.isEnergyValid()) aMuon.setCalEnergy( muon.calEnergy() );
00267
00268 if (fillCaloCompatibility_) aMuon.setCaloCompatibility( muonCaloCompatibility_.evaluate(muon) );
00269 return aMuon;
00270 }
00271
00272
00273 reco::Muon MuonIdProducer::makeMuon( const reco::MuonTrackLinks& links )
00274 {
00275 LogTrace("MuonIdentification") << "Creating a muon from a link to tracks object";
00276 reco::Muon aMuon;
00277 if ( links.trackerTrack()->pt() > ptThresholdToFillCandidateP4WithGlobalFit_ &&
00278 links.globalTrack()->pt() > ptThresholdToFillCandidateP4WithGlobalFit_ &&
00279 ( fabs(links.trackerTrack()->qoverp()-links.globalTrack()->qoverp()) <
00280 sigmaThresholdToFillCandidateP4WithGlobalFit_ * links.trackerTrack()->qoverpError() ) )
00281 aMuon = makeMuon( *(links.globalTrack()) );
00282 else
00283 aMuon = makeMuon( *(links.trackerTrack()) );
00284
00285 aMuon.setInnerTrack( links.trackerTrack() );
00286 aMuon.setOuterTrack( links.standAloneTrack() );
00287 aMuon.setGlobalTrack( links.globalTrack() );
00288
00289 if(fillGlobalTrackRefits_){
00290 if (tpfmsCollectionHandle_.isValid() && !tpfmsCollectionHandle_.failedToGet()) {
00291 reco::TrackToTrackMap::const_iterator it = tpfmsCollectionHandle_->find(links.globalTrack());
00292 if (it != tpfmsCollectionHandle_->end()) aMuon.setMuonTrack(reco::Muon::TPFMS, (it->val));
00293 }
00294 if (pickyCollectionHandle_.isValid() && !pickyCollectionHandle_.failedToGet()) {
00295 reco::TrackToTrackMap::const_iterator it = pickyCollectionHandle_->find(links.globalTrack());
00296 if (it != pickyCollectionHandle_->end()) aMuon.setMuonTrack(reco::Muon::Picky, (it->val));
00297 }
00298 if (dytCollectionHandle_.isValid() && !dytCollectionHandle_.failedToGet()) {
00299 reco::TrackToTrackMap::const_iterator it = dytCollectionHandle_->find(links.globalTrack());
00300 if (it != dytCollectionHandle_->end()) aMuon.setMuonTrack(reco::Muon::DYT, (it->val));
00301 }
00302 }
00303 return aMuon;
00304 }
00305
00306
00307 bool MuonIdProducer::isGoodTrack( const reco::Track& track )
00308 {
00309
00310 if (track.pt() < minPt_ || (track.p() < minP_ && track.p() < minPCaloMuon_)){
00311 LogTrace("MuonIdentification") << "Skipped low momentum track (Pt,P): " << track.pt() <<
00312 ", " << track.p() << " GeV";
00313 return false;
00314 }
00315
00316
00317 if ( fabs(track.eta()) > maxAbsEta_ ){
00318 LogTrace("MuonIdentification") << "Skipped track with large pseudo rapidity (Eta: " << track.eta() << " )";
00319 return false;
00320 }
00321 return true;
00322 }
00323
00324 unsigned int MuonIdProducer::chamberId( const DetId& id )
00325 {
00326 switch ( id.det() ) {
00327 case DetId::Muon:
00328 switch ( id.subdetId() ) {
00329 case MuonSubdetId::DT:
00330 {
00331 DTChamberId detId(id.rawId());
00332 return detId.rawId();
00333 }
00334 break;
00335 case MuonSubdetId::CSC:
00336 {
00337 CSCDetId detId(id.rawId());
00338 return detId.chamberId().rawId();
00339 }
00340 break;
00341 default:
00342 return 0;
00343 }
00344 default:
00345 return 0;
00346 }
00347 return 0;
00348 }
00349
00350
00351 int MuonIdProducer::overlap(const reco::Muon& muon, const reco::Track& track)
00352 {
00353 int numberOfCommonDetIds = 0;
00354 if ( ! muon.isMatchesValid() ||
00355 track.extra().isNull() ||
00356 track.extra()->recHits().isNull() ) return numberOfCommonDetIds;
00357 const std::vector<reco::MuonChamberMatch>& matches( muon.matches() );
00358 for ( std::vector<reco::MuonChamberMatch>::const_iterator match = matches.begin();
00359 match != matches.end(); ++match )
00360 {
00361 if ( match->segmentMatches.empty() ) continue;
00362 bool foundCommonDetId = false;
00363
00364 for ( TrackingRecHitRefVector::const_iterator hit = track.extra()->recHitsBegin();
00365 hit != track.extra()->recHitsEnd(); ++hit )
00366 {
00367
00368
00369
00370
00371 if ( chamberId(hit->get()->geographicalId()) == match->id.rawId() ) {
00372 foundCommonDetId = true;
00373 break;
00374 }
00375 }
00376 if ( foundCommonDetId ) {
00377 numberOfCommonDetIds++;
00378 break;
00379 }
00380 }
00381 return numberOfCommonDetIds;
00382 }
00383
00384 void MuonIdProducer::beginRun(edm::Run& iRun, const edm::EventSetup& iSetup)
00385 {
00386 edm::ESHandle<CSCGeometry> geomHandle;
00387 iSetup.get<MuonGeometryRecord>().get(geomHandle);
00388
00389 meshAlgo_->setCSCGeometry(geomHandle.product());
00390
00391 }
00392
00393 bool validateGlobalMuonPair( const reco::MuonTrackLinks& goodMuon,
00394 const reco::MuonTrackLinks& badMuon )
00395 {
00396 if ( std::min(goodMuon.globalTrack()->hitPattern().numberOfValidMuonHits(),
00397 badMuon.globalTrack()->hitPattern().numberOfValidMuonHits()) > 10 ){
00398 if ( goodMuon.globalTrack()->normalizedChi2() >
00399 badMuon.globalTrack()->normalizedChi2() )
00400 return false;
00401 else
00402 return true;
00403 }
00404 if ( goodMuon.globalTrack()->hitPattern().numberOfValidMuonHits() <
00405 badMuon.globalTrack()->hitPattern().numberOfValidMuonHits() ) return false;
00406 return true;
00407 }
00408
00409 void MuonIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00410 {
00411
00412
00413
00414 std::auto_ptr<reco::MuonCollection> outputMuons(new reco::MuonCollection);
00415 std::auto_ptr<reco::CaloMuonCollection> caloMuons( new reco::CaloMuonCollection );
00416
00417 init(iEvent, iSetup);
00418
00419 std::auto_ptr<reco::MuonTimeExtraMap> muonTimeMap(new reco::MuonTimeExtraMap());
00420 reco::MuonTimeExtraMap::Filler filler(*muonTimeMap);
00421 std::auto_ptr<reco::MuonTimeExtraMap> muonTimeMapDT(new reco::MuonTimeExtraMap());
00422 reco::MuonTimeExtraMap::Filler fillerDT(*muonTimeMapDT);
00423 std::auto_ptr<reco::MuonTimeExtraMap> muonTimeMapCSC(new reco::MuonTimeExtraMap());
00424 reco::MuonTimeExtraMap::Filler fillerCSC(*muonTimeMapCSC);
00425
00426 std::auto_ptr<reco::IsoDepositMap> trackDepMap(new reco::IsoDepositMap());
00427 reco::IsoDepositMap::Filler trackDepFiller(*trackDepMap);
00428 std::auto_ptr<reco::IsoDepositMap> ecalDepMap(new reco::IsoDepositMap());
00429 reco::IsoDepositMap::Filler ecalDepFiller(*ecalDepMap);
00430 std::auto_ptr<reco::IsoDepositMap> hcalDepMap(new reco::IsoDepositMap());
00431 reco::IsoDepositMap::Filler hcalDepFiller(*hcalDepMap);
00432 std::auto_ptr<reco::IsoDepositMap> hoDepMap(new reco::IsoDepositMap());
00433 reco::IsoDepositMap::Filler hoDepFiller(*hoDepMap);
00434 std::auto_ptr<reco::IsoDepositMap> jetDepMap(new reco::IsoDepositMap());
00435 reco::IsoDepositMap::Filler jetDepFiller(*jetDepMap);
00436
00437
00438
00439
00440 if ( muonCollectionHandle_.isValid() )
00441 for ( reco::MuonCollection::const_iterator muon = muonCollectionHandle_->begin();
00442 muon != muonCollectionHandle_->end(); ++muon )
00443 outputMuons->push_back(*muon);
00444
00445
00446 if ( linkCollectionHandle_.isValid() ){
00447 std::vector<bool> goodmuons(linkCollectionHandle_->size(),true);
00448 if ( goodmuons.size()>1 ){
00449
00450 for ( unsigned int i=0; i<linkCollectionHandle_->size()-1; ++i ){
00451 if (!checkLinks(&linkCollectionHandle_->at(i))) continue;
00452 for ( unsigned int j=i+1; j<linkCollectionHandle_->size(); ++j ){
00453 if (!checkLinks(&linkCollectionHandle_->at(j))) continue;
00454 if ( linkCollectionHandle_->at(i).trackerTrack().isNonnull() &&
00455 linkCollectionHandle_->at(i).trackerTrack() ==
00456 linkCollectionHandle_->at(j).trackerTrack() )
00457 {
00458
00459
00460
00461
00462 if ( validateGlobalMuonPair(linkCollectionHandle_->at(i),linkCollectionHandle_->at(j)) )
00463 goodmuons[j] = false;
00464 else
00465 goodmuons[i] = false;
00466 }
00467 }
00468 }
00469
00470 for ( unsigned int i=0; i<linkCollectionHandle_->size()-1; ++i ){
00471 if ( !goodmuons[i] ) continue;
00472 if (!checkLinks(&linkCollectionHandle_->at(i))) continue;
00473 for ( unsigned int j=i+1; j<linkCollectionHandle_->size(); ++j ){
00474 if ( !goodmuons[j] ) continue;
00475 if (!checkLinks(&linkCollectionHandle_->at(j))) continue;
00476 if ( linkCollectionHandle_->at(i).standAloneTrack().isNonnull() &&
00477 linkCollectionHandle_->at(i).standAloneTrack() ==
00478 linkCollectionHandle_->at(j).standAloneTrack() )
00479 {
00480 if ( validateGlobalMuonPair(linkCollectionHandle_->at(i),linkCollectionHandle_->at(j)) )
00481 goodmuons[j] = false;
00482 else
00483 goodmuons[i] = false;
00484 }
00485 }
00486 }
00487 }
00488 for ( unsigned int i=0; i<linkCollectionHandle_->size(); ++i ){
00489 if ( !goodmuons[i] ) continue;
00490 const reco::MuonTrackLinks* links = &linkCollectionHandle_->at(i);
00491 if ( ! checkLinks(links)) continue;
00492
00493 bool newMuon = true;
00494 for ( reco::MuonCollection::const_iterator muon = outputMuons->begin();
00495 muon != outputMuons->end(); ++muon )
00496 if ( muon->track() == links->trackerTrack() &&
00497 muon->standAloneMuon() == links->standAloneTrack() &&
00498 muon->combinedMuon() == links->globalTrack() )
00499 newMuon = false;
00500 if ( newMuon ) {
00501 outputMuons->push_back( makeMuon( *links ) );
00502 outputMuons->back().setType(reco::Muon::GlobalMuon | reco::Muon::StandAloneMuon);
00503 }
00504 }
00505 }
00506
00507
00508 if ( innerTrackCollectionHandle_.isValid() ) {
00509 LogTrace("MuonIdentification") << "Creating tracker muons";
00510 for ( unsigned int i = 0; i < innerTrackCollectionHandle_->size(); ++i )
00511 {
00512 const reco::Track& track = innerTrackCollectionHandle_->at(i);
00513 if ( ! isGoodTrack( track ) ) continue;
00514 bool splitTrack = false;
00515 if ( track.extra().isAvailable() &&
00516 TrackDetectorAssociator::crossedIP( track ) ) splitTrack = true;
00517 std::vector<TrackDetectorAssociator::Direction> directions;
00518 if ( splitTrack ) {
00519 directions.push_back(TrackDetectorAssociator::InsideOut);
00520 directions.push_back(TrackDetectorAssociator::OutsideIn);
00521 } else {
00522 directions.push_back(TrackDetectorAssociator::Any);
00523 }
00524 for ( std::vector<TrackDetectorAssociator::Direction>::const_iterator direction = directions.begin();
00525 direction != directions.end(); ++direction )
00526 {
00527
00528
00529 reco::Muon trackerMuon( makeMuon(iEvent, iSetup, reco::TrackRef( innerTrackCollectionHandle_, i ), reco::Muon::InnerTrack ) );
00530 trackerMuon.setType( reco::Muon::TrackerMuon );
00531 fillMuonId(iEvent, iSetup, trackerMuon, *direction);
00532
00533
00534 if ( debugWithTruthMatching_ ) {
00535
00536
00537 MuonIdTruthInfo::truthMatchMuon(iEvent, iSetup, trackerMuon);
00538 }
00539
00540
00541
00542
00543 bool newMuon = true;
00544 bool goodTrackerMuon = isGoodTrackerMuon( trackerMuon );
00545 for ( reco::MuonCollection::iterator muon = outputMuons->begin();
00546 muon != outputMuons->end(); ++muon )
00547 {
00548 if ( muon->innerTrack().get() == trackerMuon.innerTrack().get() &&
00549 cos(phiOfMuonIneteractionRegion(*muon) -
00550 phiOfMuonIneteractionRegion(trackerMuon)) > 0 )
00551 {
00552 newMuon = false;
00553 muon->setMatches( trackerMuon.matches() );
00554 if (trackerMuon.isTimeValid()) muon->setTime( trackerMuon.time() );
00555 if (trackerMuon.isEnergyValid()) muon->setCalEnergy( trackerMuon.calEnergy() );
00556 if (goodTrackerMuon) muon->setType( muon->type() | reco::Muon::TrackerMuon );
00557 LogTrace("MuonIdentification") << "Found a corresponding global muon. Set energy, matches and move on";
00558 break;
00559 }
00560 }
00561 if ( newMuon ) {
00562 if ( goodTrackerMuon ){
00563 outputMuons->push_back( trackerMuon );
00564 } else {
00565 LogTrace("MuonIdentification") << "track failed minimal number of muon matches requirement";
00566 const reco::CaloMuon& caloMuon = makeCaloMuon(trackerMuon);
00567 if ( ! caloMuon.isCaloCompatibilityValid() || caloMuon.caloCompatibility() < caloCut_ || caloMuon.p() < minPCaloMuon_) continue;
00568 caloMuons->push_back( caloMuon );
00569 }
00570 }
00571 }
00572 }
00573 }
00574
00575
00576 if ( outerTrackCollectionHandle_.isValid() ) {
00577 LogTrace("MuonIdentification") << "Looking for new muons among stand alone muon tracks";
00578 for ( unsigned int i = 0; i < outerTrackCollectionHandle_->size(); ++i )
00579 {
00580
00581 bool newMuon = true;
00582 for ( reco::MuonCollection::iterator muon = outputMuons->begin();
00583 muon != outputMuons->end(); ++muon )
00584 {
00585 if ( ! muon->standAloneMuon().isNull() ) {
00586
00587 if ( muon->standAloneMuon().get() == &(outerTrackCollectionHandle_->at(i)) ) {
00588 newMuon = false;
00589 break;
00590 }
00591 } else {
00592
00593
00594
00595
00596
00597
00598
00599 if ( overlap(*muon,outerTrackCollectionHandle_->at(i))>0 ) {
00600 LogTrace("MuonIdentification") << "Found associated tracker muon. Set a reference and move on";
00601 newMuon = false;
00602 muon->setOuterTrack( reco::TrackRef( outerTrackCollectionHandle_, i ) );
00603 muon->setType( muon->type() | reco::Muon::StandAloneMuon );
00604 break;
00605 }
00606 }
00607 }
00608 if ( newMuon ) {
00609 LogTrace("MuonIdentification") << "No associated stand alone track is found. Making a muon";
00610 outputMuons->push_back( makeMuon(iEvent, iSetup,
00611 reco::TrackRef( outerTrackCollectionHandle_, i ), reco::Muon::OuterTrack ) );
00612 outputMuons->back().setType( reco::Muon::StandAloneMuon );
00613 }
00614 }
00615 }
00616
00617 LogTrace("MuonIdentification") << "Dress up muons if it's necessary";
00618
00619 int nMuons=outputMuons->size();
00620
00621 std::vector<reco::MuonTimeExtra> dtTimeColl(nMuons);
00622 std::vector<reco::MuonTimeExtra> cscTimeColl(nMuons);
00623 std::vector<reco::MuonTimeExtra> combinedTimeColl(nMuons);
00624 std::vector<reco::IsoDeposit> trackDepColl(nMuons);
00625 std::vector<reco::IsoDeposit> ecalDepColl(nMuons);
00626 std::vector<reco::IsoDeposit> hcalDepColl(nMuons);
00627 std::vector<reco::IsoDeposit> hoDepColl(nMuons);
00628 std::vector<reco::IsoDeposit> jetDepColl(nMuons);
00629
00630
00631 unsigned int i=0;
00632 for ( reco::MuonCollection::iterator muon = outputMuons->begin(); muon != outputMuons->end(); ++muon )
00633 {
00634
00635
00636 if ( ( fillMatching_ && ! muon->isMatchesValid() ) ||
00637 ( fillEnergy_ && !muon->isEnergyValid() ) )
00638 {
00639
00640
00641 if ( muon->isStandAloneMuon() ) {
00642 if ( cos(phiOfMuonIneteractionRegion(*muon) - muon->phi()) > 0 )
00643 fillMuonId(iEvent, iSetup, *muon, TrackDetectorAssociator::InsideOut);
00644 else
00645 fillMuonId(iEvent, iSetup, *muon, TrackDetectorAssociator::OutsideIn);
00646 } else {
00647 LogTrace("MuonIdentification") << "THIS SHOULD NEVER HAPPEN";
00648 fillMuonId(iEvent, iSetup, *muon);
00649 }
00650 }
00651
00652 if (fillGlobalTrackQuality_){
00653
00654 fillGlbQuality(iEvent, iSetup, *muon);
00655 }
00656 LogDebug("MuonIdentification");
00657
00658 if (fillTrackerKink_) {
00659 fillTrackerKink(*muon);
00660 }
00661
00662
00663 if ( fillCaloCompatibility_ ) muon->setCaloCompatibility( muonCaloCompatibility_.evaluate(*muon) );
00664
00665
00666
00667 if ( fillIsolation_ ) fillMuonIsolation(iEvent, iSetup, *muon,
00668 trackDepColl[i], ecalDepColl[i], hcalDepColl[i], hoDepColl[i], jetDepColl[i]);
00669
00670
00671
00672 reco::MuonTime muonTime;
00673 reco::MuonTimeExtra dtTime;
00674 reco::MuonTimeExtra cscTime;
00675 reco::MuonTimeExtra combinedTime;
00676
00677 theTimingFiller_->fillTiming(*muon, dtTime, cscTime, combinedTime, iEvent, iSetup);
00678
00679 muonTime.nDof=combinedTime.nDof();
00680 muonTime.timeAtIpInOut=combinedTime.timeAtIpInOut();
00681 muonTime.timeAtIpInOutErr=combinedTime.timeAtIpInOutErr();
00682 muonTime.timeAtIpOutIn=combinedTime.timeAtIpOutIn();
00683 muonTime.timeAtIpOutInErr=combinedTime.timeAtIpOutInErr();
00684
00685 muon->setTime( muonTime);
00686 dtTimeColl[i] = dtTime;
00687 cscTimeColl[i] = cscTime;
00688 combinedTimeColl[i] = combinedTime;
00689
00690 i++;
00691
00692 }
00693
00694 LogTrace("MuonIdentification") << "number of muons produced: " << outputMuons->size();
00695
00696 if ( fillMatching_ ) fillArbitrationInfo( outputMuons.get() );
00697
00698 edm::OrphanHandle<reco::MuonCollection> muonHandle = iEvent.put(outputMuons);
00699
00700 filler.insert(muonHandle, combinedTimeColl.begin(), combinedTimeColl.end());
00701 filler.fill();
00702 fillerDT.insert(muonHandle, dtTimeColl.begin(), dtTimeColl.end());
00703 fillerDT.fill();
00704 fillerCSC.insert(muonHandle, cscTimeColl.begin(), cscTimeColl.end());
00705 fillerCSC.fill();
00706
00707 iEvent.put(muonTimeMap,"combined");
00708 iEvent.put(muonTimeMapDT,"dt");
00709 iEvent.put(muonTimeMapCSC,"csc");
00710
00711 if (writeIsoDeposits_ && fillIsolation_){
00712 trackDepFiller.insert(muonHandle, trackDepColl.begin(), trackDepColl.end());
00713 trackDepFiller.fill();
00714 iEvent.put(trackDepMap, trackDepositName_);
00715 ecalDepFiller.insert(muonHandle, ecalDepColl.begin(), ecalDepColl.end());
00716 ecalDepFiller.fill();
00717 iEvent.put(ecalDepMap, ecalDepositName_);
00718 hcalDepFiller.insert(muonHandle, hcalDepColl.begin(), hcalDepColl.end());
00719 hcalDepFiller.fill();
00720 iEvent.put(hcalDepMap, hcalDepositName_);
00721 hoDepFiller.insert(muonHandle, hoDepColl.begin(), hoDepColl.end());
00722 hoDepFiller.fill();
00723 iEvent.put(hoDepMap, hoDepositName_);
00724 jetDepFiller.insert(muonHandle, jetDepColl.begin(), jetDepColl.end());
00725 jetDepFiller.fill();
00726 iEvent.put(jetDepMap, jetDepositName_);
00727 }
00728
00729 iEvent.put(caloMuons);
00730 }
00731
00732
00733 bool MuonIdProducer::isGoodTrackerMuon( const reco::Muon& muon )
00734 {
00735 if(muon.track()->pt() < minPt_ || muon.track()->p() < minP_) return false;
00736 if ( addExtraSoftMuons_ &&
00737 muon.pt()<5 && fabs(muon.eta())<1.5 &&
00738 muon.numberOfMatches( reco::Muon::NoArbitration ) >= 1 ) return true;
00739 return ( muon.numberOfMatches( reco::Muon::NoArbitration ) >= minNumberOfMatches_ );
00740 }
00741
00742 void MuonIdProducer::fillMuonId(edm::Event& iEvent, const edm::EventSetup& iSetup,
00743 reco::Muon& aMuon,
00744 TrackDetectorAssociator::Direction direction)
00745 {
00746
00747 const reco::Track* track = 0;
00748 if ( ! aMuon.track().isNull() )
00749 track = aMuon.track().get();
00750 else
00751 {
00752 if ( ! aMuon.standAloneMuon().isNull() )
00753 track = aMuon.standAloneMuon().get();
00754 else
00755 throw cms::Exception("FatalError") << "Failed to fill muon id information for a muon with undefined references to tracks";
00756 }
00757
00758 TrackDetMatchInfo info = trackAssociator_.associate(iEvent, iSetup, *track, parameters_, direction);
00759
00760 if ( fillEnergy_ ) {
00761 reco::MuonEnergy muonEnergy;
00762 muonEnergy.em = info.crossedEnergy(TrackDetMatchInfo::EcalRecHits);
00763 muonEnergy.had = info.crossedEnergy(TrackDetMatchInfo::HcalRecHits);
00764 muonEnergy.ho = info.crossedEnergy(TrackDetMatchInfo::HORecHits);
00765 muonEnergy.tower = info.crossedEnergy(TrackDetMatchInfo::TowerTotal);
00766 muonEnergy.emS9 = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits,1);
00767 muonEnergy.emS25 = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits,2);
00768 muonEnergy.hadS9 = info.nXnEnergy(TrackDetMatchInfo::HcalRecHits,1);
00769 muonEnergy.hoS9 = info.nXnEnergy(TrackDetMatchInfo::HORecHits,1);
00770 muonEnergy.towerS9 = info.nXnEnergy(TrackDetMatchInfo::TowerTotal,1);
00771 muonEnergy.ecal_position = info.trkGlobPosAtEcal;
00772 muonEnergy.hcal_position = info.trkGlobPosAtHcal;
00773 if (! info.crossedEcalIds.empty() ) muonEnergy.ecal_id = info.crossedEcalIds.front();
00774 if (! info.crossedHcalIds.empty() ) muonEnergy.hcal_id = info.crossedHcalIds.front();
00775
00776 DetId emMaxId = info.findMaxDeposition(TrackDetMatchInfo::EcalRecHits,2);
00777 for(std::vector<const EcalRecHit*>::const_iterator hit=info.ecalRecHits.begin();
00778 hit!=info.ecalRecHits.end(); ++hit) {
00779 if ((*hit)->id() != emMaxId) continue;
00780 muonEnergy.emMax = (*hit)->energy();
00781 muonEnergy.ecal_time = (*hit)->time();
00782 }
00783 DetId hadMaxId = info.findMaxDeposition(TrackDetMatchInfo::HcalRecHits,1);
00784 for(std::vector<const HBHERecHit*>::const_iterator hit=info.hcalRecHits.begin();
00785 hit!=info.hcalRecHits.end(); ++hit) {
00786 if ((*hit)->id() != hadMaxId) continue;
00787 muonEnergy.hadMax = (*hit)->energy();
00788 muonEnergy.hcal_time = (*hit)->time();
00789 }
00790 aMuon.setCalEnergy( muonEnergy );
00791 }
00792 if ( ! fillMatching_ && ! aMuon.isTrackerMuon() ) return;
00793
00794
00795 std::vector<reco::MuonChamberMatch> muonChamberMatches;
00796 unsigned int nubmerOfMatchesAccordingToTrackAssociator = 0;
00797 for( std::vector<TAMuonChamberMatch>::const_iterator chamber=info.chambers.begin();
00798 chamber!=info.chambers.end(); chamber++ )
00799 {
00800 reco::MuonChamberMatch matchedChamber;
00801
00802 LocalError localError = chamber->tState.localError().positionError();
00803 matchedChamber.x = chamber->tState.localPosition().x();
00804 matchedChamber.y = chamber->tState.localPosition().y();
00805 matchedChamber.xErr = sqrt( localError.xx() );
00806 matchedChamber.yErr = sqrt( localError.yy() );
00807
00808 matchedChamber.dXdZ = chamber->tState.localDirection().z()!=0?chamber->tState.localDirection().x()/chamber->tState.localDirection().z():9999;
00809 matchedChamber.dYdZ = chamber->tState.localDirection().z()!=0?chamber->tState.localDirection().y()/chamber->tState.localDirection().z():9999;
00810
00811 AlgebraicSymMatrix55 trajectoryCovMatrix = chamber->tState.localError().matrix();
00812 matchedChamber.dXdZErr = trajectoryCovMatrix(1,1)>0?sqrt(trajectoryCovMatrix(1,1)):0;
00813 matchedChamber.dYdZErr = trajectoryCovMatrix(2,2)>0?sqrt(trajectoryCovMatrix(2,2)):0;
00814
00815 matchedChamber.edgeX = chamber->localDistanceX;
00816 matchedChamber.edgeY = chamber->localDistanceY;
00817
00818 matchedChamber.id = chamber->id;
00819 if ( ! chamber->segments.empty() ) ++nubmerOfMatchesAccordingToTrackAssociator;
00820
00821
00822 for( std::vector<TAMuonSegmentMatch>::const_iterator segment = chamber->segments.begin();
00823 segment != chamber->segments.end(); segment++ )
00824 {
00825 reco::MuonSegmentMatch matchedSegment;
00826 matchedSegment.x = segment->segmentLocalPosition.x();
00827 matchedSegment.y = segment->segmentLocalPosition.y();
00828 matchedSegment.dXdZ = segment->segmentLocalDirection.z()?segment->segmentLocalDirection.x()/segment->segmentLocalDirection.z():0;
00829 matchedSegment.dYdZ = segment->segmentLocalDirection.z()?segment->segmentLocalDirection.y()/segment->segmentLocalDirection.z():0;
00830 matchedSegment.xErr = segment->segmentLocalErrorXX>0?sqrt(segment->segmentLocalErrorXX):0;
00831 matchedSegment.yErr = segment->segmentLocalErrorYY>0?sqrt(segment->segmentLocalErrorYY):0;
00832 matchedSegment.dXdZErr = segment->segmentLocalErrorDxDz>0?sqrt(segment->segmentLocalErrorDxDz):0;
00833 matchedSegment.dYdZErr = segment->segmentLocalErrorDyDz>0?sqrt(segment->segmentLocalErrorDyDz):0;
00834 matchedSegment.t0 = segment->t0;
00835 matchedSegment.mask = 0;
00836 matchedSegment.dtSegmentRef = segment->dtSegmentRef;
00837 matchedSegment.cscSegmentRef = segment->cscSegmentRef;
00838 matchedSegment.hasZed_ = segment->hasZed;
00839 matchedSegment.hasPhi_ = segment->hasPhi;
00840
00841 bool matchedX = false;
00842 bool matchedY = false;
00843 LogTrace("MuonIdentification") << " matching local x, segment x: " << matchedSegment.x <<
00844 ", chamber x: " << matchedChamber.x << ", max: " << maxAbsDx_;
00845 LogTrace("MuonIdentification") << " matching local y, segment y: " << matchedSegment.y <<
00846 ", chamber y: " << matchedChamber.y << ", max: " << maxAbsDy_;
00847 if (matchedSegment.xErr>0 && matchedChamber.xErr>0 )
00848 LogTrace("MuonIdentification") << " xpull: " <<
00849 fabs(matchedSegment.x - matchedChamber.x)/sqrt(pow(matchedSegment.xErr,2) + pow(matchedChamber.xErr,2));
00850 if (matchedSegment.yErr>0 && matchedChamber.yErr>0 )
00851 LogTrace("MuonIdentification") << " ypull: " <<
00852 fabs(matchedSegment.y - matchedChamber.y)/sqrt(pow(matchedSegment.yErr,2) + pow(matchedChamber.yErr,2));
00853
00854 if (fabs(matchedSegment.x - matchedChamber.x) < maxAbsDx_) matchedX = true;
00855 if (fabs(matchedSegment.y - matchedChamber.y) < maxAbsDy_) matchedY = true;
00856 if (matchedSegment.xErr>0 && matchedChamber.xErr>0 &&
00857 fabs(matchedSegment.x - matchedChamber.x)/sqrt(pow(matchedSegment.xErr,2) + pow(matchedChamber.xErr,2)) < maxAbsPullX_) matchedX = true;
00858 if (matchedSegment.yErr>0 && matchedChamber.yErr>0 &&
00859 fabs(matchedSegment.y - matchedChamber.y)/sqrt(pow(matchedSegment.yErr,2) + pow(matchedChamber.yErr,2)) < maxAbsPullY_) matchedY = true;
00860 if (matchedX && matchedY) matchedChamber.segmentMatches.push_back(matchedSegment);
00861 }
00862 muonChamberMatches.push_back(matchedChamber);
00863 }
00864 aMuon.setMatches(muonChamberMatches);
00865
00866 LogTrace("MuonIdentification") << "number of muon chambers: " << aMuon.matches().size() << "\n"
00867 << "number of chambers with segments according to the associator requirements: " <<
00868 nubmerOfMatchesAccordingToTrackAssociator;
00869 LogTrace("MuonIdentification") << "number of segment matches with the producer requirements: " <<
00870 aMuon.numberOfMatches( reco::Muon::NoArbitration );
00871
00872
00873 }
00874
00875 void MuonIdProducer::fillArbitrationInfo( reco::MuonCollection* pOutputMuons )
00876 {
00877
00878
00879
00880 std::vector<std::pair<reco::MuonChamberMatch*,reco::MuonSegmentMatch*> > chamberPairs;
00881 std::vector<std::pair<reco::MuonChamberMatch*,reco::MuonSegmentMatch*> > stationPairs;
00882 std::vector<std::pair<reco::MuonChamberMatch*,reco::MuonSegmentMatch*> > arbitrationPairs;
00883
00884
00885 for( unsigned int muonIndex1 = 0; muonIndex1 < pOutputMuons->size(); ++muonIndex1 )
00886 {
00887
00888 for( std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = pOutputMuons->at(muonIndex1).matches().begin();
00889 chamberIter1 != pOutputMuons->at(muonIndex1).matches().end(); ++chamberIter1 )
00890 {
00891 if(chamberIter1->segmentMatches.empty()) continue;
00892 chamberPairs.clear();
00893
00894
00895 for( std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
00896 segmentIter1 != chamberIter1->segmentMatches.end(); ++segmentIter1 )
00897 {
00898 chamberPairs.push_back(std::make_pair(&(*chamberIter1), &(*segmentIter1)));
00899 if(!segmentIter1->isMask())
00900 {
00901 arbitrationPairs.clear();
00902 arbitrationPairs.push_back(std::make_pair(&(*chamberIter1), &(*segmentIter1)));
00903
00904
00905
00906
00907
00908 if (pOutputMuons->at(muonIndex1).isTrackerMuon()) {
00909
00910 for( unsigned int muonIndex2 = muonIndex1+1; muonIndex2 < pOutputMuons->size(); ++muonIndex2 )
00911 {
00912
00913 if (! pOutputMuons->at(muonIndex2).isTrackerMuon()) continue;
00914
00915 for( std::vector<reco::MuonChamberMatch>::iterator chamberIter2 = pOutputMuons->at(muonIndex2).matches().begin();
00916 chamberIter2 != pOutputMuons->at(muonIndex2).matches().end(); ++chamberIter2 )
00917 {
00918
00919 for( std::vector<reco::MuonSegmentMatch>::iterator segmentIter2 = chamberIter2->segmentMatches.begin();
00920 segmentIter2 != chamberIter2->segmentMatches.end(); ++segmentIter2 )
00921 {
00922 if(segmentIter2->isMask()) continue;
00923 if(fabs(segmentIter2->x - segmentIter1->x ) < 1E-3 &&
00924 fabs(segmentIter2->y - segmentIter1->y ) < 1E-3 &&
00925 fabs(segmentIter2->dXdZ - segmentIter1->dXdZ ) < 1E-3 &&
00926 fabs(segmentIter2->dYdZ - segmentIter1->dYdZ ) < 1E-3 &&
00927 fabs(segmentIter2->xErr - segmentIter1->xErr ) < 1E-3 &&
00928 fabs(segmentIter2->yErr - segmentIter1->yErr ) < 1E-3 &&
00929 fabs(segmentIter2->dXdZErr - segmentIter1->dXdZErr) < 1E-3 &&
00930 fabs(segmentIter2->dYdZErr - segmentIter1->dYdZErr) < 1E-3)
00931 arbitrationPairs.push_back(std::make_pair(&(*chamberIter2), &(*segmentIter2)));
00932 }
00933 }
00934 }
00935 }
00936
00937
00938 if(arbitrationPairs.empty()) continue;
00939 if(arbitrationPairs.size()==1) {
00940 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDRSlope);
00941 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDXSlope);
00942 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDR);
00943 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDX);
00944 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::Arbitrated);
00945 } else {
00946 sort(arbitrationPairs.begin(), arbitrationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDRSlope));
00947 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDRSlope);
00948 sort(arbitrationPairs.begin(), arbitrationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDXSlope));
00949 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDXSlope);
00950 sort(arbitrationPairs.begin(), arbitrationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDR));
00951 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDR);
00952 sort(arbitrationPairs.begin(), arbitrationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDX));
00953 arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDX);
00954 for( unsigned int it = 0; it < arbitrationPairs.size(); ++it )
00955 arbitrationPairs.at(it).second->setMask(reco::MuonSegmentMatch::Arbitrated);
00956 }
00957 }
00958
00959
00960 if( chamberIter1->id.subdetId() == MuonSubdetId::CSC && arbClean_ ) {
00961 CSCDetId cscid(chamberIter1->id);
00962 if(cscid.ring() == 4)
00963 for( std::vector<reco::MuonSegmentMatch>::iterator segmentIter2 = chamberIter1->segmentMatches.begin();
00964 segmentIter2 != chamberIter1->segmentMatches.end(); ++segmentIter2 ) {
00965 if( segmentIter1->cscSegmentRef.isNonnull() && segmentIter2->cscSegmentRef.isNonnull() )
00966 if( meshAlgo_->isDuplicateOf(segmentIter1->cscSegmentRef,segmentIter2->cscSegmentRef) &&
00967 (segmentIter2->mask & 0x1e0000) &&
00968 (segmentIter1->mask & 0x1e0000) )
00969 segmentIter2->setMask(reco::MuonSegmentMatch::BelongsToTrackByME1aClean);
00970
00971 }
00972 }
00973
00974 }
00975
00976
00977 if(chamberPairs.empty()) continue;
00978 if(chamberPairs.size()==1) {
00979 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDRSlope);
00980 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDXSlope);
00981 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDR);
00982 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDX);
00983 } else {
00984 sort(chamberPairs.begin(), chamberPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDRSlope));
00985 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDRSlope);
00986 sort(chamberPairs.begin(), chamberPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDXSlope));
00987 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDXSlope);
00988 sort(chamberPairs.begin(), chamberPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDR));
00989 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDR);
00990 sort(chamberPairs.begin(), chamberPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDX));
00991 chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDX);
00992 }
00993 }
00994
00995
00996 for( int stationIndex = 1; stationIndex < 5; ++stationIndex )
00997 for( int detectorIndex = 1; detectorIndex < 4; ++detectorIndex )
00998 {
00999 stationPairs.clear();
01000
01001
01002 for( std::vector<reco::MuonChamberMatch>::iterator chamberIter = pOutputMuons->at(muonIndex1).matches().begin();
01003 chamberIter != pOutputMuons->at(muonIndex1).matches().end(); ++chamberIter )
01004 {
01005 if(!(chamberIter->station()==stationIndex && chamberIter->detector()==detectorIndex)) continue;
01006 if(chamberIter->segmentMatches.empty()) continue;
01007
01008 for( std::vector<reco::MuonSegmentMatch>::iterator segmentIter = chamberIter->segmentMatches.begin();
01009 segmentIter != chamberIter->segmentMatches.end(); ++segmentIter )
01010 stationPairs.push_back(std::make_pair(&(*chamberIter), &(*segmentIter)));
01011 }
01012
01013 if(stationPairs.empty()) continue;
01014 if(stationPairs.size()==1) {
01015 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDRSlope);
01016 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDXSlope);
01017 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDR);
01018 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDX);
01019 } else {
01020 sort(stationPairs.begin(), stationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDRSlope));
01021 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDRSlope);
01022 sort(stationPairs.begin(), stationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDXSlope));
01023 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDXSlope);
01024 sort(stationPairs.begin(), stationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDR));
01025 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDR);
01026 sort(stationPairs.begin(), stationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDX));
01027 stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDX);
01028 }
01029 }
01030
01031 }
01032
01033 if(arbClean_) {
01034
01035 meshAlgo_->clearMesh();
01036 meshAlgo_->runMesh(pOutputMuons);
01037 }
01038 }
01039
01040 void MuonIdProducer::fillMuonIsolation(edm::Event& iEvent, const edm::EventSetup& iSetup, reco::Muon& aMuon,
01041 reco::IsoDeposit& trackDep, reco::IsoDeposit& ecalDep, reco::IsoDeposit& hcalDep, reco::IsoDeposit& hoDep,
01042 reco::IsoDeposit& jetDep)
01043 {
01044 reco::MuonIsolation isoR03, isoR05;
01045 const reco::Track* track = 0;
01046 if ( ! aMuon.track().isNull() )
01047 track = aMuon.track().get();
01048 else
01049 {
01050 if ( ! aMuon.standAloneMuon().isNull() )
01051 track = aMuon.standAloneMuon().get();
01052 else
01053 throw cms::Exception("FatalError") << "Failed to compute muon isolation information for a muon with undefined references to tracks";
01054 }
01055
01056
01057 reco::IsoDeposit depTrk = muIsoExtractorTrack_->deposit(iEvent, iSetup, *track );
01058 std::vector<reco::IsoDeposit> caloDeps = muIsoExtractorCalo_->deposits(iEvent, iSetup, *track);
01059 reco::IsoDeposit depJet = muIsoExtractorJet_->deposit(iEvent, iSetup, *track );
01060
01061 if(caloDeps.size()!=3) {
01062 LogTrace("MuonIdentification") << "Failed to fill vector of calorimeter isolation deposits!";
01063 return;
01064 }
01065
01066 reco::IsoDeposit depEcal = caloDeps.at(0);
01067 reco::IsoDeposit depHcal = caloDeps.at(1);
01068 reco::IsoDeposit depHo = caloDeps.at(2);
01069
01070 trackDep = depTrk;
01071 ecalDep = depEcal;
01072 hcalDep = depHcal;
01073 hoDep = depHo;
01074 jetDep = depJet;
01075
01076 isoR03.sumPt = depTrk.depositWithin(0.3);
01077 isoR03.emEt = depEcal.depositWithin(0.3);
01078 isoR03.hadEt = depHcal.depositWithin(0.3);
01079 isoR03.hoEt = depHo.depositWithin(0.3);
01080 isoR03.nTracks = depTrk.depositAndCountWithin(0.3).second;
01081 isoR03.nJets = depJet.depositAndCountWithin(0.3).second;
01082 isoR03.trackerVetoPt = depTrk.candEnergy();
01083 isoR03.emVetoEt = depEcal.candEnergy();
01084 isoR03.hadVetoEt = depHcal.candEnergy();
01085 isoR03.hoVetoEt = depHo.candEnergy();
01086
01087 isoR05.sumPt = depTrk.depositWithin(0.5);
01088 isoR05.emEt = depEcal.depositWithin(0.5);
01089 isoR05.hadEt = depHcal.depositWithin(0.5);
01090 isoR05.hoEt = depHo.depositWithin(0.5);
01091 isoR05.nTracks = depTrk.depositAndCountWithin(0.5).second;
01092 isoR05.nJets = depJet.depositAndCountWithin(0.5).second;
01093 isoR05.trackerVetoPt = depTrk.candEnergy();
01094 isoR05.emVetoEt = depEcal.candEnergy();
01095 isoR05.hadVetoEt = depHcal.candEnergy();
01096 isoR05.hoVetoEt = depHo.candEnergy();
01097
01098
01099 aMuon.setIsolation(isoR03, isoR05);
01100
01101 }
01102
01103 reco::Muon MuonIdProducer::makeMuon( const reco::Track& track )
01104 {
01105
01106 double energy = sqrt(track.p() * track.p() + 0.011163691);
01107 math::XYZTLorentzVector p4(track.px(),
01108 track.py(),
01109 track.pz(),
01110 energy);
01111 return reco::Muon( track.charge(), p4, track.vertex() );
01112 }
01113
01114 double MuonIdProducer::sectorPhi( const DetId& id )
01115 {
01116 double phi = 0;
01117 if( id.subdetId() == MuonSubdetId::DT ) {
01118 DTChamberId muonId(id.rawId());
01119 if ( muonId.sector() <= 12 )
01120 phi = (muonId.sector()-1)/6.*M_PI;
01121 if ( muonId.sector() == 13 ) phi = 3/6.*M_PI;
01122 if ( muonId.sector() == 14 ) phi = 9/6.*M_PI;
01123 }
01124 if( id.subdetId() == MuonSubdetId::CSC ) {
01125 CSCDetId muonId(id.rawId());
01126 phi = M_PI/4+(muonId.triggerSector()-1)/3.*M_PI;
01127 }
01128 if ( phi > M_PI ) phi -= 2*M_PI;
01129 return phi;
01130 }
01131
01132 double MuonIdProducer::phiOfMuonIneteractionRegion( const reco::Muon& muon ) const
01133 {
01134 if ( muon.isStandAloneMuon() ) return muon.standAloneMuon()->innerPosition().phi();
01135
01136 if ( muon.matches().empty() ){
01137 if ( muon.innerTrack().isAvailable() &&
01138 muon.innerTrack()->extra().isAvailable() )
01139 return muon.innerTrack()->outerPosition().phi();
01140 else
01141 return muon.phi();
01142 }
01143 return sectorPhi(muon.matches().at(0).id);
01144 }
01145
01146 void MuonIdProducer::fillGlbQuality(edm::Event& iEvent, const edm::EventSetup& iSetup, reco::Muon& aMuon)
01147 {
01148 edm::Handle<edm::ValueMap<reco::MuonQuality> > glbQualH;
01149 iEvent.getByLabel(globalTrackQualityInputTag_, glbQualH);
01150
01151 if(aMuon.isGlobalMuon() && glbQualH.isValid() && !glbQualH.failedToGet()) {
01152 aMuon.setCombinedQuality((*glbQualH)[aMuon.combinedMuon()]);
01153 }
01154
01155 LogDebug("MuonIdentification") << "tkChiVal " << aMuon.combinedQuality().trkRelChi2;
01156
01157 }
01158
01159 void MuonIdProducer::fillTrackerKink( reco::Muon& aMuon ) {
01160
01161 if (aMuon.innerTrack().isNull()) return;
01162
01163 reco::MuonQuality quality = (aMuon.isQualityValid() ? aMuon.combinedQuality() : reco::MuonQuality());
01164
01165 bool filled = trackerKinkFinder_->fillTrkKink(quality, *aMuon.innerTrack());
01166
01167 if (filled || aMuon.isQualityValid()) aMuon.setCombinedQuality(quality);
01168 }
01169
01170 bool MuonIdProducer::checkLinks(const reco::MuonTrackLinks* links) const {
01171 bool trackBAD = links->trackerTrack().isNull();
01172 bool staBAD = links->standAloneTrack().isNull();
01173 bool glbBAD = links->globalTrack().isNull();
01174 if (trackBAD || staBAD || glbBAD )
01175 {
01176 edm::LogWarning("muonIDbadLinks") << "Global muon links to constituent tracks are invalid: trkBad "
01177 <<trackBAD <<" standaloneBad "<<staBAD<<" globalBad "<<glbBAD
01178 <<". There should be no such object. Muon is skipped.";
01179 return false;
01180 }
01181 return true;
01182 }
1.7.3