FTLClusterizer/src/MTDThresholdClusterizer.cc

0001 // Our own includes
0002 #include "RecoLocalFastTime/FTLClusterizer/interface/MTDArrayBuffer.h"
0003 #include "RecoLocalFastTime/FTLClusterizer/interface/MTDThresholdClusterizer.h"
0004 #include "RecoLocalFastTime/FTLClusterizer/interface/BTLRecHitsErrorEstimatorIM.h"
0005
0006 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0007 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0008
0009 #include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h"
0010 #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
0011 #include "Geometry/MTDCommonData/interface/MTDTopologyMode.h"
0012
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015 // STL
0016 #include <stack>
0017 #include <vector>
0018 #include <iostream>
0019 #include <atomic>
0020 #include <cmath>
0021 using namespace std;
0022
0023 //----------------------------------------------------------------------------
0024 //! Constructor:
0025 //----------------------------------------------------------------------------
0026 MTDThresholdClusterizer::MTDThresholdClusterizer(edm::ParameterSet const& conf)
0027     :  // Get energy thresholds
0028       theHitThreshold(conf.getParameter<double>("HitThreshold")),
0029       theSeedThreshold(conf.getParameter<double>("SeedThreshold")),
0030       theClusterThreshold(conf.getParameter<double>("ClusterThreshold")),
0031       theTimeThreshold(conf.getParameter<double>("TimeThreshold")),
0032       thePositionThreshold(conf.getParameter<double>("PositionThreshold")),
0033       theNumOfRows(0),
0034       theNumOfCols(0),
0035       theCurrentId(0),
0036       theBuffer(theNumOfRows, theNumOfCols),
0037       bufferAlreadySet(false) {}
0038 /////////////////////////////////////////////////////////////////////////////
0039 MTDThresholdClusterizer::~MTDThresholdClusterizer() {}
0040
0041 // Configuration descriptions
0042 void MTDThresholdClusterizer::fillPSetDescription(edm::ParameterSetDescription& desc) {
0043   desc.add<double>("HitThreshold", 0.);
0044   desc.add<double>("SeedThreshold", 0.);
0045   desc.add<double>("ClusterThreshold", 0.);
0046   desc.add<double>("TimeThreshold", 10.);
0047   desc.add<double>("PositionThreshold", -1.0);
0048 }
0049
0050 //----------------------------------------------------------------------------
0051 //!  Prepare the Clusterizer to work on a particular DetUnit.  Re-init the
0052 //!  size of the panel/plaquette (so update nrows and ncols),
0053 //----------------------------------------------------------------------------
0054 bool MTDThresholdClusterizer::setup(const MTDGeometry* geom, const MTDTopology* topo, const DetId& id) {
0055   theCurrentId = id;
0056   //using geopraphicalId here
0057   const auto& thedet = geom->idToDet(id);
0058   if (thedet == nullptr) {
0059     throw cms::Exception("MTDThresholdClusterizer")
0060         << "GeographicalID: " << std::hex << id.rawId() << " is invalid!" << std::dec << std::endl;
0061     return false;
0062   }
0063   const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
0064   const RectangularMTDTopology& topol = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0065
0066   // Get the new sizes.
0067   unsigned int nrows = topol.nrows();     // rows in x
0068   unsigned int ncols = topol.ncolumns();  // cols in y
0069
0070   theNumOfRows = nrows;  // Set new sizes
0071   theNumOfCols = ncols;
0072
0073   LogDebug("MTDThresholdClusterizer") << "Buffer size [" << theNumOfRows + 1 << "," << theNumOfCols + 1 << "]";
0074
0075   if (nrows > theBuffer.rows() || ncols > theBuffer.columns()) {  // change only when a larger is needed
0076     // Resize the buffer
0077     theBuffer.setSize(nrows + 1, ncols + 1);  // +1 needed for MTD
0078     bufferAlreadySet = true;
0079   }
0080
0081   return true;
0082 }
0083 //----------------------------------------------------------------------------
0084 //!  \brief Cluster hits.
0085 //!  This method operates on a matrix of hits
0086 //!  and finds the largest contiguous cluster around
0087 //!  each seed hit.
0088 //!  Input and output data stored in DetSet
0089 //----------------------------------------------------------------------------
0090 void MTDThresholdClusterizer::clusterize(const FTLRecHitCollection& input,
0091                                          const MTDGeometry* geom,
0092                                          const MTDTopology* topo,
0093                                          FTLClusterCollection& output) {
0094   FTLRecHitCollection::const_iterator begin = input.begin();
0095   FTLRecHitCollection::const_iterator end = input.end();
0096
0097   // Do not bother for empty detectors
0098   if (begin == end) {
0099     edm::LogInfo("MTDThresholdClusterizer") << "No hits to clusterize";
0100     return;
0101   }
0102
0103   LogDebug("MTDThresholdClusterizer") << "Input collection " << input.size();
0104   assert(output.empty());
0105
0106   std::set<unsigned> geoIds;
0107   std::multimap<unsigned, unsigned> geoIdToIdx;
0108
0109   unsigned index = 0;
0110   for (const auto& hit : input) {
0111     MTDDetId mtdId = MTDDetId(hit.detid());
0112     if (mtdId.subDetector() != MTDDetId::FastTime) {
0113       throw cms::Exception("MTDThresholdClusterizer")
0114           << "MTDDetId: " << std::hex << mtdId.rawId() << " is invalid!" << std::dec << std::endl;
0115     }
0116
0117     if (mtdId.mtdSubDetector() == MTDDetId::BTL) {
0118       BTLDetId hitId(hit.detid());
0119       //for BTL topology gives different layout id
0120       DetId geoId = hitId.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topo->getMTDTopologyMode()));
0121       geoIdToIdx.emplace(geoId, index);
0122       geoIds.emplace(geoId);
0123       ++index;
0124     } else if (mtdId.mtdSubDetector() == MTDDetId::ETL) {
0125       ETLDetId hitId(hit.detid());
0126       DetId geoId = hitId.geographicalId();
0127       geoIdToIdx.emplace(geoId, index);
0128       geoIds.emplace(geoId);
0129       ++index;
0130     } else {
0131       throw cms::Exception("MTDThresholdClusterizer")
0132           << "MTDDetId: " << std::hex << mtdId.rawId() << " is invalid!" << std::dec << std::endl;
0133     }
0134   }
0135
0136   //cluster hits within geoIds (modules)
0137   for (unsigned id : geoIds) {
0138     //  Set up the clusterization on this DetId.
0139     if (!setup(geom, topo, DetId(id)))
0140       return;
0141
0142     auto range = geoIdToIdx.equal_range(id);
0143     LogDebug("MTDThresholdClusterizer") << "Matching Ids for " << std::hex << id << std::dec << " ["
0144                                         << range.first->second << "," << range.second->second << "]";
0145
0146     //  Copy MTDRecHits to the buffer array; select the seed hits
0147     //  on the way, and store them in theSeeds.
0148     for (auto itr = range.first; itr != range.second; ++itr) {
0149       const unsigned hitidx = itr->second;
0150       copy_to_buffer(begin + hitidx, geom, topo);
0151     }
0152
0153     FTLClusterCollection::FastFiller clustersOnDet(output, id);
0154
0155     for (unsigned int i = 0; i < theSeeds.size(); i++) {
0156       if (theBuffer.energy(theSeeds[i]) > theSeedThreshold) {  // Is this seed still valid?
0157         //  Make a cluster around this seed
0158         const FTLCluster& cluster = make_cluster(theSeeds[i]);
0159
0160         //  Check if the cluster is above threshold
0161         if (cluster.energy() > theClusterThreshold) {
0162           LogDebug("MTDThresholdClusterizer")
0163               << "putting in this cluster " << i << " #hits:" << cluster.size() << " E:" << cluster.energy()
0164               << " T:" << cluster.time() << " X:" << cluster.x() << " Y:" << cluster.y();
0165           clustersOnDet.push_back(cluster);
0166         }
0167       }
0168     }
0169
0170     // Erase the seeds.
0171     theSeeds.clear();
0172     //  Need to clean unused hits from the buffer array.
0173     for (auto itr = range.first; itr != range.second; ++itr) {
0174       const unsigned hitidx = itr->second;
0175       clear_buffer(begin + hitidx);
0176     }
0177   }
0178 }
0179
0180 //----------------------------------------------------------------------------
0181 //!  \brief Clear the internal buffer array.
0182 //!
0183 //!  MTDs which are not part of recognized clusters are NOT ERASED
0184 //!  during the cluster finding.  Erase them now.
0185 //!
0186 //----------------------------------------------------------------------------
0187 void MTDThresholdClusterizer::clear_buffer(RecHitIterator itr) { theBuffer.clear(itr->row(), itr->column()); }
0188
0189 //----------------------------------------------------------------------------
0190 //! \brief Copy FTLRecHit into the buffer, identify seeds.
0191 //----------------------------------------------------------------------------
0192 void MTDThresholdClusterizer::copy_to_buffer(RecHitIterator itr, const MTDGeometry* geom, const MTDTopology* topo) {
0193   MTDDetId mtdId = MTDDetId(itr->detid());
0194   int row = itr->row();
0195   int col = itr->column();
0196   GeomDetEnumerators::Location subDet = GeomDetEnumerators::invalidLoc;
0197   float energy = itr->energy();
0198   float time = itr->time();
0199   float timeError = itr->timeError();
0200   float position = itr->position();
0201   // position is the longitudinal offset that should be added into local x for bars in phi geometry
0202   LocalError local_error(0, 0, 0);
0203   GlobalPoint global_point(0, 0, 0);
0204   if (mtdId.mtdSubDetector() == MTDDetId::BTL) {
0205     subDet = GeomDetEnumerators::barrel;
0206     BTLDetId id = itr->id();
0207     DetId geoId = id.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topo->getMTDTopologyMode()));
0208     const auto& det = geom->idToDet(geoId);
0209     const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(det->topology());
0210     const RectangularMTDTopology& topol = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0211     MeasurementPoint mp(row, col);
0212     LocalPoint lp_ctr = topol.localPosition(mp);
0213     LocalPoint lp(lp_ctr.x() + position + topol.pitch().first * 0.5f, lp_ctr.y(), lp_ctr.z());
0214     // local coordinates of BTL module locates RecHits on the left edge of the bar (-9.2, -3.067, 3.067)
0215     // (position + topol.pitch().first/2.0) is the distance from the left edge to the Hit point
0216     global_point = det->toGlobal(lp);
0217     BTLRecHitsErrorEstimatorIM btlError(det, lp);
0218     local_error = btlError.localError();
0219   } else if (mtdId.mtdSubDetector() == MTDDetId::ETL) {
0220     subDet = GeomDetEnumerators::endcap;
0221     ETLDetId id = itr->id();
0222     DetId geoId = id.geographicalId();
0223     const auto& det = geom->idToDet(geoId);
0224     const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(det->topology());
0225     const RectangularMTDTopology& topol = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0226
0227     MeasurementPoint mp(row, col);
0228     LocalPoint lp = topol.localPosition(mp);
0229     global_point = det->toGlobal(lp);
0230   }
0231
0232   LogDebug("MTDThresholdClusterizer") << "ROW " << row << " COL " << col << " ENERGY " << energy << " TIME " << time;
0233   if (energy > theHitThreshold) {
0234     theBuffer.set(row, col, subDet, energy, time, timeError, local_error, global_point);
0235     if (energy > theSeedThreshold)
0236       theSeeds.push_back(FTLCluster::FTLHitPos(row, col));
0237     //sort seeds?
0238   }
0239 }
0240
0241 //----------------------------------------------------------------------------
0242 //!  \brief The actual clustering algorithm: group the neighboring hits around the seed.
0243 //----------------------------------------------------------------------------
0244 FTLCluster MTDThresholdClusterizer::make_cluster(const FTLCluster::FTLHitPos& hit) {
0245   //First we acquire the seeds for the clusters
0246
0247   GeomDetEnumerators::Location seed_subdet = theBuffer.subDet(hit.row(), hit.col());
0248   float seed_energy = theBuffer.energy(hit.row(), hit.col());
0249   float seed_time = theBuffer.time(hit.row(), hit.col());
0250   float seed_time_error = theBuffer.time_error(hit.row(), hit.col());
0251   auto const seedPoint = theBuffer.global_point(hit.row(), hit.col());
0252   double seed_error_xx = theBuffer.local_error(hit.row(), hit.col()).xx();
0253   double seed_error_yy = theBuffer.local_error(hit.row(), hit.col()).yy();
0254   theBuffer.clear(hit);
0255
0256   AccretionCluster acluster;
0257   acluster.add(hit, seed_energy, seed_time, seed_time_error);
0258
0259   bool stopClus = false;
0260   //Here we search all hits adjacent to all hits in the cluster.
0261   while (!acluster.empty() && !stopClus) {
0262     //This is the standard algorithm to find and add a hit
0263     auto curInd = acluster.top();
0264     acluster.pop();
0265     for (auto c = std::max(0, int(acluster.y[curInd] - 1));
0266          c < std::min(int(acluster.y[curInd] + 2), int(theBuffer.columns())) && !stopClus;
0267          ++c) {
0268       for (auto r = std::max(0, int(acluster.x[curInd] - 1));
0269            r < std::min(int(acluster.x[curInd] + 2), int(theBuffer.rows())) && !stopClus;
0270            ++r) {
0271         if (theBuffer.energy(r, c) > theHitThreshold) {
0272           if (std::abs(theBuffer.time(r, c) - seed_time) >
0273               theTimeThreshold *
0274                   sqrt(theBuffer.time_error(r, c) * theBuffer.time_error(r, c) + seed_time_error * seed_time_error))
0275             continue;
0276           if ((seed_subdet == GeomDetEnumerators::barrel) && (theBuffer.subDet(r, c) == GeomDetEnumerators::barrel)) {
0277             double hit_error_xx = theBuffer.local_error(r, c).xx();
0278             double hit_error_yy = theBuffer.local_error(r, c).yy();
0279             if (thePositionThreshold > 0) {
0280               if (((theBuffer.global_point(r, c) - seedPoint).mag2()) >
0281                   thePositionThreshold * thePositionThreshold *
0282                       (hit_error_xx + seed_error_xx + hit_error_yy + seed_error_yy))
0283                 continue;
0284             }
0285           }
0286           FTLCluster::FTLHitPos newhit(r, c);
0287           if (!acluster.add(newhit, theBuffer.energy(r, c), theBuffer.time(r, c), theBuffer.time_error(r, c))) {
0288             stopClus = true;
0289             break;
0290           }
0291           theBuffer.clear(newhit);
0292         }
0293       }
0294     }
0295   }  // while accretion
0296
0297   FTLCluster cluster(theCurrentId,
0298                      acluster.isize,
0299                      acluster.energy.data(),
0300                      acluster.time.data(),
0301                      acluster.timeError.data(),
0302                      acluster.x.data(),
0303                      acluster.y.data(),
0304                      acluster.xmin,
0305                      acluster.ymin);
0306   return cluster;
0307 }