Merge pull request #34985 from pieterdavid/addSiStripCalCosmicsNano_1…

…21pre1

Add SiStripCalCosmics ALCANANO

CalibTracker/SiStripCommon/BuildFile.xml

@@ -2,8 +2,11 @@
 <use name="DataFormats/Common"/>
 <use name="DataFormats/FEDRawData"/>
 <use name="DataFormats/SiStripDetId"/>
+<use name="DataFormats/TrackReco"/>
+<use name="DataFormats/TrackerRecHit2D"/>
 <use name="DataFormats/Scalers"/>
 <use name="DataFormats/OnlineMetaData"/>
+<use name="DataFormats/NanoAOD"/>
 <use name="FWCore/Framework"/>
 <use name="FWCore/MessageLogger"/>
 <use name="FWCore/Utilities"/>

CalibTracker/SiStripCommon/interface/SiStripOnTrackClusterTableProducerBase.h

@@ -0,0 +1,56 @@
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include "DataFormats/NanoAOD/interface/FlatTable.h"
+#include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
+#include "DataFormats/TrackReco/interface/Track.h"
+#include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
+
+class SiStripOnTrackClusterTableProducerBase : public edm::stream::EDProducer<> {
+public:
+  explicit SiStripOnTrackClusterTableProducerBase(const edm::ParameterSet& params)
+      : m_name(params.getParameter<std::string>("name")),
+        m_doc(params.existsAs<std::string>("doc") ? params.getParameter<std::string>("doc") : ""),
+        m_extension(params.existsAs<bool>("extension") ? params.getParameter<bool>("extension") : true),
+        m_tracks_token(consumes<edm::View<reco::Track>>(params.getParameter<edm::InputTag>("Tracks"))),
+        m_association_token(consumes<TrajTrackAssociationCollection>(params.getParameter<edm::InputTag>("Tracks"))) {
+    produces<nanoaod::FlatTable>();
+  }
+  ~SiStripOnTrackClusterTableProducerBase() override;
+
+  void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) final;
+
+  struct OnTrackCluster {
+    uint32_t det;
+    const SiStripCluster* cluster;
+    const Trajectory* traj;
+    const reco::Track* track;
+    const TrajectoryMeasurement& measurement;
+    OnTrackCluster(uint32_t detId,
+                   const SiStripCluster* stripCluster,
+                   const Trajectory* trajectory,
+                   const reco::Track* track_,
+                   const TrajectoryMeasurement& measurement_)
+        : det{detId}, cluster{stripCluster}, traj{trajectory}, track{track_}, measurement{measurement_} {}
+  };
+
+  virtual void fillTable(const std::vector<OnTrackCluster>& clusters,
+                         const edm::View<reco::Track>& tracks,
+                         nanoaod::FlatTable* table,
+                         const edm::EventSetup& iSetup) = 0;
+
+  template <typename VALUES>
+  static void addColumn(nanoaod::FlatTable* table, const std::string& name, VALUES&& values, const std::string& doc) {
+    using value_type = typename std::remove_reference<VALUES>::type::value_type;
+    table->template addColumn<value_type>(name, values, doc);
+  }
+
+private:
+  const std::string m_name;
+  const std::string m_doc;
+  bool m_extension;
+
+  const edm::EDGetTokenT<edm::View<reco::Track>> m_tracks_token;
+  const edm::EDGetTokenT<TrajTrackAssociationCollection> m_association_token;
+};

CalibTracker/SiStripCommon/plugins/BuildFile.xml

@@ -9,6 +9,7 @@
 <use name="CondFormats/RunInfo"/>
 <use name="RecoLocalTracker/SiStripClusterizer"/>
 <use name="SimDataFormats/TrackingAnalysis"/>
+<use name="DataFormats/NanoAOD"/>
 <library file="*.cc" name="CalibTrackerSiStripCommonPlugins">
   <flags EDM_PLUGIN="1"/>
 </library>

CalibTracker/SiStripCommon/plugins/SealModules.cc

@@ -20,3 +20,7 @@ DEFINE_FWK_MODULE(ShallowSimhitClustersProducer);
 DEFINE_FWK_MODULE(ShallowTracksProducer);
 DEFINE_FWK_MODULE(ShallowSimTracksProducer);
 DEFINE_FWK_MODULE(ShallowGainCalibration);
+
+#include "PhysicsTools/NanoAOD/interface/SimpleFlatTableProducer.h"
+typedef SimpleFlatTableProducer<reco::Track> SimpleTrackFlatTableProducer;
+DEFINE_FWK_MODULE(SimpleTrackFlatTableProducer);

