SiStripCommon/plugins/ShallowTrackClustersProducer.cc

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);          //link: on trk cluster --> general cluster info
0028   produces<std::vector<int>>(Prefix + "onTrkClusterIdx" + Suffix);     //link: general cluster info --> on track cluster
0029   produces<std::vector<int>>(Prefix + "onTrkClustersBegin" + Suffix);  //link: track --> onTrkInfo (range)
0030   produces<std::vector<int>>(Prefix + "onTrkClustersEnd" + Suffix);    //link: track --> onTrkInfo (range)
0031   produces<std::vector<int>>(Prefix + "trackindex" + Suffix);          //link: on trk cluster --> track index
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   //links
0078   auto clusterIdx = std::make_unique<std::vector<int>>();  //link: on trk cluster --> general cluster info
0079   auto onTrkClusterIdx =
0080       std::make_unique<std::vector<int>>(size, -1);  //link: general cluster info --> on track cluster
0081   auto onTrkClustersBegin = std::make_unique<std::vector<int>>(tracks->size(), -1);  //link: track --> on trk cluster
0082   auto onTrkClustersEnd = std::make_unique<std::vector<int>>(tracks->size(), -1);    //link: track --> on trk cluster
0083   auto trackindex = std::make_unique<std::vector<int>>();                            //link: on track cluster --> track
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;  //cluster idx --> on trk cluster idx (multiple)
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++) {  //loop over possible Hit options (1D, 2D)
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         //find if cluster was already assigned to a previous track
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 {  //otherwise store this
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;  //link: general cluster info --> on track cluster
0224         clusterIdx->push_back(i);                      //link: on trk cluster --> general cluster info
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       }  //for(unsigned h=0; h<2; h++) { //loop over possible Hit options (1D, 2D)
0266     }  //for(auto const& measurement : traj->measurements() )
0267
0268     onTrkClustersBegin->at(trk_idx) = trk_strt_idx;
0269     onTrkClustersEnd->at(trk_idx) = ontrk_cluster_idx;
0270
0271   }  //for(TrajTrackAssociationCollection::const_iterator association = associations->begin();
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 }