MCMisalignmentScaler.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    Alignment/TrackerAlignment
// Class:      MCMisalignmentScaler
//
//
// Original Author:  Gregor Mittag
//         Created:  Tue, 10 Oct 2017 15:49:17 GMT
//
//

// system include files
#include <memory>
#include <unordered_map>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
#include "CalibTracker/Records/interface/SiStripQualityRcd.h"

#include "CondCore/DBOutputService/interface/PoolDBOutputService.h"

#include "CondFormats/Alignment/interface/Alignments.h"
#include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
#include "CondFormats/DataRecord/interface/SiPixelQualityRcd.h"
#include "CondFormats/GeometryObjects/interface/PTrackerParameters.h"
#include "CondFormats/SiPixelObjects/interface/SiPixelQuality.h"

#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"

#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeomBuilderFromGeometricDet.h"
#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "Geometry/Records/interface/TrackerTopologyRcd.h"
#include "Geometry/Records/interface/PTrackerParametersRcd.h"

//
// class declaration
//

class MCMisalignmentScaler : public edm::one::EDAnalyzer<> {
public:
  explicit MCMisalignmentScaler(const edm::ParameterSet&);
  ~MCMisalignmentScaler() override = default;

  static void fillDescriptions(edm::ConfigurationDescriptions&);

private:
  using ScalerMap = std::unordered_map<unsigned int, std::unordered_map<int, double> >;

  void analyze(const edm::Event&, const edm::EventSetup&) override;
  ScalerMap decodeSubDetectors(const edm::VParameterSet&);

  // ----------member data ---------------------------
  const ScalerMap scalers_;
  const bool pullBadModulesToIdeal_;
  const double outlierPullToIdealCut_;
  bool firstEvent_{true};
};

//
// constructors and destructor
//
MCMisalignmentScaler::MCMisalignmentScaler(const edm::ParameterSet& iConfig)
    : scalers_{decodeSubDetectors(iConfig.getParameter<edm::VParameterSet>("scalers"))},
      pullBadModulesToIdeal_{iConfig.getUntrackedParameter<bool>("pullBadModulesToIdeal")},
      outlierPullToIdealCut_{iConfig.getUntrackedParameter<double>("outlierPullToIdealCut")} {}

//
// member functions
//

