DTRecHitQuality.cc

Go to the documentation of this file.

/*
 *  See header file for a description of this class.
 *
 *  \author G. Cerminara - INFN Torino
 */

#include <iostream>
#include <map>

#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "Geometry/DTGeometry/interface/DTGeometry.h"
#include "Geometry/DTGeometry/interface/DTLayer.h"
#include "Geometry/Records/interface/MuonGeometryRecord.h"
#include "Validation/DTRecHits/interface/DTHitQualityUtils.h"

#include "DTRecHitQuality.h"
#include "Histograms.h"

using namespace std;
using namespace edm;

namespace dtrechit {
  struct Histograms {
    std::unique_ptr<HRes1DHit> hRes_S1RPhi;          // RecHits, 1. step, RPh
    std::unique_ptr<HRes1DHit> hRes_S2RPhi;          // RecHits, 2. step, RPhi
    std::unique_ptr<HRes1DHit> hRes_S3RPhi;          // RecHits, 3. step, RPhi
    std::unique_ptr<HRes1DHit> hRes_S1RZ;            // RecHits, 1. step, RZ
    std::unique_ptr<HRes1DHit> hRes_S2RZ;            // RecHits, 2. step, RZ
    std::unique_ptr<HRes1DHit> hRes_S3RZ;            // RecHits, 3. step, RZ
    std::unique_ptr<HRes1DHit> hRes_S1RZ_W0;         // RecHits, 1. step, RZ, wheel 0
    std::unique_ptr<HRes1DHit> hRes_S2RZ_W0;         // RecHits, 2. step, RZ, wheel 0
    std::unique_ptr<HRes1DHit> hRes_S3RZ_W0;         // RecHits, 3. step, RZ, wheel 0
    std::unique_ptr<HRes1DHit> hRes_S1RZ_W1;         // RecHits, 1. step, RZ, wheel +-1
    std::unique_ptr<HRes1DHit> hRes_S2RZ_W1;         // RecHits, 2. step, RZ, wheel +-1
    std::unique_ptr<HRes1DHit> hRes_S3RZ_W1;         // RecHits, 3. step, RZ, wheel +-1
    std::unique_ptr<HRes1DHit> hRes_S1RZ_W2;         // RecHits, 1. step, RZ, wheel +-2
    std::unique_ptr<HRes1DHit> hRes_S2RZ_W2;         // RecHits, 2. step, RZ, wheel +-2
    std::unique_ptr<HRes1DHit> hRes_S3RZ_W2;         // RecHits, 3. step, RZ, wheel +-2
    std::unique_ptr<HRes1DHit> hRes_S1RPhi_W0;       // RecHits, 1. step, RPhi, wheel 0
    std::unique_ptr<HRes1DHit> hRes_S2RPhi_W0;       // RecHits, 2. step, RPhi, wheel 0
    std::unique_ptr<HRes1DHit> hRes_S3RPhi_W0;       // RecHits, 3. step, RPhi, wheel 0
    std::unique_ptr<HRes1DHit> hRes_S1RPhi_W1;       // RecHits, 1. step, RPhi, wheel +-1
    std::unique_ptr<HRes1DHit> hRes_S2RPhi_W1;       // RecHits, 2. step, RPhi, wheel +-1
    std::unique_ptr<HRes1DHit> hRes_S3RPhi_W1;       // RecHits, 3. step, RPhi, wheel +-1
    std::unique_ptr<HRes1DHit> hRes_S1RPhi_W2;       // RecHits, 1. step, RPhi, wheel +-2
    std::unique_ptr<HRes1DHit> hRes_S2RPhi_W2;       // RecHits, 2. step, RPhi, wheel +-2
    std::unique_ptr<HRes1DHit> hRes_S3RPhi_W2;       // RecHits, 3. step, RPhi, wheel +-2
    std::unique_ptr<HRes1DHit> hRes_S3RPhiWS[3][4];  // RecHits, 3. step, by wheel/station
    std::unique_ptr<HRes1DHit> hRes_S3RZWS[3][4];    // RecHits, 3. step, by wheel/station

