CMSSW Reference Manual

00001 // -*- C++ -*-
00002 //
00003 // Package:    MuonIdentification
00004 // Class:      MuonIdProducer
00005 // 
00006 //
00007 // Original Author:  Dmytro Kovalskyi
00008 // $Id: MuonIdProducer.cc,v 1.29.2.2 2009/01/16 03:54:34 slava77 Exp $
00009 //
00010 //
00011 
00012 
00013 // system include files
00014 #include <memory>
00015 
00016 // user include files
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/MuonTime.h"
00029 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
00030 
00031 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
00032 #include "Utilities/Timing/interface/TimerStack.h"
00033 
00034 #include <boost/regex.hpp>
00035 #include "RecoMuon/MuonIdentification/plugins/MuonIdProducer.h"
00036 #include "RecoMuon/MuonIdentification/interface/MuonIdTruthInfo.h"
00037 #include "RecoMuon/MuonIdentification/interface/MuonArbitrationMethods.h"
00038 
00039 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositExtractorFactory.h"
00040 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00041 
00042 #include <algorithm>
00043 
00044 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00045 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
00046 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
00047 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
00048 
00049 MuonIdProducer::MuonIdProducer(const edm::ParameterSet& iConfig):
00050 muIsoExtractorCalo_(0),muIsoExtractorTrack_(0),muIsoExtractorJet_(0)
00051 {
00052    produces<reco::MuonCollection>();
00053    
00054    minPt_                   = iConfig.getParameter<double>("minPt");
00055    minP_                    = iConfig.getParameter<double>("minP");
00056    minNumberOfMatches_      = iConfig.getParameter<int>("minNumberOfMatches");
00057    addExtraSoftMuons_       = iConfig.getParameter<bool>("addExtraSoftMuons");
00058    maxAbsEta_               = iConfig.getParameter<double>("maxAbsEta");
00059    maxAbsDx_                = iConfig.getParameter<double>("maxAbsDx");
00060    maxAbsPullX_             = iConfig.getParameter<double>("maxAbsPullX");
00061    maxAbsDy_                = iConfig.getParameter<double>("maxAbsDy");
00062    maxAbsPullY_             = iConfig.getParameter<double>("maxAbsPullY");
00063    fillCaloCompatibility_   = iConfig.getParameter<bool>("fillCaloCompatibility");
00064    fillEnergy_              = iConfig.getParameter<bool>("fillEnergy");
00065    fillMatching_            = iConfig.getParameter<bool>("fillMatching");
00066    fillIsolation_           = iConfig.getParameter<bool>("fillIsolation");
00067    trackPtThresholdToFillCandidateP4WithGlobalFit_ = iConfig.getParameter<double>("trackPtThresholdToFillCandidateP4WithGlobalFit");
00068    
00069    // Load TrackDetectorAssociator parameters
00070    edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
00071    parameters_.loadParameters( parameters );
00072 
00073    // Load parameters for the TimingExtractor
00074    edm::ParameterSet timingParameters = iConfig.getParameter<edm::ParameterSet>("timingParameters");
00075    theTimingExtractor_ = new MuonTimingExtractor(timingParameters);
00076    
00077    if (fillCaloCompatibility_){
00078       // Load MuonCaloCompatibility parameters
00079       parameters = iConfig.getParameter<edm::ParameterSet>("MuonCaloCompatibility");
00080       muonCaloCompatibility_.configure( parameters );
00081    }
00082 
00083    if (fillIsolation_){
00084       // Load MuIsoExtractor parameters
00085       edm::ParameterSet caloExtractorPSet = iConfig.getParameter<edm::ParameterSet>("CaloExtractorPSet");
00086       std::string caloExtractorName = caloExtractorPSet.getParameter<std::string>("ComponentName");
00087       muIsoExtractorCalo_ = IsoDepositExtractorFactory::get()->create( caloExtractorName, caloExtractorPSet);
00088 
00089       edm::ParameterSet trackExtractorPSet = iConfig.getParameter<edm::ParameterSet>("TrackExtractorPSet");
00090       std::string trackExtractorName = trackExtractorPSet.getParameter<std::string>("ComponentName");
00091       muIsoExtractorTrack_ = IsoDepositExtractorFactory::get()->create( trackExtractorName, trackExtractorPSet);
00092 
00093       edm::ParameterSet jetExtractorPSet = iConfig.getParameter<edm::ParameterSet>("JetExtractorPSet");
00094       std::string jetExtractorName = jetExtractorPSet.getParameter<std::string>("ComponentName");
00095       muIsoExtractorJet_ = IsoDepositExtractorFactory::get()->create( jetExtractorName, jetExtractorPSet);
00096    }
00097    
00098    inputCollectionLabels_ = iConfig.getParameter<std::vector<edm::InputTag> >("inputCollectionLabels");
00099    inputCollectionTypes_  = iConfig.getParameter<std::vector<std::string> >("inputCollectionTypes");
00100    if (inputCollectionLabels_.size() != inputCollectionTypes_.size()) 
00101      throw cms::Exception("ConfigurationError") << "Number of input collection labels is different from number of types. " <<
00102      "For each collection label there should be exactly one collection type specified.";
00103    if (inputCollectionLabels_.size()>4 ||inputCollectionLabels_.empty()) 
00104      throw cms::Exception("ConfigurationError") << "Number of input collections should be from 1 to 4.";
00105    
00106    debugWithTruthMatching_    = iConfig.getParameter<bool>("debugWithTruthMatching");
00107    if (debugWithTruthMatching_) edm::LogWarning("MuonIdentification") 
00108      << "========================================================================\n" 
00109      << "Debugging mode with truth matching is turned on!!! Make sure you understand what you are doing!\n"
00110      << "========================================================================\n";
00111 }
00112 
00113 
00114 MuonIdProducer::~MuonIdProducer()
00115 {
00116    if (muIsoExtractorCalo_) delete muIsoExtractorCalo_;
00117    if (muIsoExtractorTrack_) delete muIsoExtractorTrack_;
00118    if (muIsoExtractorJet_) delete muIsoExtractorJet_;
00119    if (theTimingExtractor_) delete theTimingExtractor_;
00120    // TimingReport::current()->dump(std::cout);
00121 }
00122 
00123 void MuonIdProducer::init(edm::Event& iEvent, const edm::EventSetup& iSetup)
00124 {
00125    // TimerStack timers;
00126    // timers.push("MuonIdProducer::produce::init");
00127    
00128    innerTrackCollectionHandle_.clear();
00129    outerTrackCollectionHandle_.clear();
00130    linkCollectionHandle_.clear();
00131    muonCollectionHandle_.clear();
00132    
00133    // timers.push("MuonIdProducer::produce::init::getPropagator");
00134    edm::ESHandle<Propagator> propagator;
00135    iSetup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagator);
00136    trackAssociator_.setPropagator(propagator.product());
00137    
00138    // timers.pop_and_push("MuonIdProducer::produce::init::getInputCollections");
00139    for ( unsigned int i = 0; i < inputCollectionLabels_.size(); ++i ) {
00140       if ( inputCollectionTypes_[i] == "inner tracks" ) {
00141          iEvent.getByLabel(inputCollectionLabels_[i], innerTrackCollectionHandle_);
00142          if (! innerTrackCollectionHandle_.isValid()) 
00143            throw cms::Exception("FatalError") << "Failed to get input track collection with label: " << inputCollectionLabels_[i];
00144          LogTrace("MuonIdentification") << "Number of input inner tracks: " << innerTrackCollectionHandle_->size();
00145          continue;
00146       }
00147       if ( inputCollectionTypes_[i] == "outer tracks" ) {
00148          iEvent.getByLabel(inputCollectionLabels_[i], outerTrackCollectionHandle_);
00149          if (! outerTrackCollectionHandle_.isValid()) 
00150            throw cms::Exception("FatalError") << "Failed to get input track collection with label: " << inputCollectionLabels_[i];
00151          LogTrace("MuonIdentification") << "Number of input outer tracks: " << outerTrackCollectionHandle_->size();
00152          continue;
00153       }
00154       if ( inputCollectionTypes_[i] == "links" ) {
00155          iEvent.getByLabel(inputCollectionLabels_[i], linkCollectionHandle_);
00156          if (! linkCollectionHandle_.isValid()) 
00157            throw cms::Exception("FatalError") << "Failed to get input link collection with label: " << inputCollectionLabels_[i];
00158          LogTrace("MuonIdentification") << "Number of input links: " << linkCollectionHandle_->size();
00159          continue;
00160       }
00161       if ( inputCollectionTypes_[i] == "muons" ) {
00162          iEvent.getByLabel(inputCollectionLabels_[i], muonCollectionHandle_);
00163          if (! muonCollectionHandle_.isValid()) 
00164            throw cms::Exception("FatalError") << "Failed to get input muon collection with label: " << inputCollectionLabels_[i];
00165          LogTrace("MuonIdentification") << "Number of input muons: " << muonCollectionHandle_->size();
00166          continue;
00167       }
00168       throw cms::Exception("FatalError") << "Unknown input collection type: " << inputCollectionTypes_[i];
00169    }
00170 }
00171 
00172 reco::Muon MuonIdProducer::makeMuon(edm::Event& iEvent, const edm::EventSetup& iSetup, 
00173                                      const reco::TrackRef& track, MuonIdProducer::TrackType type)
00174 {
00175    LogTrace("MuonIdentification") << "Creating a muon from a track " << track.get()->pt() << 
00176      " Pt (GeV), eta: " << track.get()->eta();
00177    reco::Muon aMuon( makeMuon( *(track.get()) ) );
00178    switch (type) {
00179     case InnerTrack: 
00180       aMuon.setInnerTrack( track );
00181       break;
00182     case OuterTrack:
00183       aMuon.setOuterTrack( track );
00184       break;
00185     case CombinedTrack:
00186       aMuon.setGlobalTrack( track );
00187       break;
00188    }
00189    return aMuon;
00190 }
00191 
00192 reco::Muon MuonIdProducer::makeMuon( const reco::MuonTrackLinks& links )
00193 {
00194    LogTrace("MuonIdentification") << "Creating a muon from a link to tracks object";
00195    reco::Muon aMuon;
00196    if ( links.trackerTrack()->pt() > trackPtThresholdToFillCandidateP4WithGlobalFit_ )
00197      aMuon = makeMuon( *(links.globalTrack()) );
00198    else
00199      aMuon = makeMuon( *(links.trackerTrack()) );
00200      
00201    aMuon.setInnerTrack( links.trackerTrack() );
00202    aMuon.setOuterTrack( links.standAloneTrack() );
00203    aMuon.setGlobalTrack( links.globalTrack() );
00204    return aMuon;
00205 }
00206 
00207 bool MuonIdProducer::isGoodTrack( const reco::Track& track )
00208 {
00209    // Pt and absolute momentum requirement
00210    if (track.pt() < minPt_ && track.p() < minP_){ 
00211       LogTrace("MuonIdentification") << "Skipped low momentum track (Pt,P): " << track.pt() <<
00212         ", " << track.p() << " GeV";
00213       return false;
00214    }
00215         
00216    // Eta requirement
00217    if ( fabs(track.eta()) > maxAbsEta_ ){
00218       LogTrace("MuonIdentification") << "Skipped track with large pseudo rapidity (Eta: " << track.eta() << " )";
00219       return false;
00220    }
00221    return true;
00222 }
00223    
00224 unsigned int MuonIdProducer::chamberId( const DetId& id )
00225 {
00226    switch ( id.det() ) {
00227     case DetId::Muon:
00228       switch ( id.subdetId() ) {
00229        case MuonSubdetId::DT:
00230            { 
00231               DTChamberId detId(id.rawId());
00232               return detId.rawId();
00233            }
00234          break;
00235        case MuonSubdetId::CSC:
00236            {
00237               CSCDetId detId(id.rawId());
00238               return detId.chamberId().rawId();
00239            }
00240          break;
00241        default:
00242          return 0;
00243       }
00244     default:
00245       return 0;
00246    }
00247    return 0;
00248 }
00249 
00250 
00251 int MuonIdProducer::overlap(const reco::Muon& muon, const reco::Track& track)
00252 {
00253    int numberOfCommonDetIds = 0;
00254    if ( ! muon.isMatchesValid() || 
00255         track.extra().isNull() ||
00256         track.extra()->recHits().isNull() ) return numberOfCommonDetIds;
00257    const std::vector<reco::MuonChamberMatch>& matches( muon.matches() );
00258    for ( std::vector<reco::MuonChamberMatch>::const_iterator match = matches.begin();
00259          match != matches.end(); ++match ) 
00260      {
00261         if ( match->segmentMatches.empty() ) continue;
00262         bool foundCommonDetId = false;
00263         
00264         for ( TrackingRecHitRefVector::const_iterator hit = track.extra()->recHitsBegin();
00265               hit != track.extra()->recHitsEnd(); ++hit )
00266           {
00267              // LogTrace("MuonIdentification") << "hit DetId: " << std::hex << hit->get()->geographicalId().rawId() <<
00268              //  "\t hit chamber DetId: " << getChamberId(hit->get()->geographicalId()) <<
00269              //  "\t segment DetId: " << match->id.rawId() << std::dec;
00270              
00271              if ( chamberId(hit->get()->geographicalId()) == match->id.rawId() ) {
00272                 foundCommonDetId = true;
00273                 break;
00274              }
00275           }
00276         if ( foundCommonDetId ) {
00277            numberOfCommonDetIds++;
00278            break;
00279         }
00280      }
00281    return numberOfCommonDetIds;
00282 }
00283 
00284 
00285 void MuonIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00286 {
00287    using namespace edm;
00288    
00289    // TimerStack timers;
00290    // timers.push("MuonIdProducer::produce");
00291 
00292    std::auto_ptr<reco::MuonCollection> outputMuons(new reco::MuonCollection);
00293    init(iEvent, iSetup);
00294 
00295    // loop over input collections
00296    
00297    // muons first
00298    if ( muonCollectionHandle_.isValid() )
00299      for ( reco::MuonCollection::const_iterator muon = muonCollectionHandle_->begin();
00300            muon !=  muonCollectionHandle_->end(); ++muon )
00301        outputMuons->push_back(*muon);
00302    
00303    // links second ( assume global muon type )
00304    if ( linkCollectionHandle_.isValid() )
00305      for ( reco::MuonTrackLinksCollection::const_iterator links = linkCollectionHandle_->begin();
00306            links != linkCollectionHandle_->end(); ++links )
00307        {
00308           if ( links->trackerTrack().isNull() ||
00309                links->standAloneTrack().isNull() ||
00310                links->globalTrack().isNull() ) 
00311             {
00312                edm::LogWarning("MuonIdentification") << "Global muon links to constituent tracks are invalid. There should be no such object. Muon is skipped.";
00313                continue;
00314             }
00315           // check if this muon is already in the list
00316           bool newMuon = true;
00317           for ( reco::MuonCollection::const_iterator muon = outputMuons->begin();
00318                 muon !=  outputMuons->end(); ++muon )
00319              if ( muon->track() == links->trackerTrack() &&
00320                   muon->standAloneMuon() == links->standAloneTrack() &&
00321                   muon->combinedMuon() == links->globalTrack() )
00322               newMuon = false;
00323           if ( newMuon ) {
00324              outputMuons->push_back( makeMuon( *links ) );
00325              outputMuons->back().setType(reco::Muon::GlobalMuon | reco::Muon::StandAloneMuon);
00326           }
00327        }
00328    
00329    // tracker muon is next
00330    if ( innerTrackCollectionHandle_.isValid() ) {
00331       LogTrace("MuonIdentification") << "Creating tracker muons";
00332       for ( unsigned int i = 0; i < innerTrackCollectionHandle_->size(); ++i )
00333         {
00334            const reco::Track& track = innerTrackCollectionHandle_->at(i);
00335            if ( ! isGoodTrack( track ) ) continue;
00336            bool splitTrack = false;
00337            if ( track.extra().isAvailable() && 
00338                 TrackDetectorAssociator::crossedIP( track ) ) splitTrack = true;
00339            std::vector<TrackDetectorAssociator::Direction> directions;
00340            if ( splitTrack ) {
00341               directions.push_back(TrackDetectorAssociator::InsideOut);
00342               directions.push_back(TrackDetectorAssociator::OutsideIn);
00343            } else {
00344               directions.push_back(TrackDetectorAssociator::Any);
00345            }
00346            for ( std::vector<TrackDetectorAssociator::Direction>::const_iterator direction = directions.begin();
00347                  direction != directions.end(); ++direction )
00348              {
00349                 // make muon
00350                 // timers.push("MuonIdProducer::produce::fillMuonId");
00351                 reco::Muon trackerMuon( makeMuon(iEvent, iSetup, reco::TrackRef( innerTrackCollectionHandle_, i ), InnerTrack ) );
00352                 trackerMuon.setType( reco::Muon::TrackerMuon );
00353                 fillMuonId(iEvent, iSetup, trackerMuon, *direction);
00354                 if ( ! isGoodTrackerMuon( trackerMuon ) ){
00355                    LogTrace("MuonIdentification") << "track failed minimal number of muon matches requirement";
00356                    continue;
00357                 }
00358                 // timers.pop();
00359           
00360                 if ( debugWithTruthMatching_ ) {
00361                    // add MC hits to a list of matched segments. 
00362                    // Since it's debugging mode - code is slow
00363                    MuonIdTruthInfo::truthMatchMuon(iEvent, iSetup, trackerMuon);
00364                 }
00365           
00366                 // check if this muon is already in the list
00367                 // have to check where muon hits are really located
00368                 // to match properly
00369                 bool newMuon = true;
00370                 for ( reco::MuonCollection::iterator muon = outputMuons->begin();
00371                       muon !=  outputMuons->end(); ++muon )
00372                   {
00373                      if ( muon->innerTrack().get() == trackerMuon.innerTrack().get() &&
00374                           cos(phiOfMuonIneteractionRegion(*muon) - 
00375                               phiOfMuonIneteractionRegion(trackerMuon)) > 0 )
00376                        {
00377                           newMuon = false;
00378                           muon->setMatches( trackerMuon.matches() );
00379                           if (trackerMuon.isTimeValid()) muon->setTime( trackerMuon.time() );
00380                           if (trackerMuon.isEnergyValid()) muon->setCalEnergy( trackerMuon.calEnergy() );
00381                           muon->setType( muon->type() | reco::Muon::TrackerMuon );
00382                           LogTrace("MuonIdentification") << "Found a corresponding global muon. Set energy, matches and move on";
00383                           break;
00384                        }
00385                   }
00386                 if ( newMuon ) outputMuons->push_back( trackerMuon );
00387              }
00388         }
00389    }
00390    
00391    // and at last the stand alone muons
00392    if ( outerTrackCollectionHandle_.isValid() ) {
00393       LogTrace("MuonIdentification") << "Looking for new muons among stand alone muon tracks";
00394       for ( unsigned int i = 0; i < outerTrackCollectionHandle_->size(); ++i )
00395         {
00396            // check if this muon is already in the list of global muons
00397            bool newMuon = true;
00398            for ( reco::MuonCollection::iterator muon = outputMuons->begin();
00399                  muon !=  outputMuons->end(); ++muon ) 
00400              {
00401                 if ( ! muon->standAloneMuon().isNull() ) {
00402                    // global muon
00403                    if ( muon->standAloneMuon().get() ==  &(outerTrackCollectionHandle_->at(i)) ) {
00404                       newMuon = false;
00405                       break;
00406                    }
00407                 } else {
00408                    // tracker muon - no direct links to the standalone muon
00409                    // since we have only a few real muons in an event, matching 
00410                    // the stand alone muon to the tracker muon by DetIds should 
00411                    // be good enough for association. At the end it's up to a 
00412                    // user to redefine the association and what it means. Here 
00413                    // we would like to avoid obvious double counting and we 
00414                    // tolerate a potential miss association
00415                    if ( overlap(*muon,outerTrackCollectionHandle_->at(i))>0 ) {
00416                       LogTrace("MuonIdentification") << "Found associated tracker muon. Set a reference and move on";
00417                       newMuon = false;
00418                       muon->setOuterTrack( reco::TrackRef( outerTrackCollectionHandle_, i ) );
00419                       muon->setType( muon->type() | reco::Muon::StandAloneMuon );
00420                       break;
00421                    }
00422                 }
00423              }
00424            if ( newMuon ) {
00425               LogTrace("MuonIdentification") << "No associated stand alone track is found. Making a muon";
00426               outputMuons->push_back( makeMuon(iEvent, iSetup, 
00427                                                reco::TrackRef( outerTrackCollectionHandle_, i ), OuterTrack ) );
00428               outputMuons->back().setType( reco::Muon::StandAloneMuon );
00429            }
00430         }
00431    }
00432    
00433    LogTrace("MuonIdentification") << "Dress up muons if it's necessary";
00434    // Fill various information
00435    for ( reco::MuonCollection::iterator muon = outputMuons->begin(); muon != outputMuons->end(); ++muon )
00436      {
00437         // Fill muonID
00438         // timers.push("MuonIdProducer::produce::fillMuonId");
00439         if ( ( fillMatching_ && ! muon->isMatchesValid() ) || 
00440              ( fillEnergy_ && !muon->isEnergyValid() ) )
00441           {
00442              // predict direction based on the muon interaction region location 
00443              // if it's available
00444              if ( muon->isStandAloneMuon() ) {
00445                 if ( cos(phiOfMuonIneteractionRegion(*muon) - muon->phi()) > 0 )
00446                   fillMuonId(iEvent, iSetup, *muon, TrackDetectorAssociator::InsideOut);
00447                 else
00448                   fillMuonId(iEvent, iSetup, *muon, TrackDetectorAssociator::OutsideIn);
00449              } else {
00450                 LogTrace("MuonIdentification") << "THIS SHOULD NEVER HAPPEN";
00451                 fillMuonId(iEvent, iSetup, *muon);
00452              }
00453           }
00454         
00455         // timers.push("MuonIdProducer::produce::fillCaloCompatibility");
00456         if ( fillCaloCompatibility_ ) muon->setCaloCompatibility( muonCaloCompatibility_.evaluate(*muon) );
00457         // timers.pop();
00458         
00459         // timers.push("MuonIdProducer::produce::fillIsolation");
00460         if ( fillIsolation_ ) fillMuonIsolation(iEvent, iSetup, *muon);
00461         // timers.pop();
00462 
00463         // fill timing information
00464         reco::MuonTime muonTime;
00465 
00466         if ( ! muon->standAloneMuon().isNull() )        
00467           muonTime = theTimingExtractor_->fillTiming(iEvent,iSetup,muon->standAloneMuon());
00468 
00469         if (muonTime.nStations) {
00470           LogTrace("MuonIdentification") << "Global 1/beta: " << muonTime.inverseBeta << " +/- " << muonTime.inverseBetaErr<<std::endl;
00471           LogTrace("MuonIdentification") << "  Free 1/beta: " << muonTime.freeInverseBeta << " +/- " << muonTime.freeInverseBetaErr<<std::endl;
00472           LogTrace("MuonIdentification") << "  Vertex time (in-out): " << muonTime.timeAtIpInOut << " +/- " << muonTime.timeAtIpInOutErr
00473                                          << "  # of points: " << muonTime.nStations <<std::endl;
00474           LogTrace("MuonIdentification") << "  Vertex time (out-in): " << muonTime.timeAtIpOutIn << " +/- " << muonTime.timeAtIpOutInErr<<std::endl;
00475           LogTrace("MuonIdentification") << "  direction: "   << muonTime.direction() << std::endl;
00476         
00477           muon->setTime(muonTime);
00478         }  
00479         
00480      }
00481         
00482    LogTrace("MuonIdentification") << "number of muons produced: " << outputMuons->size();
00483    // timers.push("MuonIdProducer::produce::fillArbitration");
00484    if ( fillMatching_ ) fillArbitrationInfo( outputMuons.get() );
00485    // timers.pop();
00486    iEvent.put(outputMuons);
00487 }
00488 
00489 void MuonIdProducer::fillTime(edm::Event& iEvent, const edm::EventSetup& iSetup,
00490                               reco::Muon& muon)
00491 {
00492    using namespace edm;
00493    
00494    ESHandle<GlobalTrackingGeometry> theTrackingGeometry;
00495    iSetup.get<GlobalTrackingGeometryRecord>().get(theTrackingGeometry);
00496 
00497    std::vector <double> distance;
00498    std::vector <double> freeT0;
00499 
00500    // loop over chambers to collect timing information assuming it was filled
00501    // earlier
00502    const std::vector<reco::MuonChamberMatch>& chambers = muon.matches();
00503         
00504    for( std::vector<reco::MuonChamberMatch>::const_iterator chamber=chambers.begin(); 
00505         chamber!=chambers.end(); ++chamber ) {
00506 
00507       const GeomDet* geomDet = theTrackingGeometry->idToDet(chamber->id());
00508 
00509       // loop over segments
00510       for( std::vector<reco::MuonSegmentMatch>::const_iterator segment=chamber->segmentMatches.begin(); 
00511            segment!=chamber->segmentMatches.end(); ++segment ) {
00512 
00513          // if we have no t0 measurement in this segment - leave
00514          if (fabs(segment->t0)<1e-6) {
00515             LogTrace("MuonIdentification") << "have no t0 measurement in this segment, leave";
00516             break;
00517          }
00518          
00519          LocalPoint segmInChamber(segment->x,segment->y,0);
00520          double dist = geomDet->toGlobal(segmInChamber).mag();
00521          
00522          distance.push_back(dist);
00523          // full time offset
00524          freeT0.push_back(segment->t0+dist/30.);
00525          LogTrace("MuonIdentification") << "Segment t0: " << segment->t0 << " local 1/beta: " << 1.+segment->t0/dist*30. << std::endl;
00526          
00527          break;
00528       }
00529           
00530    }
00531 
00532         reco::MuonTime muonTime;
00533         double invBeta = 0.;
00534         double invBeta2 = 0.;
00535         double localInvBeta;
00536         double inOutVertexTime = 0.;
00537         double inOutVertexTime2 = 0;
00538         double outInVertexTime = 0.;
00539         double outInVertexTime2 = 0;
00540         int npoints = freeT0.size();
00541 
00542         muonTime.nStations=npoints;
00543         
00544         if (npoints) {
00545         
00546            for (int i=0;i<npoints;i++) {
00547               localInvBeta = freeT0[i]/distance[i]*30.;
00548               invBeta     += localInvBeta;
00549               invBeta2    += localInvBeta*localInvBeta;
00550               inOutVertexTime  += freeT0[i] - distance[i]/30.;
00551               inOutVertexTime2 += (freeT0[i]-distance[i]/30.)*(freeT0[i]-distance[i]/30.);
00552               outInVertexTime  += freeT0[i] + distance[i]/30.;
00553               outInVertexTime2 += (freeT0[i]+distance[i]/30.)*(freeT0[i]+distance[i]/30.);
00554            }
00555         
00556            invBeta/=npoints;
00557            muonTime.inverseBeta=invBeta;                        
00558            muonTime.inverseBetaErr=sqrt(invBeta2/npoints-invBeta*invBeta);  
00559           
00560            inOutVertexTime /= npoints;
00561            outInVertexTime /= npoints;
00562            inOutVertexTime2 /= npoints;
00563            outInVertexTime2 /= npoints;
00564            
00565            muonTime.timeAtIpInOut = inOutVertexTime;
00566            muonTime.timeAtIpInOutErr = sqrt(inOutVertexTime2 - inOutVertexTime*inOutVertexTime);
00567            
00568            muonTime.timeAtIpOutIn = outInVertexTime;
00569            muonTime.timeAtIpOutInErr = sqrt(outInVertexTime2 - outInVertexTime*outInVertexTime);
00570           
00571           // do the unconstrained caclucation, if we have at least two points
00572           if (npoints>1) {
00573             double s=0.,sx=0.,sy=0.,x,y;
00574             double sxx=0.,sxy=0.;
00575             
00576             for (int i=0; i<npoints; i++) {
00577               x=distance[i]/30.;
00578               y=freeT0[i];
00579               sy+=y;
00580               sxy+=x*y;
00581               s+=1.;
00582               sx+=x;
00583               sxx+=x*x;
00584 //              LogTrace("MuonIdentification") << " FIT: x=" << x << " y= " << y << std::endl;
00585             }
00586 
00587             double d = s*sxx-sx*sx;
00588             
00589 //            muonTime.freeTime = (sxx*sy- sx*sxy)/d;
00590 //            muonTime.freeTimeErr = sqrt(s/d);
00591             muonTime.freeInverseBeta = (s*sxy - sx*sy)/d;         
00592             muonTime.freeInverseBetaErr = sqrt(sxx/d);    
00593           }
00594         } 
00595 
00596         LogTrace("MuonIdentification") << "Global 1/beta: " << muonTime.inverseBeta << " +/- " << muonTime.inverseBetaErr
00597                                        << "  # of points: " << muonTime.nStations <<std::endl;
00598         LogTrace("MuonIdentification") << "  Free 1/beta: " << muonTime.freeInverseBeta << " +/- " << muonTime.freeInverseBetaErr<<std::endl;
00599         LogTrace("MuonIdentification") << "  Vertex time (in-out): " << muonTime.timeAtIpInOut << " +/- " << muonTime.timeAtIpInOutErr<<std::endl;
00600         LogTrace("MuonIdentification") << "  Vertex time (out-in): " << muonTime.timeAtIpOutIn << " +/- " << muonTime.timeAtIpOutInErr<<std::endl;
00601         LogTrace("MuonIdentification") << "  direction: "   << muonTime.direction() << std::endl;
00602                                        
00603         muon.setTime(muonTime);                                       
00604 }
00605 
00606 // End Timing part      
00607 
00608 bool MuonIdProducer::isGoodTrackerMuon( const reco::Muon& muon )
00609 {
00610    if ( addExtraSoftMuons_ && 
00611         muon.pt()<5 && fabs(muon.eta())<1.5 && 
00612         muon.numberOfMatches( reco::Muon::NoArbitration ) >= 1 ) return true;
00613    return ( muon.numberOfMatches( reco::Muon::NoArbitration ) >= minNumberOfMatches_ );
00614 }
00615 
00616 void MuonIdProducer::fillMuonId(edm::Event& iEvent, const edm::EventSetup& iSetup,
00617                                 reco::Muon& aMuon, 
00618                                 TrackDetectorAssociator::Direction direction)
00619 {
00620    // perform track - detector association
00621    const reco::Track* track = 0;
00622    if ( ! aMuon.track().isNull() )
00623      track = aMuon.track().get();
00624    else 
00625      {
00626         if ( ! aMuon.standAloneMuon().isNull() )
00627           track = aMuon.standAloneMuon().get();
00628         else
00629           throw cms::Exception("FatalError") << "Failed to fill muon id information for a muon with undefined references to tracks"; 
00630      }
00631    
00632    TrackDetMatchInfo info = trackAssociator_.associate(iEvent, iSetup, *track, parameters_, direction);
00633    
00634    if ( fillEnergy_ ) {
00635       reco::MuonEnergy muonEnergy;
00636       muonEnergy.em  = info.crossedEnergy(TrackDetMatchInfo::EcalRecHits);
00637       muonEnergy.had = info.crossedEnergy(TrackDetMatchInfo::HcalRecHits);
00638       muonEnergy.ho  = info.crossedEnergy(TrackDetMatchInfo::HORecHits);
00639       muonEnergy.emS9  = info.nXnEnergy(TrackDetMatchInfo::EcalRecHits,1); // 3x3 energy
00640       muonEnergy.hadS9 = info.nXnEnergy(TrackDetMatchInfo::HcalRecHits,1); // 3x3 energy
00641       muonEnergy.hoS9  = info.nXnEnergy(TrackDetMatchInfo::HORecHits,1);   // 3x3 energy
00642       aMuon.setCalEnergy( muonEnergy );
00643    }
00644    if ( ! fillMatching_ ) return;
00645    
00646    // fill muon match info
00647    std::vector<reco::MuonChamberMatch> muonChamberMatches;
00648    unsigned int nubmerOfMatchesAccordingToTrackAssociator = 0;
00649    for( std::vector<TAMuonChamberMatch>::const_iterator chamber=info.chambers.begin();
00650         chamber!=info.chambers.end(); chamber++ )
00651      {
00652         reco::MuonChamberMatch matchedChamber;
00653         
00654         LocalError localError = chamber->tState.localError().positionError();
00655         matchedChamber.x = chamber->tState.localPosition().x();
00656         matchedChamber.y = chamber->tState.localPosition().y();
00657         matchedChamber.xErr = sqrt( localError.xx() );
00658         matchedChamber.yErr = sqrt( localError.yy() );
00659                                                                                                                                                             
00660         matchedChamber.dXdZ = chamber->tState.localDirection().z()!=0?chamber->tState.localDirection().x()/chamber->tState.localDirection().z():9999;
00661         matchedChamber.dYdZ = chamber->tState.localDirection().z()!=0?chamber->tState.localDirection().y()/chamber->tState.localDirection().z():9999;
00662         // DANGEROUS - compiler cannot guaranty parameters ordering
00663         AlgebraicSymMatrix55 trajectoryCovMatrix = chamber->tState.localError().matrix();
00664         matchedChamber.dXdZErr = trajectoryCovMatrix(1,1)>0?sqrt(trajectoryCovMatrix(1,1)):0;
00665         matchedChamber.dYdZErr = trajectoryCovMatrix(2,2)>0?sqrt(trajectoryCovMatrix(2,2)):0;
00666         
00667         matchedChamber.edgeX = chamber->localDistanceX;
00668         matchedChamber.edgeY = chamber->localDistanceY;
00669         
00670         matchedChamber.id = chamber->id;
00671         if ( ! chamber->segments.empty() ) ++nubmerOfMatchesAccordingToTrackAssociator;
00672         
00673         // fill segments
00674         for( std::vector<TAMuonSegmentMatch>::const_iterator segment = chamber->segments.begin();
00675              segment != chamber->segments.end(); segment++ ) 
00676           {
00677              reco::MuonSegmentMatch matchedSegment;
00678              matchedSegment.x = segment->hasPhi?segment->segmentLocalPosition.x():0;
00679              matchedSegment.y = segment->hasZed?segment->segmentLocalPosition.y():0;
00680              matchedSegment.dXdZ = segment->segmentLocalDirection.z()?segment->segmentLocalDirection.x()/segment->segmentLocalDirection.z():0;
00681              matchedSegment.dYdZ = segment->segmentLocalDirection.z()?segment->segmentLocalDirection.y()/segment->segmentLocalDirection.z():0;
00682              matchedSegment.xErr = segment->hasPhi&&segment->segmentLocalErrorXX>0?sqrt(segment->segmentLocalErrorXX):0;
00683              matchedSegment.yErr = segment->hasZed&&segment->segmentLocalErrorYY>0?sqrt(segment->segmentLocalErrorYY):0;
00684              matchedSegment.dXdZErr = segment->segmentLocalErrorDxDz>0?sqrt(segment->segmentLocalErrorDxDz):0;
00685              matchedSegment.dYdZErr = segment->segmentLocalErrorDyDz>0?sqrt(segment->segmentLocalErrorDyDz):0;
00686              matchedSegment.t0 = segment->t0;
00687              matchedSegment.mask = 0;
00688              // test segment
00689              bool matchedX = false;
00690              bool matchedY = false;
00691              LogTrace("MuonIdentification") << " matching local x, segment x: " << matchedSegment.x << 
00692                ", chamber x: " << matchedChamber.x << ", max: " << maxAbsDx_;
00693              LogTrace("MuonIdentification") << " matching local y, segment y: " << matchedSegment.y <<
00694                ", chamber y: " << matchedChamber.y << ", max: " << maxAbsDy_;
00695              if (matchedSegment.xErr>0 && matchedChamber.xErr>0 )
00696                LogTrace("MuonIdentification") << " xpull: " << 
00697                fabs(matchedSegment.x - matchedChamber.x)/sqrt(pow(matchedSegment.xErr,2) + pow(matchedChamber.xErr,2));
00698              if (matchedSegment.yErr>0 && matchedChamber.yErr>0 )
00699                LogTrace("MuonIdentification") << " ypull: " << 
00700                fabs(matchedSegment.y - matchedChamber.y)/sqrt(pow(matchedSegment.yErr,2) + pow(matchedChamber.yErr,2));
00701              
00702              if (fabs(matchedSegment.x - matchedChamber.x) < maxAbsDx_) matchedX = true;
00703              if (fabs(matchedSegment.y - matchedChamber.y) < maxAbsDy_) matchedY = true;
00704              if (matchedSegment.xErr>0 && matchedChamber.xErr>0 && 
00705                  fabs(matchedSegment.x - matchedChamber.x)/sqrt(pow(matchedSegment.xErr,2) + pow(matchedChamber.xErr,2)) < maxAbsPullX_) matchedX = true;
00706              if (matchedSegment.yErr>0 && matchedChamber.yErr>0 && 
00707                  fabs(matchedSegment.y - matchedChamber.y)/sqrt(pow(matchedSegment.yErr,2) + pow(matchedChamber.yErr,2)) < maxAbsPullY_) matchedY = true;
00708              if (matchedX && matchedY) matchedChamber.segmentMatches.push_back(matchedSegment);
00709           }
00710         muonChamberMatches.push_back(matchedChamber);
00711      }
00712    aMuon.setMatches(muonChamberMatches);
00713 
00714    LogTrace("MuonIdentification") << "number of muon chambers: " << aMuon.matches().size() << "\n" 
00715      << "number of chambers with segments according to the associator requirements: " << 
00716      nubmerOfMatchesAccordingToTrackAssociator;
00717    LogTrace("MuonIdentification") << "number of segment matches with the producer requirements: " << 
00718      aMuon.numberOfMatches( reco::Muon::NoArbitration );
00719    
00720    // fillTime( iEvent, iSetup, aMuon );
00721 }
00722 
00723 void MuonIdProducer::fillArbitrationInfo( reco::MuonCollection* pOutputMuons )
00724 {
00725    //
00726    // apply segment flags
00727    //
00728    std::vector<std::pair<reco::MuonChamberMatch*,reco::MuonSegmentMatch*> > chamberPairs;     // for chamber segment sorting
00729    std::vector<std::pair<reco::MuonChamberMatch*,reco::MuonSegmentMatch*> > stationPairs;     // for station segment sorting
00730    std::vector<std::pair<reco::MuonChamberMatch*,reco::MuonSegmentMatch*> > arbitrationPairs; // for muon segment arbitration
00731 
00732    // muonIndex1
00733    for( unsigned int muonIndex1 = 0; muonIndex1 < pOutputMuons->size(); ++muonIndex1 )
00734    {
00735       // chamberIter1
00736       for( std::vector<reco::MuonChamberMatch>::iterator chamberIter1 = pOutputMuons->at(muonIndex1).matches().begin();
00737             chamberIter1 != pOutputMuons->at(muonIndex1).matches().end(); ++chamberIter1 )
00738       {
00739          if(chamberIter1->segmentMatches.empty()) continue;
00740          chamberPairs.clear();
00741 
00742          // segmentIter1
00743          for( std::vector<reco::MuonSegmentMatch>::iterator segmentIter1 = chamberIter1->segmentMatches.begin();
00744                segmentIter1 != chamberIter1->segmentMatches.end(); ++segmentIter1 )
00745          {
00746             chamberPairs.push_back(std::make_pair(&(*chamberIter1), &(*segmentIter1)));
00747             if(!segmentIter1->isMask()) // has not yet been arbitrated
00748             {
00749                arbitrationPairs.clear();
00750                arbitrationPairs.push_back(std::make_pair(&(*chamberIter1), &(*segmentIter1)));
00751 
00752                // find identical segments with which to arbitrate
00753                // muonIndex2
00754                for( unsigned int muonIndex2 = muonIndex1+1; muonIndex2 < pOutputMuons->size(); ++muonIndex2 )
00755                {
00756                   // chamberIter2
00757                   for( std::vector<reco::MuonChamberMatch>::iterator chamberIter2 = pOutputMuons->at(muonIndex2).matches().begin();
00758                         chamberIter2 != pOutputMuons->at(muonIndex2).matches().end(); ++chamberIter2 )
00759                   {
00760                      // segmentIter2
00761                      for( std::vector<reco::MuonSegmentMatch>::iterator segmentIter2 = chamberIter2->segmentMatches.begin();
00762                            segmentIter2 != chamberIter2->segmentMatches.end(); ++segmentIter2 )
00763                      {
00764                         if(segmentIter2->isMask()) continue; // has already been arbitrated
00765                         if(fabs(segmentIter2->x       - segmentIter1->x      ) < 1E-3 &&
00766                            fabs(segmentIter2->y       - segmentIter1->y      ) < 1E-3 &&
00767                            fabs(segmentIter2->dXdZ    - segmentIter1->dXdZ   ) < 1E-3 &&
00768                            fabs(segmentIter2->dYdZ    - segmentIter1->dYdZ   ) < 1E-3 &&
00769                            fabs(segmentIter2->xErr    - segmentIter1->xErr   ) < 1E-3 &&
00770                            fabs(segmentIter2->yErr    - segmentIter1->yErr   ) < 1E-3 &&
00771                            fabs(segmentIter2->dXdZErr - segmentIter1->dXdZErr) < 1E-3 &&
00772                            fabs(segmentIter2->dYdZErr - segmentIter1->dYdZErr) < 1E-3)
00773                            arbitrationPairs.push_back(std::make_pair(&(*chamberIter2), &(*segmentIter2)));
00774                      } // segmentIter2
00775                   } // chamberIter2
00776                } // muonIndex2
00777 
00778                // arbitration segment sort
00779                if(arbitrationPairs.empty()) continue; // this should never happen
00780                if(arbitrationPairs.size()==1) {
00781                   arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDRSlope);
00782                   arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDXSlope);
00783                   arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDR);
00784                   arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDX);
00785                   arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::Arbitrated);
00786                } else {
00787                   sort(arbitrationPairs.begin(), arbitrationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDRSlope));
00788                   arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDRSlope);
00789                   sort(arbitrationPairs.begin(), arbitrationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDXSlope));
00790                   arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDXSlope);
00791                   sort(arbitrationPairs.begin(), arbitrationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDR));
00792                   arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDR);
00793                   sort(arbitrationPairs.begin(), arbitrationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BelongsToTrackByDX));
00794                   arbitrationPairs.front().second->setMask(reco::MuonSegmentMatch::BelongsToTrackByDX);
00795                   for( unsigned int it = 0; it < arbitrationPairs.size(); ++it )
00796                      arbitrationPairs.at(it).second->setMask(reco::MuonSegmentMatch::Arbitrated);
00797                }
00798             }
00799          } // segmentIter1
00800 
00801          // chamber segment sort
00802          if(chamberPairs.empty()) continue; // this should never happen
00803          if(chamberPairs.size()==1) {
00804             chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDRSlope);
00805             chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDXSlope);
00806             chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDR);
00807             chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDX);
00808          } else {
00809             sort(chamberPairs.begin(), chamberPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDRSlope));
00810             chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDRSlope);
00811             sort(chamberPairs.begin(), chamberPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDXSlope));
00812             chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDXSlope);
00813             sort(chamberPairs.begin(), chamberPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDR));
00814             chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDR);
00815             sort(chamberPairs.begin(), chamberPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInChamberByDX));
00816             chamberPairs.front().second->setMask(reco::MuonSegmentMatch::BestInChamberByDX);
00817          }
00818       } // chamberIter1
00819 
00820       // station segment sort
00821       for( int stationIndex = 1; stationIndex < 5; ++stationIndex )
00822          for( int detectorIndex = 1; detectorIndex < 4; ++detectorIndex )
00823          {
00824             stationPairs.clear();
00825 
00826             // chamberIter
00827             for( std::vector<reco::MuonChamberMatch>::iterator chamberIter = pOutputMuons->at(muonIndex1).matches().begin();
00828                   chamberIter != pOutputMuons->at(muonIndex1).matches().end(); ++chamberIter )
00829             {
00830                if(!(chamberIter->station()==stationIndex && chamberIter->detector()==detectorIndex)) continue;
00831                if(chamberIter->segmentMatches.empty()) continue;
00832 
00833                for( std::vector<reco::MuonSegmentMatch>::iterator segmentIter = chamberIter->segmentMatches.begin();
00834                      segmentIter != chamberIter->segmentMatches.end(); ++segmentIter )
00835                   stationPairs.push_back(std::make_pair(&(*chamberIter), &(*segmentIter)));
00836             } // chamberIter
00837 
00838             if(stationPairs.empty()) continue; // this may very well happen
00839             if(stationPairs.size()==1) {
00840                stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDRSlope);
00841                stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDXSlope);
00842                stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDR);
00843                stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDX);
00844             } else {
00845                sort(stationPairs.begin(), stationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDRSlope));
00846                stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDRSlope);
00847                sort(stationPairs.begin(), stationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDXSlope));
00848                stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDXSlope);
00849                sort(stationPairs.begin(), stationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDR));
00850                stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDR);
00851                sort(stationPairs.begin(), stationPairs.end(), SortMuonSegmentMatches(reco::MuonSegmentMatch::BestInStationByDX));
00852                stationPairs.front().second->setMask(reco::MuonSegmentMatch::BestInStationByDX);
00853             }
00854          }
00855 
00856    } // muonIndex1
00857 }
00858 
00859 void MuonIdProducer::fillMuonIsolation(edm::Event& iEvent, const edm::EventSetup& iSetup, reco::Muon& aMuon)
00860 {
00861    reco::MuonIsolation isoR03, isoR05;
00862    const reco::Track* track = 0;
00863    if ( ! aMuon.track().isNull() )
00864      track = aMuon.track().get();
00865    else 
00866      {
00867         if ( ! aMuon.standAloneMuon().isNull() )
00868           track = aMuon.standAloneMuon().get();
00869         else
00870           throw cms::Exception("FatalError") << "Failed to compute muon isolation information for a muon with undefined references to tracks"; 
00871      }
00872 
00873    // get deposits
00874    reco::IsoDeposit depTrk = muIsoExtractorTrack_->deposit(iEvent, iSetup, *track );
00875    std::vector<reco::IsoDeposit> caloDeps = muIsoExtractorCalo_->deposits(iEvent, iSetup, *track);
00876    reco::IsoDeposit depJet = muIsoExtractorJet_->deposit(iEvent, iSetup, *track );
00877 
00878    if(caloDeps.size()!=3) {
00879       LogTrace("MuonIdentification") << "Failed to fill vector of calorimeter isolation deposits!";
00880       return;
00881    }
00882 
00883    reco::IsoDeposit depEcal = caloDeps.at(0);
00884    reco::IsoDeposit depHcal = caloDeps.at(1);
00885    reco::IsoDeposit depHo   = caloDeps.at(2);
00886 
00887    isoR03.sumPt     = depTrk.depositWithin(0.3);
00888    isoR03.emEt      = depEcal.depositWithin(0.3);
00889    isoR03.hadEt     = depHcal.depositWithin(0.3);
00890    isoR03.hoEt      = depHo.depositWithin(0.3);
00891    isoR03.nTracks   = depTrk.depositAndCountWithin(0.3).second;
00892    isoR03.nJets     = depJet.depositAndCountWithin(0.3).second;
00893 
00894    isoR05.sumPt     = depTrk.depositWithin(0.5);
00895    isoR05.emEt      = depEcal.depositWithin(0.5);
00896    isoR05.hadEt     = depHcal.depositWithin(0.5);
00897    isoR05.hoEt      = depHo.depositWithin(0.5);
00898    isoR05.nTracks   = depTrk.depositAndCountWithin(0.5).second;
00899    isoR05.nJets     = depJet.depositAndCountWithin(0.5).second;
00900 
00901    aMuon.setIsolation(isoR03, isoR05);
00902 }
00903 
00904 reco::Muon MuonIdProducer::makeMuon( const reco::Track& track )
00905 {
00906    //FIXME: E = sqrt(p^2 + m^2), where m == 0.105658369(9)GeV 
00907    double energy = sqrt(track.p() * track.p() + 0.011163691);
00908    math::XYZTLorentzVector p4(track.px(),
00909                               track.py(),
00910                               track.pz(),
00911                               energy);
00912    return reco::Muon( track.charge(), p4, track.vertex() );
00913 }
00914 
00915 double MuonIdProducer::sectorPhi( const DetId& id )
00916 {
00917    double phi = 0;
00918    if( id.subdetId() ==  MuonSubdetId::DT ) {    // DT
00919       DTChamberId muonId(id.rawId());
00920       if ( muonId.sector() <= 12 )
00921         phi = (muonId.sector()-1)/6.*M_PI;
00922       if ( muonId.sector() == 13 ) phi = 3/6.*M_PI;
00923       if ( muonId.sector() == 14 ) phi = 9/6.*M_PI;
00924    }
00925    if( id.subdetId() == MuonSubdetId::CSC ) {    // CSC
00926       CSCDetId muonId(id.rawId());
00927       phi = M_PI/4+(muonId.triggerSector()-1)/3.*M_PI;
00928    }
00929    if ( phi > M_PI ) phi -= 2*M_PI; 
00930    return phi; 
00931 }
00932 
00933 double MuonIdProducer::phiOfMuonIneteractionRegion( const reco::Muon& muon ) const
00934 {
00935    if ( muon.isStandAloneMuon() ) return muon.standAloneMuon()->innerPosition().phi();
00936    // the rest is tracker muon only
00937    return sectorPhi(muon.matches().at(0).id);
00938 }
00939