// ------------ method called for each event  ------------
void MCMisalignmentScaler::analyze(const edm::Event&, const edm::EventSetup& iSetup) {
  if (!firstEvent_)
    return;
  firstEvent_ = false;

  // get handle on bad modules
  edm::ESHandle<SiPixelQuality> pixelModules;
  iSetup.get<SiPixelQualityRcd>().get(pixelModules);
  edm::ESHandle<SiStripQuality> stripModules;
  iSetup.get<SiStripQualityRcd>().get(stripModules);

  // get the tracker geometry
  edm::ESHandle<GeometricDet> geometricDet;
  iSetup.get<IdealGeometryRecord>().get(geometricDet);
  edm::ESHandle<PTrackerParameters> ptp;
  iSetup.get<PTrackerParametersRcd>().get(ptp);
  edm::ESHandle<TrackerTopology> tTopoHandle;
  iSetup.get<TrackerTopologyRcd>().get(tTopoHandle);
  const auto* const topology = tTopoHandle.product();
  TrackerGeomBuilderFromGeometricDet trackerBuilder;
  auto tracker = std::unique_ptr<TrackerGeometry>{trackerBuilder.build(&(*geometricDet), *ptp, topology)};

  auto dets = tracker->dets();
  std::sort(dets.begin(), dets.end(), [](const auto& a, const auto& b) {
    return a->geographicalId().rawId() < b->geographicalId().rawId();
  });

  // get the input alignment
  edm::ESHandle<Alignments> alignments;
  iSetup.get<TrackerAlignmentRcd>().get(alignments);

  if (dets.size() != alignments->m_align.size()) {
    throw cms::Exception("GeometryMismatch") << "Size mismatch between alignments (size=" << alignments->m_align.size()
                                             << ") and ideal geometry (size=" << dets.size() << ")";
  }

  Alignments rescaledAlignments{};
  {
    auto outlierCounter{0};
    auto ideal = dets.cbegin();
    const auto& ideal_end = dets.cend();
    auto misaligned = alignments->m_align.cbegin();
    for (; ideal != ideal_end; ++ideal, ++misaligned) {
      if ((*ideal)->geographicalId().rawId() != misaligned->rawId()) {
        throw cms::Exception("GeometryMismatch") << "Order differs between Dets in alignments ideal geometry.";
      }

      // determine scale factor
      const auto& subDetId = (*ideal)->geographicalId().subdetId();
      auto side = topology->side((*ideal)->geographicalId());
      if (side == 0) {
        switch (subDetId) {
          case PixelSubdetector::PixelBarrel:
            side = 1;  // both sides are treated identical -> pick one of them
            break;
          case StripSubdetector::TIB:
            side = topology->tibSide((*ideal)->geographicalId());
            break;
          case StripSubdetector::TOB:
            side = topology->tobSide((*ideal)->geographicalId());
            break;
          default:
            break;
        }
      }
      auto scaleFactor = scalers_.find(subDetId)->second.find(side)->second;

      if (pullBadModulesToIdeal_ &&
          (pixelModules->IsModuleBad(misaligned->rawId()) || stripModules->IsModuleBad(misaligned->rawId()))) {
        scaleFactor = 0.0;
      }

      auto x_diff = misaligned->translation().x() - (*ideal)->position().x();
      auto y_diff = misaligned->translation().y() - (*ideal)->position().y();
      auto z_diff = misaligned->translation().z() - (*ideal)->position().z();

      auto xx_diff = misaligned->rotation().xx() - (*ideal)->rotation().xx();
      auto xy_diff = misaligned->rotation().xy() - (*ideal)->rotation().xy();
      auto xz_diff = misaligned->rotation().xz() - (*ideal)->rotation().xz();
      auto yx_diff = misaligned->rotation().yx() - (*ideal)->rotation().yx();
      auto yy_diff = misaligned->rotation().yy() - (*ideal)->rotation().yy();
      auto yz_diff = misaligned->rotation().yz() - (*ideal)->rotation().yz();
      auto zx_diff = misaligned->rotation().zx() - (*ideal)->rotation().zx();
      auto zy_diff = misaligned->rotation().zy() - (*ideal)->rotation().zy();
      auto zz_diff = misaligned->rotation().zz() - (*ideal)->rotation().zz();

      if (outlierPullToIdealCut_ > 0.0 &&
          (x_diff * x_diff + y_diff * y_diff + z_diff * z_diff) > outlierPullToIdealCut_ * outlierPullToIdealCut_) {
        ++outlierCounter;
        edm::LogInfo("Alignment") << outlierCounter << ") Outlier found in subdetector " << subDetId
                                  << ":  delta x: " << x_diff << ",  delta y: " << y_diff << ",  delta z: " << z_diff
                                  << ",  delta xx: " << xx_diff << ",  delta xy: " << xy_diff
                                  << ",  delta xz: " << xz_diff << ",  delta yx: " << yx_diff
                                  << ",  delta yx: " << yy_diff << ",  delta yy: " << yz_diff
                                  << ",  delta zz: " << zx_diff << ",  delta zy: " << zy_diff
                                  << ",  delta zz: " << zz_diff << "\n";
        scaleFactor = 0.0;
      }

      const AlignTransform::Translation rescaledTranslation{(*ideal)->position().x() + scaleFactor * x_diff,
                                                            (*ideal)->position().y() + scaleFactor * y_diff,
                                                            (*ideal)->position().z() + scaleFactor * z_diff};

      const AlignTransform::Rotation rescaledRotation{
          CLHEP::HepRep3x3{(*ideal)->rotation().xx() + scaleFactor * xx_diff,
                           (*ideal)->rotation().xy() + scaleFactor * xy_diff,
                           (*ideal)->rotation().xz() + scaleFactor * xz_diff,
                           (*ideal)->rotation().yx() + scaleFactor * yx_diff,
                           (*ideal)->rotation().yy() + scaleFactor * yy_diff,
                           (*ideal)->rotation().yz() + scaleFactor * yz_diff,
                           (*ideal)->rotation().zx() + scaleFactor * zx_diff,
                           (*ideal)->rotation().zy() + scaleFactor * zy_diff,
                           (*ideal)->rotation().zz() + scaleFactor * zz_diff}};

      const AlignTransform rescaledTransform{rescaledTranslation, rescaledRotation, misaligned->rawId()};
      rescaledAlignments.m_align.emplace_back(std::move(rescaledTransform));
    }
  }

  edm::Service<cond::service::PoolDBOutputService> poolDb;
  if (!poolDb.isAvailable()) {
    throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
  }
  edm::LogInfo("Alignment") << "Writing rescaled tracker-alignment record.";
  const auto& since = cond::timeTypeSpecs[cond::runnumber].beginValue;
  poolDb->writeOne(&rescaledAlignments, since, "TrackerAlignmentRcd");
}