    std::unique_ptr<HEff1DHit> hEff_S1RPhi;          // RecHits, 1. step, RPhi
    std::unique_ptr<HEff1DHit> hEff_S2RPhi;          // RecHits, 2. step, RPhi
    std::unique_ptr<HEff1DHit> hEff_S3RPhi;          // RecHits, 3. step, RPhi
    std::unique_ptr<HEff1DHit> hEff_S1RZ;            // RecHits, 1. step, RZ
    std::unique_ptr<HEff1DHit> hEff_S2RZ;            // RecHits, 2. step, RZ
    std::unique_ptr<HEff1DHit> hEff_S3RZ;            // RecHits, 3. step, RZ
    std::unique_ptr<HEff1DHit> hEff_S1RZ_W0;         // RecHits, 1. step, RZ, wheel 0
    std::unique_ptr<HEff1DHit> hEff_S2RZ_W0;         // RecHits, 2. step, RZ, wheel 0
    std::unique_ptr<HEff1DHit> hEff_S3RZ_W0;         // RecHits, 3. step, RZ, wheel 0
    std::unique_ptr<HEff1DHit> hEff_S1RZ_W1;         // RecHits, 1. step, RZ, wheel +-1
    std::unique_ptr<HEff1DHit> hEff_S2RZ_W1;         // RecHits, 2. step, RZ, wheel +-1
    std::unique_ptr<HEff1DHit> hEff_S3RZ_W1;         // RecHits, 3. step, RZ, wheel +-1
    std::unique_ptr<HEff1DHit> hEff_S1RZ_W2;         // RecHits, 1. step, RZ, wheel +-2
    std::unique_ptr<HEff1DHit> hEff_S2RZ_W2;         // RecHits, 2. step, RZ, wheel +-2
    std::unique_ptr<HEff1DHit> hEff_S3RZ_W2;         // RecHits, 3. step, RZ, wheel +-2
    std::unique_ptr<HEff1DHit> hEff_S1RPhiWS[3][4];  // RecHits, 3. step, by wheel/station
    std::unique_ptr<HEff1DHit> hEff_S3RPhiWS[3][4];  // RecHits, 3. step, by wheel/station
    std::unique_ptr<HEff1DHit> hEff_S1RZWS[3][4];    // RecHits, 3. step, by wheel/station
    std::unique_ptr<HEff1DHit> hEff_S3RZWS[3][4];    // RecHits, 3. step, by wheel/station
  };
}  // namespace dtrechit

using namespace dtrechit;

// In phi SLs, The dependency on X and angle is specular in positive
// and negative wheels. Since positive and negative wheels are filled
// together into the same plots, it is useful to mirror negative wheels
// so that the actual dependency can be observerd instead of an artificially
// simmetrized one.
// Set mirrorMinusWheels to avoid this.
namespace {
  constexpr bool mirrorMinusWheels = true;
}

// Constructor
DTRecHitQuality::DTRecHitQuality(const ParameterSet &pset) {
  // Get the debug parameter for verbose output
  debug_ = pset.getUntrackedParameter<bool>("debug");
  // the name of the simhit collection
  simHitLabel_ = pset.getUntrackedParameter<InputTag>("simHitLabel");
  simHitToken_ = consumes<PSimHitContainer>(pset.getUntrackedParameter<InputTag>("simHitLabel"));
  // the name of the 1D rec hit collection
  recHitLabel_ = pset.getUntrackedParameter<InputTag>("recHitLabel");
  recHitToken_ = consumes<DTRecHitCollection>(pset.getUntrackedParameter<InputTag>("recHitLabel"));
  // the name of the 2D rec hit collection
  segment2DLabel_ = pset.getUntrackedParameter<InputTag>("segment2DLabel");
  segment2DToken_ = consumes<DTRecSegment2DCollection>(pset.getUntrackedParameter<InputTag>("segment2DLabel"));
  // the name of the 4D rec hit collection
  segment4DLabel_ = pset.getUntrackedParameter<InputTag>("segment4DLabel");
  segment4DToken_ = consumes<DTRecSegment4DCollection>(pset.getUntrackedParameter<InputTag>("segment4DLabel"));
  // Switches for analysis at various steps
  doStep1_ = pset.getUntrackedParameter<bool>("doStep1", false);
  doStep2_ = pset.getUntrackedParameter<bool>("doStep2", false);
  doStep3_ = pset.getUntrackedParameter<bool>("doStep3", false);
  doall_ = pset.getUntrackedParameter<bool>("doall", false);
  local_ = pset.getUntrackedParameter<bool>("local", true);
}

