SiStripCommon/plugins/SiStripGainCalibTableProducer.cc

0001 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0002 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0003
0004 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0005 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0006 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0007 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0008 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0009
0010 #include "CalibFormats/SiStripObjects/interface/SiStripGain.h"
0011 #include "CalibTracker/Records/interface/SiStripGainRcd.h"
0012
0013 #include "CalibTracker/SiStripCommon/interface/SiStripOnTrackClusterTableProducerBase.h"
0014
0015 class SiStripGainCalibTableProducer : public SiStripOnTrackClusterTableProducerBase {
0016 public:
0017   explicit SiStripGainCalibTableProducer(const edm::ParameterSet& params)
0018       : SiStripOnTrackClusterTableProducerBase(params), m_tkGeomToken{esConsumes<>()}, m_gainToken{esConsumes<>()} {}
0019
0020   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0021     edm::ParameterSetDescription desc;
0022     desc.add<std::string>("name", "cluster");
0023     desc.add<std::string>("doc", "");
0024     desc.add<bool>("extension", false);
0025     desc.add<edm::InputTag>("Tracks", edm::InputTag{"generalTracks"});
0026     descriptions.add("siStripGainCalibTable", desc);
0027   }
0028
0029   void fillTable(const std::vector<OnTrackCluster>& clusters,
0030                  const edm::View<reco::Track>& tracks,
0031                  nanoaod::FlatTable* table,
0032                  const edm::EventSetup& iSetup) final;
0033
0034 private:
0035   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> m_tkGeomToken;
0036   const edm::ESGetToken<SiStripGain, SiStripGainRcd> m_gainToken;
0037
0038   std::map<DetId, double> m_thicknessMap;
0039   double thickness(DetId id, const TrackerGeometry* tGeom);
0040 };
0041
0042 namespace {
0043   bool isFarFromBorder(const TrajectoryStateOnSurface& trajState, uint32_t detId, const TrackerGeometry* tGeom) {
0044     const auto gdu = tGeom->idToDetUnit(detId);
0045     if ((!dynamic_cast<const StripGeomDetUnit*>(gdu)) && (!dynamic_cast<const PixelGeomDetUnit*>(gdu))) {
0046       edm::LogWarning("SiStripGainCalibTableProducer")
0047           << "DetId " << detId << " does not seem to belong to the tracker";
0048       return false;
0049     }
0050     const auto plane = gdu->surface();
0051     const auto trapBounds = dynamic_cast<const TrapezoidalPlaneBounds*>(&plane.bounds());
0052     const auto rectBounds = dynamic_cast<const RectangularPlaneBounds*>(&plane.bounds());
0053
0054     static constexpr double distFromBorder = 1.0;
0055     double halfLength = 0.;
0056     if (trapBounds) {
0057       halfLength = trapBounds->parameters()[3];
0058     } else if (rectBounds) {
0059       halfLength = .5 * gdu->surface().bounds().length();
0060     } else {
0061       return false;
0062     }
0063
0064     const auto pos = trajState.localPosition();
0065     const auto posError = trajState.localError().positionError();
0066     if (std::abs(pos.y()) + posError.yy() >= (halfLength - distFromBorder))
0067       return false;
0068
0069     return true;
0070   }
0071 }  // namespace
0072
0073 double SiStripGainCalibTableProducer::thickness(DetId id, const TrackerGeometry* tGeom) {
0074   const auto it = m_thicknessMap.find(id);
0075   if (m_thicknessMap.end() != it) {
0076     return it->second;
0077   } else {
0078     double detThickness = 1.;
0079     const auto gdu = tGeom->idToDetUnit(id);
0080     const auto isPixel = (dynamic_cast<const PixelGeomDetUnit*>(gdu) != nullptr);
0081     const auto isStrip = (dynamic_cast<const StripGeomDetUnit*>(gdu) != nullptr);
0082     if (!isPixel && !isStrip) {
0083       edm::LogWarning("SiStripGainCalibTableProducer")
0084           << "DetId " << id.rawId() << " doesn't seem to belong to the Tracker";
0085     } else {
0086       detThickness = gdu->surface().bounds().thickness();
0087     }
0088     m_thicknessMap[id] = detThickness;
0089     return detThickness;
0090   }
0091 }
0092
0093 void SiStripGainCalibTableProducer::fillTable(const std::vector<OnTrackCluster>& clusters,
0094                                               const edm::View<reco::Track>& tracks,
0095                                               nanoaod::FlatTable* table,
0096                                               const edm::EventSetup& iSetup) {
0097   edm::ESHandle<TrackerGeometry> tGeom = iSetup.getHandle(m_tkGeomToken);
0098   edm::ESHandle<SiStripGain> stripGains = iSetup.getHandle(m_gainToken);
0099
0100   std::vector<double> c_localdirx;
0101   std::vector<double> c_localdiry;
0102   std::vector<double> c_localdirz;
0103   std::vector<uint16_t> c_firststrip;
0104   std::vector<uint16_t> c_nstrips;
0105   std::vector<bool> c_saturation;
0106   std::vector<bool> c_overlapping;
0107   std::vector<bool> c_farfromedge;
0108   std::vector<int> c_charge;
0109   std::vector<double> c_path;
0110   // NOTE only very few types are supported by NanoAOD, but more could be added (to discuss with XPOG / core software)
0111   // NOTE removed amplitude vector, I don't think it was used anywhere
0112   std::vector<float> c_gainused;      // NOTE was double
0113   std::vector<float> c_gainusedTick;  // NOTE was double
0114   for (const auto clus : clusters) {
0115     const auto& ampls = clus.cluster->amplitudes();
0116     const int firstStrip = clus.cluster->firstStrip();
0117     const int nStrips = ampls.size();
0118     double prevGain = -1;
0119     double prevGainTick = -1;
0120     if (stripGains.isValid()) {
0121       prevGain = stripGains->getApvGain(firstStrip / 128, stripGains->getRange(clus.det, 1), 1);
0122       prevGainTick = stripGains->getApvGain(firstStrip / 128, stripGains->getRange(clus.det, 0), 1);
0123     }
0124     const unsigned int charge = clus.cluster->charge();
0125     const bool saturation = std::any_of(ampls.begin(), ampls.end(), [](uint8_t amp) { return amp >= 254; });
0126
0127     const bool overlapping = (((firstStrip % 128) == 0) || ((firstStrip / 128) != ((firstStrip + nStrips) / 128)) ||
0128                               (((firstStrip + nStrips) % 128) == 127));
0129     const auto& trajState = clus.measurement.updatedState();
0130     const auto trackDir = trajState.localDirection();
0131     const auto cosine = trackDir.z() / trackDir.mag();
0132     const auto path = (10. * thickness(clus.det, tGeom.product())) / std::abs(cosine);
0133     const auto farFromEdge = isFarFromBorder(trajState, clus.det, tGeom.product());
0134     c_localdirx.push_back(trackDir.x());
0135     c_localdiry.push_back(trackDir.y());
0136     c_localdirz.push_back(trackDir.z());
0137     c_firststrip.push_back(firstStrip);
0138     c_nstrips.push_back(nStrips);
0139     c_saturation.push_back(saturation);
0140     c_overlapping.push_back(overlapping);
0141     c_farfromedge.push_back(farFromEdge);
0142     c_charge.push_back(charge);
0143     c_path.push_back(path);
0144     c_gainused.push_back(prevGain);
0145     c_gainusedTick.push_back(prevGainTick);
0146   }
0147   // addColumn(table, "localdirx", c_localdirx, "<doc>");
0148   // addColumn(table, "localdiry", c_localdiry, "<doc>");
0149   // addColumn(table, "localdirz", c_localdirz, "<doc>");
0150   // addColumn(table, "firststrip", c_firststrip, "<doc>");
0151   // addColumn(table, "nstrips", c_nstrips, "<doc>");
0152   addColumn(table, "saturation", c_saturation, "<doc>");
0153   addColumn(table, "overlapping", c_overlapping, "<doc>");
0154   addColumn(table, "farfromedge", c_farfromedge, "<doc>");
0155   addColumn(table, "charge", c_charge, "<doc>");
0156   // addColumn(table, "path", c_path, "<doc>");
0157   // ExtendedCalibTree: also charge/path
0158   addColumn(table, "gainused", c_gainused, "<doc>");
0159   addColumn(table, "gainusedTick", c_gainusedTick, "<doc>");
0160 }
0161
0162 #include "FWCore/PluginManager/interface/ModuleDef.h"
0163 #include "FWCore/Framework/interface/MakerMacros.h"
0164 DEFINE_FWK_MODULE(SiStripGainCalibTableProducer);