MCMisalignmentScaler::ScalerMap MCMisalignmentScaler::decodeSubDetectors(const edm::VParameterSet& psets) {
  // initialize scaler map
  ScalerMap subDetMap;
  for (unsigned int subDetId = 1; subDetId <= 6; ++subDetId) {
    subDetMap[subDetId][1] = 1.0;
    subDetMap[subDetId][2] = 1.0;
  }

  // apply scale factors from configuration
  for (const auto& pset : psets) {
    const auto& name = pset.getUntrackedParameter<std::string>("subDetector");
    const auto& factor = pset.getUntrackedParameter<double>("factor");

    std::vector<int> sides;
    if (name.find("-") != std::string::npos)
      sides.push_back(1);
    if (name.find("+") != std::string::npos)
      sides.push_back(2);
    if (sides.empty()) {  // -> use both sides
      sides.push_back(1);
      sides.push_back(2);
    }

    if (name.find("Tracker") != std::string::npos) {
      for (unsigned int subDetId = 1; subDetId <= 6; ++subDetId) {
        for (const auto& side : sides)
          subDetMap[subDetId][side] *= factor;
      }
      if (sides.size() == 1) {
        // if only one side to be scaled
        // -> scale also the other side for PXB (subdetid = 1)
        subDetMap[PixelSubdetector::PixelBarrel][std::abs(sides[0] - 2) + 1] *= factor;
      }
    } else if (name.find("PXB") != std::string::npos) {
      // ignore sides for PXB
      subDetMap[PixelSubdetector::PixelBarrel][1] *= factor;
      subDetMap[PixelSubdetector::PixelBarrel][2] *= factor;
    } else if (name.find("PXF") != std::string::npos) {
      for (const auto& side : sides)
        subDetMap[PixelSubdetector::PixelEndcap][side] *= factor;
    } else if (name.find("TIB") != std::string::npos) {
      for (const auto& side : sides)
        subDetMap[StripSubdetector::TIB][side] *= factor;
    } else if (name.find("TOB") != std::string::npos) {
      for (const auto& side : sides)
        subDetMap[StripSubdetector::TOB][side] *= factor;
    } else if (name.find("TID") != std::string::npos) {
      for (const auto& side : sides)
        subDetMap[StripSubdetector::TID][side] *= factor;
    } else if (name.find("TEC") != std::string::npos) {
      for (const auto& side : sides)
        subDetMap[StripSubdetector::TEC][side] *= factor;
    } else {
      throw cms::Exception("BadConfig") << "@SUB=MCMisalignmentScaler::decodeSubDetectors\n"
                                        << "Unknown tracker subdetector: " << name
                                        << "\nSupported options: Tracker, PXB, PXF, TIB, TOB, TID, TEC "
                                        << "(possibly decorated with '+' or '-')";
    }
  }

  std::stringstream logInfo;
  logInfo << "MC misalignment scale factors:\n";
  for (const auto& subdet : subDetMap) {
    logInfo << "  Subdet " << subdet.first << "\n";
    for (const auto& side : subdet.second) {
      logInfo << "    side " << side.first << ": " << side.second << "\n";
    }
    logInfo << "\n";
  }
  edm::LogInfo("Alignment") << logInfo.str();

  return subDetMap;
}

// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
void MCMisalignmentScaler::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  desc.setComment(
      "Creates rescaled MC misalignment scenario. "
      "PoolDBOutputService must be set up for 'TrackerAlignmentRcd'.");
  edm::ParameterSetDescription descScaler;
  descScaler.setComment(
      "ParameterSet specifying the tracker part to be scaled "
      "by a given factor.");
  descScaler.addUntracked<std::string>("subDetector", "Tracker");
  descScaler.addUntracked<double>("factor", 1.0);
  desc.addVPSet("scalers", descScaler, std::vector<edm::ParameterSet>(1));
  desc.addUntracked<bool>("pullBadModulesToIdeal", false);
  desc.addUntracked<double>("outlierPullToIdealCut", -1.0);
  descriptions.add("mcMisalignmentScaler", desc);
}

//define this as a plug-in
DEFINE_FWK_MODULE(MCMisalignmentScaler);