void DTRecHitQuality::bookHistograms(DQMStore::ConcurrentBooker &booker,
                                     edm::Run const &run,
                                     edm::EventSetup const &setup,
                                     Histograms &histograms) const {
  if (doall_ && doStep1_) {
    histograms.hRes_S1RPhi = std::make_unique<HRes1DHit>("S1RPhi", booker, true, local_);  // RecHits, 1. step, RPhi
    histograms.hRes_S1RPhi_W0 =
        std::make_unique<HRes1DHit>("S1RPhi_W0", booker, true, local_);  // RecHits, 1. step, RZ, wheel 0
    histograms.hRes_S1RPhi_W1 =
        std::make_unique<HRes1DHit>("S1RPhi_W1", booker, true, local_);  // RecHits, 1. step, RZ, wheel +-1
    histograms.hRes_S1RPhi_W2 =
        std::make_unique<HRes1DHit>("S1RPhi_W2", booker, true, local_);  // RecHits, 1. step, RZ, wheel +-2
    histograms.hRes_S1RZ = std::make_unique<HRes1DHit>("S1RZ", booker, true, local_);  // RecHits, 1. step, RZ
    histograms.hRes_S1RZ_W0 =
        std::make_unique<HRes1DHit>("S1RZ_W0", booker, true, local_);  // RecHits, 1. step, RZ, wheel 0
    histograms.hRes_S1RZ_W1 =
        std::make_unique<HRes1DHit>("S1RZ_W1", booker, true, local_);  // RecHits, 1. step, RZ, wheel +-1
    histograms.hRes_S1RZ_W2 =
        std::make_unique<HRes1DHit>("S1RZ_W2", booker, true, local_);          // RecHits, 1. step, RZ, wheel +-2
    histograms.hEff_S1RPhi = std::make_unique<HEff1DHit>("S1RPhi", booker);    // RecHits, 1. step, RPhi
    histograms.hEff_S1RZ = std::make_unique<HEff1DHit>("S1RZ", booker);        // RecHits, 1. step, RZ
    histograms.hEff_S1RZ_W0 = std::make_unique<HEff1DHit>("S1RZ_W0", booker);  // RecHits, 1. step, RZ, wheel 0
    histograms.hEff_S1RZ_W1 = std::make_unique<HEff1DHit>("S1RZ_W1", booker);  // RecHits, 1. step, RZ, wheel +-1
    histograms.hEff_S1RZ_W2 = std::make_unique<HEff1DHit>("S1RZ_W2", booker);  // RecHits, 1. step, RZ, wheel +-2
  }
  if (doall_ && doStep2_) {
    histograms.hRes_S2RPhi = std::make_unique<HRes1DHit>("S2RPhi", booker, true, local_);  // RecHits, 2. step, RPhi
    histograms.hRes_S2RPhi_W0 =
        std::make_unique<HRes1DHit>("S2RPhi_W0", booker, true, local_);  // RecHits, 2. step, RPhi, wheel 0
    histograms.hRes_S2RPhi_W1 =
        std::make_unique<HRes1DHit>("S2RPhi_W1", booker, true, local_);  // RecHits, 2. step, RPhi, wheel +-1
    histograms.hRes_S2RPhi_W2 =
        std::make_unique<HRes1DHit>("S2RPhi_W2", booker, true, local_);  // RecHits, 2. step, RPhi, wheel +-2
    histograms.hRes_S2RZ = std::make_unique<HRes1DHit>("S2RZ", booker, true, local_);  // RecHits, 2. step, RZ
    histograms.hRes_S2RZ_W0 =
        std::make_unique<HRes1DHit>("S2RZ_W0", booker, true, local_);  // RecHits, 2. step, RZ, wheel 0
    histograms.hRes_S2RZ_W1 =
        std::make_unique<HRes1DHit>("S2RZ_W1", booker, true, local_);  // RecHits, 2. step, RZ, wheel +-1
    histograms.hRes_S2RZ_W2 =
        std::make_unique<HRes1DHit>("S2RZ_W2", booker, true, local_);          // RecHits, 2. step, RZ, wheel +-2
    histograms.hEff_S2RPhi = std::make_unique<HEff1DHit>("S2RPhi", booker);    // RecHits, 2. step, RPhi
    histograms.hEff_S2RZ_W0 = std::make_unique<HEff1DHit>("S2RZ_W0", booker);  // RecHits, 2. step, RZ, wheel 0
    histograms.hEff_S2RZ_W1 = std::make_unique<HEff1DHit>("S2RZ_W1", booker);  // RecHits, 2. step, RZ, wheel +-1
    histograms.hEff_S2RZ_W2 = std::make_unique<HEff1DHit>("S2RZ_W2", booker);  // RecHits, 2. step, RZ, wheel +-2
    histograms.hEff_S2RZ = std::make_unique<HEff1DHit>("S2RZ", booker);        // RecHits, 2. step, RZ
  }
  if (doStep3_) {
    histograms.hRes_S3RPhi = std::make_unique<HRes1DHit>("S3RPhi", booker, doall_, local_);  // RecHits, 3. step, RPhi
    histograms.hRes_S3RPhi_W0 =
        std::make_unique<HRes1DHit>("S3RPhi_W0", booker, doall_, local_);  // RecHits, 3. step, RPhi, wheel 0
    histograms.hRes_S3RPhi_W1 = std::make_unique<HRes1DHit>("S3RPhi_W1",
                                                            booker,
                                                            doall_,
                                                            local_);  // RecHits, 3. step, RPhi, wheel +-1
    histograms.hRes_S3RPhi_W2 = std::make_unique<HRes1DHit>("S3RPhi_W2",
                                                            booker,
                                                            doall_,
                                                            local_);  // RecHits, 3. step, RPhi, wheel +-2
    histograms.hRes_S3RZ = std::make_unique<HRes1DHit>("S3RZ", booker, doall_, local_);  // RecHits, 3. step, RZ
    histograms.hRes_S3RZ_W0 =
        std::make_unique<HRes1DHit>("S3RZ_W0", booker, doall_, local_);  // RecHits, 3. step, RZ, wheel 0
    histograms.hRes_S3RZ_W1 =
        std::make_unique<HRes1DHit>("S3RZ_W1", booker, doall_, local_);  // RecHits, 3. step, RZ, wheel +-1
    histograms.hRes_S3RZ_W2 =
        std::make_unique<HRes1DHit>("S3RZ_W2", booker, doall_, local_);  // RecHits, 3. step, RZ, wheel +-2

    if (local_) {
      // Plots with finer granularity, not to be included in DQM
      TString name1 = "RPhi_W";
      TString name2 = "RZ_W";
      for (long w = 0; w <= 2; ++w) {
        for (long s = 1; s <= 4; ++s) {
          histograms.hRes_S3RPhiWS[w][s - 1] =
              std::make_unique<HRes1DHit>(("S3" + name1 + w + "_St" + s).Data(), booker, doall_, local_);
          histograms.hEff_S1RPhiWS[w][s - 1] =
              std::make_unique<HEff1DHit>(("S1" + name1 + w + "_St" + s).Data(), booker);
          histograms.hEff_S3RPhiWS[w][s - 1] =
              std::make_unique<HEff1DHit>(("S3" + name1 + w + "_St" + s).Data(), booker);
          if (s != 4) {
            histograms.hRes_S3RZWS[w][s - 1] =
                std::make_unique<HRes1DHit>(("S3" + name2 + w + "_St" + s).Data(), booker, doall_, local_);
            histograms.hEff_S1RZWS[w][s - 1] =
                std::make_unique<HEff1DHit>(("S1" + name2 + w + "_St" + s).Data(), booker);
            histograms.hEff_S3RZWS[w][s - 1] =
                std::make_unique<HEff1DHit>(("S3" + name2 + w + "_St" + s).Data(), booker);
          }
        }
      }
    }

    if (doall_) {
      histograms.hEff_S3RPhi = std::make_unique<HEff1DHit>("S3RPhi", booker);    // RecHits, 3. step, RPhi
      histograms.hEff_S3RZ = std::make_unique<HEff1DHit>("S3RZ", booker);        // RecHits, 3. step, RZ
      histograms.hEff_S3RZ_W0 = std::make_unique<HEff1DHit>("S3RZ_W0", booker);  // RecHits, 3. step, RZ, wheel 0
      histograms.hEff_S3RZ_W1 = std::make_unique<HEff1DHit>("S3RZ_W1", booker);  // RecHits, 3. step, RZ, wheel +-1
      histograms.hEff_S3RZ_W2 = std::make_unique<HEff1DHit>("S3RZ_W2", booker);  // RecHits, 3. step, RZ, wheel +-2
    }
  }
}

