File indexing completed on 2024-04-06 11:59:46
0001 #include "ShallowGainCalibration.h"
0002
0003 using namespace edm;
0004 using namespace reco;
0005 using namespace std;
0006
0007 ShallowGainCalibration::ShallowGainCalibration(const edm::ParameterSet& iConfig)
0008 : tracks_token_(consumes<edm::View<reco::Track>>(iConfig.getParameter<edm::InputTag>("Tracks"))),
0009 association_token_(consumes<TrajTrackAssociationCollection>(iConfig.getParameter<edm::InputTag>("Tracks"))),
0010 trackerGeometry_token_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()),
0011 gain_token_(esConsumes<SiStripGain, SiStripGainRcd>()),
0012 tkGeom_token_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()),
0013 Suffix(iConfig.getParameter<std::string>("Suffix")),
0014 Prefix(iConfig.getParameter<std::string>("Prefix")) {
0015 produces<std::vector<int>>(Prefix + "trackindex" + Suffix);
0016 produces<std::vector<unsigned int>>(Prefix + "rawid" + Suffix);
0017 produces<std::vector<double>>(Prefix + "localdirx" + Suffix);
0018 produces<std::vector<double>>(Prefix + "localdiry" + Suffix);
0019 produces<std::vector<double>>(Prefix + "localdirz" + Suffix);
0020 produces<std::vector<unsigned short>>(Prefix + "firststrip" + Suffix);
0021 produces<std::vector<unsigned short>>(Prefix + "nstrips" + Suffix);
0022 produces<std::vector<bool>>(Prefix + "saturation" + Suffix);
0023 produces<std::vector<bool>>(Prefix + "overlapping" + Suffix);
0024 produces<std::vector<bool>>(Prefix + "farfromedge" + Suffix);
0025 produces<std::vector<unsigned int>>(Prefix + "charge" + Suffix);
0026 produces<std::vector<double>>(Prefix + "path" + Suffix);
0027 #ifdef ExtendedCALIBTree
0028 produces<std::vector<double>>(Prefix + "chargeoverpath" + Suffix);
0029 #endif
0030 produces<std::vector<unsigned char>>(Prefix + "amplitude" + Suffix);
0031 produces<std::vector<double>>(Prefix + "gainused" + Suffix);
0032 produces<std::vector<double>>(Prefix + "gainusedTick" + Suffix);
0033 }
0034
0035 void ShallowGainCalibration::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0036 auto trackindex = std::make_unique<std::vector<int>>();
0037 auto rawid = std::make_unique<std::vector<unsigned int>>();
0038 auto localdirx = std::make_unique<std::vector<double>>();
0039 auto localdiry = std::make_unique<std::vector<double>>();
0040 auto localdirz = std::make_unique<std::vector<double>>();
0041 auto firststrip = std::make_unique<std::vector<unsigned short>>();
0042 auto nstrips = std::make_unique<std::vector<unsigned short>>();
0043 auto saturation = std::make_unique<std::vector<bool>>();
0044 auto overlapping = std::make_unique<std::vector<bool>>();
0045 auto farfromedge = std::make_unique<std::vector<bool>>();
0046 auto charge = std::make_unique<std::vector<unsigned int>>();
0047 auto path = std::make_unique<std::vector<double>>();
0048 #ifdef ExtendedCALIBTree
0049 auto chargeoverpath = std::make_unique<std::vector<double>>();
0050 #endif
0051 auto amplitude = std::make_unique<std::vector<unsigned char>>();
0052 auto gainused = std::make_unique<std::vector<double>>();
0053 auto gainusedTick = std::make_unique<std::vector<double>>();
0054
0055 m_tracker = &iSetup.getData(trackerGeometry_token_);
0056 edm::ESHandle<SiStripGain> gainHandle = iSetup.getHandle(gain_token_);
0057 edm::Handle<edm::View<reco::Track>> tracks;
0058 iEvent.getByToken(tracks_token_, tracks);
0059 edm::Handle<TrajTrackAssociationCollection> associations;
0060 iEvent.getByToken(association_token_, associations);
0061
0062 for (TrajTrackAssociationCollection::const_iterator association = associations->begin();
0063 association != associations->end();
0064 association++) {
0065 const Trajectory* traj = association->key.get();
0066 const reco::Track* track = association->val.get();
0067
0068 vector<TrajectoryMeasurement> measurements = traj->measurements();
0069 for (vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin();
0070 measurement_it != measurements.end();
0071 measurement_it++) {
0072 TrajectoryStateOnSurface trajState = measurement_it->updatedState();
0073 if (!trajState.isValid())
0074 continue;
0075
0076 const TrackingRecHit* hit = (*measurement_it->recHit()).hit();
0077 const SiStripRecHit1D* sistripsimple1dhit = dynamic_cast<const SiStripRecHit1D*>(hit);
0078 const SiStripRecHit2D* sistripsimplehit = dynamic_cast<const SiStripRecHit2D*>(hit);
0079 const SiStripMatchedRecHit2D* sistripmatchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
0080 const SiPixelRecHit* sipixelhit = dynamic_cast<const SiPixelRecHit*>(hit);
0081
0082 const SiPixelCluster* PixelCluster = nullptr;
0083 const SiStripCluster* StripCluster = nullptr;
0084 uint32_t DetId = 0;
0085
0086 for (unsigned int h = 0; h < 2; h++) {
0087 if (!sistripmatchedhit && h == 1) {
0088 continue;
0089 } else if (sistripmatchedhit && h == 0) {
0090 StripCluster = &sistripmatchedhit->monoCluster();
0091 DetId = sistripmatchedhit->monoId();
0092 } else if (sistripmatchedhit && h == 1) {
0093 StripCluster = &sistripmatchedhit->stereoCluster();
0094 ;
0095 DetId = sistripmatchedhit->stereoId();
0096 } else if (sistripsimplehit) {
0097 StripCluster = (sistripsimplehit->cluster()).get();
0098 DetId = sistripsimplehit->geographicalId().rawId();
0099 } else if (sistripsimple1dhit) {
0100 StripCluster = (sistripsimple1dhit->cluster()).get();
0101 DetId = sistripsimple1dhit->geographicalId().rawId();
0102 } else if (sipixelhit) {
0103 PixelCluster = (sipixelhit->cluster()).get();
0104 DetId = sipixelhit->geographicalId().rawId();
0105 } else {
0106 continue;
0107 }
0108
0109 LocalVector trackDirection = trajState.localDirection();
0110 double cosine = trackDirection.z() / trackDirection.mag();
0111 bool Saturation = false;
0112 bool Overlapping = false;
0113 unsigned int Charge = 0;
0114 double Path = (10.0 * thickness(DetId)) / fabs(cosine);
0115 double PrevGain = -1;
0116 double PrevGainTick = -1;
0117 int FirstStrip = 0;
0118 int NStrips = 0;
0119
0120 if (StripCluster) {
0121 const auto& Ampls = StripCluster->amplitudes();
0122 FirstStrip = StripCluster->firstStrip();
0123 NStrips = Ampls.size();
0124 int APVId = FirstStrip / 128;
0125
0126 if (gainHandle.isValid()) {
0127 PrevGain = gainHandle->getApvGain(APVId, gainHandle->getRange(DetId, 1), 1);
0128 PrevGainTick = gainHandle->getApvGain(APVId, gainHandle->getRange(DetId, 0), 1);
0129 }
0130
0131 for (unsigned int a = 0; a < Ampls.size(); a++) {
0132 Charge += Ampls[a];
0133 if (Ampls[a] >= 254)
0134 Saturation = true;
0135 amplitude->push_back(Ampls[a]);
0136 }
0137
0138 if (FirstStrip == 0)
0139 Overlapping = true;
0140 if (FirstStrip == 128)
0141 Overlapping = true;
0142 if (FirstStrip == 256)
0143 Overlapping = true;
0144 if (FirstStrip == 384)
0145 Overlapping = true;
0146 if (FirstStrip == 512)
0147 Overlapping = true;
0148 if (FirstStrip == 640)
0149 Overlapping = true;
0150
0151 if (FirstStrip <= 127 && FirstStrip + Ampls.size() > 127)
0152 Overlapping = true;
0153 if (FirstStrip <= 255 && FirstStrip + Ampls.size() > 255)
0154 Overlapping = true;
0155 if (FirstStrip <= 383 && FirstStrip + Ampls.size() > 383)
0156 Overlapping = true;
0157 if (FirstStrip <= 511 && FirstStrip + Ampls.size() > 511)
0158 Overlapping = true;
0159 if (FirstStrip <= 639 && FirstStrip + Ampls.size() > 639)
0160 Overlapping = true;
0161
0162 if (FirstStrip + Ampls.size() == 127)
0163 Overlapping = true;
0164 if (FirstStrip + Ampls.size() == 255)
0165 Overlapping = true;
0166 if (FirstStrip + Ampls.size() == 383)
0167 Overlapping = true;
0168 if (FirstStrip + Ampls.size() == 511)
0169 Overlapping = true;
0170 if (FirstStrip + Ampls.size() == 639)
0171 Overlapping = true;
0172 if (FirstStrip + Ampls.size() == 767)
0173 Overlapping = true;
0174 } else if (PixelCluster) {
0175 const auto& Ampls = PixelCluster->pixelADC();
0176 int FirstRow = PixelCluster->minPixelRow();
0177 int FirstCol = PixelCluster->minPixelCol();
0178 FirstStrip = ((FirstRow / 80) << 3 | (FirstCol / 52)) * 128;
0179 NStrips = 0;
0180 Saturation = false;
0181 Overlapping = false;
0182
0183 for (unsigned int a = 0; a < Ampls.size(); a++) {
0184 Charge += Ampls[a];
0185 if (Ampls[a] >= 254)
0186 Saturation = true;
0187 }
0188 }
0189 #ifdef ExtendedCALIBTree
0190 double ChargeOverPath = (double)Charge / Path;
0191 #endif
0192
0193 trackindex->push_back(shallow::findTrackIndex(tracks, track));
0194 rawid->push_back(DetId);
0195 localdirx->push_back(trackDirection.x());
0196 localdiry->push_back(trackDirection.y());
0197 localdirz->push_back(trackDirection.z());
0198 firststrip->push_back(FirstStrip);
0199 nstrips->push_back(NStrips);
0200 saturation->push_back(Saturation);
0201 overlapping->push_back(Overlapping);
0202 farfromedge->push_back(StripCluster ? isFarFromBorder(&trajState, DetId, &iSetup) : true);
0203 charge->push_back(Charge);
0204 path->push_back(Path);
0205 #ifdef ExtendedCALIBTree
0206 chargeoverpath->push_back(ChargeOverPath);
0207 #endif
0208 gainused->push_back(PrevGain);
0209 gainusedTick->push_back(PrevGainTick);
0210 }
0211 }
0212 }
0213
0214 iEvent.put(std::move(trackindex), Prefix + "trackindex" + Suffix);
0215 iEvent.put(std::move(rawid), Prefix + "rawid" + Suffix);
0216 iEvent.put(std::move(localdirx), Prefix + "localdirx" + Suffix);
0217 iEvent.put(std::move(localdiry), Prefix + "localdiry" + Suffix);
0218 iEvent.put(std::move(localdirz), Prefix + "localdirz" + Suffix);
0219 iEvent.put(std::move(firststrip), Prefix + "firststrip" + Suffix);
0220 iEvent.put(std::move(nstrips), Prefix + "nstrips" + Suffix);
0221 iEvent.put(std::move(saturation), Prefix + "saturation" + Suffix);
0222 iEvent.put(std::move(overlapping), Prefix + "overlapping" + Suffix);
0223 iEvent.put(std::move(farfromedge), Prefix + "farfromedge" + Suffix);
0224 iEvent.put(std::move(charge), Prefix + "charge" + Suffix);
0225 iEvent.put(std::move(path), Prefix + "path" + Suffix);
0226 #ifdef ExtendedCALIBTree
0227 iEvent.put(std::move(chargeoverpath), Prefix + "chargeoverpath" + Suffix);
0228 #endif
0229 iEvent.put(std::move(amplitude), Prefix + "amplitude" + Suffix);
0230 iEvent.put(std::move(gainused), Prefix + "gainused" + Suffix);
0231 iEvent.put(std::move(gainusedTick), Prefix + "gainusedTick" + Suffix);
0232 }
0233
0234 bool ShallowGainCalibration::isFarFromBorder(TrajectoryStateOnSurface* trajState,
0235 const uint32_t detid,
0236 const edm::EventSetup* iSetup) {
0237 edm::ESHandle<TrackerGeometry> tkGeom = iSetup->getHandle(tkGeom_token_);
0238
0239 LocalPoint HitLocalPos = trajState->localPosition();
0240 LocalError HitLocalError = trajState->localError().positionError();
0241
0242 const GeomDetUnit* it = tkGeom->idToDetUnit(DetId(detid));
0243 if (dynamic_cast<const StripGeomDetUnit*>(it) == nullptr && dynamic_cast<const PixelGeomDetUnit*>(it) == nullptr) {
0244 throw cms::Exception("Logic Error") << "\t\t this detID doesn't seem to belong to the Tracker";
0245 }
0246
0247 const BoundPlane plane = it->surface();
0248 const TrapezoidalPlaneBounds* trapezoidalBounds(dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
0249 const RectangularPlaneBounds* rectangularBounds(dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
0250
0251 double DistFromBorder = 1.0;
0252 double HalfLength = it->surface().bounds().length() / 2.0;
0253
0254 if (trapezoidalBounds) {
0255 std::array<const float, 4> const& parameters = (*trapezoidalBounds).parameters();
0256 HalfLength = parameters[3];
0257 } else if (rectangularBounds) {
0258 HalfLength = it->surface().bounds().length() / 2.0;
0259 } else {
0260 return false;
0261 }
0262
0263 if (fabs(HitLocalPos.y()) + HitLocalError.yy() >= (HalfLength - DistFromBorder))
0264 return false;
0265
0266 return true;
0267 }
0268
0269 double ShallowGainCalibration::thickness(DetId id) {
0270 map<DetId, double>::iterator th = m_thicknessMap.find(id);
0271 if (th != m_thicknessMap.end())
0272 return (*th).second;
0273 else {
0274 double detThickness = 1.;
0275
0276 const GeomDetUnit* it = m_tracker->idToDetUnit(DetId(id));
0277 bool isPixel = dynamic_cast<const PixelGeomDetUnit*>(it) != nullptr;
0278 bool isStrip = dynamic_cast<const StripGeomDetUnit*>(it) != nullptr;
0279 if (!isPixel && !isStrip) {
0280 throw cms::Exception("Logic Error") << "\t\t this detID doesn't seem to belong to the Tracker";
0281 } else {
0282 detThickness = it->surface().bounds().thickness();
0283 }
0284
0285 m_thicknessMap[id] = detThickness;
0286 return detThickness;
0287 }
0288 }
