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, ncols);
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 +/- DT:" << cluster.time() << " +/- " << cluster.timeError() << " X:" << cluster.x()
0165               << " Y:" << cluster.y() << " xpos +/- err " << cluster.getClusterPosX() << " +/- "
0166               << cluster.getClusterErrorX();
0167           clustersOnDet.push_back(cluster);
0168         }
0169       }
0170     }
0171
0172     // Erase the seeds.
0173     theSeeds.clear();
0174     //  Need to clean unused hits from the buffer array.
0175     for (auto itr = range.first; itr != range.second; ++itr) {
0176       const unsigned hitidx = itr->second;
0177       clear_buffer(begin + hitidx);
0178     }
0179   }
0180 }
0181
0182 //----------------------------------------------------------------------------
0183 //!  \brief Clear the internal buffer array.
0184 //!
0185 //!  MTDs which are not part of recognized clusters are NOT ERASED
0186 //!  during the cluster finding.  Erase them now.
0187 //!
0188 //----------------------------------------------------------------------------
0189 void MTDThresholdClusterizer::clear_buffer(RecHitIterator itr) { theBuffer.clear(itr->row(), itr->column()); }
0190
0191 //----------------------------------------------------------------------------
0192 //! \brief Copy FTLRecHit into the buffer, identify seeds.
0193 //----------------------------------------------------------------------------
0194 void MTDThresholdClusterizer::copy_to_buffer(RecHitIterator itr, const MTDGeometry* geom, const MTDTopology* topo) {
0195   MTDDetId mtdId = MTDDetId(itr->detid());
0196   int row = itr->row();
0197   int col = itr->column();
0198   GeomDetEnumerators::Location subDet = GeomDetEnumerators::invalidLoc;
0199   float energy = itr->energy();
0200   float time = itr->time();
0201   float timeError = itr->timeError();
0202   float position = itr->position();
0203   float xpos = 0.f;
0204   // position is the longitudinal offset that should be added into local x for bars in phi geometry
0205   LocalPoint local_point(0, 0, 0);
0206   LocalError local_error(0, 0, 0);
0207   if (mtdId.mtdSubDetector() == MTDDetId::BTL) {
0208     subDet = GeomDetEnumerators::barrel;
0209     BTLDetId id = itr->id();
0210     DetId geoId = id.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topo->getMTDTopologyMode()));
0211     const auto& det = geom->idToDet(geoId);
0212     const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(det->topology());
0213     const RectangularMTDTopology& topol = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0214
0215     LocalPoint lp_pixel(position, 0, 0);
0216     local_point = topol.pixelToModuleLocalPoint(lp_pixel, row, col);
0217     BTLRecHitsErrorEstimatorIM btlError(det, local_point);
0218     local_error = btlError.localError();
0219     xpos = local_point.x();
0220   } else if (mtdId.mtdSubDetector() == MTDDetId::ETL) {
0221     subDet = GeomDetEnumerators::endcap;
0222     ETLDetId id = itr->id();
0223     DetId geoId = id.geographicalId();
0224     const auto& det = geom->idToDet(geoId);
0225     const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(det->topology());
0226     const RectangularMTDTopology& topol = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0227
0228     LocalPoint lp_pixel(0, 0, 0);
0229     local_point = topol.pixelToModuleLocalPoint(lp_pixel, row, col);
0230   }
0231
0232   LogDebug("MTDThresholdClusterizer") << "DetId " << mtdId.rawId() << " subd " << mtdId.mtdSubDetector() << " row/col "
0233                                       << row << " / " << col << " energy " << energy << " time " << time
0234                                       << " time error " << timeError << " local_point " << local_point.x() << " "
0235                                       << local_point.y() << " " << local_point.z() << " local error "
0236                                       << std::sqrt(local_error.xx()) << " " << std::sqrt(local_error.yy()) << " xpos "
0237                                       << xpos;
0238   if (energy > theHitThreshold) {
0239     theBuffer.set(row, col, subDet, energy, time, timeError, local_error, local_point, xpos);
0240     if (energy > theSeedThreshold)
0241       theSeeds.push_back(FTLCluster::FTLHitPos(row, col));
0242     //sort seeds?
0243   }
0244 }
0245
0246 //----------------------------------------------------------------------------
0247 //!  \brief The actual clustering algorithm: group the neighboring hits around the seed.
0248 //----------------------------------------------------------------------------
0249 FTLCluster MTDThresholdClusterizer::make_cluster(const FTLCluster::FTLHitPos& hit) {
0250   //First we acquire the seeds for the clusters
0251
0252   GeomDetEnumerators::Location seed_subdet = theBuffer.subDet(hit.row(), hit.col());
0253   float seed_energy = theBuffer.energy(hit.row(), hit.col());
0254   float seed_time = theBuffer.time(hit.row(), hit.col());
0255   float seed_time_error = theBuffer.time_error(hit.row(), hit.col());
0256   auto const seedPoint = theBuffer.local_point(hit.row(), hit.col());
0257   double seed_error_xx = theBuffer.local_error(hit.row(), hit.col()).xx();
0258   double seed_error_yy = theBuffer.local_error(hit.row(), hit.col()).yy();
0259   double seed_xpos = theBuffer.xpos(hit.row(), hit.col());
0260   theBuffer.clear(hit);
0261
0262   AccretionCluster acluster;
0263   acluster.add(hit, seed_energy, seed_time, seed_time_error);
0264
0265   // for BTL position along crystals add auxiliary vectors
0266   std::array<float, AccretionCluster::MAXSIZE> pixel_x{{-99999.}};
0267   std::array<float, AccretionCluster::MAXSIZE> pixel_errx2{{-99999.}};
0268   if (seed_subdet == GeomDetEnumerators::barrel) {
0269     pixel_x[acluster.top()] = seed_xpos;
0270     pixel_errx2[acluster.top()] = seed_error_xx;
0271   }
0272
0273   bool stopClus = false;
0274   //Here we search all hits adjacent to all hits in the cluster.
0275   while (!acluster.empty() && !stopClus) {
0276     //This is the standard algorithm to find and add a hit
0277     auto curInd = acluster.top();
0278     acluster.pop();
0279     for (auto c = std::max(0, int(acluster.y[curInd] - 1));
0280          c < std::min(int(acluster.y[curInd] + 2), int(theBuffer.columns())) && !stopClus;
0281          ++c) {
0282       for (auto r = std::max(0, int(acluster.x[curInd] - 1));
0283            r < std::min(int(acluster.x[curInd] + 2), int(theBuffer.rows())) && !stopClus;
0284            ++r) {
0285         LogDebug("MTDThresholdClusterizer")
0286             << "Clustering subdet " << seed_subdet << " around " << curInd << " X,Y = " << acluster.x[curInd] << ","
0287             << acluster.y[curInd] << " r,c = " << r << "," << c << " energy,time = " << theBuffer.energy(r, c) << " "
0288             << theBuffer.time(r, c);
0289         if (theBuffer.energy(r, c) > theHitThreshold) {
0290           if (std::abs(theBuffer.time(r, c) - seed_time) >
0291               theTimeThreshold *
0292                   sqrt(theBuffer.time_error(r, c) * theBuffer.time_error(r, c) + seed_time_error * seed_time_error))
0293             continue;
0294           if ((seed_subdet == GeomDetEnumerators::barrel) && (theBuffer.subDet(r, c) == GeomDetEnumerators::barrel) &&
0295               (thePositionThreshold > 0)) {
0296             double hit_error_xx = theBuffer.local_error(r, c).xx();
0297             double hit_error_yy = theBuffer.local_error(r, c).yy();
0298             if (((theBuffer.local_point(r, c) - seedPoint).mag2()) >
0299                 thePositionThreshold * thePositionThreshold *
0300                     (hit_error_xx + seed_error_xx + hit_error_yy + seed_error_yy))
0301               continue;
0302           }
0303           FTLCluster::FTLHitPos newhit(r, c);
0304           if (!acluster.add(newhit, theBuffer.energy(r, c), theBuffer.time(r, c), theBuffer.time_error(r, c))) {
0305             stopClus = true;
0306             break;
0307           }
0308           if (theBuffer.subDet(r, c) == GeomDetEnumerators::barrel) {
0309             pixel_x[acluster.top()] = theBuffer.xpos(r, c);
0310             pixel_errx2[acluster.top()] = theBuffer.local_error(r, c).xx();
0311           }
0312           theBuffer.clear(newhit);
0313         }
0314       }
0315     }
0316   }  // while accretion
0317
0318   FTLCluster cluster(theCurrentId,
0319                      acluster.isize,
0320                      acluster.energy.data(),
0321                      acluster.time.data(),
0322                      acluster.timeError.data(),
0323                      acluster.x.data(),
0324                      acluster.y.data(),
0325                      acluster.xmin,
0326                      acluster.ymin);
0327
0328   // For BTL compute the optimal position along crystal and uncertainty on it in absolute length units
0329
0330   if (seed_subdet == GeomDetEnumerators::barrel) {
0331     float sumW(0.f), sumXW(0.f), sumXW2(0.f);
0332     for (unsigned int index = 0; index < acluster.top(); index++) {
0333       sumW += acluster.energy[index];
0334       sumXW += acluster.energy[index] * pixel_x[index];
0335       sumXW2 += acluster.energy[index] * acluster.energy[index] * pixel_errx2[index];
0336     }
0337     cluster.setClusterPosX(sumXW / sumW);
0338     cluster.setClusterErrorX(std::sqrt(sumXW2) / sumW);
0339   }
0340
0341   return cluster;
0342 }