// The real analysis
void DTRecHitQuality::dqmAnalyze(edm::Event const &event,
                                 edm::EventSetup const &setup,
                                 Histograms const &histograms) const {
  if (debug_) {
    cout << "--- [DTRecHitQuality] Analysing Event: #Run: " << event.id().run() << " #Event: " << event.id().event()
         << endl;
  }

  // Get the DT Geometry
  ESHandle<DTGeometry> dtGeom;
  setup.get<MuonGeometryRecord>().get(dtGeom);

  // Get the SimHit collection from the event
  Handle<PSimHitContainer> simHits;
  event.getByToken(simHitToken_, simHits);

  // Map simhits per wire
  map<DTWireId, PSimHitContainer> simHitsPerWire = DTHitQualityUtils::mapSimHitsPerWire(*(simHits.product()));

  //=======================================================================================
  // RecHit analysis at Step 1
  if (doStep1_ && doall_) {
    if (debug_) {
      cout << "  -- DTRecHit S1: begin analysis:" << endl;
    }
    // Get the rechit collection from the event
    Handle<DTRecHitCollection> dtRecHits;
    event.getByToken(recHitToken_, dtRecHits);

    if (!dtRecHits.isValid()) {
      if (debug_) {
        cout << "[DTRecHitQuality]**Warning: no 1DRechits with label: " << recHitLabel_ << " in this event, skipping!"
             << endl;
      }
      return;
    }

    // Map rechits per wire
    auto const &recHitsPerWire = map1DRecHitsPerWire(dtRecHits.product());
    compute(dtGeom.product(), simHitsPerWire, recHitsPerWire, histograms, 1);
  }

  //=======================================================================================
  // RecHit analysis at Step 2
  if (doStep2_ && doall_) {
    if (debug_) {
      cout << "  -- DTRecHit S2: begin analysis:" << endl;
    }

    // Get the 2D rechits from the event
    Handle<DTRecSegment2DCollection> segment2Ds;
    event.getByToken(segment2DToken_, segment2Ds);

    if (!segment2Ds.isValid()) {
      if (debug_) {
        cout << "[DTRecHitQuality]**Warning: no 2DSegments with label: " << segment2DLabel_
             << " in this event, skipping!" << endl;
      }

    } else {
      // Map rechits per wire
      auto const &recHitsPerWire = map1DRecHitsPerWire(segment2Ds.product());
      compute(dtGeom.product(), simHitsPerWire, recHitsPerWire, histograms, 2);
    }
  }

  //=======================================================================================
  // RecHit analysis at Step 3
  if (doStep3_) {
    if (debug_) {
      cout << "  -- DTRecHit S3: begin analysis:" << endl;
    }

    // Get the 4D rechits from the event
    Handle<DTRecSegment4DCollection> segment4Ds;
    event.getByToken(segment4DToken_, segment4Ds);

    if (!segment4Ds.isValid()) {
      if (debug_) {
        cout << "[DTRecHitQuality]**Warning: no 4D Segments with label: " << segment4DLabel_
             << " in this event, skipping!" << endl;
      }
      return;
    }

    // Map rechits per wire
    auto const &recHitsPerWire = map1DRecHitsPerWire(segment4Ds.product());
    compute(dtGeom.product(), simHitsPerWire, recHitsPerWire, histograms, 3);
  }
}

