File indexing completed on 2024-04-06 11:56:11
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorPluginFactory.h"
0012 #include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorBase.h"
0013 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFromTrack.h"
0014
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/Utilities/interface/InputTag.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/Framework/interface/EventSetup.h"
0019 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0020 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0021
0022 #include <sstream>
0023 #include "TProfile.h"
0024 #include "TH2F.h"
0025 #include "TH1F.h"
0026
0027 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0028 #include "TrackingTools/Records/interface/DetIdAssociatorRecord.h"
0029 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0030 #include "TrackingTools/TrackAssociator/interface/DetIdAssociator.h"
0031 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0032 #include "MagneticField/Engine/interface/MagneticField.h"
0033
0034 class AlignmentMonitorMuonVsCurvature : public AlignmentMonitorBase {
0035 public:
0036 AlignmentMonitorMuonVsCurvature(const edm::ParameterSet &cfg, edm::ConsumesCollector iC);
0037 ~AlignmentMonitorMuonVsCurvature() override {}
0038
0039 void book() override;
0040
0041 void event(const edm::Event &iEvent,
0042 const edm::EventSetup &iSetup,
0043 const ConstTrajTrackPairCollection &iTrajTracks) override;
0044 void processMuonResidualsFromTrack(MuonResidualsFromTrack &mrft, const Trajectory *traj = nullptr);
0045
0046 private:
0047
0048 const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> m_esTokenGBTGeom;
0049 const edm::ESGetToken<DetIdAssociator, DetIdAssociatorRecord> m_esTokenDetId;
0050 const edm::ESGetToken<Propagator, TrackingComponentsRecord> m_esTokenProp;
0051 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> m_esTokenMF;
0052 const MuonResidualsFromTrack::BuilderToken m_esTokenBuilder;
0053
0054
0055 const edm::InputTag m_muonCollectionTag;
0056 const double m_minTrackPt;
0057 const double m_minTrackP;
0058 const int m_minTrackerHits;
0059 const double m_maxTrackerRedChi2;
0060 const bool m_allowTIDTEC;
0061 const bool m_minNCrossedChambers;
0062 const double m_maxDxy;
0063 const int m_minDT13Hits;
0064 const int m_minDT2Hits;
0065 const int m_minCSCHits;
0066 const int m_layer;
0067 const std::string m_propagator;
0068 const bool m_doDT;
0069 const bool m_doCSC;
0070
0071 const edm::EDGetTokenT<reco::BeamSpot> bsToken_;
0072 const edm::EDGetTokenT<reco::MuonCollection> muonToken_;
0073
0074 enum { kDeltaX = 0, kDeltaDxDz, kNumComponents };
0075
0076
0077 TH2F *th2f_wheel_st_sector[5][4][14][kNumComponents];
0078 TProfile *tprofile_wheel_st_sector[5][4][14][kNumComponents];
0079
0080
0081 TH2F *th2f_st_ring_chamber[8][3][36][kNumComponents];
0082 TProfile *tprofile_st_ring_chamber[8][3][36][kNumComponents];
0083
0084 TH1F *th1f_trackerRedChi2;
0085 TH1F *th1f_trackerRedChi2Diff;
0086 };
0087
0088 AlignmentMonitorMuonVsCurvature::AlignmentMonitorMuonVsCurvature(const edm::ParameterSet &cfg,
0089 edm::ConsumesCollector iC)
0090 : AlignmentMonitorBase(cfg, iC, "AlignmentMonitorMuonVsCurvature"),
0091 m_esTokenGBTGeom(iC.esConsumes()),
0092 m_esTokenDetId(iC.esConsumes(edm::ESInputTag("", "MuonDetIdAssociator"))),
0093 m_esTokenProp(iC.esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))),
0094 m_esTokenMF(iC.esConsumes()),
0095 m_esTokenBuilder(iC.esConsumes(MuonResidualsFromTrack::builderESInputTag())),
0096 m_muonCollectionTag(cfg.getParameter<edm::InputTag>("muonCollectionTag")),
0097 m_minTrackPt(cfg.getParameter<double>("minTrackPt")),
0098 m_minTrackP(cfg.getParameter<double>("minTrackP")),
0099 m_minTrackerHits(cfg.getParameter<int>("minTrackerHits")),
0100 m_maxTrackerRedChi2(cfg.getParameter<double>("maxTrackerRedChi2")),
0101 m_allowTIDTEC(cfg.getParameter<bool>("allowTIDTEC")),
0102 m_minNCrossedChambers(cfg.getParameter<int>("minNCrossedChambers")),
0103 m_maxDxy(cfg.getParameter<double>("maxDxy")),
0104 m_minDT13Hits(cfg.getParameter<int>("minDT13Hits")),
0105 m_minDT2Hits(cfg.getParameter<int>("minDT2Hits")),
0106 m_minCSCHits(cfg.getParameter<int>("minCSCHits")),
0107 m_layer(cfg.getParameter<int>("layer")),
0108 m_propagator(cfg.getParameter<std::string>("propagator")),
0109 m_doDT(cfg.getParameter<bool>("doDT")),
0110 m_doCSC(cfg.getParameter<bool>("doCSC")),
0111 bsToken_(iC.consumes<reco::BeamSpot>(m_beamSpotTag)),
0112 muonToken_(iC.consumes<reco::MuonCollection>(m_muonCollectionTag)) {}
0113
0114 void AlignmentMonitorMuonVsCurvature::book() {
0115
0116 std::string wheelname[5] = {"wheelm2_", "wheelm1_", "wheelz_", "wheelp1_", "wheelp2_"};
0117 if (m_doDT)
0118 for (int wheel = -2; wheel <= 2; wheel++)
0119 for (int station = 1; station <= 4; station++)
0120 for (int sector = 1; sector <= 14; sector++) {
0121 if (station != 4 && sector > 12)
0122 continue;
0123
0124 char stationname[20];
0125 sprintf(stationname, "st%d_", station);
0126
0127 char sectorname[20];
0128 sprintf(sectorname, "sector%02d_", sector);
0129
0130 for (int component = 0; component < kNumComponents; component++) {
0131 std::stringstream th2f_name, tprofile_name;
0132 th2f_name << "th2f_" << wheelname[wheel + 2] << stationname << sectorname;
0133 tprofile_name << "tprofile_" << wheelname[wheel + 2] << stationname << sectorname;
0134
0135 double yminmax = 50., xminmax = 0.05;
0136 if (m_minTrackPt > 0.)
0137 xminmax = 1. / m_minTrackPt;
0138 int ynbins = 50;
0139 if (component == kDeltaX) {
0140 th2f_name << "deltax";
0141 tprofile_name << "deltax";
0142 } else if (component == kDeltaDxDz) {
0143 th2f_name << "deltadxdz";
0144 tprofile_name << "deltadxdz";
0145 }
0146
0147 th2f_wheel_st_sector[wheel + 2][station - 1][sector - 1][component] =
0148 book2D("/iterN/", th2f_name.str(), "", 30, -xminmax, xminmax, ynbins, -yminmax, yminmax);
0149 tprofile_wheel_st_sector[wheel + 2][station - 1][sector - 1][component] =
0150 bookProfile("/iterN/", tprofile_name.str(), "", 30, -xminmax, xminmax);
0151 }
0152 }
0153
0154
0155 std::string stname[8] = {"Ep_S1_", "Ep_S2_", "Ep_S3_", "Ep_S4_", "Em_S1_", "Em_S2_", "Em_S3_", "Em_S4_"};
0156 if (m_doCSC)
0157 for (int station = 0; station < 8; station++)
0158 for (int ring = 1; ring <= 3; ring++)
0159 for (int chamber = 1; chamber <= 36; chamber++) {
0160 int st = station % 4 + 1;
0161 if (st > 1 && ring > 2)
0162 continue;
0163 if (st > 1 && ring == 1 && chamber > 18)
0164 continue;
0165
0166 char ringname[20];
0167 sprintf(ringname, "R%d_", ring);
0168
0169 char chname[20];
0170 sprintf(chname, "C%02d_", chamber);
0171
0172 for (int component = 0; component < kNumComponents; component++) {
0173 std::stringstream componentname;
0174 double yminmax = 50., xminmax = 0.05;
0175 if (m_minTrackPt > 0.)
0176 xminmax = 1. / m_minTrackPt;
0177 if (ring == 1)
0178 xminmax *= 0.5;
0179 if (component == kDeltaX) {
0180 componentname << "deltax";
0181 } else if (component == kDeltaDxDz) {
0182 componentname << "deltadxdz";
0183 }
0184
0185 std::stringstream th2f_name, tprofile_name;
0186 th2f_name << "th2f_" << stname[station] << ringname << chname << componentname.str();
0187 tprofile_name << "tprofile_" << stname[station] << ringname << chname << componentname.str();
0188
0189 th2f_st_ring_chamber[station][ring - 1][chamber - 1][component] =
0190 book2D("/iterN/", th2f_name.str(), "", 30, -xminmax, xminmax, 100, -yminmax, yminmax);
0191 tprofile_st_ring_chamber[station][ring - 1][chamber - 1][component] =
0192 bookProfile("/iterN/", tprofile_name.str(), "", 30, -xminmax, xminmax);
0193 }
0194 }
0195
0196 th1f_trackerRedChi2 = book1D("/iterN/", "trackerRedChi2", "Refit tracker reduced chi^2", 100, 0., 30.);
0197 th1f_trackerRedChi2Diff =
0198 book1D("/iterN/", "trackerRedChi2Diff", "Fit-minus-refit tracker reduced chi^2", 100, -5., 5.);
0199 }
0200
0201 void AlignmentMonitorMuonVsCurvature::event(const edm::Event &iEvent,
0202 const edm::EventSetup &iSetup,
0203 const ConstTrajTrackPairCollection &trajtracks) {
0204 const edm::Handle<reco::BeamSpot> &beamSpot = iEvent.getHandle(bsToken_);
0205
0206 const GlobalTrackingGeometry *globalGeometry = &iSetup.getData(m_esTokenGBTGeom);
0207 const DetIdAssociator *muonDetIdAssociator_ = &iSetup.getData(m_esTokenDetId);
0208 const Propagator *prop = &iSetup.getData(m_esTokenProp);
0209 const MagneticField *magneticField = &iSetup.getData(m_esTokenMF);
0210 auto builder = iSetup.getHandle(m_esTokenBuilder);
0211
0212 if (m_muonCollectionTag.label().empty())
0213 {
0214 for (ConstTrajTrackPairCollection::const_iterator trajtrack = trajtracks.begin(); trajtrack != trajtracks.end();
0215 ++trajtrack) {
0216 const Trajectory *traj = (*trajtrack).first;
0217 const reco::Track *track = (*trajtrack).second;
0218
0219 if (track->pt() > m_minTrackPt && track->p() > m_minTrackP && fabs(track->dxy(beamSpot->position())) < m_maxDxy) {
0220 MuonResidualsFromTrack muonResidualsFromTrack(
0221 builder, magneticField, globalGeometry, muonDetIdAssociator_, prop, traj, track, pNavigator(), 1000.);
0222 processMuonResidualsFromTrack(muonResidualsFromTrack, traj);
0223 }
0224 }
0225 } else {
0226 const edm::Handle<reco::MuonCollection> &muons = iEvent.getHandle(muonToken_);
0227
0228 for (reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
0229 if (!(muon->isTrackerMuon() && muon->innerTrack().isNonnull()))
0230 continue;
0231
0232 if (m_minTrackPt < muon->pt() && m_minTrackP < muon->p() &&
0233 fabs(muon->innerTrack()->dxy(beamSpot->position())) < m_maxDxy) {
0234 MuonResidualsFromTrack muonResidualsFromTrack(globalGeometry, &(*muon), pNavigator(), 100.);
0235 processMuonResidualsFromTrack(muonResidualsFromTrack);
0236 }
0237 }
0238 }
0239 }
0240
0241 void AlignmentMonitorMuonVsCurvature::processMuonResidualsFromTrack(MuonResidualsFromTrack &mrft,
0242 const Trajectory *traj) {
0243 if (mrft.trackerNumHits() < m_minTrackerHits)
0244 return;
0245 if (!m_allowTIDTEC && mrft.contains_TIDTEC())
0246 return;
0247
0248 int nMuChambers = 0;
0249 std::vector<DetId> chamberIds = mrft.chamberIds();
0250 for (unsigned ch = 0; ch < chamberIds.size(); ch++)
0251 if (chamberIds[ch].det() == DetId::Muon)
0252 nMuChambers++;
0253 if (nMuChambers < m_minNCrossedChambers)
0254 return;
0255
0256 th1f_trackerRedChi2->Fill(mrft.trackerRedChi2());
0257 th1f_trackerRedChi2Diff->Fill(mrft.getTrack()->normalizedChi2() - mrft.trackerRedChi2());
0258
0259 if (mrft.normalizedChi2() > m_maxTrackerRedChi2)
0260 return;
0261
0262 double qoverpt = mrft.getTrack()->charge() / mrft.getTrack()->pt();
0263 double qoverpz = 0.;
0264 if (fabs(mrft.getTrack()->pz()) > 0.01)
0265 qoverpz = mrft.getTrack()->charge() / fabs(mrft.getTrack()->pz());
0266
0267 for (std::vector<DetId>::const_iterator chamberId = chamberIds.begin(); chamberId != chamberIds.end(); ++chamberId) {
0268 if (chamberId->det() != DetId::Muon)
0269 continue;
0270
0271 if (m_doDT && chamberId->subdetId() == MuonSubdetId::DT) {
0272 DTChamberId dtid(chamberId->rawId());
0273 MuonChamberResidual *dt13 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT13);
0274
0275 if (dt13 != nullptr && dt13->numHits() >= m_minDT13Hits) {
0276 int wheel = dtid.wheel() + 2;
0277 int station = dtid.station() - 1;
0278 int sector = dtid.sector() - 1;
0279
0280 double resid_x = 10. * dt13->global_residual();
0281 double resid_dxdz = 1000. * dt13->global_resslope();
0282
0283 if (fabs(resid_x) < 100. && fabs(resid_dxdz) < 100.) {
0284 th2f_wheel_st_sector[wheel][station][sector][kDeltaX]->Fill(qoverpt, resid_x);
0285 tprofile_wheel_st_sector[wheel][station][sector][kDeltaX]->Fill(qoverpt, resid_x);
0286 th2f_wheel_st_sector[wheel][station][sector][kDeltaDxDz]->Fill(qoverpt, resid_dxdz);
0287 tprofile_wheel_st_sector[wheel][station][sector][kDeltaDxDz]->Fill(qoverpt, resid_dxdz);
0288 }
0289 }
0290 }
0291
0292 if (m_doCSC && chamberId->subdetId() == MuonSubdetId::CSC) {
0293 CSCDetId cscid(chamberId->rawId());
0294 MuonChamberResidual *csc = mrft.chamberResidual(*chamberId, MuonChamberResidual::kCSC);
0295
0296 if (csc != nullptr && csc->numHits() >= m_minCSCHits) {
0297 int station = 4 * cscid.endcap() + cscid.station() - 5;
0298 int ring = cscid.ring() - 1;
0299 if (cscid.station() == 1 && cscid.ring() == 4)
0300 ring = 0;
0301 int chamber = cscid.chamber() - 1;
0302
0303 double resid_x = 10. * csc->global_residual();
0304 double resid_dxdz = 1000. * csc->global_resslope();
0305
0306 if (fabs(resid_x) < 100. && fabs(resid_dxdz) < 100.) {
0307 th2f_st_ring_chamber[station][ring][chamber][kDeltaX]->Fill(qoverpz, resid_x);
0308 tprofile_st_ring_chamber[station][ring][chamber][kDeltaX]->Fill(qoverpz, resid_x);
0309 th2f_st_ring_chamber[station][ring][chamber][kDeltaDxDz]->Fill(qoverpz, resid_dxdz);
0310 tprofile_st_ring_chamber[station][ring][chamber][kDeltaDxDz]->Fill(qoverpz, resid_dxdz);
0311 }
0312 }
0313 }
0314
0315 }
0316 }
0317
0318 DEFINE_EDM_PLUGIN(AlignmentMonitorPluginFactory, AlignmentMonitorMuonVsCurvature, "AlignmentMonitorMuonVsCurvature");