CreateTrackerAlignmentRcds.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    Alignment/TrackerAlignment
// Class:      CreateIdealTkAlRecords
//
//
// Original Author:  Gregor Mittag
//         Created:  Tue, 26 Apr 2016 09:45:13 GMT
//
//


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

// 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 "CondCore/DBOutputService/interface/PoolDBOutputService.h"

#include "CondFormats/Alignment/interface/Alignments.h"
#include "CondFormats/Alignment/interface/AlignTransform.h"
#include "CondFormats/Alignment/interface/AlignmentErrorsExtended.h"
#include "CondFormats/Alignment/interface/AlignTransformError.h"
#include "CondFormats/Alignment/interface/AlignTransformErrorExtended.h"
#include "CondFormats/Alignment/interface/AlignmentSurfaceDeformations.h"
#include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
#include "CondFormats/AlignmentRecord/interface/TrackerAlignmentErrorExtendedRcd.h"
#include "CondFormats/AlignmentRecord/interface/TrackerSurfaceDeformationRcd.h"
#include "CondFormats/GeometryObjects/interface/PTrackerParameters.h"

#include "DataFormats/TrackerCommon/interface/TrackerTopology.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"

#include "CLHEP/Vector/RotationInterfaces.h"

//
// class declaration
//

class CreateIdealTkAlRecords : public edm::one::EDAnalyzer<>  {
public:
  explicit CreateIdealTkAlRecords(const edm::ParameterSet&);
  ~CreateIdealTkAlRecords();

  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
  static std::string toString(const GeomDetEnumerators::SubDetector&);
  static GeomDetEnumerators::SubDetector toSubDetector(const std::string& sub);
  static std::vector<GeomDetEnumerators::SubDetector>
  toSubDetectors(const std::vector<std::string>& subs);

private:
  virtual void analyze(const edm::Event&, const edm::EventSetup&) override;
  void clearAlignmentInfos();
  std::unique_ptr<TrackerGeometry> retrieveGeometry(const edm::EventSetup&);
  void addAlignmentInfo(const GeomDet&);
  void alignToGT(const edm::EventSetup&);
  void writeToDB();

  // ----------member data ---------------------------
  const std::vector<GeomDetEnumerators::SubDetector> skipSubDetectors_;
  const bool alignToGlobalTag_;
  const bool createReferenceRcd_;
  bool firstEvent_;
  Alignments alignments_;
  AlignmentErrorsExtended alignmentErrors_;
  AlignmentSurfaceDeformations alignmentSurfaceDeformations_;
  std::vector<uint32_t> rawIDs_;
  std::vector<GeomDetEnumerators::SubDetector> subDets_;
};


//
// constructors and destructor
//
CreateIdealTkAlRecords::CreateIdealTkAlRecords(const edm::ParameterSet& iConfig) :
  skipSubDetectors_(toSubDetectors(iConfig.getUntrackedParameter<std::vector<std::string> >("skipSubDetectors"))),
  alignToGlobalTag_(iConfig.getUntrackedParameter<bool>("alignToGlobalTag")),
  createReferenceRcd_(iConfig.getUntrackedParameter<bool>("createReferenceRcd")),
  firstEvent_(true)
{
}


CreateIdealTkAlRecords::~CreateIdealTkAlRecords()
{
}


//
// member functions
//

// ------------ method called for each event  ------------
void
CreateIdealTkAlRecords::analyze(const edm::Event&, const edm::EventSetup& iSetup)
{
  if (firstEvent_) {
    clearAlignmentInfos();
    const auto tracker = retrieveGeometry(iSetup);

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

    for (const auto& det: dets) addAlignmentInfo(*det);
    if (alignToGlobalTag_ && !createReferenceRcd_) alignToGT(iSetup);
    writeToDB();
    firstEvent_ = false;
  }
}