CalibTracker/SiStripCommon/plugins/SiStripLorentzAngleRunInfoTableProducer.cc

@@ -0,0 +1,107 @@
+#include "FWCore/Framework/interface/global/EDProducer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+#include "DataFormats/NanoAOD/interface/FlatTable.h"
+
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+
+#include "MagneticField/Engine/interface/MagneticField.h"
+#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
+#include "CondFormats/SiStripObjects/interface/SiStripLorentzAngle.h"
+#include "CalibTracker/Records/interface/SiStripDependentRecords.h"
+#include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
+#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
+#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
+
+#include "CalibTracker/SiStripCommon/interface/ShallowTools.h"
+
+class SiStripLorentzAngleRunInfoTableProducer : public edm::global::EDProducer<edm::BeginRunProducer> {
+public:
+  explicit SiStripLorentzAngleRunInfoTableProducer(const edm::ParameterSet& params)
+      : m_name{params.getParameter<std::string>("name")},
+        m_magFieldName{params.getParameter<std::string>("magFieldName")},
+        m_doc{params.existsAs<std::string>("doc") ? params.getParameter<std::string>("doc") : ""},
+        m_tkGeomToken{esConsumes<edm::Transition::BeginRun>()},
+        m_magFieldToken{esConsumes<edm::Transition::BeginRun>()},
+        m_lorentzAngleToken{esConsumes<edm::Transition::BeginRun>()} {
+    produces<nanoaod::FlatTable, edm::Transition::BeginRun>();
+    produces<nanoaod::FlatTable, edm::Transition::BeginRun>("magField");
+  }
+
+  void produce(edm::StreamID, edm::Event&, edm::EventSetup const&) const override {}
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+    edm::ParameterSetDescription desc;
+    desc.add<std::string>("name", "Det");
+    desc.add<std::string>("magFieldName", "magField");
+    desc.add<std::string>("doc", "Run info for the Lorentz angle measurement");
+    descriptions.add("siStripLorentzAngleRunInfoTable", desc);
+  }
+
+  void globalBeginRunProduce(edm::Run& iRun, edm::EventSetup const& iSetup) const override;
+
+private:
+  const std::string m_name, m_magFieldName;
+  const std::string m_doc;
+  edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> m_tkGeomToken;
+  edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> m_magFieldToken;
+  edm::ESGetToken<SiStripLorentzAngle, SiStripLorentzAngleDepRcd> m_lorentzAngleToken;
+};
+
+namespace {
+  template <typename VALUES>
+  void addColumn(nanoaod::FlatTable* table, const std::string& name, VALUES&& values, const std::string& doc) {
+    using value_type = typename std::remove_reference<VALUES>::type::value_type;
+    table->template addColumn<value_type>(name, values, doc);
+  }
+}  // namespace
+
+void SiStripLorentzAngleRunInfoTableProducer::globalBeginRunProduce(edm::Run& iRun,
+                                                                    edm::EventSetup const& iSetup) const {
+  const auto& tkGeom = iSetup.getData(m_tkGeomToken);
+  const auto& magField = iSetup.getData(m_magFieldToken);
+  const auto& lorentzAngle = iSetup.getData(m_lorentzAngleToken);
+  std::vector<uint32_t> c_rawid;
+  std::vector<float> c_globalZofunitlocalY, c_localB, c_BdotY, c_driftx, c_drifty, c_driftz, c_lorentzAngle;
+
+  auto dets = tkGeom.detsTIB();
+  dets.insert(dets.end(), tkGeom.detsTID().begin(), tkGeom.detsTID().end());
+  dets.insert(dets.end(), tkGeom.detsTOB().begin(), tkGeom.detsTOB().end());
+  dets.insert(dets.end(), tkGeom.detsTEC().begin(), tkGeom.detsTEC().end());
+  for (auto det : dets) {
+    auto detid = det->geographicalId().rawId();
+    const StripGeomDetUnit* stripDet = dynamic_cast<const StripGeomDetUnit*>(tkGeom.idToDet(det->geographicalId()));
+    if (stripDet) {
+      c_rawid.push_back(detid);
+      c_globalZofunitlocalY.push_back(stripDet->toGlobal(LocalVector(0, 1, 0)).z());
+      const auto locB = magField.inTesla(stripDet->surface().position());
+      c_localB.push_back(locB.mag());
+      c_BdotY.push_back(stripDet->surface().toLocal(locB).y());
+      const auto drift = shallow::drift(stripDet, magField, lorentzAngle);
+      c_driftx.push_back(drift.x());
+      c_drifty.push_back(drift.y());
+      c_driftz.push_back(drift.z());
+      c_lorentzAngle.push_back(lorentzAngle.getLorentzAngle(detid));
+    }
+  }
+  auto out = std::make_unique<nanoaod::FlatTable>(c_rawid.size(), m_name, false, false);
+  addColumn(out.get(), "rawid", c_rawid, "DetId");
+  addColumn(out.get(), "globalZofunitlocalY", c_globalZofunitlocalY, "z component of a local unit vector along y");
+  addColumn(out.get(), "localB", c_localB, "Local magnitude of the magnetic field");
+  addColumn(out.get(), "BdotY", c_BdotY, "Magnetic field projection on the local y axis");
+  addColumn(out.get(), "driftx", c_driftx, "x component of the drift vector");
+  addColumn(out.get(), "drifty", c_drifty, "y component of the drift vector");
+  addColumn(out.get(), "driftz", c_driftz, "z component of the drift vector");
+  addColumn(out.get(), "lorentzAngle", c_lorentzAngle, "Lorentz angle from database");
+  iRun.put(std::move(out));
+
+  auto out2 = std::make_unique<nanoaod::FlatTable>(1, m_magFieldName, true, false);
+  out2->addColumnValue<float>(
+      "origin", magField.inTesla(GlobalPoint(0, 0, 0)).z(), "z-component of the magnetic field at (0,0,0) in Tesla");
+  iRun.put(std::move(out2), "magField");
+}
+
+#include "FWCore/PluginManager/interface/ModuleDef.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+DEFINE_FWK_MODULE(SiStripLorentzAngleRunInfoTableProducer);

