File indexing completed on 2024-04-06 11:59:47
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 }
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
0111
0112 std::vector<float> c_gainused;
0113 std::vector<float> c_gainusedTick;
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
0148
0149
0150
0151
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
0157
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);