CommonAlignmentProducer/plugins/AlignmentMuonHIPTrajectorySelector.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    AlignmentMuonHIPTrajectorySelector
0004 // Class:      AlignmentMuonHIPTrajectorySelector
0005 //
0006 /**\class AlignmentMuonHIPTrajectorySelector AlignmentMuonHIPTrajectorySelector.cc Alignment/CommonAlignmentProducer/plugins/AlignmentMuonHIPTrajectorySelector.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Jim Pivarski
0015 //         Created:  Wed Feb 20 10:56:46 CST 2008
0016 // $Id: AlignmentMuonHIPTrajectorySelector.cc,v 1.4 2010/01/06 15:26:10 mussgill Exp $
0017 //
0018 //
0019
0020 // system include files
0021 #include <memory>
0022 #include <map>
0023 #include <fstream>
0024
0025 // user include files
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/stream/EDProducer.h"
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0032 #include "FWCore/Utilities/interface/InputTag.h"
0033 #include "FWCore/Framework/interface/Event.h"
0034 #include "FWCore/Framework/interface/EventSetup.h"
0035
0036 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0037 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0038 #include "DataFormats/TrackReco/interface/Track.h"
0039 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0040 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0041 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0042 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0043 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0044 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0045 #include "DataFormats/DetId/interface/DetId.h"
0046 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0047 #include "FWCore/ServiceRegistry/interface/Service.h"
0048 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0049 #include "TH1F.h"
0050
0051 //
0052 // class decleration
0053 //
0054
0055 class AlignmentMuonHIPTrajectorySelector : public edm::stream::EDProducer<> {
0056 public:
0057   explicit AlignmentMuonHIPTrajectorySelector(const edm::ParameterSet&);
0058   ~AlignmentMuonHIPTrajectorySelector() override = default;
0059
0060   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0061
0062 private:
0063   void produce(edm::Event&, const edm::EventSetup&) override;
0064
0065   // ---------- member data --------------------------------
0066   const edm::InputTag m_input;
0067   const double m_minPt;
0068   const double m_maxTrackerForwardRedChi2;
0069   const int m_minTrackerDOF;
0070   const double m_maxMuonResidual;
0071
0072   const bool m_hists;
0073   const edm::EDGetTokenT<TrajTrackAssociationCollection> mapToken_;
0074   TH1F *m_pt, *m_tracker_forwardredchi2, *m_tracker_dof;
0075   TH1F *m_resid_before, *m_resid_after;
0076 };
0077
0078 //
0079 // constants, enums and typedefs
0080 //
0081
0082 //
0083 // static data member definitions
0084 //
0085
0086 //
0087 // constructors and destructor
0088 //
0089 AlignmentMuonHIPTrajectorySelector::AlignmentMuonHIPTrajectorySelector(const edm::ParameterSet& iConfig)
0090     : m_input(iConfig.getParameter<edm::InputTag>("input")),
0091       m_minPt(iConfig.getParameter<double>("minPt")),
0092       m_maxTrackerForwardRedChi2(iConfig.getParameter<double>("maxTrackerForwardRedChi2")),
0093       m_minTrackerDOF(iConfig.getParameter<int>("minTrackerDOF")),
0094       m_maxMuonResidual(iConfig.getParameter<double>("maxMuonResidual")),
0095       m_hists(iConfig.getParameter<bool>("hists")),
0096       mapToken_(consumes<TrajTrackAssociationCollection>(m_input)),
0097       m_pt(nullptr),
0098       m_tracker_forwardredchi2(nullptr),
0099       m_tracker_dof(nullptr) {
0100   produces<TrajTrackAssociationCollection>();
0101
0102   if (m_hists) {
0103     edm::Service<TFileService> fs;
0104     m_pt = fs->make<TH1F>("pt", "Transverse momentum (GeV)", 100, 0., 100.);
0105     m_tracker_forwardredchi2 =
0106         fs->make<TH1F>("trackerForwardRedChi2", "forward-biased reduced chi2 in tracker", 100, 0., 5.);
0107     m_tracker_dof = fs->make<TH1F>("trackerDOF", "DOF in tracker", 61, -0.5, 60.5);
0108     m_resid_before = fs->make<TH1F>("residBefore", "muon residuals before cut (cm)", 100, -20, 20);
0109     m_resid_after = fs->make<TH1F>("residAfter", "muon residuals after cut (cm)", 100, -20, 20);
0110   }
0111 }
0112
0113 //
0114 // member functions
0115 //
0116 void AlignmentMuonHIPTrajectorySelector::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0117   edm::ParameterSetDescription desc;
0118   desc.setUnknown();
0119   descriptions.addDefault(desc);
0120 }
0121
0122 // ------------ method called to produce the data  ------------
0123 void AlignmentMuonHIPTrajectorySelector::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0124   // input
0125   const edm::Handle<TrajTrackAssociationCollection>& originalTrajTrackMap = iEvent.getHandle(mapToken_);
0126
0127   // output
0128   auto newTrajTrackMap = std::make_unique<TrajTrackAssociationCollection>();
0129
0130   TrajectoryStateCombiner tsoscomb;
0131
0132   for (TrajTrackAssociationCollection::const_iterator iPair = originalTrajTrackMap->begin();
0133        iPair != originalTrajTrackMap->end();
0134        ++iPair) {
0135     if (m_hists) {
0136       m_pt->Fill((*(*iPair).val).pt());
0137     }
0138
0139     if ((*(*iPair).val).pt() > m_minPt) {
0140       std::vector<TrajectoryMeasurement> measurements = (*(*iPair).key).measurements();
0141
0142       bool has_bad_residual = false;
0143
0144       double tracker_forwardchi2 = 0.;
0145       double tracker_dof = 0.;
0146       for (std::vector<TrajectoryMeasurement>::const_iterator im = measurements.begin(); im != measurements.end();
0147            ++im) {
0148         const TrajectoryMeasurement& meas = *im;
0149         auto hit = &(*meas.recHit());
0150         const DetId id = hit->geographicalId();
0151
0152         if (hit->isValid() && id.det() == DetId::Tracker) {
0153           if (hit->dimension() == 1) {
0154             double residual = meas.forwardPredictedState().localPosition().x() - hit->localPosition().x();
0155             double error2 =
0156                 meas.forwardPredictedState().localError().positionError().xx() + hit->localPositionError().xx();
0157
0158             tracker_forwardchi2 += residual * residual / error2;
0159             tracker_dof += 1.;
0160           } else if (hit->dimension() == 2) {
0161             double residualx = meas.forwardPredictedState().localPosition().x() - hit->localPosition().x();
0162             double residualy = meas.forwardPredictedState().localPosition().y() - hit->localPosition().y();
0163             double errorxx2 =
0164                 meas.forwardPredictedState().localError().positionError().xx() + hit->localPositionError().xx();
0165             double errorxy2 =
0166                 meas.forwardPredictedState().localError().positionError().xy() + hit->localPositionError().xy();
0167             double erroryy2 =
0168                 meas.forwardPredictedState().localError().positionError().yy() + hit->localPositionError().yy();
0169
0170             tracker_forwardchi2 +=
0171                 (residualx * residualx + residualy * residualy) / (errorxx2 + 2. * errorxy2 + erroryy2);
0172             tracker_dof += 2.;
0173           }
0174         }  // end if a tracker hit
0175
0176         if (hit->isValid() && id.det() == DetId::Muon &&
0177             (id.subdetId() == MuonSubdetId::DT || id.subdetId() == MuonSubdetId::CSC)) {
0178           TrajectoryStateOnSurface tsosc =
0179               tsoscomb.combine(meas.forwardPredictedState(), meas.backwardPredictedState());
0180           double residual = tsosc.localPosition().x() - hit->localPosition().x();
0181           m_resid_before->Fill(residual);
0182           if (fabs(residual) > m_maxMuonResidual) {
0183             has_bad_residual = true;
0184           }
0185         }  // end if a muon hit
0186       }
0187       tracker_dof -= 5.;
0188       double tracker_forwardredchi2 = tracker_forwardchi2 / tracker_dof;
0189
0190       if (m_hists) {
0191         m_tracker_forwardredchi2->Fill(tracker_forwardredchi2);
0192         m_tracker_dof->Fill(tracker_dof);
0193
0194         for (std::vector<TrajectoryMeasurement>::const_iterator im = measurements.begin(); im != measurements.end();
0195              ++im) {
0196           const TrajectoryMeasurement& meas = *im;
0197           auto hit = &(*meas.recHit());
0198           const DetId id = hit->geographicalId();
0199
0200           if (!has_bad_residual) {
0201             if (hit->isValid() && id.det() == DetId::Muon &&
0202                 (id.subdetId() == MuonSubdetId::DT || id.subdetId() == MuonSubdetId::CSC)) {
0203               TrajectoryStateOnSurface tsosc =
0204                   tsoscomb.combine(meas.forwardPredictedState(), meas.backwardPredictedState());
0205               double residual = tsosc.localPosition().x() - hit->localPosition().x();
0206               m_resid_after->Fill(residual);
0207             }
0208           }  // end if residuals pass cut
0209         }  // end second loop over hits
0210       }  // end if filling histograms
0211
0212       if (tracker_forwardredchi2 < m_maxTrackerForwardRedChi2 && tracker_dof >= m_minTrackerDOF && !has_bad_residual) {
0213         newTrajTrackMap->insert((*iPair).key, (*iPair).val);
0214       }  // end if passes tracker cuts
0215     }  // end if passes pT cut
0216   }  // end loop over original trajTrackMap
0217
0218   // put it in the Event
0219   iEvent.put(std::move(newTrajTrackMap));
0220 }
0221
0222 //define this as a plug-in
0223 DEFINE_FWK_MODULE(AlignmentMuonHIPTrajectorySelector);