CMS 3D CMS Logo
CSCDQM_StripClusterFinder.cc
Go to the documentation of this file.
1 
2 #include "CSCDQM_StripClusterFinder.h"
3 
4 namespace cscdqm {
5 
6  StripClusterFinder::StripClusterFinder(int l, int s, int cf, int st) {
7  //
8  // Options
9  //
10  // fOpt = new CalibOptions();
11  LayerNmb = l;
12  TimeSliceNmb = s;
13  StripNmb = cf * 16;
14  }
15  void StripClusterFinder::DoAction(int LayerId, float* Cathodes) {
16  int TimeId, StripId;
17  this->LId = LayerId;
18  MEStripClusters.clear();
19  StripClusterFitData PulseHeightMapTMP;
20 
21  thePulseHeightMap.clear();
22 
23  // fill
24  //===========================================================================
25 
26  for (StripId = 0; StripId < StripNmb; StripId++) {
27  for (TimeId = 0; TimeId < TimeSliceNmb; TimeId++) {
28  PulseHeightMapTMP.height_[TimeId] = *(Cathodes + StripNmb * (TimeSliceNmb * LayerId + TimeId) + StripId);
29  }
30  PulseHeightMapTMP.bx_ = 0;
31  PulseHeightMapTMP.channel_ = StripId; // was StripId
32  thePulseHeightMap.push_back(PulseHeightMapTMP);
33  }
34  sort(thePulseHeightMap.begin(), thePulseHeightMap.end(), Sort());
35  //===========================================================================
36 
37  if (thePulseHeightMap.empty())
38  return;
39  SearchMax();
40  SearchBorders();
41  Match();
42  RefindMax();
43  /*
44  int val;
45  for(i=0;i<MEStripClusters.size();i++){
46  val=MEStripClusters[i].LFTBNDStrip;
47  MEStripClusters[i].LFTBNDStrip=thePulseHeightMap[val].channel_;
48  val=MEStripClusters[i].IRTBNDStrip;
49  MEStripClusters[i].IRTBNDStrip=thePulseHeightMap[val].channel_;
50  for(j=0;j<MEStripClusters[i].localMax.size();j++){
51  val=MEStripClusters[i].localMax[j].Strip;
52  MEStripClusters[i].localMax[j].Strip=thePulseHeightMap[val].channel_;
53  }
54  }
55  */
56 
57  float sumstrip;
58  float sumtime;
59  float sumheight;
60 
61  for (i = 0; i < MEStripClusters.size(); i++) {
62  MEStripClusters[i].ClusterPulseMapHeight.clear();
63  for (j = 0; j < thePulseHeightMap.size(); j++) {
64  if (thePulseHeightMap[j].channel_ >= MEStripClusters[i].LFTBNDStrip &&
65  thePulseHeightMap[j].channel_ <= MEStripClusters[i].IRTBNDStrip)
66  MEStripClusters[i].ClusterPulseMapHeight.push_back(thePulseHeightMap[j]);
67  }
68  sumstrip = 0;
69  sumtime = 0;
70  sumheight = 0;
71  for (uint32_t k = 0; k < MEStripClusters[i].ClusterPulseMapHeight.size(); k++) {
72  for (int l = 0; l < 16; l++) {
73  sumstrip += MEStripClusters[i].ClusterPulseMapHeight[k].height_[l] *
74  (MEStripClusters[i].ClusterPulseMapHeight[k].channel_ + 1);
75  sumtime += MEStripClusters[i].ClusterPulseMapHeight[k].height_[l] * (l + 1);
76  sumheight += MEStripClusters[i].ClusterPulseMapHeight[k].height_[l];
77  }
78  }
79  if (sumheight) {
80  MEStripClusters[i].Mean[0] = sumstrip / sumheight;
81  MEStripClusters[i].Mean[1] = sumtime / sumheight;
82  }
83  }
84  // printClusters();
85  return;
86  }
87 
88  void StripClusterFinder::SearchMax(void) {
89  StripCluster tmpCluster;
90  for (i = 1; i < (thePulseHeightMap.size() - 1); i++) {
91  if (thePulseHeightMap[i].channel_ == 63 || thePulseHeightMap[i].channel_ == 64)
92  continue;
93  for (j = 1; j < 15; j++) {
94  if (thePulseHeightMap[i].height_[j] > thePulseHeightMap[i - 1].height_[j] &&
95  thePulseHeightMap[i].height_[j] > thePulseHeightMap[i + 1].height_[j] &&
96  thePulseHeightMap[i].height_[j] > thePulseHeightMap[i].height_[j - 1] &&
97  thePulseHeightMap[i].height_[j] > thePulseHeightMap[i].height_[j + 1] &&
98  thePulseHeightMap[i].height_[j] > thePulseHeightMap[i - 1].height_[j - 1] &&
99  thePulseHeightMap[i].height_[j] > thePulseHeightMap[i - 1].height_[j + 1] &&
100  thePulseHeightMap[i].height_[j] > thePulseHeightMap[i + 1].height_[j - 1] &&
101  thePulseHeightMap[i].height_[j] > thePulseHeightMap[i + 1].height_[j + 1]) {
102  tmpCluster.localMax.clear();
103  localMaxTMP.Strip = i;
104  localMaxTMP.Time = j;
105  tmpCluster.localMax.push_back(localMaxTMP);
106  tmpCluster.LayerId = this->LId;
107  tmpCluster.LFTBNDTime = -100;
108  tmpCluster.LFTBNDStrip = -100;
109  tmpCluster.IRTBNDTime = -100;
110  tmpCluster.IRTBNDStrip = -100;
111  MEStripClusters.push_back(tmpCluster);
112  }
113  }
114  }
115  return;
116  }
117  void StripClusterFinder::SearchBorders(void) {
118  uint32_t iS, iT, iL, jL, iR, jR;
119 
120  // SEARCHING PARAMETERS OF THE CLASTERS (LEFT DOWN & RIGHT UP)
121 
122  for (i = 0; i < MEStripClusters.size(); i++) {
123  if (MEStripClusters[i].localMax.empty()) {
124  std::cout << "!!!Warning Cluster has'nt local Maxima" << std::endl;
125  continue;
126  }
127  iS = MEStripClusters[i].localMax[0].Strip;
128  iT = MEStripClusters[i].localMax[0].Time;
129  // LEFT DOWN
130  // strip
131  MEStripClusters[i].LFTBNDStrip = 0;
132  for (iL = iS - 1; iL > 0; iL--) {
133  if (thePulseHeightMap[iL].channel_ == 64) {
134  MEStripClusters[i].LFTBNDStrip = iL;
135  break;
136  }
137  if (thePulseHeightMap[iL].height_[iT] == 0.) {
138  MEStripClusters[i].LFTBNDStrip = iL + 1;
139  break;
140  }
141  }
142  //time
143  MEStripClusters[i].LFTBNDTime = 0;
144  for (jL = iT - 1; jL > 0; jL--) {
145  if (thePulseHeightMap[iS].height_[jL] == 0.) {
146  MEStripClusters[i].LFTBNDTime = jL + 1;
147  break;
148  }
149  }
150  // RIGHT UP
151  //strip
152  MEStripClusters[i].IRTBNDStrip = thePulseHeightMap.size() - 1;
153  for (iR = iS + 1; iR < thePulseHeightMap.size(); iR++) {
154  if (thePulseHeightMap[iR].channel_ == 63) {
155  MEStripClusters[i].IRTBNDStrip = iR;
156  break;
157  }
158  if (thePulseHeightMap[iR].height_[iT] == 0.) {
159  MEStripClusters[i].IRTBNDStrip = iR - 1;
160  break;
161  }
162  }
163  //time
164  MEStripClusters[i].IRTBNDTime = 15;
165  for (jR = iT + 1; jR < 16; jR++) {
166  if (thePulseHeightMap[iS].height_[jR] == 0.) {
167  MEStripClusters[i].IRTBNDTime = jR - 1;
168  break;
169  }
170  }
171  }
172  return;
173  }
174 
175  void StripClusterFinder::Match(void) {
176  // MATCHING THE OVERLAPING CLASTERS
177  bool find2match;
178  find2match = true;
179  do {
180  find2match = FindAndMatch();
181  } while (find2match);
182 
183  return;
184  }
185 
186  bool StripClusterFinder::FindAndMatch(void) {
187  // Find clusters to match
188  icstart = 0;
189  for (ic1 = icstart; ic1 < MEStripClusters.size(); ic1++) {
190  IC1MIN = MEStripClusters[ic1].LFTBNDStrip;
191  IC1MAX = MEStripClusters[ic1].IRTBNDStrip;
192  JC1MIN = MEStripClusters[ic1].LFTBNDTime;
193  JC1MAX = MEStripClusters[ic1].IRTBNDTime;
194  for (ic2 = ic1 + 1; ic2 < MEStripClusters.size(); ic2++) {
195  IC2MIN = MEStripClusters[ic2].LFTBNDStrip;
196  IC2MAX = MEStripClusters[ic2].IRTBNDStrip;
197  JC2MIN = MEStripClusters[ic2].LFTBNDTime;
198  JC2MAX = MEStripClusters[ic2].IRTBNDTime;
199  if ((IC2MIN >= IC1MIN && IC2MIN <= IC1MAX && JC2MIN >= JC1MIN && JC2MIN <= JC1MAX) ||
200  (IC2MIN >= IC1MIN && IC2MIN <= IC1MAX && JC2MAX >= JC1MIN && JC2MAX <= JC1MAX) ||
201  (IC2MAX >= IC1MIN && IC2MAX <= IC1MAX && JC2MIN >= JC1MIN && JC2MIN <= JC1MAX) ||
202  (IC2MAX >= IC1MIN && IC2MAX <= IC1MAX && JC2MAX >= JC1MIN && JC2MAX <= JC1MAX)) {
203  KillCluster();
204  return true;
205  } else {
206  if ((IC1MIN >= IC2MIN && IC1MIN <= IC2MAX && JC1MIN >= JC2MIN && JC1MIN <= JC2MAX) ||
207  (IC1MIN >= IC2MIN && IC1MIN <= IC2MAX && JC1MAX >= JC2MIN && JC1MAX <= JC2MAX) ||
208  (IC1MAX >= IC2MIN && IC1MAX <= IC2MAX && JC1MIN >= JC2MIN && JC1MIN <= JC2MAX) ||
209  (IC1MAX >= IC2MIN && IC1MAX <= IC2MAX && JC1MAX >= JC2MIN && JC1MAX <= JC2MAX)) {
210  KillCluster();
211  return true;
212  }
213  }
214  }
215  }
216  return false;
217  }
218  void StripClusterFinder::KillCluster(void) {
219  // Match Clusters and kill one of clusters.
220  if (IC1MIN < IC2MIN)
221  MEStripClusters[ic1].LFTBNDStrip = IC1MIN;
222  else
223  MEStripClusters[ic1].LFTBNDStrip = IC2MIN;
224  if (JC1MIN < JC2MIN)
225  MEStripClusters[ic1].LFTBNDTime = JC1MIN;
226  else
227  MEStripClusters[ic1].LFTBNDTime = JC2MIN;
228  if (IC1MAX > IC2MAX)
229  MEStripClusters[ic1].IRTBNDStrip = IC1MAX;
230  else
231  MEStripClusters[ic1].IRTBNDStrip = IC2MAX;
232  if (JC1MAX > JC2MAX)
233  MEStripClusters[ic1].IRTBNDTime = JC1MAX;
234  else
235  MEStripClusters[ic1].IRTBNDTime = JC2MAX;
236 
237  MEStripClusters.erase(MEStripClusters.begin() + ic2);
238  icstart = ic1;
239 
240  return;
241  }
242  void StripClusterFinder::RefindMax(void) {
243  int iLS, iRS, iLT, iRT;
244  int iS, jT;
245  int ilocal;
246  float GlobalMax;
247  bool Erased;
248  // SEARCHING EXTREMUMS IN THE CLASTERS
249 
250  for (i = 0; i < MEStripClusters.size(); i++) {
251  MEStripClusters[i].localMax.clear();
252  ilocal = 0;
253  iLS = MEStripClusters[i].LFTBNDStrip;
254  iRS = MEStripClusters[i].IRTBNDStrip;
255  iLT = MEStripClusters[i].LFTBNDTime;
256  iRT = MEStripClusters[i].IRTBNDTime;
257 
258  /*
259  for(iS=iLS+1;iS<=iRS-1;iS++){
260  for(jT=iLT+1;jT<=iRT-1;jT++){
261  */
262  for (iS = iLS; iS <= iRS; iS++) {
263  if (thePulseHeightMap[iS].channel_ == 63 || thePulseHeightMap[iS].channel_ == 64)
264  continue;
265  for (jT = iLT; jT <= iRT; jT++) {
266  if (iS == 0 || jT == 0 || iS == 79 || jT == 7)
267  continue;
268  if (thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS - 1].height_[jT] &&
269  thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS + 1].height_[jT] &&
270  thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS].height_[jT - 1] &&
271  thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS].height_[jT + 1] &&
272  thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS - 1].height_[jT - 1] &&
273  thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS - 1].height_[jT + 1] &&
274  thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS + 1].height_[jT - 1] &&
275  thePulseHeightMap[iS].height_[jT] > thePulseHeightMap[iS + 1].height_[jT + 1]) {
276  localMaxTMP.Strip = iS;
277  localMaxTMP.Time = jT;
278  MEStripClusters[i].localMax.push_back(localMaxTMP);
279  ilocal++;
280  }
281  }
282  }
283  // kill local maximums rellated to noise, maximums with pulse height less then 10% of Global max of clust.
284  //fing Global Max
285  GlobalMax = 0;
286  if (!MEStripClusters[i].localMax.empty()) {
287  //std::cout << "Cluster: " << i << " Number of local maximums before erase: "
288  // << MEStripClusters[i].localMax.size() << std::endl;
289  for (j = 0; j < MEStripClusters[i].localMax.size(); j++) {
290  iS = MEStripClusters[i].localMax[j].Strip;
291  jT = MEStripClusters[i].localMax[j].Time;
292  /*
293  std::cout << "Current Max:"
294  << " " << iS
295  << " " << jT
296  << " " << thePulseHeightMap[iS].height_[jT] << std::endl;
297  */
298  if (thePulseHeightMap[iS].height_[jT] > GlobalMax)
299  GlobalMax = thePulseHeightMap[iS].height_[jT];
300  }
301  GlobalMax = (float)(GlobalMax / 10.);
302  //erase noise localMaximums
303  do {
304  Erased = false;
305  for (j = 0; j < MEStripClusters[i].localMax.size(); j++) {
306  iS = MEStripClusters[i].localMax[j].Strip;
307  jT = MEStripClusters[i].localMax[j].Time;
308  if (thePulseHeightMap[iS].height_[jT] < GlobalMax) {
309  MEStripClusters[i].localMax.erase(MEStripClusters[i].localMax.begin() + j);
310  Erased = true;
311  break;
312  }
313  }
314  } while (Erased);
315 
316  //debug outs
317  //std::cout << "Cluster: " << i << " Number of local maximums: "
318  // << MEStripClusters[i].localMax.size() << std::endl;
319  /*
320  for(j=0;j<MEStripClusters[i].localMax.size();j++){
321  iS=MEStripClusters[i].localMax[j].Strip;
322  jT=MEStripClusters[i].localMax[j].Time;
323  std::cout << "Local Max: " << j << " Strip: " << iS << " Time: " << jT
324  << " Height: " << thePulseHeightMap[iS].height_[jT]
325  << " Cut Value: " << GlobalMax << std::endl;
326  }
327  */
328  }
329  }
330  return;
331  }
332  void StripClusterFinder::printClusters(void) {
333  int iS, jT;
334  std::cout << "====================================================================" << std::endl;
335  std::cout << "debug information from StripClusterFinder" << std::endl;
336  for (i = 0; i < MEStripClusters.size(); i++) {
337  if (MEStripClusters[i].localMax.empty())
338  continue;
339  std::cout << " Cluster: " << i + 1 << " Number of local Maximums " << MEStripClusters[i].localMax.size()
340  << std::endl;
341  for (j = 0; j < MEStripClusters[i].localMax.size(); j++) {
342  iS = MEStripClusters[i].localMax[j].Strip;
343  jT = MEStripClusters[i].localMax[j].Time;
344 
345  // std::cout << "Local Max: " << j << " Strip: " << iS << " Time: " << jT << std::endl;
346  for (uint32_t k = 0; k < MEStripClusters[i].ClusterPulseMapHeight.size(); k++) {
347  if (MEStripClusters[i].ClusterPulseMapHeight[k].channel_ == iS)
348  std::cout << "Local Max: " << j + 1 << " Strip: " << iS + 1 << " Time: " << jT + 1
349  << " Height: " << MEStripClusters[i].ClusterPulseMapHeight[k].height_[jT] << std::endl;
350  }
351  }
352  for (uint32_t k = 0; k < MEStripClusters[i].ClusterPulseMapHeight.size(); k++) {
353  std::cout << "Strip: " << MEStripClusters[i].ClusterPulseMapHeight[k].channel_ + 1;
354  for (int l = 0; l < 16; l++)
355  std::cout << " " << MEStripClusters[i].ClusterPulseMapHeight[k].height_[l];
356  std::cout << std::endl;
357  }
358 
359  std::cout << " Left Top corner strip: " << MEStripClusters[i].LFTBNDStrip + 1 << " "
360  << " time: " << MEStripClusters[i].LFTBNDTime + 1 << std::endl;
361  std::cout << " Right Bottom corner strip: " << MEStripClusters[i].IRTBNDStrip + 1 << " "
362  << " time: " << MEStripClusters[i].IRTBNDTime + 1 << std::endl;
363  }
364  std::cout << "======================================================================" << std::endl;
365  return;
366  }
367  bool StripClusterFinder::Sort::operator()(const StripClusterFitData& a, const StripClusterFitData& b) const {
368  return a.channel_ < b.channel_;
369  }
370 
371 } // namespace cscdqm
cscdqm::StripCluster::localMax
std::vector< ClusterLocalMax > localMax
Definition: CSCDQM_StripCluster.h:20
cscdqm::StripClusterFitData
Strip Cluster Fit Data Object.
Definition: CSCDQM_StripClusterFitData.h:12
CSCDQM_StripClusterFinder.h
cscdqm::StripClusterFinder::StripClusterFinder
StripClusterFinder(int l, int s, int cf, int st)
Definition: CSCDQM_StripClusterFinder.cc:6
cscdqm::StripClusterFinder::Sort::operator()
bool operator()(const StripClusterFitData &a, const StripClusterFitData &b) const
Definition: CSCDQM_StripClusterFinder.cc:367
cscdqm::StripClusterFinder::thePulseHeightMap
std::vector< StripClusterFitData > thePulseHeightMap
Definition: CSCDQM_StripClusterFinder.h:28
cscdqm::StripClusterFinder::DoAction
void DoAction(int layerId, float *cathodes)
Definition: CSCDQM_StripClusterFinder.cc:15
cscdqm::StripClusterFinder::MEStripClusters
std::vector< StripCluster > MEStripClusters
Definition: CSCDQM_StripClusterFinder.h:35
b
double b
Definition: hdecay.h:118
a
double a
Definition: hdecay.h:119
gather_cfg.cout
cout
Definition: gather_cfg.py:144
Generated for CMSSW Reference Manual by doxygen 1.8.11