File indexing completed on 2024-09-07 04:35:09
0001 #include "ShallowTrackClustersProducer.h"
0002
0003 #include "CalibTracker/SiStripCommon/interface/ShallowTools.h"
0004
0005 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0006 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0007 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0008 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0009 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0010 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0011 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015
0016 #include <map>
0017
0018 ShallowTrackClustersProducer::ShallowTrackClustersProducer(const edm::ParameterSet& iConfig)
0019 : tracks_token_(consumes<edm::View<reco::Track>>(iConfig.getParameter<edm::InputTag>("Tracks"))),
0020 association_token_(consumes<TrajTrackAssociationCollection>(iConfig.getParameter<edm::InputTag>("Tracks"))),
0021 clusters_token_(consumes<edmNew::DetSetVector<SiStripCluster>>(iConfig.getParameter<edm::InputTag>("Clusters"))),
0022 geomToken_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()),
0023 magFieldToken_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0024 laToken_(esConsumes<SiStripLorentzAngle, SiStripLorentzAngleDepRcd>()),
0025 Suffix(iConfig.getParameter<std::string>("Suffix")),
0026 Prefix(iConfig.getParameter<std::string>("Prefix")) {
0027 produces<std::vector<int>>(Prefix + "clusterIdx" + Suffix);
0028 produces<std::vector<int>>(Prefix + "onTrkClusterIdx" + Suffix);
0029 produces<std::vector<int>>(Prefix + "onTrkClustersBegin" + Suffix);
0030 produces<std::vector<int>>(Prefix + "onTrkClustersEnd" + Suffix);
0031 produces<std::vector<int>>(Prefix + "trackindex" + Suffix);
0032
0033 produces<std::vector<unsigned int>>(Prefix + "trackmulti" + Suffix);
0034 produces<std::vector<float>>(Prefix + "localtheta" + Suffix);
0035 produces<std::vector<float>>(Prefix + "localphi" + Suffix);
0036 produces<std::vector<float>>(Prefix + "localpitch" + Suffix);
0037 produces<std::vector<float>>(Prefix + "localx" + Suffix);
0038 produces<std::vector<float>>(Prefix + "localy" + Suffix);
0039 produces<std::vector<float>>(Prefix + "localz" + Suffix);
0040 produces<std::vector<float>>(Prefix + "strip" + Suffix);
0041 produces<std::vector<float>>(Prefix + "globaltheta" + Suffix);
0042 produces<std::vector<float>>(Prefix + "globalphi" + Suffix);
0043 produces<std::vector<float>>(Prefix + "globalx" + Suffix);
0044 produces<std::vector<float>>(Prefix + "globaly" + Suffix);
0045 produces<std::vector<float>>(Prefix + "globalz" + Suffix);
0046 produces<std::vector<float>>(Prefix + "insidistance" + Suffix);
0047 produces<std::vector<float>>(Prefix + "covered" + Suffix);
0048 produces<std::vector<float>>(Prefix + "projwidth" + Suffix);
0049 produces<std::vector<float>>(Prefix + "BdotY" + Suffix);
0050
0051 produces<std::vector<float>>(Prefix + "rhlocalx" + Suffix);
0052 produces<std::vector<float>>(Prefix + "rhlocaly" + Suffix);
0053 produces<std::vector<float>>(Prefix + "rhlocalxerr" + Suffix);
0054 produces<std::vector<float>>(Prefix + "rhlocalyerr" + Suffix);
0055 produces<std::vector<float>>(Prefix + "rhglobalx" + Suffix);
0056 produces<std::vector<float>>(Prefix + "rhglobaly" + Suffix);
0057 produces<std::vector<float>>(Prefix + "rhglobalz" + Suffix);
0058 produces<std::vector<float>>(Prefix + "rhstrip" + Suffix);
0059 produces<std::vector<float>>(Prefix + "rhmerr" + Suffix);
0060
0061 produces<std::vector<float>>(Prefix + "ubstrip" + Suffix);
0062 produces<std::vector<float>>(Prefix + "ubmerr" + Suffix);
0063
0064 produces<std::vector<float>>(Prefix + "driftx" + Suffix);
0065 produces<std::vector<float>>(Prefix + "drifty" + Suffix);
0066 produces<std::vector<float>>(Prefix + "driftz" + Suffix);
0067 produces<std::vector<float>>(Prefix + "globalZofunitlocalY" + Suffix);
0068 }
0069
0070 void ShallowTrackClustersProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0071 shallow::CLUSTERMAP clustermap = shallow::make_cluster_map(iEvent, clusters_token_);
0072 edm::Handle<edm::View<reco::Track>> tracks;
0073 iEvent.getByToken(tracks_token_, tracks);
0074
0075 int size = clustermap.size();
0076
0077
0078 auto clusterIdx = std::make_unique<std::vector<int>>();
0079 auto onTrkClusterIdx =
0080 std::make_unique<std::vector<int>>(size, -1);
0081 auto onTrkClustersBegin = std::make_unique<std::vector<int>>(tracks->size(), -1);
0082 auto onTrkClustersEnd = std::make_unique<std::vector<int>>(tracks->size(), -1);
0083 auto trackindex = std::make_unique<std::vector<int>>();
0084 clusterIdx->reserve(size);
0085 trackindex->reserve(size);
0086
0087 auto trackmulti = std::make_unique<std::vector<unsigned int>>();
0088 trackmulti->reserve(size);
0089 auto localtheta = std::make_unique<std::vector<float>>();
0090 localtheta->reserve(size);
0091 auto localphi = std::make_unique<std::vector<float>>();
0092 localphi->reserve(size);
0093 auto localpitch = std::make_unique<std::vector<float>>();
0094 localpitch->reserve(size);
0095 auto localx = std::make_unique<std::vector<float>>();
0096 localx->reserve(size);
0097 auto localy = std::make_unique<std::vector<float>>();
0098 localy->reserve(size);
0099 auto localz = std::make_unique<std::vector<float>>();
0100 localz->reserve(size);
0101 auto strip = std::make_unique<std::vector<float>>();
0102 strip->reserve(size);
0103 auto globaltheta = std::make_unique<std::vector<float>>();
0104 globaltheta->reserve(size);
0105 auto globalphi = std::make_unique<std::vector<float>>();
0106 globalphi->reserve(size);
0107 auto globalx = std::make_unique<std::vector<float>>();
0108 globalx->reserve(size);
0109 auto globaly = std::make_unique<std::vector<float>>();
0110 globaly->reserve(size);
0111 auto globalz = std::make_unique<std::vector<float>>();
0112 globalz->reserve(size);
0113 auto insidistance = std::make_unique<std::vector<float>>();
0114 insidistance->reserve(size);
0115 auto projwidth = std::make_unique<std::vector<float>>();
0116 projwidth->reserve(size);
0117 auto BdotY = std::make_unique<std::vector<float>>();
0118 BdotY->reserve(size);
0119 auto covered = std::make_unique<std::vector<float>>();
0120 covered->reserve(size);
0121 auto rhlocalx = std::make_unique<std::vector<float>>();
0122 rhlocalx->reserve(size);
0123 auto rhlocaly = std::make_unique<std::vector<float>>();
0124 rhlocaly->reserve(size);
0125 auto rhlocalxerr = std::make_unique<std::vector<float>>();
0126 rhlocalxerr->reserve(size);
0127 auto rhlocalyerr = std::make_unique<std::vector<float>>();
0128 rhlocalyerr->reserve(size);
0129 auto rhglobalx = std::make_unique<std::vector<float>>();
0130 rhglobalx->reserve(size);
0131 auto rhglobaly = std::make_unique<std::vector<float>>();
0132 rhglobaly->reserve(size);
0133 auto rhglobalz = std::make_unique<std::vector<float>>();
0134 rhglobalz->reserve(size);
0135 auto rhstrip = std::make_unique<std::vector<float>>();
0136 rhstrip->reserve(size);
0137 auto rhmerr = std::make_unique<std::vector<float>>();
0138 rhmerr->reserve(size);
0139 auto ubstrip = std::make_unique<std::vector<float>>();
0140 ubstrip->reserve(size);
0141 auto ubmerr = std::make_unique<std::vector<float>>();
0142 ubmerr->reserve(size);
0143 auto driftx = std::make_unique<std::vector<float>>();
0144 driftx->reserve(size);
0145 auto drifty = std::make_unique<std::vector<float>>();
0146 drifty->reserve(size);
0147 auto driftz = std::make_unique<std::vector<float>>();
0148 driftz->reserve(size);
0149 auto globalZofunitlocalY = std::make_unique<std::vector<float>>();
0150 globalZofunitlocalY->reserve(size);
0151
0152 edm::ESHandle<TrackerGeometry> theTrackerGeometry = iSetup.getHandle(geomToken_);
0153 edm::ESHandle<MagneticField> magfield = iSetup.getHandle(magFieldToken_);
0154 edm::ESHandle<SiStripLorentzAngle> SiStripLorentzAngle = iSetup.getHandle(laToken_);
0155
0156 edm::Handle<TrajTrackAssociationCollection> associations;
0157 iEvent.getByToken(association_token_, associations);
0158
0159 TrajectoryStateCombiner combiner;
0160
0161 size_t ontrk_cluster_idx = 0;
0162 std::map<size_t, std::vector<size_t>> mapping;
0163
0164 for (TrajTrackAssociationCollection::const_iterator association = associations->begin();
0165 association != associations->end();
0166 association++) {
0167 const Trajectory* traj = association->key.get();
0168 const reco::Track* track = association->val.get();
0169 int trk_idx = shallow::findTrackIndex(tracks, track);
0170 size_t trk_strt_idx = ontrk_cluster_idx;
0171
0172 for (auto const& measurement : traj->measurements()) {
0173 const TrajectoryStateOnSurface& tsos = measurement.updatedState();
0174 const TrajectoryStateOnSurface unbiased =
0175 combiner(measurement.forwardPredictedState(), measurement.backwardPredictedState());
0176
0177 const TrackingRecHit* hit = measurement.recHit()->hit();
0178 const SiStripRecHit1D* hit1D = dynamic_cast<const SiStripRecHit1D*>(hit);
0179 const SiStripRecHit2D* hit2D = dynamic_cast<const SiStripRecHit2D*>(hit);
0180 const SiStripMatchedRecHit2D* matchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
0181
0182 for (unsigned h = 0; h < 2; h++) {
0183 const SiStripCluster* cluster_ptr;
0184 if (!matchedhit && h == 1)
0185 continue;
0186 else if (matchedhit && h == 0)
0187 cluster_ptr = &matchedhit->monoCluster();
0188 else if (matchedhit && h == 1)
0189 cluster_ptr = &matchedhit->stereoCluster();
0190 else if (hit2D)
0191 cluster_ptr = (hit2D->cluster()).get();
0192 else if (hit1D)
0193 cluster_ptr = (hit1D->cluster()).get();
0194 else
0195 continue;
0196
0197 shallow::CLUSTERMAP::const_iterator cluster =
0198 clustermap.find(std::make_pair(hit->geographicalId().rawId(), cluster_ptr->firstStrip()));
0199 if (cluster == clustermap.end())
0200 throw cms::Exception("Logic Error") << "Cluster not found: this could be a configuration error" << std::endl;
0201
0202 unsigned i = cluster->second;
0203
0204
0205 auto already_visited = mapping.find(i);
0206 int nassociations = 1;
0207 if (already_visited != mapping.end()) {
0208 nassociations += already_visited->second.size();
0209 for (size_t idx : already_visited->second) {
0210 trackmulti->at(idx)++;
0211 }
0212 already_visited->second.push_back(ontrk_cluster_idx);
0213 } else {
0214 std::vector<size_t> single = {ontrk_cluster_idx};
0215 mapping.insert(std::make_pair(i, single));
0216 }
0217
0218 const StripGeomDetUnit* theStripDet =
0219 dynamic_cast<const StripGeomDetUnit*>(theTrackerGeometry->idToDet(hit->geographicalId()));
0220 LocalVector drift = shallow::drift(theStripDet, *magfield, *SiStripLorentzAngle);
0221
0222 if (nassociations == 1)
0223 onTrkClusterIdx->at(i) = ontrk_cluster_idx;
0224 clusterIdx->push_back(i);
0225 trackmulti->push_back(nassociations);
0226 trackindex->push_back(trk_idx);
0227 localtheta->push_back((theStripDet->toLocal(tsos.globalDirection())).theta());
0228 localphi->push_back((theStripDet->toLocal(tsos.globalDirection())).phi());
0229 localpitch->push_back(
0230 (theStripDet->specificTopology()).localPitch(theStripDet->toLocal(tsos.globalPosition())));
0231 localx->push_back((theStripDet->toLocal(tsos.globalPosition())).x());
0232 localy->push_back((theStripDet->toLocal(tsos.globalPosition())).y());
0233 localz->push_back((theStripDet->toLocal(tsos.globalPosition())).z());
0234 strip->push_back((theStripDet->specificTopology()).strip(theStripDet->toLocal(tsos.globalPosition())));
0235 globaltheta->push_back(tsos.globalDirection().theta());
0236 globalphi->push_back(tsos.globalDirection().phi());
0237 globalx->push_back(tsos.globalPosition().x());
0238 globaly->push_back(tsos.globalPosition().y());
0239 globalz->push_back(tsos.globalPosition().z());
0240 insidistance->push_back(1. / fabs(cos(localtheta->at(ontrk_cluster_idx))));
0241 projwidth->push_back(tan(localtheta->at(ontrk_cluster_idx)) * cos(localphi->at(ontrk_cluster_idx)));
0242 BdotY->push_back((theStripDet->surface()).toLocal(magfield->inTesla(theStripDet->surface().position())).y());
0243 covered->push_back(drift.z() / localpitch->at(ontrk_cluster_idx) *
0244 fabs(projwidth->at(ontrk_cluster_idx) - drift.x() / drift.z()));
0245 rhlocalx->push_back(hit->localPosition().x());
0246 rhlocaly->push_back(hit->localPosition().y());
0247 rhlocalxerr->push_back(sqrt(hit->localPositionError().xx()));
0248 rhlocalyerr->push_back(sqrt(hit->localPositionError().yy()));
0249 rhglobalx->push_back(theStripDet->toGlobal(hit->localPosition()).x());
0250 rhglobaly->push_back(theStripDet->toGlobal(hit->localPosition()).y());
0251 rhglobalz->push_back(theStripDet->toGlobal(hit->localPosition()).z());
0252 rhstrip->push_back(theStripDet->specificTopology().strip(hit->localPosition()));
0253 rhmerr->push_back(sqrt(
0254 theStripDet->specificTopology().measurementError(hit->localPosition(), hit->localPositionError()).uu()));
0255 ubstrip->push_back(theStripDet->specificTopology().strip(unbiased.localPosition()));
0256 ubmerr->push_back(sqrt(theStripDet->specificTopology()
0257 .measurementError(unbiased.localPosition(), unbiased.localError().positionError())
0258 .uu()));
0259 driftx->push_back(drift.x());
0260 drifty->push_back(drift.y());
0261 driftz->push_back(drift.z());
0262 globalZofunitlocalY->push_back((theStripDet->toGlobal(LocalVector(0, 1, 0))).z());
0263
0264 ontrk_cluster_idx++;
0265 }
0266 }
0267
0268 onTrkClustersBegin->at(trk_idx) = trk_strt_idx;
0269 onTrkClustersEnd->at(trk_idx) = ontrk_cluster_idx;
0270
0271 }
0272
0273 iEvent.put(std::move(clusterIdx), Prefix + "clusterIdx" + Suffix);
0274 iEvent.put(std::move(onTrkClusterIdx), Prefix + "onTrkClusterIdx" + Suffix);
0275 iEvent.put(std::move(onTrkClustersBegin), Prefix + "onTrkClustersBegin" + Suffix);
0276 iEvent.put(std::move(onTrkClustersEnd), Prefix + "onTrkClustersEnd" + Suffix);
0277 iEvent.put(std::move(trackindex), Prefix + "trackindex" + Suffix);
0278
0279 iEvent.put(std::move(trackmulti), Prefix + "trackmulti" + Suffix);
0280 iEvent.put(std::move(localtheta), Prefix + "localtheta" + Suffix);
0281 iEvent.put(std::move(localphi), Prefix + "localphi" + Suffix);
0282 iEvent.put(std::move(localpitch), Prefix + "localpitch" + Suffix);
0283 iEvent.put(std::move(localx), Prefix + "localx" + Suffix);
0284 iEvent.put(std::move(localy), Prefix + "localy" + Suffix);
0285 iEvent.put(std::move(localz), Prefix + "localz" + Suffix);
0286 iEvent.put(std::move(strip), Prefix + "strip" + Suffix);
0287 iEvent.put(std::move(globaltheta), Prefix + "globaltheta" + Suffix);
0288 iEvent.put(std::move(globalphi), Prefix + "globalphi" + Suffix);
0289 iEvent.put(std::move(globalx), Prefix + "globalx" + Suffix);
0290 iEvent.put(std::move(globaly), Prefix + "globaly" + Suffix);
0291 iEvent.put(std::move(globalz), Prefix + "globalz" + Suffix);
0292 iEvent.put(std::move(insidistance), Prefix + "insidistance" + Suffix);
0293 iEvent.put(std::move(covered), Prefix + "covered" + Suffix);
0294 iEvent.put(std::move(projwidth), Prefix + "projwidth" + Suffix);
0295 iEvent.put(std::move(BdotY), Prefix + "BdotY" + Suffix);
0296 iEvent.put(std::move(rhlocalx), Prefix + "rhlocalx" + Suffix);
0297 iEvent.put(std::move(rhlocaly), Prefix + "rhlocaly" + Suffix);
0298 iEvent.put(std::move(rhlocalxerr), Prefix + "rhlocalxerr" + Suffix);
0299 iEvent.put(std::move(rhlocalyerr), Prefix + "rhlocalyerr" + Suffix);
0300 iEvent.put(std::move(rhglobalx), Prefix + "rhglobalx" + Suffix);
0301 iEvent.put(std::move(rhglobaly), Prefix + "rhglobaly" + Suffix);
0302 iEvent.put(std::move(rhglobalz), Prefix + "rhglobalz" + Suffix);
0303 iEvent.put(std::move(rhstrip), Prefix + "rhstrip" + Suffix);
0304 iEvent.put(std::move(rhmerr), Prefix + "rhmerr" + Suffix);
0305 iEvent.put(std::move(ubstrip), Prefix + "ubstrip" + Suffix);
0306 iEvent.put(std::move(ubmerr), Prefix + "ubmerr" + Suffix);
0307 iEvent.put(std::move(driftx), Prefix + "driftx" + Suffix);
0308 iEvent.put(std::move(drifty), Prefix + "drifty" + Suffix);
0309 iEvent.put(std::move(driftz), Prefix + "driftz" + Suffix);
0310 iEvent.put(std::move(globalZofunitlocalY), Prefix + "globalZofunitlocalY" + Suffix);
0311 }