File indexing completed on 2021-09-16 03:23:10
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021 #include <memory>
0022 #include <map>
0023 #include <fstream>
0024
0025
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
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
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
0076
0077
0078
0079
0080
0081
0082
0083
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
0112
0113
0114
0115 void AlignmentMuonHIPTrajectorySelector::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0116
0117 edm::Handle<TrajTrackAssociationCollection> originalTrajTrackMap;
0118 iEvent.getByLabel(m_input, originalTrajTrackMap);
0119
0120
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 }
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 }
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 }
0202 }
0203 }
0204
0205 if (tracker_forwardredchi2 < m_maxTrackerForwardRedChi2 && tracker_dof >= m_minTrackerDOF && !has_bad_residual) {
0206 newTrajTrackMap->insert((*iPair).key, (*iPair).val);
0207 }
0208 }
0209 }
0210
0211
0212 iEvent.put(std::move(newTrajTrackMap));
0213 }
0214
0215
0216 DEFINE_FWK_MODULE(AlignmentMuonHIPTrajectorySelector);