Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-20 10:41:38

0001 /*
0002  * Package:     CommonAlignmentProducer
0003  * Class  :     AlignmentMonitorMuonVsCurvature
0004  *
0005  * Original Author:  Jim Pivarski
0006  *         Created:  Fri Feb 19 21:45:02 CET 2010
0007  *
0008  * $Id:$
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   // es token
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   // parameters
0055   edm::InputTag m_muonCollectionTag;
0056   double m_minTrackPt;
0057   double m_minTrackP;
0058   int m_minTrackerHits;
0059   double m_maxTrackerRedChi2;
0060   bool m_allowTIDTEC;
0061   bool m_minNCrossedChambers;
0062   double m_maxDxy;
0063   int m_minDT13Hits;
0064   int m_minDT2Hits;
0065   int m_minCSCHits;
0066   int m_layer;
0067   std::string m_propagator;
0068   bool m_doDT;
0069   bool m_doCSC;
0070 
0071   enum { kDeltaX = 0, kDeltaDxDz, kNumComponents };
0072 
0073   // DT [wheel][station][sector]
0074   TH2F *th2f_wheel_st_sector[5][4][14][kNumComponents];
0075   TProfile *tprofile_wheel_st_sector[5][4][14][kNumComponents];
0076 
0077   // CSC [endcap*station][ring][chamber]
0078   TH2F *th2f_st_ring_chamber[8][3][36][kNumComponents];
0079   TProfile *tprofile_st_ring_chamber[8][3][36][kNumComponents];
0080 
0081   TH1F *th1f_trackerRedChi2;
0082   TH1F *th1f_trackerRedChi2Diff;
0083 };
0084 
0085 AlignmentMonitorMuonVsCurvature::AlignmentMonitorMuonVsCurvature(const edm::ParameterSet &cfg,
0086                                                                  edm::ConsumesCollector iC)
0087     : AlignmentMonitorBase(cfg, iC, "AlignmentMonitorMuonVsCurvature"),
0088       m_esTokenGBTGeom(iC.esConsumes()),
0089       m_esTokenDetId(iC.esConsumes(edm::ESInputTag("", "MuonDetIdAssociator"))),
0090       m_esTokenProp(iC.esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))),
0091       m_esTokenMF(iC.esConsumes()),
0092       m_esTokenBuilder(iC.esConsumes(MuonResidualsFromTrack::builderESInputTag())),
0093       m_muonCollectionTag(cfg.getParameter<edm::InputTag>("muonCollectionTag")),
0094       m_minTrackPt(cfg.getParameter<double>("minTrackPt")),
0095       m_minTrackP(cfg.getParameter<double>("minTrackP")),
0096       m_minTrackerHits(cfg.getParameter<int>("minTrackerHits")),
0097       m_maxTrackerRedChi2(cfg.getParameter<double>("maxTrackerRedChi2")),
0098       m_allowTIDTEC(cfg.getParameter<bool>("allowTIDTEC")),
0099       m_minNCrossedChambers(cfg.getParameter<int>("minNCrossedChambers")),
0100       m_maxDxy(cfg.getParameter<double>("maxDxy")),
0101       m_minDT13Hits(cfg.getParameter<int>("minDT13Hits")),
0102       m_minDT2Hits(cfg.getParameter<int>("minDT2Hits")),
0103       m_minCSCHits(cfg.getParameter<int>("minCSCHits")),
0104       m_layer(cfg.getParameter<int>("layer")),
0105       m_propagator(cfg.getParameter<std::string>("propagator")),
0106       m_doDT(cfg.getParameter<bool>("doDT")),
0107       m_doCSC(cfg.getParameter<bool>("doCSC")) {}
0108 
0109 void AlignmentMonitorMuonVsCurvature::book() {
0110   // DT
0111   std::string wheelname[5] = {"wheelm2_", "wheelm1_", "wheelz_", "wheelp1_", "wheelp2_"};
0112   if (m_doDT)
0113     for (int wheel = -2; wheel <= 2; wheel++)
0114       for (int station = 1; station <= 4; station++)
0115         for (int sector = 1; sector <= 14; sector++) {
0116           if (station != 4 && sector > 12)
0117             continue;
0118 
0119           char stationname[20];
0120           sprintf(stationname, "st%d_", station);
0121 
0122           char sectorname[20];
0123           sprintf(sectorname, "sector%02d_", sector);
0124 
0125           for (int component = 0; component < kNumComponents; component++) {
0126             std::stringstream th2f_name, tprofile_name;
0127             th2f_name << "th2f_" << wheelname[wheel + 2] << stationname << sectorname;
0128             tprofile_name << "tprofile_" << wheelname[wheel + 2] << stationname << sectorname;
0129 
0130             double yminmax = 50., xminmax = 0.05;
0131             if (m_minTrackPt > 0.)
0132               xminmax = 1. / m_minTrackPt;
0133             int ynbins = 50;
0134             if (component == kDeltaX) {
0135               th2f_name << "deltax";
0136               tprofile_name << "deltax";
0137             } else if (component == kDeltaDxDz) {
0138               th2f_name << "deltadxdz";
0139               tprofile_name << "deltadxdz";
0140             }
0141 
0142             th2f_wheel_st_sector[wheel + 2][station - 1][sector - 1][component] =
0143                 book2D("/iterN/", th2f_name.str(), "", 30, -xminmax, xminmax, ynbins, -yminmax, yminmax);
0144             tprofile_wheel_st_sector[wheel + 2][station - 1][sector - 1][component] =
0145                 bookProfile("/iterN/", tprofile_name.str(), "", 30, -xminmax, xminmax);
0146           }
0147         }
0148 
0149   // CSC
0150   std::string stname[8] = {"Ep_S1_", "Ep_S2_", "Ep_S3_", "Ep_S4_", "Em_S1_", "Em_S2_", "Em_S3_", "Em_S4_"};
0151   if (m_doCSC)
0152     for (int station = 0; station < 8; station++)
0153       for (int ring = 1; ring <= 3; ring++)
0154         for (int chamber = 1; chamber <= 36; chamber++) {
0155           int st = station % 4 + 1;
0156           if (st > 1 && ring > 2)
0157             continue;  // only station 1 has more then 2 rings
0158           if (st > 1 && ring == 1 && chamber > 18)
0159             continue;  // ring 1 stations 1,2,3 have 18 chambers
0160 
0161           char ringname[20];
0162           sprintf(ringname, "R%d_", ring);
0163 
0164           char chname[20];
0165           sprintf(chname, "C%02d_", chamber);
0166 
0167           for (int component = 0; component < kNumComponents; component++) {
0168             std::stringstream componentname;
0169             double yminmax = 50., xminmax = 0.05;
0170             if (m_minTrackPt > 0.)
0171               xminmax = 1. / m_minTrackPt;
0172             if (ring == 1)
0173               xminmax *= 0.5;
0174             if (component == kDeltaX) {
0175               componentname << "deltax";
0176             } else if (component == kDeltaDxDz) {
0177               componentname << "deltadxdz";
0178             }
0179 
0180             std::stringstream th2f_name, tprofile_name;
0181             th2f_name << "th2f_" << stname[station] << ringname << chname << componentname.str();
0182             tprofile_name << "tprofile_" << stname[station] << ringname << chname << componentname.str();
0183 
0184             th2f_st_ring_chamber[station][ring - 1][chamber - 1][component] =
0185                 book2D("/iterN/", th2f_name.str(), "", 30, -xminmax, xminmax, 100, -yminmax, yminmax);
0186             tprofile_st_ring_chamber[station][ring - 1][chamber - 1][component] =
0187                 bookProfile("/iterN/", tprofile_name.str(), "", 30, -xminmax, xminmax);
0188           }
0189         }
0190 
0191   th1f_trackerRedChi2 = book1D("/iterN/", "trackerRedChi2", "Refit tracker reduced chi^2", 100, 0., 30.);
0192   th1f_trackerRedChi2Diff =
0193       book1D("/iterN/", "trackerRedChi2Diff", "Fit-minus-refit tracker reduced chi^2", 100, -5., 5.);
0194 }
0195 
0196 void AlignmentMonitorMuonVsCurvature::event(const edm::Event &iEvent,
0197                                             const edm::EventSetup &iSetup,
0198                                             const ConstTrajTrackPairCollection &trajtracks) {
0199   edm::Handle<reco::BeamSpot> beamSpot;
0200   iEvent.getByLabel(m_beamSpotTag, beamSpot);
0201 
0202   const GlobalTrackingGeometry *globalGeometry = &iSetup.getData(m_esTokenGBTGeom);
0203   const DetIdAssociator *muonDetIdAssociator_ = &iSetup.getData(m_esTokenDetId);
0204   const Propagator *prop = &iSetup.getData(m_esTokenProp);
0205   const MagneticField *magneticField = &iSetup.getData(m_esTokenMF);
0206   auto builder = iSetup.getHandle(m_esTokenBuilder);
0207 
0208   if (m_muonCollectionTag.label().empty())  // use trajectories
0209   {
0210     for (ConstTrajTrackPairCollection::const_iterator trajtrack = trajtracks.begin(); trajtrack != trajtracks.end();
0211          ++trajtrack) {
0212       const Trajectory *traj = (*trajtrack).first;
0213       const reco::Track *track = (*trajtrack).second;
0214 
0215       if (track->pt() > m_minTrackPt && track->p() > m_minTrackP && fabs(track->dxy(beamSpot->position())) < m_maxDxy) {
0216         MuonResidualsFromTrack muonResidualsFromTrack(
0217             builder, magneticField, globalGeometry, muonDetIdAssociator_, prop, traj, track, pNavigator(), 1000.);
0218         processMuonResidualsFromTrack(muonResidualsFromTrack, traj);
0219       }  // end if track pT is within range
0220     }    // end loop over tracks
0221   } else {
0222     edm::Handle<reco::MuonCollection> muons;
0223     iEvent.getByLabel(m_muonCollectionTag, muons);
0224 
0225     for (reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
0226       if (!(muon->isTrackerMuon() && muon->innerTrack().isNonnull()))
0227         continue;
0228 
0229       if (m_minTrackPt < muon->pt() && m_minTrackP < muon->p() &&
0230           fabs(muon->innerTrack()->dxy(beamSpot->position())) < m_maxDxy) {
0231         MuonResidualsFromTrack muonResidualsFromTrack(globalGeometry, &(*muon), pNavigator(), 100.);
0232         processMuonResidualsFromTrack(muonResidualsFromTrack);
0233       }
0234     }
0235   }
0236 }
0237 
0238 void AlignmentMonitorMuonVsCurvature::processMuonResidualsFromTrack(MuonResidualsFromTrack &mrft,
0239                                                                     const Trajectory *traj) {
0240   if (mrft.trackerNumHits() < m_minTrackerHits)
0241     return;
0242   if (!m_allowTIDTEC && mrft.contains_TIDTEC())
0243     return;
0244 
0245   int nMuChambers = 0;
0246   std::vector<DetId> chamberIds = mrft.chamberIds();
0247   for (unsigned ch = 0; ch < chamberIds.size(); ch++)
0248     if (chamberIds[ch].det() == DetId::Muon)
0249       nMuChambers++;
0250   if (nMuChambers < m_minNCrossedChambers)
0251     return;
0252 
0253   th1f_trackerRedChi2->Fill(mrft.trackerRedChi2());
0254   th1f_trackerRedChi2Diff->Fill(mrft.getTrack()->normalizedChi2() - mrft.trackerRedChi2());
0255 
0256   if (mrft.normalizedChi2() > m_maxTrackerRedChi2)
0257     return;
0258 
0259   double qoverpt = mrft.getTrack()->charge() / mrft.getTrack()->pt();
0260   double qoverpz = 0.;
0261   if (fabs(mrft.getTrack()->pz()) > 0.01)
0262     qoverpz = mrft.getTrack()->charge() / fabs(mrft.getTrack()->pz());
0263 
0264   for (std::vector<DetId>::const_iterator chamberId = chamberIds.begin(); chamberId != chamberIds.end(); ++chamberId) {
0265     if (chamberId->det() != DetId::Muon)
0266       continue;
0267 
0268     if (m_doDT && chamberId->subdetId() == MuonSubdetId::DT) {
0269       DTChamberId dtid(chamberId->rawId());
0270       MuonChamberResidual *dt13 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT13);
0271 
0272       if (dt13 != nullptr && dt13->numHits() >= m_minDT13Hits) {
0273         int wheel = dtid.wheel() + 2;
0274         int station = dtid.station() - 1;
0275         int sector = dtid.sector() - 1;
0276 
0277         double resid_x = 10. * dt13->global_residual();
0278         double resid_dxdz = 1000. * dt13->global_resslope();
0279 
0280         if (fabs(resid_x) < 100. && fabs(resid_dxdz) < 100.) {
0281           th2f_wheel_st_sector[wheel][station][sector][kDeltaX]->Fill(qoverpt, resid_x);
0282           tprofile_wheel_st_sector[wheel][station][sector][kDeltaX]->Fill(qoverpt, resid_x);
0283           th2f_wheel_st_sector[wheel][station][sector][kDeltaDxDz]->Fill(qoverpt, resid_dxdz);
0284           tprofile_wheel_st_sector[wheel][station][sector][kDeltaDxDz]->Fill(qoverpt, resid_dxdz);
0285         }
0286       }  // if it's a good segment
0287     }    // if DT
0288 
0289     if (m_doCSC && chamberId->subdetId() == MuonSubdetId::CSC) {
0290       CSCDetId cscid(chamberId->rawId());
0291       MuonChamberResidual *csc = mrft.chamberResidual(*chamberId, MuonChamberResidual::kCSC);
0292 
0293       if (csc != nullptr && csc->numHits() >= m_minCSCHits) {
0294         int station = 4 * cscid.endcap() + cscid.station() - 5;
0295         int ring = cscid.ring() - 1;
0296         if (cscid.station() == 1 && cscid.ring() == 4)
0297           ring = 0;  // join ME1/a to ME1/b
0298         int chamber = cscid.chamber() - 1;
0299 
0300         double resid_x = 10. * csc->global_residual();
0301         double resid_dxdz = 1000. * csc->global_resslope();
0302 
0303         if (fabs(resid_x) < 100. && fabs(resid_dxdz) < 100.) {
0304           th2f_st_ring_chamber[station][ring][chamber][kDeltaX]->Fill(qoverpz, resid_x);
0305           tprofile_st_ring_chamber[station][ring][chamber][kDeltaX]->Fill(qoverpz, resid_x);
0306           th2f_st_ring_chamber[station][ring][chamber][kDeltaDxDz]->Fill(qoverpz, resid_dxdz);
0307           tprofile_st_ring_chamber[station][ring][chamber][kDeltaDxDz]->Fill(qoverpz, resid_dxdz);
0308         }
0309       }  // if it's a good segment
0310     }    // if CSC
0311 
0312   }  // end loop over chamberIds
0313 }
0314 
0315 DEFINE_EDM_PLUGIN(AlignmentMonitorPluginFactory, AlignmentMonitorMuonVsCurvature, "AlignmentMonitorMuonVsCurvature");