Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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