DTSegment/src/DTClusterer.cc

0001 /******* \class DTClusterer *******
0002  *
0003  * Description:
0004  *
0005  *  detailed description
0006  *
0007  * \author : Stefano Lacaprara - INFN LNL <stefano.lacaprara@pd.infn.it>
0008  *
0009  * Modification:
0010  *
0011  *********************************/
0012
0013 /* This Class Header */
0014 #include "RecoLocalMuon/DTSegment/src/DTClusterer.h"
0015
0016 /* Collaborating Class Header */
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/ESHandle.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include "FWCore/Utilities/interface/InputTag.h"
0021 using namespace edm;
0022
0023 #include "Geometry/DTGeometry/interface/DTLayer.h"
0024 #include "Geometry/DTGeometry/interface/DTSuperLayer.h"
0025 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0026 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0027
0028 #include "DataFormats/DTRecHit/interface/DTRecHit1DPair.h"
0029 #include "DataFormats/DTRecHit/interface/DTRecClusterCollection.h"
0030 #include "DataFormats/DTRecHit/interface/DTRangeMapAccessor.h"
0031
0032 /* C++ Headers */
0033 #include <iostream>
0034 using namespace std;
0035
0036 /* ====================================================================== */
0037
0038 /* Constructor */
0039 DTClusterer::DTClusterer(const edm::ParameterSet& pset) {
0040   // Set verbose output
0041   debug = pset.getUntrackedParameter<bool>("debug");
0042
0043   // the name of the 1D rec hits collection
0044   recHits1DToken_ = consumes<DTRecHitCollection>(pset.getParameter<InputTag>("recHits1DLabel"));
0045   // min number of hits to build a cluster
0046   theMinHits = pset.getParameter<unsigned int>("minHits");
0047   // min number of hits to build a cluster
0048   theMinLayers = pset.getParameter<unsigned int>("minLayers");
0049
0050   dtGeomToken_ = esConsumes();
0051   if (debug)
0052     cout << "[DTClusterer] Constructor called" << endl;
0053
0054   produces<DTRecClusterCollection>();
0055 }
0056
0057 /* Destructor */
0058 DTClusterer::~DTClusterer() {}
0059
0060 /* Operations */
0061 void DTClusterer::produce(edm::StreamID, edm::Event& event, const edm::EventSetup& setup) const {
0062   if (debug)
0063     cout << "[DTClusterer] produce called" << endl;
0064   // Get the DT Geometry
0065   ESHandle<DTGeometry> dtGeom = setup.getHandle(dtGeomToken_);
0066
0067   // Get the 1D rechits from the event
0068   Handle<DTRecHitCollection> allHits;
0069   event.getByToken(recHits1DToken_, allHits);
0070
0071   // Create the pointer to the collection which will store the rechits
0072   auto clusters = std::make_unique<DTRecClusterCollection>();
0073
0074   // Iterate through all hit collections ordered by LayerId
0075   DTRecHitCollection::id_iterator dtLayerIt;
0076   DTSuperLayerId oldSlId;
0077   for (dtLayerIt = allHits->id_begin(); dtLayerIt != allHits->id_end(); ++dtLayerIt) {
0078     // The layerId
0079     DTLayerId layerId = (*dtLayerIt);
0080     const DTSuperLayerId SLId = layerId.superlayerId();
0081     if (SLId == oldSlId)
0082       continue;  // I'm on the same SL as before
0083     oldSlId = SLId;
0084
0085     if (debug)
0086       cout << "Reconstructing the clusters in " << SLId << endl;
0087
0088     const DTSuperLayer* sl = dtGeom->superLayer(SLId);
0089
0090     // Get all the rec hit in the same superLayer in which layerId relies
0091     DTRecHitCollection::range range = allHits->get(DTRangeMapAccessor::layersBySuperLayer(SLId));
0092
0093     // Fill the vector with the 1D RecHit
0094     vector<DTRecHit1DPair> pairs(range.first, range.second);
0095     if (debug)
0096       cout << "Number of 1D-RecHit pairs " << pairs.size() << endl;
0097     vector<DTSLRecCluster> clus = buildClusters(sl, pairs);
0098     if (debug)
0099       cout << "Number of clusters build " << clus.size() << endl;
0100     if (!clus.empty())
0101       clusters->put(sl->id(), clus.begin(), clus.end());
0102   }
0103   event.put(std::move(clusters));
0104 }
0105
0106 vector<DTSLRecCluster> DTClusterer::buildClusters(const DTSuperLayer* sl, vector<DTRecHit1DPair>& pairs) const {
0107   // create a vector of hits with wire position in SL frame
0108   vector<pair<float, DTRecHit1DPair> > hits = initHits(sl, pairs);
0109
0110   vector<DTSLRecCluster> result;
0111   // loop over pairs
0112   vector<DTRecHit1DPair> adiacentPairs;
0113   float lastPos = hits.front().first;
0114   const float cellWidth = 4.2;  // cm
0115   float sum = 0.;
0116   float sum2 = 0.;
0117
0118   for (vector<pair<float, DTRecHit1DPair> >::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
0119     if (debug)
0120       cout << "Hit: " << (*hit).first << " lastPos: " << lastPos << endl;
0121     // start from first hits
0122     // two cells are adiacente if their position is closer than cell width
0123     if (abs((*hit).first - lastPos) > cellWidth) {
0124       if (adiacentPairs.size() >= theMinHits && differentLayers(adiacentPairs) >= theMinLayers) {
0125         // if not, build the cluster with so far collection hits and go on
0126         float mean = sum / adiacentPairs.size();
0127         float err2 = sum2 / adiacentPairs.size() - mean * mean;
0128         DTSLRecCluster cluster(sl->id(), LocalPoint(mean, 0., 0.), LocalError(err2, 0., 0.), adiacentPairs);
0129         if (debug)
0130           cout << "Cluster " << cluster << endl;
0131         result.push_back(cluster);
0132       }
0133       // clean the vector
0134       adiacentPairs.clear();
0135       sum = 0.;
0136       sum2 = 0.;
0137     }
0138     // if adiacente, add them to a vector
0139     adiacentPairs.push_back((*hit).second);
0140     if (debug)
0141       cout << "adiacentPairs " << adiacentPairs.size() << endl;
0142     sum += (*hit).first;
0143     sum2 += (*hit).first * (*hit).first;
0144
0145     lastPos = (*hit).first;
0146   }
0147   // build the last cluster
0148   if (adiacentPairs.size() >= theMinHits && differentLayers(adiacentPairs) >= theMinLayers) {
0149     float mean = sum / adiacentPairs.size();
0150     float err2 = sum2 / adiacentPairs.size() - mean * mean;
0151     DTSLRecCluster cluster(sl->id(), LocalPoint(mean, 0., 0.), LocalError(err2, 0., 0.), adiacentPairs);
0152     if (debug)
0153       cout << "Cluster " << cluster << endl;
0154     result.push_back(cluster);
0155   }
0156
0157   return result;
0158 }
0159
0160 vector<pair<float, DTRecHit1DPair> > DTClusterer::initHits(const DTSuperLayer* sl,
0161                                                            vector<DTRecHit1DPair>& pairs) const {
0162   vector<pair<float, DTRecHit1DPair> > result;
0163   for (vector<DTRecHit1DPair>::const_iterator pair = pairs.begin(); pair != pairs.end(); ++pair) {
0164     // get wire
0165     DTWireId wid = (*pair).wireId();
0166     // get Layer
0167     const DTLayer* lay = sl->layer(wid.layer());
0168     // get wire position in SL (only x)
0169     LocalPoint posInLayer(lay->specificTopology().wirePosition(wid.wire()), 0., 0.);
0170     LocalPoint posInSL = sl->toLocal(lay->toGlobal(posInLayer));
0171     // put the pair into result
0172     result.push_back(make_pair(posInSL.x(), *pair));
0173   }
0174   // sorted by x
0175   sort(result.begin(), result.end(), sortClusterByX());
0176
0177   return result;
0178 }
0179
0180 unsigned int DTClusterer::differentLayers(vector<DTRecHit1DPair>& hits) const {
0181   // Count the number of different layers
0182   int layers = 0;
0183   unsigned int result = 0;
0184   for (vector<DTRecHit1DPair>::const_iterator hit = hits.begin(); hit != hits.end(); ++hit) {
0185     int pos = (1 << ((*hit).wireId().layer() - 1));
0186     if (!(pos & layers)) {
0187       result++;
0188       layers |= pos;
0189     }
0190   }
0191   return result;
0192 }