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