Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-09-07 04:34:28

0001 /*
0002  * Package:     CommonAlignmentProducer
0003  * Class  :     AlignmentMonitorSegmentDifferences
0004  *
0005  * Original Author:  Jim Pivarski
0006  *         Created:  Mon Nov 12 13:30:14 CST 2007
0007  *
0008  * $Id: AlignmentMonitorSegmentDifferences.cc,v 1.5 2011/04/15 23:09:38 khotilov Exp $
0009  */
0010 
0011 #include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorPluginFactory.h"
0012 #include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorBase.h"
0013 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFromTrack.h"
0014 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsPositionFitter.h"
0015 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsAngleFitter.h"
0016 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsTwoBin.h"
0017 
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/Utilities/interface/InputTag.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0023 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0024 
0025 #include <sstream>
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 AlignmentMonitorSegmentDifferences : public AlignmentMonitorBase {
0035 public:
0036   AlignmentMonitorSegmentDifferences(const edm::ParameterSet &cfg, edm::ConsumesCollector iC);
0037   ~AlignmentMonitorSegmentDifferences() 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);
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   const edm::InputTag m_muonCollectionTag;
0056   const double m_minTrackPt;
0057   const double m_minTrackP;
0058   const double m_maxDxy;
0059   const int m_minTrackerHits;
0060   const double m_maxTrackerRedChi2;
0061   const bool m_allowTIDTEC;
0062   const bool m_minNCrossedChambers;
0063   const int m_minDT13Hits;
0064   const int m_minDT2Hits;
0065   const int m_minCSCHits;
0066   const bool m_doDT;
0067   const bool m_doCSC;
0068 
0069   const edm::EDGetTokenT<reco::BeamSpot> bsToken_;
0070   const edm::EDGetTokenT<reco::MuonCollection> muonToken_;
0071 
0072   // wheel, sector, stationdiff
0073   TProfile *m_dt13_resid[5][12][3];
0074   TProfile *m_dt13_slope[5][12][3];
0075   TProfile *m_dt2_resid[5][12][2];
0076   TProfile *m_dt2_slope[5][12][2];
0077   TH1F *m_posdt13_resid[5][12][3];
0078   TH1F *m_posdt13_slope[5][12][3];
0079   TH1F *m_posdt2_resid[5][12][2];
0080   TH1F *m_posdt2_slope[5][12][2];
0081   TH1F *m_negdt13_resid[5][12][3];
0082   TH1F *m_negdt13_slope[5][12][3];
0083   TH1F *m_negdt2_resid[5][12][2];
0084   TH1F *m_negdt2_slope[5][12][2];
0085 
0086   // endcap, chamber, stationdiff
0087   TProfile *m_cscouter_resid[2][36][2];
0088   TProfile *m_cscouter_slope[2][36][2];
0089   TProfile *m_cscinner_resid[2][18][3];
0090   TProfile *m_cscinner_slope[2][18][3];
0091   TH1F *m_poscscouter_resid[2][36][2];
0092   TH1F *m_poscscouter_slope[2][36][2];
0093   TH1F *m_poscscinner_resid[2][18][3];
0094   TH1F *m_poscscinner_slope[2][18][3];
0095   TH1F *m_negcscouter_resid[2][36][2];
0096   TH1F *m_negcscouter_slope[2][36][2];
0097   TH1F *m_negcscinner_resid[2][18][3];
0098   TH1F *m_negcscinner_slope[2][18][3];
0099 
0100   // cross-system segdiffs: // [endcap][dtsector]
0101   TH1F *m_x_pos_dt1_csc1_resid[2][12];  // tracks going through DT W+-2 St1 and CSC St1 Ring3
0102   TH1F *m_x_pos_dt1_csc2_resid[2][12];  // tracks going through DT W+-2 St1 and CSC St2 Ring2
0103   TH1F *m_x_pos_dt2_csc1_resid[2][12];  // tracks going through DT W+-2 St2 and CSC St1 Ring3
0104   TH1F *m_x_neg_dt1_csc1_resid[2][12];
0105   TH1F *m_x_neg_dt1_csc2_resid[2][12];
0106   TH1F *m_x_neg_dt2_csc1_resid[2][12];
0107 };
0108 
0109 //
0110 // constants, enums and typedefs
0111 //
0112 
0113 //
0114 // static data member definitions
0115 //
0116 
0117 //
0118 // member functions
0119 //
0120 
0121 AlignmentMonitorSegmentDifferences::AlignmentMonitorSegmentDifferences(const edm::ParameterSet &cfg,
0122                                                                        edm::ConsumesCollector iC)
0123     : AlignmentMonitorBase(cfg, iC, "AlignmentMonitorSegmentDifferences"),
0124       m_esTokenGBTGeom(iC.esConsumes()),
0125       m_esTokenDetId(iC.esConsumes(edm::ESInputTag("", "MuonDetIdAssociator"))),
0126       m_esTokenProp(iC.esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))),
0127       m_esTokenMF(iC.esConsumes()),
0128       m_esTokenBuilder(iC.esConsumes(MuonResidualsFromTrack::builderESInputTag())),
0129       m_muonCollectionTag(cfg.getParameter<edm::InputTag>("muonCollectionTag")),
0130       m_minTrackPt(cfg.getParameter<double>("minTrackPt")),
0131       m_minTrackP(cfg.getParameter<double>("minTrackP")),
0132       m_maxDxy(cfg.getParameter<double>("maxDxy")),
0133       m_minTrackerHits(cfg.getParameter<int>("minTrackerHits")),
0134       m_maxTrackerRedChi2(cfg.getParameter<double>("maxTrackerRedChi2")),
0135       m_allowTIDTEC(cfg.getParameter<bool>("allowTIDTEC")),
0136       m_minNCrossedChambers(cfg.getParameter<int>("minNCrossedChambers")),
0137       m_minDT13Hits(cfg.getParameter<int>("minDT13Hits")),
0138       m_minDT2Hits(cfg.getParameter<int>("minDT2Hits")),
0139       m_minCSCHits(cfg.getParameter<int>("minCSCHits")),
0140       m_doDT(cfg.getParameter<bool>("doDT")),
0141       m_doCSC(cfg.getParameter<bool>("doCSC")),
0142       bsToken_(iC.consumes<reco::BeamSpot>(m_beamSpotTag)),
0143       muonToken_(iC.consumes<reco::MuonCollection>(m_muonCollectionTag)) {}
0144 
0145 void AlignmentMonitorSegmentDifferences::book() {
0146   char name[225], pos[228], neg[228];
0147 
0148   double max_curv = 1. / m_minTrackPt;
0149 
0150   if (m_doDT)
0151     for (int wheel = -2; wheel <= +2; wheel++) {
0152       char wheel_label[][2] = {"A", "B", "C", "D", "E"};
0153       for (int sector = 1; sector <= 12; sector++) {
0154         char wheel_sector[50];
0155         sprintf(wheel_sector, "%s_%02d", wheel_label[wheel + 2], sector);
0156 
0157         int nb = 100;
0158         double wnd = 25.;
0159 
0160         sprintf(name, "dt13_resid_%s_12", wheel_sector);
0161         sprintf(pos, "pos%s", name);
0162         sprintf(neg, "neg%s", name);
0163         m_dt13_resid[wheel + 2][sector - 1][0] =
0164             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0165         m_posdt13_resid[wheel + 2][sector - 1][0] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0166         m_negdt13_resid[wheel + 2][sector - 1][0] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0167 
0168         sprintf(name, "dt13_resid_%s_23", wheel_sector);
0169         sprintf(pos, "pos%s", name);
0170         sprintf(neg, "neg%s", name);
0171         m_dt13_resid[wheel + 2][sector - 1][1] =
0172             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0173         m_posdt13_resid[wheel + 2][sector - 1][1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0174         m_negdt13_resid[wheel + 2][sector - 1][1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0175 
0176         sprintf(name, "dt13_resid_%s_34", wheel_sector);
0177         sprintf(pos, "pos%s", name);
0178         sprintf(neg, "neg%s", name);
0179         m_dt13_resid[wheel + 2][sector - 1][2] =
0180             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0181         m_posdt13_resid[wheel + 2][sector - 1][2] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0182         m_negdt13_resid[wheel + 2][sector - 1][2] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0183 
0184         sprintf(name, "dt2_resid_%s_12", wheel_sector);
0185         sprintf(pos, "pos%s", name);
0186         sprintf(neg, "neg%s", name);
0187         m_dt2_resid[wheel + 2][sector - 1][0] =
0188             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -200., 200., " ");
0189         m_posdt2_resid[wheel + 2][sector - 1][0] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0190         m_negdt2_resid[wheel + 2][sector - 1][0] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0191 
0192         sprintf(name, "dt2_resid_%s_23", wheel_sector);
0193         sprintf(pos, "pos%s", name);
0194         sprintf(neg, "neg%s", name);
0195         m_dt2_resid[wheel + 2][sector - 1][1] =
0196             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -200., 200., " ");
0197         m_posdt2_resid[wheel + 2][sector - 1][1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0198         m_negdt2_resid[wheel + 2][sector - 1][1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0199 
0200         sprintf(name, "dt13_slope_%s_12", wheel_sector);
0201         sprintf(pos, "pos%s", name);
0202         sprintf(neg, "neg%s", name);
0203         m_dt13_slope[wheel + 2][sector - 1][0] =
0204             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0205         m_posdt13_slope[wheel + 2][sector - 1][0] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0206         m_negdt13_slope[wheel + 2][sector - 1][0] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0207 
0208         sprintf(name, "dt13_slope_%s_23", wheel_sector);
0209         sprintf(pos, "pos%s", name);
0210         sprintf(neg, "neg%s", name);
0211         m_dt13_slope[wheel + 2][sector - 1][1] =
0212             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0213         m_posdt13_slope[wheel + 2][sector - 1][1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0214         m_negdt13_slope[wheel + 2][sector - 1][1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0215 
0216         sprintf(name, "dt13_slope_%s_34", wheel_sector);
0217         sprintf(pos, "pos%s", name);
0218         sprintf(neg, "neg%s", name);
0219         m_dt13_slope[wheel + 2][sector - 1][2] =
0220             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0221         m_posdt13_slope[wheel + 2][sector - 1][2] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0222         m_negdt13_slope[wheel + 2][sector - 1][2] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0223 
0224         sprintf(name, "dt2_slope_%s_12", wheel_sector);
0225         sprintf(pos, "pos%s", name);
0226         sprintf(neg, "neg%s", name);
0227         m_dt2_slope[wheel + 2][sector - 1][0] =
0228             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -1000., 1000., " ");
0229         m_posdt2_slope[wheel + 2][sector - 1][0] = book1D("/iterN/", pos, pos, nb, -100., 100.);
0230         m_negdt2_slope[wheel + 2][sector - 1][0] = book1D("/iterN/", neg, neg, nb, -100., 100.);
0231 
0232         sprintf(name, "dt2_slope_%s_23", wheel_sector);
0233         sprintf(pos, "pos%s", name);
0234         sprintf(neg, "neg%s", name);
0235         m_dt2_slope[wheel + 2][sector - 1][1] =
0236             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -1000., 1000., " ");
0237         m_posdt2_slope[wheel + 2][sector - 1][1] = book1D("/iterN/", pos, pos, nb, -100., 100.);
0238         m_negdt2_slope[wheel + 2][sector - 1][1] = book1D("/iterN/", neg, neg, nb, -100., 100.);
0239       }
0240     }
0241 
0242   if (m_doCSC)
0243     for (int endcap = 1; endcap <= 2; endcap++) {
0244       std::string endcapletter;
0245       if (endcap == 1)
0246         endcapletter = "p";
0247       else if (endcap == 2)
0248         endcapletter = "m";
0249 
0250       for (int chamber = 1; chamber <= 36; chamber++) {
0251         char ec_chamber[50];
0252         sprintf(ec_chamber, "%s_%02d", endcapletter.c_str(), chamber);
0253 
0254         int nb = 100;
0255         double wnd = 60.;
0256 
0257         sprintf(name, "cscouter_resid_%s_12", ec_chamber);
0258         sprintf(pos, "pos%s", name);
0259         sprintf(neg, "neg%s", name);
0260         m_cscouter_resid[endcap - 1][chamber - 1][0] =
0261             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0262         m_poscscouter_resid[endcap - 1][chamber - 1][0] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0263         m_negcscouter_resid[endcap - 1][chamber - 1][0] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0264 
0265         sprintf(name, "cscouter_resid_%s_23", ec_chamber);
0266         sprintf(pos, "pos%s", name);
0267         sprintf(neg, "neg%s", name);
0268         m_cscouter_resid[endcap - 1][chamber - 1][1] =
0269             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0270         m_poscscouter_resid[endcap - 1][chamber - 1][1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0271         m_negcscouter_resid[endcap - 1][chamber - 1][1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0272 
0273         sprintf(name, "cscouter_slope_%s_12", ec_chamber);
0274         sprintf(pos, "pos%s", name);
0275         sprintf(neg, "neg%s", name);
0276         m_cscouter_slope[endcap - 1][chamber - 1][0] =
0277             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0278         m_poscscouter_slope[endcap - 1][chamber - 1][0] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0279         m_negcscouter_slope[endcap - 1][chamber - 1][0] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0280 
0281         sprintf(name, "cscouter_slope_%s_23", ec_chamber);
0282         sprintf(pos, "pos%s", name);
0283         sprintf(neg, "neg%s", name);
0284         m_cscouter_slope[endcap - 1][chamber - 1][1] =
0285             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0286         m_poscscouter_slope[endcap - 1][chamber - 1][1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0287         m_negcscouter_slope[endcap - 1][chamber - 1][1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0288       }
0289 
0290       for (int chamber = 1; chamber <= 18; chamber++) {
0291         char ec_chamber[50];
0292         sprintf(ec_chamber, "%s_%02d", endcapletter.c_str(), chamber);
0293 
0294         int nb = 100;
0295         double wnd = 40.;
0296 
0297         sprintf(name, "cscinner_resid_%s_12", ec_chamber);
0298         sprintf(pos, "pos%s", name);
0299         sprintf(neg, "neg%s", name);
0300         m_cscinner_resid[endcap - 1][chamber - 1][0] =
0301             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0302         m_poscscinner_resid[endcap - 1][chamber - 1][0] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0303         m_negcscinner_resid[endcap - 1][chamber - 1][0] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0304 
0305         sprintf(name, "cscinner_resid_%s_23", ec_chamber);
0306         sprintf(pos, "pos%s", name);
0307         sprintf(neg, "neg%s", name);
0308         m_cscinner_resid[endcap - 1][chamber - 1][1] =
0309             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0310         m_poscscinner_resid[endcap - 1][chamber - 1][1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0311         m_negcscinner_resid[endcap - 1][chamber - 1][1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0312 
0313         sprintf(name, "cscinner_resid_%s_34", ec_chamber);
0314         sprintf(pos, "pos%s", name);
0315         sprintf(neg, "neg%s", name);
0316         m_cscinner_resid[endcap - 1][chamber - 1][2] =
0317             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0318         m_poscscinner_resid[endcap - 1][chamber - 1][2] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0319         m_negcscinner_resid[endcap - 1][chamber - 1][2] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0320 
0321         sprintf(name, "cscinner_slope_%s_12", ec_chamber);
0322         sprintf(pos, "pos%s", name);
0323         sprintf(neg, "neg%s", name);
0324         m_cscinner_slope[endcap - 1][chamber - 1][0] =
0325             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0326         m_poscscinner_slope[endcap - 1][chamber - 1][0] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0327         m_negcscinner_slope[endcap - 1][chamber - 1][0] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0328 
0329         sprintf(name, "cscinner_slope_%s_23", ec_chamber);
0330         sprintf(pos, "pos%s", name);
0331         sprintf(neg, "neg%s", name);
0332         m_cscinner_slope[endcap - 1][chamber - 1][1] =
0333             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0334         m_poscscinner_slope[endcap - 1][chamber - 1][1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0335         m_negcscinner_slope[endcap - 1][chamber - 1][1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0336 
0337         sprintf(name, "cscinner_slope_%s_34", ec_chamber);
0338         sprintf(pos, "pos%s", name);
0339         sprintf(neg, "neg%s", name);
0340         m_cscinner_slope[endcap - 1][chamber - 1][2] =
0341             bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
0342         m_poscscinner_slope[endcap - 1][chamber - 1][2] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0343         m_negcscinner_slope[endcap - 1][chamber - 1][2] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0344       }
0345     }
0346 
0347   // cross-system
0348   for (int e = 1; e <= 2; e++)
0349     for (int s = 1; s <= 12; s++) {
0350       char endcap_sector[50];
0351       if (e == 1)
0352         sprintf(endcap_sector, "Wp2S%02d", s);
0353       if (e == 2)
0354         sprintf(endcap_sector, "Wm2S%02d", s);
0355 
0356       int nb = 200;
0357       double wnd = 100.;
0358 
0359       sprintf(pos, "pos_x_dt1_csc1_%s", endcap_sector);
0360       sprintf(neg, "neg_x_dt1_csc1_%s", endcap_sector);
0361       m_x_pos_dt1_csc1_resid[e - 1][s - 1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0362       m_x_neg_dt1_csc1_resid[e - 1][s - 1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0363 
0364       sprintf(pos, "pos_x_dt1_csc2_%s", endcap_sector);
0365       sprintf(neg, "neg_x_dt1_csc2_%s", endcap_sector);
0366       m_x_pos_dt1_csc2_resid[e - 1][s - 1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0367       m_x_neg_dt1_csc2_resid[e - 1][s - 1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0368 
0369       sprintf(pos, "pos_x_dt2_csc1_%s", endcap_sector);
0370       sprintf(neg, "neg_x_dt2_csc1_%s", endcap_sector);
0371       m_x_pos_dt2_csc1_resid[e - 1][s - 1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
0372       m_x_neg_dt2_csc1_resid[e - 1][s - 1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
0373     }
0374 }
0375 
0376 void AlignmentMonitorSegmentDifferences::event(const edm::Event &iEvent,
0377                                                const edm::EventSetup &iSetup,
0378                                                const ConstTrajTrackPairCollection &trajtracks) {
0379   const edm::Handle<reco::BeamSpot> &beamSpot = iEvent.getHandle(bsToken_);
0380 
0381   const GlobalTrackingGeometry *globalGeometry = &iSetup.getData(m_esTokenGBTGeom);
0382   const DetIdAssociator *muonDetIdAssociator_ = &iSetup.getData(m_esTokenDetId);
0383   const Propagator *prop = &iSetup.getData(m_esTokenProp);
0384   const MagneticField *magneticField = &iSetup.getData(m_esTokenMF);
0385   auto builder = iSetup.getHandle(m_esTokenBuilder);
0386 
0387   if (m_muonCollectionTag.label().empty())  // use trajectories
0388   {
0389     for (ConstTrajTrackPairCollection::const_iterator trajtrack = trajtracks.begin(); trajtrack != trajtracks.end();
0390          ++trajtrack) {
0391       const Trajectory *traj = (*trajtrack).first;
0392       const reco::Track *track = (*trajtrack).second;
0393 
0394       if (track->pt() > m_minTrackPt && track->p() > m_minTrackP && fabs(track->dxy(beamSpot->position())) < m_maxDxy) {
0395         MuonResidualsFromTrack muonResidualsFromTrack(
0396             builder, magneticField, globalGeometry, muonDetIdAssociator_, prop, traj, track, pNavigator(), 1000.);
0397         processMuonResidualsFromTrack(muonResidualsFromTrack);
0398       }
0399     }  // end loop over tracks
0400   } else {
0401     const edm::Handle<reco::MuonCollection> &muons = iEvent.getHandle(muonToken_);
0402 
0403     for (reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
0404       if (!(muon->isTrackerMuon() && muon->innerTrack().isNonnull()))
0405         continue;
0406 
0407       if (m_minTrackPt < muon->pt() && m_minTrackP < muon->p() &&
0408           fabs(muon->innerTrack()->dxy(beamSpot->position())) < m_maxDxy) {
0409         MuonResidualsFromTrack muonResidualsFromTrack(globalGeometry, &(*muon), pNavigator(), 100.);
0410         processMuonResidualsFromTrack(muonResidualsFromTrack);
0411       }
0412     }
0413   }
0414 }
0415 
0416 void AlignmentMonitorSegmentDifferences::processMuonResidualsFromTrack(MuonResidualsFromTrack &mrft) {
0417   if (mrft.trackerNumHits() < m_minTrackerHits)
0418     return;
0419   if (!m_allowTIDTEC && mrft.contains_TIDTEC())
0420     return;
0421   if (mrft.normalizedChi2() > m_maxTrackerRedChi2)
0422     return;
0423 
0424   int nMuChambers = 0;
0425   std::vector<DetId> chamberIds = mrft.chamberIds();
0426   for (unsigned ch = 0; ch < chamberIds.size(); ch++)
0427     if (chamberIds[ch].det() == DetId::Muon)
0428       nMuChambers++;
0429   if (nMuChambers < m_minNCrossedChambers)
0430     return;
0431 
0432   double qoverpt = (mrft.getTrack()->charge() > 0 ? 1. : -1.) / mrft.getTrack()->pt();
0433   double qoverpz = 0.;
0434   if (fabs(mrft.getTrack()->pz()) > 0.01)
0435     qoverpz = mrft.getTrack()->charge() / fabs(mrft.getTrack()->pz());
0436 
0437   for (std::vector<DetId>::const_iterator chamberId = chamberIds.begin(); chamberId != chamberIds.end(); ++chamberId) {
0438     if (chamberId->det() != DetId::Muon)
0439       continue;
0440 
0441     // **************** DT ****************
0442     if (m_doDT && chamberId->subdetId() == MuonSubdetId::DT) {
0443       MuonChamberResidual *dt13 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT13);
0444       MuonChamberResidual *dt2 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT2);
0445 
0446       if (dt13 != nullptr && dt13->numHits() >= m_minDT13Hits) {
0447         DTChamberId thisid(chamberId->rawId());
0448         for (std::vector<DetId>::const_iterator otherId = chamberIds.begin(); otherId != chamberIds.end(); ++otherId) {
0449           if (otherId->det() == DetId::Muon && otherId->subdetId() == MuonSubdetId::DT) {
0450             DTChamberId thatid(otherId->rawId());
0451             if (thisid.rawId() != thatid.rawId() && thisid.wheel() == thatid.wheel() &&
0452                 thisid.sector() == thatid.sector()) {
0453               MuonChamberResidual *dt13other = mrft.chamberResidual(*otherId, MuonChamberResidual::kDT13);
0454               if (dt13other != nullptr && dt13other->numHits() >= m_minDT13Hits) {
0455                 double slopediff = 1000. * (dt13->global_resslope() - dt13other->global_resslope());
0456                 //double length = dt13->chamberAlignable()->surface().toGlobal(align::LocalPoint(0,0,0)).perp() -
0457                 //                dt13other->chamberAlignable()->surface().toGlobal(align::LocalPoint(0,0,0)).perp();
0458                 //double residdiff = 10. * (dt13->global_residual() + length*dt13->global_resslope() - dt13other->global_residual());
0459                 double residdiff = 10. * (dt13->global_residual() - dt13other->global_residual());
0460 
0461                 int st = 0;
0462                 if (thatid.station() - thisid.station() == 1)
0463                   st = thisid.station();
0464                 if (st > 0) {
0465                   m_dt13_resid[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(qoverpt, residdiff);
0466                   m_dt13_slope[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(qoverpt, slopediff);
0467                   if (qoverpt > 0) {
0468                     m_posdt13_resid[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(residdiff);
0469                     m_posdt13_slope[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(slopediff);
0470                   } else {
0471                     m_negdt13_resid[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(residdiff);
0472                     m_negdt13_slope[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(slopediff);
0473                   }
0474                 }
0475               }  // end other numhits
0476             }  // end this near other
0477           }  // end other is DT
0478 
0479           // cross-system: other is CSC
0480           // only do it for DT stubs in W+-2 St1&2:
0481           if (!(abs(thisid.wheel()) == 2 && (thisid.station() == 1 || thisid.station() == 2)))
0482             continue;
0483           if (otherId->det() == DetId::Muon && otherId->subdetId() == MuonSubdetId::CSC) {
0484             CSCDetId thatid(otherId->rawId());
0485             //only do it for CSC stubs in St1R3 or St2R2:
0486             if (!((thatid.station() == 1 && thatid.ring() == 3) || (thatid.station() == 2 && thatid.ring() == 2)))
0487               continue;
0488 
0489             MuonChamberResidual *cscother = mrft.chamberResidual(*otherId, MuonChamberResidual::kCSC);
0490             if (cscother != nullptr && cscother->numHits() >= m_minCSCHits) {
0491               // scale to adjust the csc residual size to be comparabe to dt's one
0492               double csc_scale = dt13->chamberAlignable()
0493                                      ->surface()
0494                                      .toGlobal(align::LocalPoint(dt13->trackx(), dt13->tracky(), 0))
0495                                      .perp() /
0496                                  cscother->chamberAlignable()
0497                                      ->surface()
0498                                      .toGlobal(align::LocalPoint(cscother->trackx(), cscother->tracky(), 0))
0499                                      .perp();
0500               double residdiff = 10. * (dt13->global_residual() - cscother->global_residual() * csc_scale);
0501               if (thisid.station() == 1 && thatid.station() == 1) {
0502                 if (qoverpt > 0)
0503                   m_x_pos_dt1_csc1_resid[thatid.endcap() - 1][thisid.sector() - 1]->Fill(residdiff);
0504                 else
0505                   m_x_neg_dt1_csc1_resid[thatid.endcap() - 1][thisid.sector() - 1]->Fill(residdiff);
0506               } else if (thisid.station() == 1 && thatid.station() == 2) {
0507                 if (qoverpt > 0)
0508                   m_x_pos_dt1_csc2_resid[thatid.endcap() - 1][thisid.sector() - 1]->Fill(residdiff);
0509                 else
0510                   m_x_neg_dt1_csc2_resid[thatid.endcap() - 1][thisid.sector() - 1]->Fill(residdiff);
0511               } else if (thisid.station() == 2 && thatid.station() == 1) {
0512                 if (qoverpt > 0)
0513                   m_x_pos_dt2_csc1_resid[thatid.endcap() - 1][thisid.sector() - 1]->Fill(residdiff);
0514                 else
0515                   m_x_neg_dt2_csc1_resid[thatid.endcap() - 1][thisid.sector() - 1]->Fill(residdiff);
0516               }
0517             }
0518           }  // end other is CSC
0519         }  // end loop over other
0520       }  // end if DT13
0521 
0522       // z-direction
0523       if (dt2 != nullptr && dt2->numHits() >= m_minDT2Hits && (dt2->chi2() / double(dt2->ndof())) < 2.0) {
0524         DTChamberId thisid(chamberId->rawId());
0525         for (std::vector<DetId>::const_iterator otherId = chamberIds.begin(); otherId != chamberIds.end(); ++otherId) {
0526           if (otherId->det() == DetId::Muon && otherId->subdetId() == MuonSubdetId::DT) {
0527             DTChamberId thatid(otherId->rawId());
0528             if (thisid.rawId() != thatid.rawId() && thisid.wheel() == thatid.wheel() &&
0529                 thisid.sector() == thatid.sector()) {
0530               MuonChamberResidual *dt2other = mrft.chamberResidual(*otherId, MuonChamberResidual::kDT2);
0531               if (dt2other != nullptr && dt2other->numHits() >= m_minDT2Hits) {
0532                 double slopediff = 1000. * (dt2->global_resslope() - dt2other->global_resslope());
0533                 //double length = dt2->chamberAlignable()->surface().toGlobal(align::LocalPoint(0,0,0)).perp() -
0534                 //                dt2other->chamberAlignable()->surface().toGlobal(align::LocalPoint(0,0,0)).perp();
0535                 //double residdiff = 10. * (dt2->global_residual() + length*dt2->global_resslope() - dt2other->global_residual());
0536                 double residdiff = 10. * (dt2->global_residual() - dt2other->global_residual());
0537 
0538                 int st = 0;
0539                 if (thatid.station() - thisid.station() == 1)
0540                   st = thisid.station();
0541                 if (st > 0) {
0542                   m_dt2_resid[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(qoverpt, residdiff);
0543                   m_dt2_slope[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(qoverpt, slopediff);
0544                   if (qoverpt > 0) {
0545                     m_posdt2_resid[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(residdiff);
0546                     m_posdt2_slope[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(slopediff);
0547                   } else {
0548                     m_negdt2_resid[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(residdiff);
0549                     m_negdt2_slope[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(slopediff);
0550                   }
0551                 }
0552               }  // end other numhits
0553             }  // end this near other
0554           }  // end other is DT
0555         }  // end loop over other
0556       }  // end if DT2
0557     }  // end if DT
0558 
0559     // **************** CSC ****************
0560     else if (m_doCSC && chamberId->subdetId() == MuonSubdetId::CSC) {
0561       MuonChamberResidual *csc = mrft.chamberResidual(*chamberId, MuonChamberResidual::kCSC);
0562       if (csc->numHits() >= m_minCSCHits) {
0563         CSCDetId thisid(chamberId->rawId());
0564         for (std::vector<DetId>::const_iterator otherId = chamberIds.begin(); otherId != chamberIds.end(); ++otherId) {
0565           if (otherId->det() == DetId::Muon && otherId->subdetId() == MuonSubdetId::CSC) {
0566             CSCDetId thatid(otherId->rawId());
0567             if (thisid.rawId() != thatid.rawId() && thisid.endcap() == thatid.endcap()) {
0568               MuonChamberResidual *cscother = mrft.chamberResidual(*otherId, MuonChamberResidual::kCSC);
0569               if (cscother != nullptr && cscother->numHits() >= m_minCSCHits) {
0570                 double slopediff = 1000. * (csc->global_resslope() - cscother->global_resslope());
0571                 //double length = csc->chamberAlignable()->surface().toGlobal(align::LocalPoint(0,0,0)).z() -
0572                 //                cscother->chamberAlignable()->surface().toGlobal(align::LocalPoint(0,0,0)).z();
0573                 //double residdiff = 10. * (csc->global_residual() + length*csc->global_resslope() - cscother->global_residual());
0574                 double residdiff = 10. * (csc->global_residual() - cscother->global_residual());
0575 
0576                 int thischamber = thisid.chamber();
0577                 int thisring = thisid.ring();
0578                 if (thisid.station() == 1 && (thisring == 1 || thisring == 4)) {
0579                   thischamber = (thischamber - 1) / 2 + 1;
0580                   thisring = 1;
0581                 }
0582 
0583                 if (thisring == thatid.ring() && thischamber == thatid.chamber()) {
0584                   bool inner = (thisring == 1);
0585                   bool outer = (thisring == 2);
0586                   int st = 0;
0587                   if (thatid.station() - thisid.station() == 1 && (inner || thisid.station() < 3))
0588                     st = thisid.station();
0589 
0590                   if (outer && st > 0) {
0591                     m_cscouter_resid[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(qoverpz, residdiff);
0592                     m_cscouter_slope[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(qoverpz, slopediff);
0593                     if (qoverpz > 0) {
0594                       m_poscscouter_resid[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(residdiff);
0595                       m_poscscouter_slope[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(slopediff);
0596                     } else {
0597                       m_negcscouter_resid[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(residdiff);
0598                       m_negcscouter_slope[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(slopediff);
0599                     }
0600                   }
0601                   if (inner && st > 0) {
0602                     m_cscinner_resid[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(qoverpz, residdiff);
0603                     m_cscinner_slope[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(qoverpz, slopediff);
0604                     if (qoverpz > 0) {
0605                       m_poscscinner_resid[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(residdiff);
0606                       m_poscscinner_slope[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(slopediff);
0607                     } else {
0608                       m_negcscinner_resid[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(residdiff);
0609                       m_negcscinner_slope[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(slopediff);
0610                     }
0611                   }
0612                 }  // end of same ring&chamber
0613               }  // end other min numhits
0614             }  // end this near other
0615           }  // end other is CSC
0616         }  // end loop over other
0617 
0618       }  // end if this min numhits
0619     }  // end if CSC
0620 
0621   }  // end loop over chamberIds
0622 }
0623 
0624 DEFINE_EDM_PLUGIN(AlignmentMonitorPluginFactory,
0625                   AlignmentMonitorSegmentDifferences,
0626                   "AlignmentMonitorSegmentDifferences");