AlignmentMonitorSegmentDifferences.cc

Go to the documentation of this file.

/*
 * Package:     CommonAlignmentProducer
 * Class  :     AlignmentMonitorSegmentDifferences
 *
 * Original Author:  Jim Pivarski
 *         Created:  Mon Nov 12 13:30:14 CST 2007
 *
 * $Id: AlignmentMonitorSegmentDifferences.cc,v 1.5 2011/04/15 23:09:38 khotilov Exp $
 */

#include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorPluginFactory.h"
#include "Alignment/CommonAlignmentMonitor/interface/AlignmentMonitorBase.h"
#include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFromTrack.h"
#include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsPositionFitter.h"
#include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsAngleFitter.h"
#include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsTwoBin.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"

#include <sstream>

#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
#include "TrackingTools/Records/interface/DetIdAssociatorRecord.h"
#include "TrackingTools/GeomPropagators/interface/Propagator.h"
#include "TrackingTools/TrackAssociator/interface/DetIdAssociator.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "MagneticField/Engine/interface/MagneticField.h"

class AlignmentMonitorSegmentDifferences : public AlignmentMonitorBase {
public:
  AlignmentMonitorSegmentDifferences(const edm::ParameterSet &cfg, edm::ConsumesCollector iC);
  ~AlignmentMonitorSegmentDifferences() override {}

  void book() override;

  void event(const edm::Event &iEvent,
             const edm::EventSetup &iSetup,
             const ConstTrajTrackPairCollection &iTrajTracks) override;
  void processMuonResidualsFromTrack(MuonResidualsFromTrack &mrft);

private:
  // es token
  const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> m_esTokenGBTGeom;
  const edm::ESGetToken<DetIdAssociator, DetIdAssociatorRecord> m_esTokenDetId;
  const edm::ESGetToken<Propagator, TrackingComponentsRecord> m_esTokenProp;
  const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> m_esTokenMF;
  const MuonResidualsFromTrack::BuilderToken m_esTokenBuilder;

  // parameters
  edm::InputTag m_muonCollectionTag;
  double m_minTrackPt;
  double m_minTrackP;
  double m_maxDxy;
  int m_minTrackerHits;
  double m_maxTrackerRedChi2;
  bool m_allowTIDTEC;
  bool m_minNCrossedChambers;
  int m_minDT13Hits;
  int m_minDT2Hits;
  int m_minCSCHits;
  bool m_doDT;
  bool m_doCSC;

  // wheel, sector, stationdiff
  TProfile *m_dt13_resid[5][12][3];
  TProfile *m_dt13_slope[5][12][3];
  TProfile *m_dt2_resid[5][12][2];
  TProfile *m_dt2_slope[5][12][2];
  TH1F *m_posdt13_resid[5][12][3];
  TH1F *m_posdt13_slope[5][12][3];
  TH1F *m_posdt2_resid[5][12][2];
  TH1F *m_posdt2_slope[5][12][2];
  TH1F *m_negdt13_resid[5][12][3];
  TH1F *m_negdt13_slope[5][12][3];
  TH1F *m_negdt2_resid[5][12][2];
  TH1F *m_negdt2_slope[5][12][2];

  // endcap, chamber, stationdiff
  TProfile *m_cscouter_resid[2][36][2];
  TProfile *m_cscouter_slope[2][36][2];
  TProfile *m_cscinner_resid[2][18][3];
  TProfile *m_cscinner_slope[2][18][3];
  TH1F *m_poscscouter_resid[2][36][2];
  TH1F *m_poscscouter_slope[2][36][2];
  TH1F *m_poscscinner_resid[2][18][3];
  TH1F *m_poscscinner_slope[2][18][3];
  TH1F *m_negcscouter_resid[2][36][2];
  TH1F *m_negcscouter_slope[2][36][2];
  TH1F *m_negcscinner_resid[2][18][3];
  TH1F *m_negcscinner_slope[2][18][3];