CalibTracker/SiStripCommon/plugins/SiStripPositionCorrectionsTableProducer.cc

@@ -0,0 +1,75 @@
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+
+#include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
+#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
+#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
+
+#include "RecoLocalTracker/SiStripClusterizer/interface/SiStripClusterInfo.h"
+
+#include "CalibTracker/SiStripCommon/interface/SiStripOnTrackClusterTableProducerBase.h"
+
+class SiStripPositionCorrectionsTableProducer : public SiStripOnTrackClusterTableProducerBase {
+public:
+  explicit SiStripPositionCorrectionsTableProducer(const edm::ParameterSet& params)
+      : SiStripOnTrackClusterTableProducerBase(params),
+        m_clusterInfo(consumesCollector()),
+        m_tkGeomToken{esConsumes<>()} {}
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+    edm::ParameterSetDescription desc;
+    desc.add<std::string>("name", "cluster");
+    desc.add<std::string>("doc", "On-track cluster properties for Lorentz angle and backplane correction measurement");
+    desc.add<bool>("extension", false);
+    desc.add<edm::InputTag>("Tracks", edm::InputTag{"generalTracks"});
+    descriptions.add("siStripPositionCorrectionsTable", desc);
+  }
+
+  void fillTable(const std::vector<OnTrackCluster>& clusters,
+                 const edm::View<reco::Track>& tracks,
+                 nanoaod::FlatTable* table,
+                 const edm::EventSetup& iSetup) final;
+
+private:
+  SiStripClusterInfo m_clusterInfo;
+  edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> m_tkGeomToken;
+};
+
+void SiStripPositionCorrectionsTableProducer::fillTable(const std::vector<OnTrackCluster>& clusters,
+                                                        const edm::View<reco::Track>& tracks,
+                                                        nanoaod::FlatTable* table,
+                                                        const edm::EventSetup& iSetup) {
+  const auto& tkGeom = iSetup.getData(m_tkGeomToken);
+  std::vector<uint32_t> c_nstrips;
+  std::vector<float> c_barycenter, c_variance, c_localdirx, c_localdiry, c_localdirz, c_localx, c_rhlocalx,
+      c_rhlocalxerr;
+  for (const auto clus : clusters) {
+    c_nstrips.push_back(clus.cluster->amplitudes().size());
+    m_clusterInfo.setCluster(*clus.cluster, clus.det);
+    c_variance.push_back(m_clusterInfo.variance());
+    const auto& trajState = clus.measurement.updatedState();
+    const auto trackDir = trajState.localDirection();
+    c_localdirx.push_back(trackDir.x());
+    c_localdiry.push_back(trackDir.y());
+    c_localdirz.push_back(trackDir.z());
+    const auto hit = clus.measurement.recHit()->hit();
+    const auto stripDet = dynamic_cast<const StripGeomDetUnit*>(tkGeom.idToDet(hit->geographicalId()));
+    c_barycenter.push_back(stripDet->specificTopology().localPosition(clus.cluster->barycenter()).x());
+    c_localx.push_back(stripDet->toLocal(trajState.globalPosition()).x());
+    c_rhlocalx.push_back(hit->localPosition().x());
+    c_rhlocalxerr.push_back(hit->localPositionError().xx());
+  }
+  addColumn(table, "nstrips", c_nstrips, "cluster width");
+  addColumn(table, "variance", c_variance, "Cluster variance");
+  addColumn(table, "localdirx", c_localdirx, "x component of the local track direction");
+  addColumn(table, "localdiry", c_localdiry, "y component of the local track direction");
+  addColumn(table, "localdirz", c_localdirz, "z component of the local track direction");
+  addColumn(table, "barycenter", c_barycenter, "Cluster barycenter (local x without corrections)");
+  addColumn(table, "localx", c_localx, "Track local x");
+  addColumn(table, "rhlocalx", c_rhlocalx, "RecHit local x");
+  addColumn(table, "rhlocalxerr", c_rhlocalxerr, "RecHit local x uncertainty");
+}
+
+#include "FWCore/PluginManager/interface/ModuleDef.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+DEFINE_FWK_MODULE(SiStripPositionCorrectionsTableProducer);