std::string
CreateIdealTkAlRecords::toString(const GeomDetEnumerators::SubDetector& sub)
{
  switch (sub) {
  case GeomDetEnumerators::PixelBarrel: return "PixelBarrel";
  case GeomDetEnumerators::PixelEndcap: return "PixelEndcap";
  case GeomDetEnumerators::TIB:         return "TIB";
  case GeomDetEnumerators::TOB:         return "TOB";
  case GeomDetEnumerators::TID:         return "TID";
  case GeomDetEnumerators::TEC:         return "TEC";
  case GeomDetEnumerators::CSC:         return "CSC";
  case GeomDetEnumerators::DT:          return "DT";
  case GeomDetEnumerators::RPCBarrel:   return "RPCBarrel";
  case GeomDetEnumerators::RPCEndcap:   return "RPCEndcap";
  case GeomDetEnumerators::GEM:         return "GEM";
  case GeomDetEnumerators::ME0:         return "ME0";
  case GeomDetEnumerators::P2OTB:       return "P2OTB";
  case GeomDetEnumerators::P2OTEC:      return "P2OTEC";
  case GeomDetEnumerators::P1PXB:       return "P1PXB";
  case GeomDetEnumerators::P1PXEC:      return "P1PXEC";
  case GeomDetEnumerators::P2PXB:       return "P2PXB";
  case GeomDetEnumerators::P2PXEC:      return "P2PXEC";
  case GeomDetEnumerators::invalidDet:  return "invalidDet";
  default:
    throw cms::Exception("UnknownSubdetector");
  }
}


GeomDetEnumerators::SubDetector
CreateIdealTkAlRecords::toSubDetector(const std::string& sub)
{
  if (sub == "PixelBarrel")      return GeomDetEnumerators::PixelBarrel;
  else if (sub == "PixelEndcap") return GeomDetEnumerators::PixelEndcap;
  else if (sub == "TIB")         return GeomDetEnumerators::TIB;
  else if (sub == "TOB")         return GeomDetEnumerators::TOB;
  else if (sub == "TID")         return GeomDetEnumerators::TID;
  else if (sub == "TEC")         return GeomDetEnumerators::TEC;
  else if (sub == "CSC")         return GeomDetEnumerators::CSC;
  else if (sub == "DT")          return GeomDetEnumerators::DT;
  else if (sub == "RPCBarrel")   return GeomDetEnumerators::RPCBarrel;
  else if (sub == "RPCEndcap")   return GeomDetEnumerators::RPCEndcap;
  else if (sub == "GEM")         return GeomDetEnumerators::GEM;
  else if (sub == "ME0")         return GeomDetEnumerators::ME0;
  else if (sub == "P2OTB")       return GeomDetEnumerators::P2OTB;
  else if (sub == "P2OTEC")      return GeomDetEnumerators::P2OTEC;
  else if (sub == "P1PXB")       return GeomDetEnumerators::P1PXB;
  else if (sub == "P1PXEC")      return GeomDetEnumerators::P1PXEC;
  else if (sub == "P2PXB")       return GeomDetEnumerators::P2PXB;
  else if (sub == "P2PXEC")      return GeomDetEnumerators::P2PXEC;
  else if (sub == "invalidDet")  return GeomDetEnumerators::invalidDet;
  else throw cms::Exception("UnknownSubdetector") << sub;
}


std::vector<GeomDetEnumerators::SubDetector>
CreateIdealTkAlRecords::toSubDetectors(const std::vector<std::string>& subs)
{
  std::vector<GeomDetEnumerators::SubDetector> result;
  for (const auto& sub: subs) result.emplace_back(toSubDetector(sub));
  return result;
}


void
CreateIdealTkAlRecords::clearAlignmentInfos()
{
  alignments_.clear();
  alignmentErrors_.clear();
  alignmentSurfaceDeformations_ = AlignmentSurfaceDeformations{};
  rawIDs_.clear();
}


std::unique_ptr<TrackerGeometry>
CreateIdealTkAlRecords::retrieveGeometry(const edm::EventSetup& iSetup)
{
  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 tTopo = tTopoHandle.product();

  TrackerGeomBuilderFromGeometricDet trackerBuilder;

  return std::unique_ptr<TrackerGeometry> {
    trackerBuilder.build(&(*geometricDet), *ptp, tTopo )};
}