  // cross-system segdiffs: // [endcap][dtsector]
  TH1F *m_x_pos_dt1_csc1_resid[2][12];  // tracks going through DT W+-2 St1 and CSC St1 Ring3
  TH1F *m_x_pos_dt1_csc2_resid[2][12];  // tracks going through DT W+-2 St1 and CSC St2 Ring2
  TH1F *m_x_pos_dt2_csc1_resid[2][12];  // tracks going through DT W+-2 St2 and CSC St1 Ring3
  TH1F *m_x_neg_dt1_csc1_resid[2][12];
  TH1F *m_x_neg_dt1_csc2_resid[2][12];
  TH1F *m_x_neg_dt2_csc1_resid[2][12];
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// member functions
//

AlignmentMonitorSegmentDifferences::AlignmentMonitorSegmentDifferences(const edm::ParameterSet &cfg,
                                                                       edm::ConsumesCollector iC)
    : AlignmentMonitorBase(cfg, iC, "AlignmentMonitorSegmentDifferences"),
      m_esTokenGBTGeom(iC.esConsumes()),
      m_esTokenDetId(iC.esConsumes(edm::ESInputTag("", "MuonDetIdAssociator"))),
      m_esTokenProp(iC.esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))),
      m_esTokenMF(iC.esConsumes()),
      m_esTokenBuilder(iC.esConsumes(MuonResidualsFromTrack::builderESInputTag())),
      m_muonCollectionTag(cfg.getParameter<edm::InputTag>("muonCollectionTag")),
      m_minTrackPt(cfg.getParameter<double>("minTrackPt")),
      m_minTrackP(cfg.getParameter<double>("minTrackP")),
      m_maxDxy(cfg.getParameter<double>("maxDxy")),
      m_minTrackerHits(cfg.getParameter<int>("minTrackerHits")),
      m_maxTrackerRedChi2(cfg.getParameter<double>("maxTrackerRedChi2")),
      m_allowTIDTEC(cfg.getParameter<bool>("allowTIDTEC")),
      m_minNCrossedChambers(cfg.getParameter<int>("minNCrossedChambers")),
      m_minDT13Hits(cfg.getParameter<int>("minDT13Hits")),
      m_minDT2Hits(cfg.getParameter<int>("minDT2Hits")),
      m_minCSCHits(cfg.getParameter<int>("minCSCHits")),
      m_doDT(cfg.getParameter<bool>("doDT")),
      m_doCSC(cfg.getParameter<bool>("doCSC")) {}