// Return a map between DTRecHit1DPair and wireId
map<DTWireId, vector<DTRecHit1DPair>> DTRecHitQuality::map1DRecHitsPerWire(
    const DTRecHitCollection *dt1DRecHitPairs) const {
  map<DTWireId, vector<DTRecHit1DPair>> ret;

  for (const auto &dt1DRecHitPair : *dt1DRecHitPairs) {
    ret[dt1DRecHitPair.wireId()].push_back(dt1DRecHitPair);
  }

  return ret;
}

// Return a map between DTRecHit1D at S2 and wireId
map<DTWireId, vector<DTRecHit1D>> DTRecHitQuality::map1DRecHitsPerWire(
    const DTRecSegment2DCollection *segment2Ds) const {
  map<DTWireId, vector<DTRecHit1D>> ret;

  // Loop over all 2D segments
  for (const auto &segment2D : *segment2Ds) {
    vector<DTRecHit1D> component1DHits = segment2D.specificRecHits();
    // Loop over all component 1D hits
    for (auto &component1DHit : component1DHits) {
      ret[component1DHit.wireId()].push_back(component1DHit);
    }
  }
  return ret;
}

// Return a map between DTRecHit1D at S3 and wireId
map<DTWireId, std::vector<DTRecHit1D>> DTRecHitQuality::map1DRecHitsPerWire(
    const DTRecSegment4DCollection *segment4Ds) const {
  map<DTWireId, vector<DTRecHit1D>> ret;
  // Loop over all 4D segments
  for (const auto &segment4D : *segment4Ds) {
    // Get component 2D segments
    vector<const TrackingRecHit *> segment2Ds = segment4D.recHits();
    // Loop over 2D segments:
    for (auto &segment2D : segment2Ds) {
      // Get 1D component rechits
      vector<const TrackingRecHit *> hits = segment2D->recHits();
      // Loop over them
      for (auto &hit : hits) {
        const auto *hit1D = dynamic_cast<const DTRecHit1D *>(hit);
        ret[hit1D->wireId()].push_back(*hit1D);
      }
    }
  }

  return ret;
}