CalibTracker/SiStripCommon/src/SiStripOnTrackClusterTableProducerBase.cc

@@ -0,0 +1,68 @@
+#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
+#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
+#include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
+
+#include "CalibTracker/SiStripCommon/interface/SiStripOnTrackClusterTableProducerBase.h"
+
+SiStripOnTrackClusterTableProducerBase::~SiStripOnTrackClusterTableProducerBase() {}
+
+namespace {
+  int findTrackIndex(const edm::View<reco::Track>& tracks, const reco::Track* track) {
+    for (auto iTr = tracks.begin(); iTr != tracks.end(); ++iTr) {
+      if (&(*iTr) == track) {
+        return iTr - tracks.begin();
+      }
+    }
+    return -2;
+  }
+}  // namespace
+
+void SiStripOnTrackClusterTableProducerBase::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  edm::Handle<edm::View<reco::Track>> tracks;
+  iEvent.getByToken(m_tracks_token, tracks);
+  edm::Handle<TrajTrackAssociationCollection> trajTrackAssociations;
+  iEvent.getByToken(m_association_token, trajTrackAssociations);
+
+  std::vector<OnTrackCluster> clusters{};
+
+  for (const auto& assoc : *trajTrackAssociations) {
+    const auto traj = assoc.key.get();
+    const auto track = assoc.val.get();
+
+    for (const auto& meas : traj->measurements()) {
+      const auto& trajState = meas.updatedState();
+      if (!trajState.isValid())
+        continue;
+
+      // there can be 2 (stereo module), 1 (no stereo module), or 0 (no strip hit) clusters per measurement
+      const auto trechit = meas.recHit()->hit();
+      const auto simple1d = dynamic_cast<const SiStripRecHit1D*>(trechit);
+      const auto simple = dynamic_cast<const SiStripRecHit2D*>(trechit);
+      const auto matched = dynamic_cast<const SiStripMatchedRecHit2D*>(trechit);
+      if (matched) {
+        clusters.emplace_back(matched->monoId(), &matched->monoCluster(), traj, track, meas);
+        clusters.emplace_back(matched->stereoId(), &matched->stereoCluster(), traj, track, meas);
+      } else if (simple) {
+        clusters.emplace_back(simple->geographicalId().rawId(), simple->cluster().get(), traj, track, meas);
+      } else if (simple1d) {
+        clusters.emplace_back(simple1d->geographicalId().rawId(), simple1d->cluster().get(), traj, track, meas);
+      }
+    }
+  }
+
+  auto out = std::make_unique<nanoaod::FlatTable>(clusters.size(), m_name, false, m_extension);
+  if (!m_extension) {
+    std::vector<int> c_trackindex;
+    c_trackindex.reserve(clusters.size());
+    std::vector<uint32_t> c_rawid;
+    c_rawid.reserve(clusters.size());
+    for (const auto clus : clusters) {
+      c_trackindex.push_back(findTrackIndex(*tracks, clus.track));
+      c_rawid.push_back(clus.det);
+    }
+    addColumn(out.get(), "trackindex", c_trackindex, "Track index");
+    addColumn(out.get(), "rawid", c_rawid, "DetId");
+  }
+  fillTable(clusters, *tracks, out.get(), iSetup);
+  iEvent.put(std::move(out));
+}