void
CreateIdealTkAlRecords::addAlignmentInfo(const GeomDet& det)
{
  const auto subDetector = toString(det.subDetector());
  const auto& detId = det.geographicalId().rawId();
  const auto& pos = det.position();
  const auto& rot = det.rotation();
  rawIDs_.push_back(detId);
  subDets_.push_back(det.subDetector());

  // TrackerAlignmentRcd entry
  if (createReferenceRcd_) {
    alignments_.m_align.emplace_back(AlignTransform(AlignTransform::Translation(),
                                                    AlignTransform::Rotation(),
                                                    detId));
  } else {
    const AlignTransform::Translation translation(pos.x(), pos.y(), pos.z());
    const AlignTransform::Rotation rotation(
        CLHEP::HepRep3x3(rot.xx(),rot.xy(),rot.xz(),
                         rot.yx(),rot.yy(),rot.yz(),
                         rot.zx(),rot.zy(),rot.zz()));
    const auto& eulerAngles = rotation.eulerAngles();
    LogDebug("Alignment")
      << "============================================================\n"
      << "subdetector: " << subDetector << "\n"
      << "detId:       " << detId << "\n"
      << "------------------------------------------------------------\n"
      << "     x: " << pos.x() << "\n"
      << "     y: " << pos.y() << "\n"
      << "     z: " << pos.z() << "\n"
      << "   phi: " << eulerAngles.phi() << "\n"
      << " theta: " << eulerAngles.theta() << "\n"
      << "   psi: " << eulerAngles.psi() << "\n"
      << "============================================================\n";
    alignments_.m_align.emplace_back(AlignTransform(translation, rotation, detId));
  }

  // TrackerAlignmentErrorExtendedRcd entry
  const AlignTransformError::SymMatrix zeroAPEs(6, 0);
  alignmentErrors_.m_alignError.emplace_back(AlignTransformErrorExtended(zeroAPEs, detId));
}


