Belle II Software  release-05-01-25
HoughPlaneBoolean.cc
1 //-----------------------------------------------------------------------------
2 // $Id$
3 //-----------------------------------------------------------------------------
4 // Filename : HoughPlaneBoolean.cc
5 // Section : TRG CDC
6 // Owner : Yoshihito Iwasaki
7 // Email : yoshihito.iwasaki@kek.jp
8 //-----------------------------------------------------------------------------
9 // Description : A class to represent a Hough parameter plane with simple
10 // counter.
11 //-----------------------------------------------------------------------------
12 // $Log$
13 //-----------------------------------------------------------------------------
14 
15 #define TRGCDC_SHORT_NAMES
16 
17 #include <iostream>
18 #include "trg/cdc/HoughPlaneBoolean.h"
19 
20 using namespace std;
21 
22 namespace Belle2 {
28  TRGCDCHoughPlaneBoolean::TRGCDCHoughPlaneBoolean(const std::string& name,
29  const TCHTransformation& trans,
30  unsigned nX,
31  float xMin,
32  float xMax,
33  unsigned nY,
34  float yMin,
35  float yMax)
36  : TRGCDCHoughPlaneBase(name, trans, nX, xMin, xMax, nY, yMin, yMax),
37  _n(nX * nY / 32 + 1),
38  _cell(new unsigned[_n]),
39  _nPatterns(0),
40  _patterns(0),
41  _nActive(0),
42  _reverse(0)
43  {
44  clear();
45  }
46 
48  {
49  for (unsigned i = 0; i < _nPatterns; i++)
50  delete [] _patterns[i];
51  delete [] _patterns;
52  delete [] _nActive;
53  delete [] _cell;
54  delete [] _reverse;
55  }
56 
57  void
59  float ry,
60  float targetCharge,
61  int weight)
62  {
63 
64  const HepGeom::Point3D<double> r(rx, ry, 0);
65 
66 //
67 // cout << "yMax=" << yMax() << endl;
68 //
69 
70  //...phi loop...
71  for (unsigned i = 0; i < nX(); i++) {
72  const float x0 = xSize() * float(i);
73  const HepGeom::Point3D<double> center(cos(x0), sin(x0), 0);
74  float charge = r.cross(center).z();
75  if (targetCharge != 0)
76  if (targetCharge * charge > 0)
77  continue;
78 
79  float y0 = transformation().y(rx, ry, x0);
80  const float x1 = xSize() * float(i + 1);
81  float y1 = transformation().y(rx, ry, x1);
82 
83 //
84 // cout << "x0,x1,y0,y1=" << x0 << "," << x1 << "," << y0 << "," << y1
85 // << endl;
86 //
87 
88  //...Check y position...
89  if ((y0 == 0) && (y1 == 0))
90  continue;
91  else if ((y0 > yMax()) && (y1 > yMax()))
92  continue;
93  else if ((y0 < yMin()) && (y1 < yMin()))
94  continue;
95 
96  //...Divergence here...
97  if ((y0 == 0) && (y1 != 0))
98  y0 = yMax();
99  else if ((y0 != 0) && (y1 == 0))
100  y1 = yMax();
101 
102  //...Adjust location...
103  if (y0 < yMin())
104  y0 = yMin();
105  if (y1 < yMin())
106  y1 = yMin();
107  if (y0 > yMax())
108  y0 = yMax();
109  if (y1 > yMax())
110  y1 = yMax();
111 
112  //...Location in the plane...
113  int iY0 = int((y0 - yMin()) / ySize());
114  int iY1 = int((y1 - yMin()) / ySize());
115 
116 //
117 // cout << "x0,x1,y0,y1=" << x0 << "," << x1 << "," << y0 << "," << y1
118 // << "," << iY0 << "," << iY1 << endl;
119 //
120 
121  //...Sorting...
122  if (iY0 > iY1) {
123  const int tmp = iY0;
124  iY0 = iY1;
125  iY1 = tmp;
126  }
127 
128  //...Voting...
129  for (unsigned j = (unsigned) iY0; j < (unsigned)(iY1 + 1); j++) {
130  const unsigned k = i * nY() + j;
131  const unsigned b0 = k / 32;
132  const unsigned b1 = k % 32;
133  if (weight > 0)
134  _cell[b0] |= (1 << b1);
135  else
136  _cell[b0] &= (~(1 << b1));
137  }
138  }
139  }
140 
141  void
143  float ry,
144  float targetCharge,
145  int weight)
146  {
147 
148  const HepGeom::Point3D<double> r(rx, ry, 0);
149 
150  //...phi loop...
151  for (unsigned i = 0; i < nX(); i++) {
152  const float x0 = xSize() * float(i);
153  const HepGeom::Point3D<double> phi(cos(x0), sin(x0), 0);
154  float charge = r.cross(phi).z();
155  if (targetCharge != 0)
156  if (targetCharge * charge > 0)
157  continue;
158 
159  const float y0 = transformation().y(rx, ry, x0);
160  const float x1 = xSize() * float(i + 1);
161  const float y1 = transformation().y(rx, ry, x1);
162 
163  //...Location in the plane...
164  int iY0 = int((y0 - yMin()) / ySize());
165  int iY1 = int((y1 - yMin()) / ySize());
166 
167  //...This is special implementation for Circle Hough...
168  if (transformation().diverge(rx, ry, x0, x1)) {
169  if (iY0 > 0) {
170  if (iY0 >= (int) nY()) continue;
171  iY1 = nY() - 1;
172  } else {
173  if (iY1 >= (int) nY()) continue;
174  iY0 = iY1;
175  iY1 = nY() - 1;
176  }
177  }
178 
179  //...Sorting...
180  if (iY0 > iY1) {
181  const int tmp = iY0;
182  iY0 = iY1;
183  iY1 = tmp;
184  }
185 
186  //...Both out of region ?...
187  if (iY1 < 0) continue;
188  if (iY0 >= (int) nY()) continue;
189 
190  //...In region ?...
191  if (iY0 < 0) iY0 = 0;
192  if (iY0 >= (int) nY()) iY0 = nY() - 1;
193  if (iY1 < 0) iY1 = 0;
194  if (iY1 >= (int) nY()) iY1 = nY() - 1;
195 
196  //...Voting...
197  for (unsigned j = (unsigned) iY0; j < (unsigned)(iY1 + 1); j++) {
198  const unsigned k = i * nY() + j;
199  const unsigned b0 = k / 32;
200  const unsigned b1 = k % 32;
201  if (weight > 0)
202  _cell[b0] |= (1 << b1);
203  else
204  _cell[b0] &= (~(1 << b1));
205  }
206  }
207  }
208 
209  void
211  {
212 
213  //...Check status...
214  if (_patterns[id]) {
215  std::cout << "TRGCDCHoughPlaneBoolean::registerPattern !!! "
216  << "a pattern(id=" << id << ") was already registered"
217  << std::endl;
218  return;
219  }
220  if (id >= _nPatterns) {
221  std::cout << "TRGCDCHoughPlaneBoolean::registerPattern !!! "
222  << "given id(" << id << ") is over range(" << _nPatterns << ")"
223  << std::endl;
224  return;
225  }
226  const unsigned n = nX() * nY();
227 
228  //...Check # of active cells...
229  unsigned nActive = 0;
230  for (unsigned i = 0; i < n; i++) {
231  const unsigned j = i / 32;
232  const unsigned k = i % 32;
233  if ((_cell[j] >> k) & 1) {
234  ++nActive;
235  _reverse[i].push_back(id);
236  }
237  }
238  _nActive[id] = nActive;
239 
240  //...Create an array...
241  _patterns[id] = new unsigned[nActive];
242 
243  //...Store them...
244  unsigned l = 0;
245  for (unsigned i = 0; i < n; i++) {
246  const unsigned j = i / 32;
247  const unsigned k = i % 32;
248  if ((_cell[j] >> k) & 1) {
249  _patterns[id][l] = i;
250  l++;
251  }
252  }
253  }
254 
255  void
257  {
258  if (_nPatterns) {
259  std::cout << "TRGCDCHoughPlaneBoolean::preparePatterns !!! "
260  << "patterns were already allocated" << std::endl;
261  return;
262  }
263 
264  _nPatterns = nPatterns;
265  _patterns = new unsigned * [_nPatterns];
266  _nActive = new unsigned[_nPatterns];
267  memset(_patterns, 0, _nPatterns * sizeof(unsigned*));
268  _reverse = new vector<unsigned>[nX() * nY()];
269  }
270 
271  void
272  TRGCDCHoughPlaneBoolean::vote(unsigned patternId, int weight)
273  {
274  const unsigned n = _nActive[patternId];
275  for (unsigned i = 0; i < n; i++) {
276  add(_patterns[patternId][i], weight);
277  }
278  }
279 
281 } // namespace Belle2
Belle2::TRGCDCHoughPlaneBoolean::clear
void clear(void) override
clear all entries.
Definition: HoughPlaneBoolean.h:143
Belle2::TRGCDCHoughPlaneBase::yMax
float yMax(void) const
returns max. of y.
Definition: HoughPlaneBase.h:320
Belle2::TRGCDCHoughPlaneBoolean::vote
void vote(float rx, float ry, int weight=1) override
Votes.
Definition: HoughPlaneBoolean.h:151
Belle2::TRGCDCHoughPlaneBase::id
void id(unsigned serialId, unsigned &x, unsigned &y) const
returns x and y for serialID.
Definition: HoughPlaneBase.h:484
Belle2::TRGCDCHoughPlaneBase::charge
float charge(void) const
returns charge for this plane.
Definition: HoughPlaneBase.h:230
Belle2::TRGCDCHoughPlaneBoolean::_cell
unsigned *const _cell
cell
Definition: HoughPlaneBoolean.h:104
Belle2::TRGCDCHoughPlaneBoolean::add
void add(unsigned cellId, int weight) override
Add to a cell.
Definition: HoughPlaneBoolean.h:196
Belle2::TRGCDCHoughPlaneBoolean::_patterns
unsigned ** _patterns
patterns
Definition: HoughPlaneBoolean.h:108
Belle2::TRGCDCHoughPlaneBase::nY
unsigned nY(void) const
return # of y bins.
Definition: HoughPlaneBase.h:297
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TRGCDCHoughPlaneBoolean::_nActive
unsigned * _nActive
number of active cells
Definition: HoughPlaneBoolean.h:110
Belle2::TRGCDCHoughPlaneBoolean::registerPattern
void registerPattern(unsigned id) override
registers a pattern..
Definition: HoughPlaneBoolean.cc:210
Belle2::TRGCDCHoughPlaneBase::transformation
const TRGCDCHoughTransformation & transformation(void) const
returns Hough transformation object.
Definition: HoughPlaneBase.h:523
Belle2::TRGCDCHoughPlaneBoolean::~TRGCDCHoughPlaneBoolean
virtual ~TRGCDCHoughPlaneBoolean()
Destructor.
Definition: HoughPlaneBoolean.cc:47
Belle2::TRGCDCHoughPlaneBoolean::_nPatterns
unsigned _nPatterns
number of patterns
Definition: HoughPlaneBoolean.h:106
Belle2::TRGCDCHoughPlaneBase::ySize
float ySize(void) const
returns size of y bin.
Definition: HoughPlaneBase.h:336
Belle2::TRGCDCHoughPlaneBoolean::preparePatterns
void preparePatterns(unsigned nPatterns)
allocate memory for patterns.
Definition: HoughPlaneBoolean.cc:256
Belle2::TRGCDCHoughPlaneBoolean::voteUsedInTrasan
void voteUsedInTrasan(float rx, float ry, float charge, int weight=1)
Votes with charge decision.
Definition: HoughPlaneBoolean.cc:142
Belle2::TRGCDCHoughPlaneBoolean::patternId
const std::vector< unsigned > & patternId(unsigned cellId) const
returns pattern ID which activates specified cell.
Definition: HoughPlaneBoolean.h:206
Belle2::TRGCDCHoughTransformation::y
virtual float y(float xReal, float yReal, float x) const =0
returns Y coordinate in a Hough parameter plane.
HepGeom::Point3D< double >
Belle2::TRGCDCHoughPlaneBoolean::_reverse
std::vector< unsigned > * _reverse
Pattern ID's for each cell.
Definition: HoughPlaneBoolean.h:113
Belle2::TRGCDCHoughPlaneBase::nX
unsigned nX(void) const
returns # of x bins.
Definition: HoughPlaneBase.h:251
Belle2::TRGCDCHoughPlaneBase::yMin
float yMin(void) const
returns min. of y.
Definition: HoughPlaneBase.h:304
Belle2::TRGCDCHoughPlaneBase::xSize
float xSize(void) const
returns size of x bin.
Definition: HoughPlaneBase.h:290
Belle2::TRGCDCHoughPlaneBase
A class to represent a Hough parameter plane.
Definition: HoughPlaneBase.h:34