void AlignmentMonitorSegmentDifferences::book() {
  char name[225], pos[228], neg[228];

  double max_curv = 1. / m_minTrackPt;

  if (m_doDT)
    for (int wheel = -2; wheel <= +2; wheel++) {
      char wheel_label[][2] = {"A", "B", "C", "D", "E"};
      for (int sector = 1; sector <= 12; sector++) {
        char wheel_sector[50];
        sprintf(wheel_sector, "%s_%02d", wheel_label[wheel + 2], sector);

        int nb = 100;
        double wnd = 25.;

        sprintf(name, "dt13_resid_%s_12", wheel_sector);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_dt13_resid[wheel + 2][sector - 1][0] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_posdt13_resid[wheel + 2][sector - 1][0] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negdt13_resid[wheel + 2][sector - 1][0] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "dt13_resid_%s_23", wheel_sector);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_dt13_resid[wheel + 2][sector - 1][1] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_posdt13_resid[wheel + 2][sector - 1][1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negdt13_resid[wheel + 2][sector - 1][1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "dt13_resid_%s_34", wheel_sector);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_dt13_resid[wheel + 2][sector - 1][2] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_posdt13_resid[wheel + 2][sector - 1][2] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negdt13_resid[wheel + 2][sector - 1][2] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "dt2_resid_%s_12", wheel_sector);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_dt2_resid[wheel + 2][sector - 1][0] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -200., 200., " ");
        m_posdt2_resid[wheel + 2][sector - 1][0] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negdt2_resid[wheel + 2][sector - 1][0] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "dt2_resid_%s_23", wheel_sector);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_dt2_resid[wheel + 2][sector - 1][1] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -200., 200., " ");
        m_posdt2_resid[wheel + 2][sector - 1][1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negdt2_resid[wheel + 2][sector - 1][1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "dt13_slope_%s_12", wheel_sector);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_dt13_slope[wheel + 2][sector - 1][0] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_posdt13_slope[wheel + 2][sector - 1][0] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negdt13_slope[wheel + 2][sector - 1][0] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "dt13_slope_%s_23", wheel_sector);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_dt13_slope[wheel + 2][sector - 1][1] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_posdt13_slope[wheel + 2][sector - 1][1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negdt13_slope[wheel + 2][sector - 1][1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "dt13_slope_%s_34", wheel_sector);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_dt13_slope[wheel + 2][sector - 1][2] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_posdt13_slope[wheel + 2][sector - 1][2] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negdt13_slope[wheel + 2][sector - 1][2] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "dt2_slope_%s_12", wheel_sector);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_dt2_slope[wheel + 2][sector - 1][0] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -1000., 1000., " ");
        m_posdt2_slope[wheel + 2][sector - 1][0] = book1D("/iterN/", pos, pos, nb, -100., 100.);
        m_negdt2_slope[wheel + 2][sector - 1][0] = book1D("/iterN/", neg, neg, nb, -100., 100.);

        sprintf(name, "dt2_slope_%s_23", wheel_sector);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_dt2_slope[wheel + 2][sector - 1][1] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -1000., 1000., " ");
        m_posdt2_slope[wheel + 2][sector - 1][1] = book1D("/iterN/", pos, pos, nb, -100., 100.);
        m_negdt2_slope[wheel + 2][sector - 1][1] = book1D("/iterN/", neg, neg, nb, -100., 100.);
      }
    }

  if (m_doCSC)
    for (int endcap = 1; endcap <= 2; endcap++) {
      std::string endcapletter;
      if (endcap == 1)
        endcapletter = "p";
      else if (endcap == 2)
        endcapletter = "m";

      for (int chamber = 1; chamber <= 36; chamber++) {
        char ec_chamber[50];
        sprintf(ec_chamber, "%s_%02d", endcapletter.c_str(), chamber);

        int nb = 100;
        double wnd = 60.;

        sprintf(name, "cscouter_resid_%s_12", ec_chamber);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_cscouter_resid[endcap - 1][chamber - 1][0] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_poscscouter_resid[endcap - 1][chamber - 1][0] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negcscouter_resid[endcap - 1][chamber - 1][0] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "cscouter_resid_%s_23", ec_chamber);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_cscouter_resid[endcap - 1][chamber - 1][1] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_poscscouter_resid[endcap - 1][chamber - 1][1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negcscouter_resid[endcap - 1][chamber - 1][1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "cscouter_slope_%s_12", ec_chamber);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_cscouter_slope[endcap - 1][chamber - 1][0] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_poscscouter_slope[endcap - 1][chamber - 1][0] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negcscouter_slope[endcap - 1][chamber - 1][0] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "cscouter_slope_%s_23", ec_chamber);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_cscouter_slope[endcap - 1][chamber - 1][1] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_poscscouter_slope[endcap - 1][chamber - 1][1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negcscouter_slope[endcap - 1][chamber - 1][1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
      }

      for (int chamber = 1; chamber <= 18; chamber++) {
        char ec_chamber[50];
        sprintf(ec_chamber, "%s_%02d", endcapletter.c_str(), chamber);

        int nb = 100;
        double wnd = 40.;

        sprintf(name, "cscinner_resid_%s_12", ec_chamber);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_cscinner_resid[endcap - 1][chamber - 1][0] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_poscscinner_resid[endcap - 1][chamber - 1][0] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negcscinner_resid[endcap - 1][chamber - 1][0] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "cscinner_resid_%s_23", ec_chamber);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_cscinner_resid[endcap - 1][chamber - 1][1] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_poscscinner_resid[endcap - 1][chamber - 1][1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negcscinner_resid[endcap - 1][chamber - 1][1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "cscinner_resid_%s_34", ec_chamber);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_cscinner_resid[endcap - 1][chamber - 1][2] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_poscscinner_resid[endcap - 1][chamber - 1][2] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negcscinner_resid[endcap - 1][chamber - 1][2] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "cscinner_slope_%s_12", ec_chamber);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_cscinner_slope[endcap - 1][chamber - 1][0] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_poscscinner_slope[endcap - 1][chamber - 1][0] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negcscinner_slope[endcap - 1][chamber - 1][0] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "cscinner_slope_%s_23", ec_chamber);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_cscinner_slope[endcap - 1][chamber - 1][1] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_poscscinner_slope[endcap - 1][chamber - 1][1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negcscinner_slope[endcap - 1][chamber - 1][1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

        sprintf(name, "cscinner_slope_%s_34", ec_chamber);
        sprintf(pos, "pos%s", name);
        sprintf(neg, "neg%s", name);
        m_cscinner_slope[endcap - 1][chamber - 1][2] =
            bookProfile("/iterN/", name, name, 20, -max_curv, max_curv, 1, -100., 100., " ");
        m_poscscinner_slope[endcap - 1][chamber - 1][2] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
        m_negcscinner_slope[endcap - 1][chamber - 1][2] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
      }
    }

  // cross-system
  for (int e = 1; e <= 2; e++)
    for (int s = 1; s <= 12; s++) {
      char endcap_sector[50];
      if (e == 1)
        sprintf(endcap_sector, "Wp2S%02d", s);
      if (e == 2)
        sprintf(endcap_sector, "Wm2S%02d", s);

      int nb = 200;
      double wnd = 100.;

      sprintf(pos, "pos_x_dt1_csc1_%s", endcap_sector);
      sprintf(neg, "neg_x_dt1_csc1_%s", endcap_sector);
      m_x_pos_dt1_csc1_resid[e - 1][s - 1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
      m_x_neg_dt1_csc1_resid[e - 1][s - 1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

      sprintf(pos, "pos_x_dt1_csc2_%s", endcap_sector);
      sprintf(neg, "neg_x_dt1_csc2_%s", endcap_sector);
      m_x_pos_dt1_csc2_resid[e - 1][s - 1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
      m_x_neg_dt1_csc2_resid[e - 1][s - 1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);

      sprintf(pos, "pos_x_dt2_csc1_%s", endcap_sector);
      sprintf(neg, "neg_x_dt2_csc1_%s", endcap_sector);
      m_x_pos_dt2_csc1_resid[e - 1][s - 1] = book1D("/iterN/", pos, pos, nb, -wnd, wnd);
      m_x_neg_dt2_csc1_resid[e - 1][s - 1] = book1D("/iterN/", neg, neg, nb, -wnd, wnd);
    }
}

void AlignmentMonitorSegmentDifferences::event(const edm::Event &iEvent,
                                               const edm::EventSetup &iSetup,
                                               const ConstTrajTrackPairCollection &trajtracks) {
  edm::Handle<reco::BeamSpot> beamSpot;
  iEvent.getByLabel(m_beamSpotTag, beamSpot);

  const GlobalTrackingGeometry *globalGeometry = &iSetup.getData(m_esTokenGBTGeom);
  const DetIdAssociator *muonDetIdAssociator_ = &iSetup.getData(m_esTokenDetId);
  const Propagator *prop = &iSetup.getData(m_esTokenProp);
  const MagneticField *magneticField = &iSetup.getData(m_esTokenMF);
  auto builder = iSetup.getHandle(m_esTokenBuilder);

  if (m_muonCollectionTag.label().empty())  // use trajectories
  {
    for (ConstTrajTrackPairCollection::const_iterator trajtrack = trajtracks.begin(); trajtrack != trajtracks.end();
         ++trajtrack) {
      const Trajectory *traj = (*trajtrack).first;
      const reco::Track *track = (*trajtrack).second;

      if (track->pt() > m_minTrackPt && track->p() > m_minTrackP && fabs(track->dxy(beamSpot->position())) < m_maxDxy) {
        MuonResidualsFromTrack muonResidualsFromTrack(
            builder, magneticField, globalGeometry, muonDetIdAssociator_, prop, traj, track, pNavigator(), 1000.);
        processMuonResidualsFromTrack(muonResidualsFromTrack);
      }
    }  // end loop over tracks
  } else {
    edm::Handle<reco::MuonCollection> muons;
    iEvent.getByLabel(m_muonCollectionTag, muons);

    for (reco::MuonCollection::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
      if (!(muon->isTrackerMuon() && muon->innerTrack().isNonnull()))
        continue;

      if (m_minTrackPt < muon->pt() && m_minTrackP < muon->p() &&
          fabs(muon->innerTrack()->dxy(beamSpot->position())) < m_maxDxy) {
        MuonResidualsFromTrack muonResidualsFromTrack(globalGeometry, &(*muon), pNavigator(), 100.);
        processMuonResidualsFromTrack(muonResidualsFromTrack);
      }
    }
  }
}

void AlignmentMonitorSegmentDifferences::processMuonResidualsFromTrack(MuonResidualsFromTrack &mrft) {
  if (mrft.trackerNumHits() < m_minTrackerHits)
    return;
  if (!m_allowTIDTEC && mrft.contains_TIDTEC())
    return;
  if (mrft.normalizedChi2() > m_maxTrackerRedChi2)
    return;

  int nMuChambers = 0;
  std::vector<DetId> chamberIds = mrft.chamberIds();
  for (unsigned ch = 0; ch < chamberIds.size(); ch++)
    if (chamberIds[ch].det() == DetId::Muon)
      nMuChambers++;
  if (nMuChambers < m_minNCrossedChambers)
    return;

  double qoverpt = (mrft.getTrack()->charge() > 0 ? 1. : -1.) / mrft.getTrack()->pt();
  double qoverpz = 0.;
  if (fabs(mrft.getTrack()->pz()) > 0.01)
    qoverpz = mrft.getTrack()->charge() / fabs(mrft.getTrack()->pz());

  for (std::vector<DetId>::const_iterator chamberId = chamberIds.begin(); chamberId != chamberIds.end(); ++chamberId) {
    if (chamberId->det() != DetId::Muon)
      continue;

    // **************** DT ****************
    if (m_doDT && chamberId->subdetId() == MuonSubdetId::DT) {
      MuonChamberResidual *dt13 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT13);
      MuonChamberResidual *dt2 = mrft.chamberResidual(*chamberId, MuonChamberResidual::kDT2);

      if (dt13 != nullptr && dt13->numHits() >= m_minDT13Hits) {
        DTChamberId thisid(chamberId->rawId());
        for (std::vector<DetId>::const_iterator otherId = chamberIds.begin(); otherId != chamberIds.end(); ++otherId) {
          if (otherId->det() == DetId::Muon && otherId->subdetId() == MuonSubdetId::DT) {
            DTChamberId thatid(otherId->rawId());
            if (thisid.rawId() != thatid.rawId() && thisid.wheel() == thatid.wheel() &&
                thisid.sector() == thatid.sector()) {
              MuonChamberResidual *dt13other = mrft.chamberResidual(*otherId, MuonChamberResidual::kDT13);
              if (dt13other != nullptr && dt13other->numHits() >= m_minDT13Hits) {
                double slopediff = 1000. * (dt13->global_resslope() - dt13other->global_resslope());
                //double length = dt13->chamberAlignable()->surface().toGlobal(align::LocalPoint(0,0,0)).perp() -
                //                dt13other->chamberAlignable()->surface().toGlobal(align::LocalPoint(0,0,0)).perp();
                //double residdiff = 10. * (dt13->global_residual() + length*dt13->global_resslope() - dt13other->global_residual());
                double residdiff = 10. * (dt13->global_residual() - dt13other->global_residual());

                int st = 0;
                if (thatid.station() - thisid.station() == 1)
                  st = thisid.station();
                if (st > 0) {
                  m_dt13_resid[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(qoverpt, residdiff);
                  m_dt13_slope[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(qoverpt, slopediff);
                  if (qoverpt > 0) {
                    m_posdt13_resid[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(residdiff);
                    m_posdt13_slope[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(slopediff);
                  } else {
                    m_negdt13_resid[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(residdiff);
                    m_negdt13_slope[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(slopediff);
                  }
                }
              }  // end other numhits
            }    // end this near other
          }      // end other is DT

          // cross-system: other is CSC
          // only do it for DT stubs in W+-2 St1&2:
          if (!(abs(thisid.wheel()) == 2 && (thisid.station() == 1 || thisid.station() == 2)))
            continue;
          if (otherId->det() == DetId::Muon && otherId->subdetId() == MuonSubdetId::CSC) {
            CSCDetId thatid(otherId->rawId());
            //only do it for CSC stubs in St1R3 or St2R2:
            if (!((thatid.station() == 1 && thatid.ring() == 3) || (thatid.station() == 2 && thatid.ring() == 2)))
              continue;

            MuonChamberResidual *cscother = mrft.chamberResidual(*otherId, MuonChamberResidual::kCSC);
            if (cscother != nullptr && cscother->numHits() >= m_minCSCHits) {
              // scale to adjust the csc residual size to be comparabe to dt's one
              double csc_scale = dt13->chamberAlignable()
                                     ->surface()
                                     .toGlobal(align::LocalPoint(dt13->trackx(), dt13->tracky(), 0))
                                     .perp() /
                                 cscother->chamberAlignable()
                                     ->surface()
                                     .toGlobal(align::LocalPoint(cscother->trackx(), cscother->tracky(), 0))
                                     .perp();
              double residdiff = 10. * (dt13->global_residual() - cscother->global_residual() * csc_scale);
              if (thisid.station() == 1 && thatid.station() == 1) {
                if (qoverpt > 0)
                  m_x_pos_dt1_csc1_resid[thatid.endcap() - 1][thisid.sector() - 1]->Fill(residdiff);
                else
                  m_x_neg_dt1_csc1_resid[thatid.endcap() - 1][thisid.sector() - 1]->Fill(residdiff);
              } else if (thisid.station() == 1 && thatid.station() == 2) {
                if (qoverpt > 0)
                  m_x_pos_dt1_csc2_resid[thatid.endcap() - 1][thisid.sector() - 1]->Fill(residdiff);
                else
                  m_x_neg_dt1_csc2_resid[thatid.endcap() - 1][thisid.sector() - 1]->Fill(residdiff);
              } else if (thisid.station() == 2 && thatid.station() == 1) {
                if (qoverpt > 0)
                  m_x_pos_dt2_csc1_resid[thatid.endcap() - 1][thisid.sector() - 1]->Fill(residdiff);
                else
                  m_x_neg_dt2_csc1_resid[thatid.endcap() - 1][thisid.sector() - 1]->Fill(residdiff);
              }
            }
          }  // end other is CSC
        }    // end loop over other
      }      // end if DT13

      // z-direction
      if (dt2 != nullptr && dt2->numHits() >= m_minDT2Hits && (dt2->chi2() / double(dt2->ndof())) < 2.0) {
        DTChamberId thisid(chamberId->rawId());
        for (std::vector<DetId>::const_iterator otherId = chamberIds.begin(); otherId != chamberIds.end(); ++otherId) {
          if (otherId->det() == DetId::Muon && otherId->subdetId() == MuonSubdetId::DT) {
            DTChamberId thatid(otherId->rawId());
            if (thisid.rawId() != thatid.rawId() && thisid.wheel() == thatid.wheel() &&
                thisid.sector() == thatid.sector()) {
              MuonChamberResidual *dt2other = mrft.chamberResidual(*otherId, MuonChamberResidual::kDT2);
              if (dt2other != nullptr && dt2other->numHits() >= m_minDT2Hits) {
                double slopediff = 1000. * (dt2->global_resslope() - dt2other->global_resslope());
                //double length = dt2->chamberAlignable()->surface().toGlobal(align::LocalPoint(0,0,0)).perp() -
                //                dt2other->chamberAlignable()->surface().toGlobal(align::LocalPoint(0,0,0)).perp();
                //double residdiff = 10. * (dt2->global_residual() + length*dt2->global_resslope() - dt2other->global_residual());
                double residdiff = 10. * (dt2->global_residual() - dt2other->global_residual());

                int st = 0;
                if (thatid.station() - thisid.station() == 1)
                  st = thisid.station();
                if (st > 0) {
                  m_dt2_resid[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(qoverpt, residdiff);
                  m_dt2_slope[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(qoverpt, slopediff);
                  if (qoverpt > 0) {
                    m_posdt2_resid[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(residdiff);
                    m_posdt2_slope[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(slopediff);
                  } else {
                    m_negdt2_resid[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(residdiff);
                    m_negdt2_slope[thisid.wheel() + 2][thisid.sector() - 1][st - 1]->Fill(slopediff);
                  }
                }
              }  // end other numhits
            }    // end this near other
          }      // end other is DT
        }        // end loop over other
      }          // end if DT2
    }            // end if DT

    // **************** CSC ****************
    else if (m_doCSC && chamberId->subdetId() == MuonSubdetId::CSC) {
      MuonChamberResidual *csc = mrft.chamberResidual(*chamberId, MuonChamberResidual::kCSC);
      if (csc->numHits() >= m_minCSCHits) {
        CSCDetId thisid(chamberId->rawId());
        for (std::vector<DetId>::const_iterator otherId = chamberIds.begin(); otherId != chamberIds.end(); ++otherId) {
          if (otherId->det() == DetId::Muon && otherId->subdetId() == MuonSubdetId::CSC) {
            CSCDetId thatid(otherId->rawId());
            if (thisid.rawId() != thatid.rawId() && thisid.endcap() == thatid.endcap()) {
              MuonChamberResidual *cscother = mrft.chamberResidual(*otherId, MuonChamberResidual::kCSC);
              if (cscother != nullptr && cscother->numHits() >= m_minCSCHits) {
                double slopediff = 1000. * (csc->global_resslope() - cscother->global_resslope());
                //double length = csc->chamberAlignable()->surface().toGlobal(align::LocalPoint(0,0,0)).z() -
                //                cscother->chamberAlignable()->surface().toGlobal(align::LocalPoint(0,0,0)).z();
                //double residdiff = 10. * (csc->global_residual() + length*csc->global_resslope() - cscother->global_residual());
                double residdiff = 10. * (csc->global_residual() - cscother->global_residual());

                int thischamber = thisid.chamber();
                int thisring = thisid.ring();
                if (thisid.station() == 1 && (thisring == 1 || thisring == 4)) {
                  thischamber = (thischamber - 1) / 2 + 1;
                  thisring = 1;
                }

                if (thisring == thatid.ring() && thischamber == thatid.chamber()) {
                  bool inner = (thisring == 1);
                  bool outer = (thisring == 2);
                  int st = 0;
                  if (thatid.station() - thisid.station() == 1 && (inner || thisid.station() < 3))
                    st = thisid.station();

                  if (outer && st > 0) {
                    m_cscouter_resid[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(qoverpz, residdiff);
                    m_cscouter_slope[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(qoverpz, slopediff);
                    if (qoverpz > 0) {
                      m_poscscouter_resid[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(residdiff);
                      m_poscscouter_slope[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(slopediff);
                    } else {
                      m_negcscouter_resid[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(residdiff);
                      m_negcscouter_slope[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(slopediff);
                    }
                  }
                  if (inner && st > 0) {
                    m_cscinner_resid[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(qoverpz, residdiff);
                    m_cscinner_slope[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(qoverpz, slopediff);
                    if (qoverpz > 0) {
                      m_poscscinner_resid[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(residdiff);
                      m_poscscinner_slope[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(slopediff);
                    } else {
                      m_negcscinner_resid[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(residdiff);
                      m_negcscinner_slope[thisid.endcap() - 1][thischamber - 1][st - 1]->Fill(slopediff);
                    }
                  }
                }  // end of same ring&chamber
              }    // end other min numhits
            }      // end this near other
          }        // end other is CSC
        }          // end loop over other

      }  // end if this min numhits
    }    // end if CSC

  }  // end loop over chamberIds
}

DEFINE_EDM_PLUGIN(AlignmentMonitorPluginFactory,
                  AlignmentMonitorSegmentDifferences,
                  "AlignmentMonitorSegmentDifferences");