CMS 3D CMS Logo

DTClusterer.cc
Go to the documentation of this file.
1 /******* \class DTClusterer *******
2  *
3  * Description:
4  *
5  * detailed description
6  *
7  * \author : Stefano Lacaprara - INFN LNL <stefano.lacaprara@pd.infn.it>
8  *
9  * Modification:
10  *
11  *********************************/
12 
13 /* This Class Header */
15 
16 /* Collaborating Class Header */
21 using namespace edm;
22 
27 
31 
32 /* C++ Headers */
33 #include <iostream>
34 using namespace std;
35 
36 /* ====================================================================== */
37 
38 /* Constructor */
40  // Set verbose output
41  debug = pset.getUntrackedParameter<bool>("debug");
42 
43  // the name of the 1D rec hits collection
44  recHits1DToken_ = consumes<DTRecHitCollection>(pset.getParameter<InputTag>("recHits1DLabel"));
45  // min number of hits to build a cluster
46  theMinHits = pset.getParameter<unsigned int>("minHits");
47  // min number of hits to build a cluster
48  theMinLayers = pset.getParameter<unsigned int>("minLayers");
49 
50  if(debug)
51  cout << "[DTClusterer] Constructor called" << endl;
52 
53  produces<DTRecClusterCollection>();
54 }
55 
56 /* Destructor */
58 }
59 
60 /* Operations */
62  if(debug)
63  cout << "[DTClusterer] produce called" << endl;
64  // Get the DT Geometry
65  ESHandle<DTGeometry> dtGeom;
66  setup.get<MuonGeometryRecord>().get(dtGeom);
67 
68  // Get the 1D rechits from the event
70  event.getByToken(recHits1DToken_, allHits);
71 
72  // Create the pointer to the collection which will store the rechits
73  auto clusters = std::make_unique<DTRecClusterCollection>();
74 
75  // Iterate through all hit collections ordered by LayerId
76  DTRecHitCollection::id_iterator dtLayerIt;
77  DTSuperLayerId oldSlId;
78  for (dtLayerIt = allHits->id_begin(); dtLayerIt != allHits->id_end(); ++dtLayerIt){
79  // The layerId
80  DTLayerId layerId = (*dtLayerIt);
81  const DTSuperLayerId SLId = layerId.superlayerId();
82  if (SLId==oldSlId) continue; // I'm on the same SL as before
83  oldSlId = SLId;
84 
85  if(debug) cout <<"Reconstructing the clusters in "<< SLId << endl;
86 
87  const DTSuperLayer* sl = dtGeom->superLayer(SLId);
88 
89  // Get all the rec hit in the same superLayer in which layerId relies
91  allHits->get(DTRangeMapAccessor::layersBySuperLayer(SLId));
92 
93  // Fill the vector with the 1D RecHit
94  vector<DTRecHit1DPair> pairs(range.first,range.second);
95  if(debug) cout << "Number of 1D-RecHit pairs " << pairs.size() << endl;
96  vector<DTSLRecCluster> clus=buildClusters(sl, pairs);
97  if(debug) cout << "Number of clusters build " << clus.size() << endl;
98  if (!clus.empty() )
99  clusters->put(sl->id(), clus.begin(), clus.end());
100  }
101  event.put(std::move(clusters));
102 }
103 
104 vector<DTSLRecCluster> DTClusterer::buildClusters(const DTSuperLayer* sl,
105  vector<DTRecHit1DPair>& pairs) {
106  // create a vector of hits with wire position in SL frame
107  vector<pair<float, DTRecHit1DPair> > hits = initHits(sl, pairs);
108 
109  vector<DTSLRecCluster> result;
110  // loop over pairs
111  vector<DTRecHit1DPair> adiacentPairs;
112  float lastPos=hits.front().first;
113  const float cellWidth=4.2; // cm
114  float sum=0.;
115  float sum2=0.;
116 
117  for (vector<pair<float, DTRecHit1DPair> >::const_iterator hit=hits.begin(); hit!=hits.end();
118  ++hit) {
119  if(debug) cout << "Hit: " << (*hit).first << " lastPos: " << lastPos << endl;
120  // start from first hits
121  // two cells are adiacente if their position is closer than cell width
122  if(abs((*hit).first-lastPos)>cellWidth ) {
123  if (adiacentPairs.size() >= theMinHits && differentLayers(adiacentPairs)>=theMinLayers) {
124  // if not, build the cluster with so far collection hits and go on
125  float mean=sum/adiacentPairs.size();
126  float err2=sum2/adiacentPairs.size()-mean*mean;
127  DTSLRecCluster cluster(sl->id(),
128  LocalPoint(mean,0.,0.),
129  LocalError(err2,0.,0.),
130  adiacentPairs);
131  if(debug) cout << "Cluster " << cluster << endl;
132  result.push_back(cluster);
133  }
134  // clean the vector
135  adiacentPairs.clear();
136  sum=0.;
137  sum2=0.;
138  }
139  // if adiacente, add them to a vector
140  adiacentPairs.push_back((*hit).second);
141  if(debug) cout << "adiacentPairs " << adiacentPairs.size() << endl;
142  sum+=(*hit).first;
143  sum2+=(*hit).first*(*hit).first;
144 
145  lastPos=(*hit).first;
146  }
147  // build the last cluster
148  if (adiacentPairs.size() >= theMinHits && differentLayers(adiacentPairs)>=theMinLayers) {
149  float mean=sum/adiacentPairs.size();
150  float err2=sum2/adiacentPairs.size()-mean*mean;
151  DTSLRecCluster cluster(sl->id(),
152  LocalPoint(mean,0.,0.),
153  LocalError(err2,0.,0.),
154  adiacentPairs);
155  if(debug) cout << "Cluster " << cluster << endl;
156  result.push_back(cluster);
157  }
158 
159  return result;
160 }
161 
162 vector<pair<float, DTRecHit1DPair> > DTClusterer::initHits(const DTSuperLayer* sl,
163  vector<DTRecHit1DPair>& pairs) {
164  vector<pair<float, DTRecHit1DPair> > result;
165  for (vector<DTRecHit1DPair>::const_iterator pair=pairs.begin();
166  pair!=pairs.end(); ++pair ) {
167 
168  // get wire
169  DTWireId wid = (*pair).wireId();
170  // get Layer
171  const DTLayer* lay= sl->layer(wid.layer());
172  // get wire position in SL (only x)
173  LocalPoint
174  posInLayer(lay->specificTopology().wirePosition(wid.wire()),0.,0.);
175  LocalPoint posInSL = sl->toLocal(lay->toGlobal(posInLayer));
176  // put the pair into result
177  result.push_back(make_pair(posInSL.x(), *pair));
178  }
179  // sorted by x
180  sort(result.begin(), result.end(), sortClusterByX());
181 
182  return result;
183 }
184 
185 unsigned int DTClusterer::differentLayers(vector<DTRecHit1DPair>& hits) {
186 // Count the number of different layers
187  int layers=0;
188  unsigned int result=0;
189  for (vector<DTRecHit1DPair>::const_iterator hit=hits.begin(); hit!=hits.end();
190  ++hit) {
191  int pos = (1 << ((*hit).wireId().layer() -1 ));
192  if (! (pos & layers)) {
193  result++;
194  layers |= pos;
195  }
196  }
197  return result;
198 }
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
float wirePosition(int wireNumber) const
Returns the x position in the layer of a given wire number.
Definition: DTTopology.cc:86
unsigned int differentLayers(std::vector< DTRecHit1DPair > &hits)
Definition: DTClusterer.cc:185
std::pair< const_iterator, const_iterator > range
iterator range
Definition: RangeMap.h:50
Point3DBase< Scalar, LocalTag > LocalPoint
Definition: Definitions.h:32
std::vector< LayerSetAndLayers > layers(const SeedingLayerSetsHits &sets)
Definition: LayerTriplets.cc:4
~DTClusterer() override
Definition: DTClusterer.cc:57
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:54
static std::pair< DTLayerId, DTSuperLayerIdComparator > layersBySuperLayer(DTSuperLayerId slId)
Access by SL objects written into a RangeMap by layer.
void produce(edm::Event &event, const edm::EventSetup &setup) override
Definition: DTClusterer.cc:61
const DTLayer * layer(const DTLayerId &id) const
Return the layer corresponding to the given id.
Definition: DTSuperLayer.cc:68
def setup(process, global_tag, zero_tesla=False)
Definition: GeneralSetup.py:1
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:69
int layer() const
Return the layer number.
Definition: DTLayerId.h:53
DTSuperLayerId superlayerId() const
Return the corresponding SuperLayerId.
Definition: DTLayerId.h:59
DTSuperLayerId id() const
Return the DetId of this SL.
Definition: DTSuperLayer.cc:36
const DTTopology & specificTopology() const
Definition: DTLayer.cc:42
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< DTSLRecCluster > buildClusters(const DTSuperLayer *sl, std::vector< DTRecHit1DPair > &pairs)
Definition: DTClusterer.cc:104
int wire() const
Return the wire number.
Definition: DTWireId.h:56
DTClusterer(const edm::ParameterSet &)
Definition: DTClusterer.cc:39
#define debug
Definition: HDRShower.cc:19
const T & get() const
Definition: EventSetup.h:55
HLT enums.
std::vector< std::pair< float, DTRecHit1DPair > > initHits(const DTSuperLayer *sl, std::vector< DTRecHit1DPair > &pairs)
Definition: DTClusterer.cc:162
T x() const
Definition: PV3DBase.h:62
def move(src, dest)
Definition: eostools.py:510
Definition: event.py:1
const DTSuperLayer * superLayer(const DTSuperLayerId &id) const
Return a DTSuperLayer given its id.
Definition: DTGeometry.cc:104