SiStripCommon/plugins/ShallowSimhitClustersProducer.cc

0001 #include "ShallowSimhitClustersProducer.h"
0002
0003 #include "CalibTracker/SiStripCommon/interface/ShallowTools.h"
0004 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0005 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0006 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0010
0011 ShallowSimhitClustersProducer::ShallowSimhitClustersProducer(const edm::ParameterSet& iConfig)
0012     : clusters_token_(consumes<edmNew::DetSetVector<SiStripCluster>>(iConfig.getParameter<edm::InputTag>("Clusters"))),
0013       geomToken_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()),
0014       magFieldToken_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0015       laToken_(esConsumes<SiStripLorentzAngle, SiStripLorentzAngleRcd>(
0016           edm::ESInputTag{"", iConfig.getParameter<std::string>("runningMode")})),
0017       Prefix(iConfig.getParameter<std::string>("Prefix")),
0018       runningmode_(iConfig.getParameter<std::string>("runningMode")) {
0019   std::vector<edm::InputTag> simhits_tags = iConfig.getParameter<std::vector<edm::InputTag>>("InputTags");
0020   for (const auto& itag : simhits_tags) {
0021     simhits_tokens_.push_back(consumes<std::vector<PSimHit>>(itag));
0022   }
0023
0024   produces<std::vector<unsigned>>(Prefix + "hits");
0025   produces<std::vector<float>>(Prefix + "strip");
0026   produces<std::vector<float>>(Prefix + "localtheta");
0027   produces<std::vector<float>>(Prefix + "localphi");
0028   produces<std::vector<float>>(Prefix + "localx");
0029   produces<std::vector<float>>(Prefix + "localy");
0030   produces<std::vector<float>>(Prefix + "localz");
0031   produces<std::vector<float>>(Prefix + "momentum");
0032   produces<std::vector<float>>(Prefix + "energyloss");
0033   produces<std::vector<float>>(Prefix + "time");
0034   produces<std::vector<int>>(Prefix + "particle");
0035   produces<std::vector<unsigned short>>(Prefix + "process");
0036 }
0037
0038 void ShallowSimhitClustersProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0039   shallow::CLUSTERMAP clustermap = shallow::make_cluster_map(iEvent, clusters_token_);
0040
0041   int size = clustermap.size();
0042   auto hits = std::make_unique<std::vector<unsigned>>(size, 0);
0043   auto strip = std::make_unique<std::vector<float>>(size, -100);
0044   auto localtheta = std::make_unique<std::vector<float>>(size, -100);
0045   auto localphi = std::make_unique<std::vector<float>>(size, -100);
0046   auto localx = std::make_unique<std::vector<float>>(size, -100);
0047   auto localy = std::make_unique<std::vector<float>>(size, -100);
0048   auto localz = std::make_unique<std::vector<float>>(size, -100);
0049   auto momentum = std::make_unique<std::vector<float>>(size, 0);
0050   auto energyloss = std::make_unique<std::vector<float>>(size, -1);
0051   auto time = std::make_unique<std::vector<float>>(size, -1);
0052   auto particle = std::make_unique<std::vector<int>>(size, -500);
0053   auto process = std::make_unique<std::vector<unsigned short>>(size, 0);
0054
0055   edm::ESHandle<TrackerGeometry> theTrackerGeometry = iSetup.getHandle(geomToken_);
0056   edm::ESHandle<MagneticField> magfield = iSetup.getHandle(magFieldToken_);
0057   edm::ESHandle<SiStripLorentzAngle> SiStripLorentzAngle = iSetup.getHandle(laToken_);
0058   edm::Handle<edmNew::DetSetVector<SiStripCluster>> clusters;
0059   iEvent.getByLabel("siStripClusters", "", clusters);
0060
0061   for (auto& simhit_token : simhits_tokens_) {
0062     edm::Handle<std::vector<PSimHit>> simhits;
0063     iEvent.getByToken(simhit_token, simhits);
0064     for (auto const& hit : *simhits) {
0065       const uint32_t id = hit.detUnitId();
0066       const StripGeomDetUnit* theStripDet = dynamic_cast<const StripGeomDetUnit*>(theTrackerGeometry->idToDet(id));
0067       const LocalVector drift = shallow::drift(theStripDet, *magfield, *SiStripLorentzAngle);
0068
0069       const float driftedstrip_ = theStripDet->specificTopology().strip(hit.localPosition() + 0.5 * drift);
0070       const float hitstrip_ = theStripDet->specificTopology().strip(hit.localPosition());
0071
0072       shallow::CLUSTERMAP::const_iterator cluster = match_cluster(id, driftedstrip_, clustermap, *clusters);
0073       if (cluster != clustermap.end()) {
0074         unsigned i = cluster->second;
0075         hits->at(i) += 1;
0076         if (hits->at(i) == 1) {
0077           strip->at(i) = hitstrip_;
0078           localtheta->at(i) = hit.thetaAtEntry();
0079           localphi->at(i) = hit.phiAtEntry();
0080           localx->at(i) = hit.localPosition().x();
0081           localy->at(i) = hit.localPosition().y();
0082           localz->at(i) = hit.localPosition().z();
0083           momentum->at(i) = hit.pabs();
0084           energyloss->at(i) = hit.energyLoss();
0085           time->at(i) = hit.timeOfFlight();
0086           particle->at(i) = hit.particleType();
0087           process->at(i) = hit.processType();
0088         }
0089       }
0090     }
0091   }
0092
0093   iEvent.put(std::move(hits), Prefix + "hits");
0094   iEvent.put(std::move(strip), Prefix + "strip");
0095   iEvent.put(std::move(localtheta), Prefix + "localtheta");
0096   iEvent.put(std::move(localphi), Prefix + "localphi");
0097   iEvent.put(std::move(localx), Prefix + "localx");
0098   iEvent.put(std::move(localy), Prefix + "localy");
0099   iEvent.put(std::move(localz), Prefix + "localz");
0100   iEvent.put(std::move(momentum), Prefix + "momentum");
0101   iEvent.put(std::move(energyloss), Prefix + "energyloss");
0102   iEvent.put(std::move(time), Prefix + "time");
0103   iEvent.put(std::move(particle), Prefix + "particle");
0104   iEvent.put(std::move(process), Prefix + "process");
0105 }
0106
0107 shallow::CLUSTERMAP::const_iterator ShallowSimhitClustersProducer::match_cluster(
0108     const unsigned& id,
0109     const float& strip_,
0110     const shallow::CLUSTERMAP& clustermap,
0111     const edmNew::DetSetVector<SiStripCluster>& clusters) const {
0112   shallow::CLUSTERMAP::const_iterator cluster = clustermap.end();
0113   edmNew::DetSetVector<SiStripCluster>::const_iterator clustersDetSet = clusters.find(id);
0114   if (clustersDetSet != clusters.end()) {
0115     edmNew::DetSet<SiStripCluster>::const_iterator left, right = clustersDetSet->begin();
0116     while (right != clustersDetSet->end() && strip_ > right->barycenter())
0117       right++;
0118     left = right - 1;
0119     if (right != clustersDetSet->end() && right != clustersDetSet->begin()) {
0120       unsigned firstStrip =
0121           (right->barycenter() - strip_) < (strip_ - left->barycenter()) ? right->firstStrip() : left->firstStrip();
0122       cluster = clustermap.find(std::make_pair(id, firstStrip));
0123     } else if (right != clustersDetSet->begin())
0124       cluster = clustermap.find(std::make_pair(id, left->firstStrip()));
0125     else
0126       cluster = clustermap.find(std::make_pair(id, right->firstStrip()));
0127   }
0128   return cluster;
0129 }