// Compute SimHit distance from wire (cm)
float DTRecHitQuality::simHitDistFromWire(const DTLayer *layer, const DTWireId &wireId, const PSimHit &hit) const {
  float xwire = layer->specificTopology().wirePosition(wireId.wire());
  LocalPoint entryP = hit.entryPoint();
  LocalPoint exitP = hit.exitPoint();
  float xEntry = entryP.x() - xwire;
  float xExit = exitP.x() - xwire;

  return fabs(xEntry - (entryP.z() * (xExit - xEntry)) / (exitP.z() - entryP.z()));  // FIXME: check...
}

// Compute SimHit impact angle (in direction perp to wire), in the SL RF
float DTRecHitQuality::simHitImpactAngle(const DTLayer *layer, const DTWireId &wireId, const PSimHit &hit) const {
  LocalPoint entryP = hit.entryPoint();
  LocalPoint exitP = hit.exitPoint();
  float theta = (exitP.x() - entryP.x()) / (exitP.z() - entryP.z());
  return atan(theta);
}

// Compute SimHit distance from FrontEnd
float DTRecHitQuality::simHitDistFromFE(const DTLayer *layer, const DTWireId &wireId, const PSimHit &hit) const {
  LocalPoint entryP = hit.entryPoint();
  LocalPoint exitP = hit.exitPoint();
  float wireLenght = layer->specificTopology().cellLenght();
  // FIXME: should take only wireLenght/2.;
  // moreover, pos+cellLenght/2. is shorter than the distance from FE.
  // In fact it would make more sense to make plots vs y.
  return (entryP.y() + exitP.y()) / 2. + wireLenght;
}

// Find the RecHit closest to the muon SimHit
template <typename type>
const type *DTRecHitQuality::findBestRecHit(const DTLayer *layer,
                                            const DTWireId &wireId,
                                            const vector<type> &recHits,
                                            const float simHitDist) const {
  float res = 99999;
  const type *theBestRecHit = nullptr;
  // Loop over RecHits within the cell
  for (auto recHit = recHits.begin(); recHit != recHits.end(); recHit++) {
    float distTmp = recHitDistFromWire(*recHit, layer);
    if (fabs(distTmp - simHitDist) < res) {
      res = fabs(distTmp - simHitDist);
      theBestRecHit = &(*recHit);
    }
  }  // End of loop over RecHits within the cell

  return theBestRecHit;
}

// Compute the distance from wire (cm) of a hits in a DTRecHit1DPair
float DTRecHitQuality::recHitDistFromWire(const DTRecHit1DPair &hitPair, const DTLayer *layer) const {
  // Compute the rechit distance from wire
  return fabs(hitPair.localPosition(DTEnums::Left).x() - hitPair.localPosition(DTEnums::Right).x()) / 2.;
}

// Compute the distance from wire (cm) of a hits in a DTRecHit1D
float DTRecHitQuality::recHitDistFromWire(const DTRecHit1D &recHit, const DTLayer *layer) const {
  return fabs(recHit.localPosition().x() - layer->specificTopology().wirePosition(recHit.wireId().wire()));
}