void
CreateIdealTkAlRecords::alignToGT(const edm::EventSetup& iSetup)
{
  LogDebug("Alignment") << "Aligning to global tag\n";

  edm::ESHandle<Alignments> alignments;
  iSetup.get<TrackerAlignmentRcd>().get(alignments);
  edm::ESHandle<AlignmentErrorsExtended> alignmentErrors;
  iSetup.get<TrackerAlignmentErrorExtendedRcd>().get(alignmentErrors);
  edm::ESHandle<AlignmentSurfaceDeformations> surfaceDeformations;
  iSetup.get<TrackerSurfaceDeformationRcd>().get(surfaceDeformations);

  if (alignments->m_align.size() != alignmentErrors->m_alignError.size())
    throw cms::Exception("GeometryMismatch")
      << "Size mismatch between alignments (size=" << alignments->m_align.size()
      << ") and alignment errors (size=" << alignmentErrors->m_alignError.size()
      << ")";

  std::vector<uint32_t> commonIDs;
  auto itAlignErr = alignmentErrors->m_alignError.cbegin();
  for (auto itAlign = alignments->m_align.cbegin();
       itAlign != alignments->m_align.cend();
       ++itAlign, ++itAlignErr) {
    const auto id = itAlign->rawId();
    auto found = std::find(rawIDs_.cbegin(), rawIDs_.cend(), id);
    if (found != rawIDs_.cend()) {
      if (id != itAlignErr->rawId())
        throw cms::Exception("GeometryMismatch")
          << "DetId mismatch between alignments (rawId=" << id
          << ") and alignment errors (rawId=" << itAlignErr->rawId() << ")";

      const auto index = std::distance(rawIDs_.cbegin(), found);
      if (std::find(skipSubDetectors_.begin(),
                    skipSubDetectors_.end(),
                    subDets_[index]) != skipSubDetectors_.end()) continue;

      if (alignments_.m_align[index].rawId()
          != alignmentErrors_.m_alignError[index].rawId())
        throw cms::Exception("GeometryMismatch")
          << "DetId mismatch between alignments (rawId="
          << alignments_.m_align[index].rawId()
          << ") and alignment errors (rawId="
          << alignmentErrors_.m_alignError[index].rawId() << ")";

      LogDebug("Alignment")
    << "============================================================\n"
    << "\nGeometry content (" << toString(subDets_[index]) << ", "
    << alignments_.m_align[index].rawId() << "):\n"
    << "\tx: " << alignments_.m_align[index].translation().x()
    << "\ty: " << alignments_.m_align[index].translation().y()
    << "\tz: " << alignments_.m_align[index].translation().z()
    << "\tphi: " << alignments_.m_align[index].rotation().phi()
    << "\ttheta: " << alignments_.m_align[index].rotation().theta()
    << "\tpsi: " << alignments_.m_align[index].rotation().psi()
    << "============================================================\n";
      alignments_.m_align[index] = *itAlign;
      alignmentErrors_.m_alignError[index] = *itAlignErr;
      commonIDs.push_back(id);
      LogDebug("Alignment")
    << "============================================================\n"
    << "Global tag content (" << toString(subDets_[index]) << ", "
    << alignments_.m_align[index].rawId() << "):\n"
    << "\tx: " << alignments_.m_align[index].translation().x()
    << "\ty: " << alignments_.m_align[index].translation().y()
    << "\tz: " << alignments_.m_align[index].translation().z()
    << "\tphi: " << alignments_.m_align[index].rotation().phi()
    << "\ttheta: " << alignments_.m_align[index].rotation().theta()
    << "\tpsi: " << alignments_.m_align[index].rotation().psi()
    << "============================================================\n";
    }
  }

  // - surface deformations are stored differently
  //   -> different treatment
  // - the above payloads contain also entries for ideal modules
  // - no entry is created for ideal surfaces
  //   -> size of surface deformation payload does not necessarily match the
  //      size of the other tracker alignment payload
  for (const auto& id: commonIDs) {
    // search for common raw ID in surface deformation items
    auto item = std::find_if(surfaceDeformations->items().cbegin(),
                             surfaceDeformations->items().cend(),
                             [&id](const auto& i) { return i.m_rawId == id; });
    if (item == surfaceDeformations->items().cend()) continue; // not found

    // copy surface deformation item
    const auto index = std::distance(surfaceDeformations->items().cbegin(), item);
    const auto beginEndPair = surfaceDeformations->parameters(index);
    std::vector<align::Scalar> params(beginEndPair.first, beginEndPair.second);
    alignmentSurfaceDeformations_.add(item->m_rawId,
                                      item->m_parametrizationType,
                                      params);
  }
}


void
CreateIdealTkAlRecords::writeToDB()
{
  const auto& since = cond::timeTypeSpecs[cond::runnumber].beginValue;

  edm::Service<cond::service::PoolDBOutputService> poolDb;
  if (!poolDb.isAvailable()) {
    throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
  }

  edm::LogInfo("Alignment")
    << "Writing ideal tracker-alignment records.";
  poolDb->writeOne(&alignments_, since, "TrackerAlignmentRcd");
  poolDb->writeOne(&alignmentErrors_, since, "TrackerAlignmentErrorExtendedRcd");
  poolDb->writeOne(&alignmentSurfaceDeformations_, since, "TrackerSurfaceDeformationRcd");
}


// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
void
CreateIdealTkAlRecords::fillDescriptions(edm::ConfigurationDescriptions& descriptions)
{
  edm::ParameterSetDescription desc;
  desc.setComment("Creates ideal TrackerAlignmentRcd and TrackerAlignmentErrorExtendedRcd "
                  "from the loaded tracker geometry. "
                  "PoolDBOutputService must be set up for these records.");
  desc.addUntracked<bool>("alignToGlobalTag", false);
  desc.addUntracked<std::vector<std::string> >("skipSubDetectors",
                                               std::vector<std::string>{});
  desc.addUntracked<bool>("createReferenceRcd", false);
  descriptions.add("createIdealTkAlRecords", desc);
}


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