template <typename type>
void DTRecHitQuality::compute(const DTGeometry *dtGeom,
                              const std::map<DTWireId, std::vector<PSimHit>> &simHitsPerWire,
                              const std::map<DTWireId, std::vector<type>> &recHitsPerWire,
                              Histograms const &histograms,
                              int step) const {
  // Loop over cells with a muon SimHit
  for (const auto &wireAndSHits : simHitsPerWire) {
    DTWireId wireId = wireAndSHits.first;
    int wheel = wireId.wheel();
    int sl = wireId.superLayer();

    vector<PSimHit> simHitsInCell = wireAndSHits.second;

    // Get the layer
    const DTLayer *layer = dtGeom->layer(wireId);

    // Look for a mu hit in the cell
    const PSimHit *muSimHit = DTHitQualityUtils::findMuSimHit(simHitsInCell);
    if (muSimHit == nullptr) {
      if (debug_) {
        cout << "   No mu SimHit in channel: " << wireId << ", skipping! " << endl;
      }
      continue;  // Skip this cell
    }

    // Find the distance of the simhit from the wire
    float simHitWireDist = simHitDistFromWire(layer, wireId, *muSimHit);
    // Skip simhits out of the cell
    if (simHitWireDist > 2.1) {
      if (debug_) {
        cout << "  [DTRecHitQuality]###Warning: The mu SimHit in out of the "
                "cell, skipping!"
             << endl;
      }
      continue;  // Skip this cell
    }
    GlobalPoint simHitGlobalPos = layer->toGlobal(muSimHit->localPosition());

    // find SH impact angle
    float simHitTheta = simHitImpactAngle(layer, wireId, *muSimHit);

    // find SH distance from FE
    float simHitFEDist = simHitDistFromFE(layer, wireId, *muSimHit);

    bool recHitReconstructed = false;

    // Look for RecHits in the same cell
    if (recHitsPerWire.find(wireId) == recHitsPerWire.end()) {
      // No RecHit found in this cell
      if (debug_) {
        cout << "   No RecHit found at Step: " << step << " in cell: " << wireId << endl;
      }
    } else {
      recHitReconstructed = true;
      // vector<type> recHits = (*wireAndRecHits).second;
      const vector<type> &recHits = recHitsPerWire.at(wireId);
      if (debug_) {
        cout << "   " << recHits.size() << " RecHits, Step " << step << " in channel: " << wireId << endl;
      }

      // Find the best RecHit
      const type *theBestRecHit = findBestRecHit(layer, wireId, recHits, simHitWireDist);

      float recHitWireDist = recHitDistFromWire(*theBestRecHit, layer);
      if (debug_) {
        cout << "    SimHit distance from wire: " << simHitWireDist << endl
             << "    SimHit distance from FE:   " << simHitFEDist << endl
             << "    SimHit angle in layer RF:  " << simHitTheta << endl
             << "    RecHit distance from wire: " << recHitWireDist << endl;
      }
      float recHitErr = recHitPositionError(*theBestRecHit);
      HRes1DHit *hRes = nullptr;
      HRes1DHit *hResTot = nullptr;

      // Mirror angle in phi so that + and - wheels can be plotted together
      if (mirrorMinusWheels && wheel < 0 && sl != 2) {
        simHitTheta *= -1.;
        // Note: local X, if used, would have to be mirrored as well
      }

      // Fill residuals and pulls
      // Select the histo to be filled
      if (step == 1) {
        // Step 1
        if (sl != 2) {
          hResTot = histograms.hRes_S1RPhi.get();
          if (wheel == 0) {
            hRes = histograms.hRes_S1RPhi_W0.get();
          }
          if (abs(wheel) == 1) {
            hRes = histograms.hRes_S1RPhi_W1.get();
          }
          if (abs(wheel) == 2) {
            hRes = histograms.hRes_S1RPhi_W2.get();
          }
        } else {
          hResTot = histograms.hRes_S1RZ.get();
          if (wheel == 0) {
            hRes = histograms.hRes_S1RZ_W0.get();
          }
          if (abs(wheel) == 1) {
            hRes = histograms.hRes_S1RZ_W1.get();
          }
          if (abs(wheel) == 2) {
            hRes = histograms.hRes_S1RZ_W2.get();
          }
        }

      } else if (step == 2) {
        // Step 2
        if (sl != 2) {
          hRes = histograms.hRes_S2RPhi.get();
          if (wheel == 0) {
            hRes = histograms.hRes_S2RPhi_W0.get();
          }
          if (abs(wheel) == 1) {
            hRes = histograms.hRes_S2RPhi_W1.get();
          }
          if (abs(wheel) == 2) {
            hRes = histograms.hRes_S2RPhi_W2.get();
          }
        } else {
          hResTot = histograms.hRes_S2RZ.get();
          if (wheel == 0) {
            hRes = histograms.hRes_S2RZ_W0.get();
          }
          if (abs(wheel) == 1) {
            hRes = histograms.hRes_S2RZ_W1.get();
          }
          if (abs(wheel) == 2) {
            hRes = histograms.hRes_S2RZ_W2.get();
          }
        }

      } else if (step == 3) {
        // Step 3
        if (sl != 2) {
          hResTot = histograms.hRes_S3RPhi.get();
          if (wheel == 0) {
            hRes = histograms.hRes_S3RPhi_W0.get();
          }
          if (abs(wheel) == 1) {
            hRes = histograms.hRes_S3RPhi_W1.get();
          }
          if (abs(wheel) == 2) {
            hRes = histograms.hRes_S3RPhi_W2.get();
          }
          if (local_) {
            histograms.hRes_S3RPhiWS[abs(wheel)][wireId.station() - 1]->fill(simHitWireDist,
                                                                             simHitTheta,
                                                                             simHitFEDist,
                                                                             recHitWireDist,
                                                                             simHitGlobalPos.eta(),
                                                                             simHitGlobalPos.phi(),
                                                                             recHitErr,
                                                                             wireId.station());
          }
        } else {
          hResTot = histograms.hRes_S3RZ.get();
          if (wheel == 0) {
            hRes = histograms.hRes_S3RZ_W0.get();
          }
          if (abs(wheel) == 1) {
            hRes = histograms.hRes_S3RZ_W1.get();
          }
          if (abs(wheel) == 2) {
            hRes = histograms.hRes_S3RZ_W2.get();
          }
          if (local_) {
            histograms.hRes_S3RZWS[abs(wheel)][wireId.station() - 1]->fill(simHitWireDist,
                                                                           simHitTheta,
                                                                           simHitFEDist,
                                                                           recHitWireDist,
                                                                           simHitGlobalPos.eta(),
                                                                           simHitGlobalPos.phi(),
                                                                           recHitErr,
                                                                           wireId.station());
          }
        }
      }

      // Fill
      hRes->fill(simHitWireDist,
                 simHitTheta,
                 simHitFEDist,
                 recHitWireDist,
                 simHitGlobalPos.eta(),
                 simHitGlobalPos.phi(),
                 recHitErr,
                 wireId.station());
      if (hResTot != nullptr) {
        hResTot->fill(simHitWireDist,
                      simHitTheta,
                      simHitFEDist,
                      recHitWireDist,
                      simHitGlobalPos.eta(),
                      simHitGlobalPos.phi(),
                      recHitErr,
                      wireId.station());
      }
    }

    // Fill Efficiencies
    if (doall_) {
      HEff1DHit *hEff = nullptr;
      HEff1DHit *hEffTot = nullptr;
      if (step == 1) {
        // Step 1
        if (sl != 2) {
          hEff = histograms.hEff_S1RPhi.get();
          if (local_) {
            histograms.hEff_S1RPhiWS[abs(wheel)][wireId.station() - 1]->fill(
                simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed);
          }
        } else {
          hEffTot = histograms.hEff_S1RZ.get();
          if (wheel == 0) {
            hEff = histograms.hEff_S1RZ_W0.get();
          }
          if (abs(wheel) == 1) {
            hEff = histograms.hEff_S1RZ_W1.get();
          }
          if (abs(wheel) == 2) {
            hEff = histograms.hEff_S1RZ_W2.get();
          }
          if (local_) {
            histograms.hEff_S1RZWS[abs(wheel)][wireId.station() - 1]->fill(
                simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed);
          }
        }

      } else if (step == 2) {
        // Step 2
        if (sl != 2) {
          hEff = histograms.hEff_S2RPhi.get();
        } else {
          hEffTot = histograms.hEff_S2RZ.get();
          if (wheel == 0) {
            hEff = histograms.hEff_S2RZ_W0.get();
          }
          if (abs(wheel) == 1) {
            hEff = histograms.hEff_S2RZ_W1.get();
          }
          if (abs(wheel) == 2) {
            hEff = histograms.hEff_S2RZ_W2.get();
          }
        }

      } else if (step == 3) {
        // Step 3
        if (sl != 2) {
          hEff = histograms.hEff_S3RPhi.get();
          if (local_) {
            histograms.hEff_S3RPhiWS[abs(wheel)][wireId.station() - 1]->fill(
                simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed);
          }
        } else {
          hEffTot = histograms.hEff_S3RZ.get();
          if (wheel == 0) {
            hEff = histograms.hEff_S3RZ_W0.get();
          }
          if (abs(wheel) == 1) {
            hEff = histograms.hEff_S3RZ_W1.get();
          }
          if (abs(wheel) == 2) {
            hEff = histograms.hEff_S3RZ_W2.get();
          }
          if (local_) {
            histograms.hEff_S3RZWS[abs(wheel)][wireId.station() - 1]->fill(
                simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed);
          }
        }
      }
      // Fill
      hEff->fill(simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed);
      if (hEffTot != nullptr) {
        hEffTot->fill(simHitWireDist, simHitGlobalPos.eta(), simHitGlobalPos.phi(), recHitReconstructed);
      }
    }
  }
}

// Return the error on the measured (cm) coordinate
float DTRecHitQuality::recHitPositionError(const DTRecHit1DPair &recHit) const {
  return sqrt(recHit.localPositionError(DTEnums::Left).xx());
}

// Return the error on the measured (cm) coordinate
float DTRecHitQuality::recHitPositionError(const DTRecHit1D &recHit) const {
  return sqrt(recHit.localPositionError().xx());
}

// declare this as a framework plugin
#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(DTRecHitQuality);