Belle II Software  release-05-01-25
TRGCDCHoughMapping.cc
1 #define TRG_SHORT_NAMES
2 
3 #include <iostream>
4 #include <fstream>
5 #include <vector>
6 #include <math.h>
7 #include <iomanip>
8 #include "trg/trg/Utilities.h"
9 
10 using namespace std;
11 using namespace Belle2;
12 
13 #define NAME "HoughMapping"
14 #define VERSION "version 0.01"
15 
17 struct XY {
19  double x;
21  double y;
22 };
24 struct Plane {
26  unsigned x;
28  unsigned y;
29 };
30 
31 void printHeader(ofstream& out, const string& function);
32 void superLayer(const unsigned id);
33 void printMapping(ofstream& out, vector<vector<vector<int>>>& HPcell);
34 
35 //Hough Plane Parameter
36 //Phi Range (X Axial)
37 const double PI2 = 2 * M_PI;
38 //log10(r/cm) Range (Y Axial)
39 const double minY = 1.823908740944321;
40 const double maxY = 3.204119982655926;
41 //cell number
42 const unsigned nX = 160;
43 const unsigned nY = 16;
44 
45 vector<vector<vector<int>>> HPcellM(nY, vector<vector<int>>(nX, vector<int>()));
46 ofstream outputM("minus_total.dat");
47 
48 vector<vector<vector<int>>> HPcellP(nY, vector<vector<int>>(nX, vector<int>()));
49 ofstream outputP("plus_total.dat");
50 
51 
52 //...C++ for TSIM...
53 const string fncM = "HoughMappingMinus.cc";
54 const string fncP = "HoughMappingPlus.cc";
55 ofstream outcM(fncM);
56 ofstream outcP(fncP);
57 
58 int
59 main(int, char**)
60 {
61  cout << NAME << " ... " << VERSION << endl;
62 
63  printHeader(outcM, "HoughMappingMinus");
64  printHeader(outcP, "HoughMappingPlus");
65 
66  for (unsigned isl = 0; isl < 5 ; isl++)
67  superLayer(isl);
68 
69  printMapping(outputM, HPcellM);
70  printMapping(outputP, HPcellP);
71 
72  outcM << endl << "}" << endl;
73  outcP << endl << "}" << endl;
74 
75  outcM.close();
76  outcP.close();
77 
78  //...Termination...
79  cout << "Files generated" << endl;
80  cout << " c++ for tsim firmware : " << fncM << endl;
81  cout << " c++ for tsim firmware : " << fncP << endl;
82 }
83 
84 void
85 printHeader(ofstream& out, const string& function)
86 {
87  //...Date...
88  string ts = TRGUtil::dateString();
89 
90  out << "// This file is generated by " << NAME << "(" << VERSION << ")" << endl;
91  out << "// " << ts << endl << endl;
92  out << "#define TRGCDC_SHORT_NAMES" << endl;
93  out << "#include \"trg/trg/State.h\"" << endl;
94  out << "#include \"trg/cdc/Tracker2D.h\"" << endl;
95  out << "using namespace std;" << endl;
96  out << "using namespace Belle2;" << endl;
97  out << "void" << endl;
98  out << "TCTracker2D::" << function << "(void) {" << endl;
99 
100  out << " //...TS hit map..." << endl;
101  out << " TRGState SL0_TS = _ts.subset(0, 160);" << endl;
102  out << " TRGState SL2_TS = _ts.subset(160, 192);" << endl;
103  out << " TRGState SL4_TS = _ts.subset(160 + 192, 256);" << endl;
104  out << " TRGState SL6_TS = _ts.subset(160 + 192 + 256, 320);" << endl;
105  out << " TRGState SL8_TS = _ts.subset(160 + 192 + 256 + 320, 384);" << endl;
106 
107  out << " //...Hough cells..." << endl;
108 
109  for (unsigned isl = 0; isl < 5; isl++) {
110  for (unsigned iy = 1; iy < nY + 1; iy++) {
111  out << " TRGState SL" << to_string(isl * 2) << "_row"
112  << to_string(iy) << "(" << to_string(nX) << ");" << endl;
113  }
114  }
115 }
116 
117 void
118 superLayer(const unsigned id)
119 {
120  //Radius of SL_center cell
121  double r_SL = 0;
122  //Number of TS each SL
123  int N_TS_SL = 0;
124  int SL = 2 * id;
125 
126  if (id == 0) {
127  r_SL = 19.8;
128  N_TS_SL = 160;
129  } else if (id == 1) {
130  r_SL = 40.16;
131  N_TS_SL = 192;
132  } else if (id == 2) {
133  r_SL = 62.0;
134  N_TS_SL = 256;
135  } else if (id == 3) {
136  r_SL = 83.84;
137  N_TS_SL = 320;
138  } else if (id == 4) {
139  r_SL = 105.68;
140  N_TS_SL = 384;
141  } else {
142  cout << NAME << " !!! bad super layer ID" << endl;
143  exit(-1);
144  }
145 
146  //Hough Plane database
147  vector<XY> xymatrix;
148  XY xy = {0, 0};
149 
151 
152  for (int i = 0; i < N_TS_SL; ++i) {
153  xy.x = r_SL * cos((PI2 / N_TS_SL) * i);
154  xy.y = r_SL * sin((PI2 / N_TS_SL) * i);
155  xymatrix.push_back(xy);
156  }
158 
159  const double r0 = minY;
160  const double phi0 = 0;
161  const double dr = (maxY - minY) / nY;
162  const double dphi = PI2 / nX;
163  double r1;
164  double r2;
165  double phi1;
166  double phi2;
167  double minus1;
168  double minus2;
169  double plus1;
170  double plus2;
171 
172  const string vhM = "UT3_0_SL" + to_string(id * 2) + ".vhd";
173  const string vhP = "UT3_0_SL" + to_string(id * 2) + "_p.vhd";
174  ofstream outputfM(vhM);
175  ofstream outputfP(vhP);
176 
177  //generate firware code(Minus)
178  outputfM << "library IEEE;" << endl;
179  outputfM << "use IEEE.STD_LOGIC_1164.ALL;" << endl;
180  outputfM << " " << endl;
181  outputfM << " " << endl;
182  outputfM << "entity UT3_0_SL" << SL << " is" << endl;
183  outputfM << " " << endl;
184  outputfM << "Port (" << endl;
185  for (unsigned irow = 1; irow < nY + 1; ++irow) {
186  outputfM << " SL" << SL << "_row" << left << setw(2) << to_string(irow)
187  << " : out STD_LOGIC_VECTOR (79 downto 40);" << endl;
188  }
189  outputfM << " SL" << SL << "_TS : in STD_LOGIC_VECTOR ("
190  << N_TS_SL / 2 << " downto 0));" << endl;
191  outputfM << "end UT3_0_SL" << SL << ";" << endl;
192  outputfM << " " << endl;
193  outputfM << " " << endl;
194  outputfM << "architecture Behavioral of UT3_0_SL" << SL << " is" << endl;
195  outputfM << " " << endl;
196  outputfM << "begin" << endl;
197  outputfM << " " << endl;
198 
199  //generate firmware code(Plus)
200  outputfP << "library IEEE;" << endl;
201  outputfP << "use IEEE.STD_LOGIC_1164.ALL;" << endl;
202  outputfP << " " << endl;
203  outputfP << " " << endl;
204  outputfP << "entity UT3_0_SL" << SL << "_P is" << endl;
205  outputfP << " " << endl;
206  outputfP << "Port (" << endl;
207  for (unsigned irow = 1; irow < nY + 1; ++irow) {
208  outputfP << " SL" << SL << "_row" << left << setw(2) << to_string(irow)
209  << " : out STD_LOGIC_VECTOR (39 downto 0);" << endl;
210  }
211  outputfP << " SL" << SL << "_TS : in STD_LOGIC_VECTOR ("
212  << N_TS_SL / 2 << " downto 0));" << endl;
213  outputfP << "end UT3_0_SL" << SL << "_P;" << endl;
214  outputfP << " " << endl;
215  outputfP << " " << endl;
216  outputfP << "architecture Behavioral of UT3_0_SL" << SL << "_P is" << endl;
217  outputfP << " " << endl;
218  outputfP << "begin" << endl;
219  outputfP << " " << endl;
220 
221  //vertical
222  for (unsigned k = 0 ; k < nY ; ++k) {
223  //horizontal
224  for (unsigned j = 0 ; j < nX ; ++j) {
225  double ffM = 0;
226  double ffP = 0;
227 
228  if (j > 39 && j < 80)
229  outputfM << "SL" << SL << "_row" << k + 1 << "(" << j << ")<=";
230  if (j < 40)
231  outputfP << "SL" << SL << "_row" << k + 1 << "(" << j << ")<=";
232  outcM << " SL" << SL << "_row" << k + 1 << ".set(" << j << ", ";
233  outcP << " SL" << SL << "_row" << k + 1 << ".set(" << j << ", ";
234  bool firstM = true;
235  bool firstP = true;
236 
237  phi1 = phi0 + j * dphi;
238  phi2 = phi0 + (j + 1) * dphi;
239 
240  //TS
241  for (int i = 0 ; i < N_TS_SL ; i++) {
242  // calculate r(phi) at Hough cell borders phi1 and phi2
243  r1 = ((xymatrix[i].x * xymatrix[i].x) + (xymatrix[i].y * xymatrix[i].y)) /
244  ((2 * xymatrix[i].x * cos(phi1)) + (2 * xymatrix[i].y * sin(phi1)));
245  r2 = ((xymatrix[i].x * xymatrix[i].x) + (xymatrix[i].y * xymatrix[i].y)) /
246  ((2 * xymatrix[i].x * cos(phi2)) + (2 * xymatrix[i].y * sin(phi2)));
247 
248  /* Check whether f(phi) = log(r(phi)) crosses the Hough cell
249  * defined by (phi1, phi2, log(r1), log(r2))
250  * The slope determines the charge of the track.
251  *
252  * Since f is not defined for all phi, 3 cases can occur:
253  * 1. f(phi1) and f(phi2) both defined:
254  * compare f(phi1) and f(phi2) to log(r1) and log(r2),
255  * get slope from f(phi2) - f(phi1)
256  * 2. f(phi1) defined, f(phi2) not defined (or vice-versa):
257  * compare f(phi1) to log(r1) and log(r2),
258  * slope is known
259  * 3. f(phi1) and f(phi2) both not defined: no entry
260  */
261 
262  //...minus...
263  if (r1 >= 0 && r2 >= 0 && r1 < r2) {
264  /* positive slope:
265  * crossing if f(phi1) < log(r2) and f(phi2) > log(r1)
266  */
267  minus1 = r0 + (k + 1) * dr - log10(r1);
268  minus2 = r0 + k * dr - log10(r2);
269  if (minus1 * minus2 <= 0.0) {
270  HPcellM[k][j].push_back(SL);
271  HPcellM[k][j].push_back(i);
272 
273  if (! firstM)
274  outcM << " or ";
275 
276  if (ffM != 0) {
277  if (j > 39 && j < 80)
278  outputfM << "or ";
279  }
280 
281  ffM++;
282  if (j > 39 && j < 80)
283  outputfM << "SL" << SL << "_TS(" << i << ") ";
284 
285  outcM << "SL" << SL << "_TS[" << i << "]";
286  firstM = false;
287  }
288  } else if (r2 < 0 && r1 >= 0) {
289  /* positive slope, f(phi2) = inf:
290  * crossing if f(phi1) < log(r2)
291  */
292  minus1 = r0 + (k + 1) * dr - log10(r1);
293  if (minus1 > 0) {
294  HPcellM[k][j].push_back(SL);
295  HPcellM[k][j].push_back(i);
296 
297  if (! firstM)
298  outcM << " or ";
299 
300  if (ffM != 0) {
301  if (j > 39 && j < 80)
302  outputfM << "or ";
303  }
304  ffM++;
305  if (j > 39 && j < 80)
306  outputfM << "SL" << SL << "_TS(" << i << ") ";
307 
308  outcM << "SL" << SL << "_TS[" << i << "]";
309  firstM = false;
310  }
311  }
312 
313  //plus
314  if (r1 >= 0 && r2 >= 0 && r2 < r1) {
315  /* negative slope:
316  * crossing if f(phi2) < log(r2) and f(phi1) > log(r1)
317  */
318  plus1 = r0 + (k + 1) * dr - log10(r2);
319  plus2 = r0 + k * dr - log10(r1);
320  if (plus1 * plus2 <= 0.0) {
321  HPcellP[k][j].push_back(SL);
322  HPcellP[k][j].push_back(i);
323  if (! firstP)
324  outcP << " or ";
325 
326  if (ffP != 0) {
327  if (j < 40)
328  outputfP << "or ";
329  }
330  ffP++;
331  if (j < 40)
332  outputfP << "SL" << SL << "_TS(" << i << ") ";
333 
334  outcP << "SL" << SL << "_TS[" << i << "]";
335  firstP = false;
336  }
337  } else if (r1 < 0 && r2 >= 0) {
338  /* negative slope, f(phi1) = inf:
339  * crossing if f(phi2) < log(r2)
340  */
341  plus1 = r0 + (k + 1) * dr - log10(r2);
342  if (plus1 > 0) {
343  HPcellP[k][j].push_back(SL);
344  HPcellP[k][j].push_back(i);
345  if (! firstP)
346  outcP << " or ";
347 
348  if (ffP != 0) {
349  if (j < 40)
350  outputfP << "or ";
351  }
352  ffP++;
353  if (j < 40)
354  outputfP << "SL" << SL << "_TS(" << i << ") ";
355 
356  outcP << "SL" << SL << "_TS[" << i << "]";
357  firstP = false;
358  }
359  }
360  }
361 
362  if (j > 39 && j < 80)
363  outputfM << ";" << endl;
364  if (j < 40)
365  outputfP << ";" << endl;
366  outcM << ");" << endl;
367  outcP << ");" << endl;
368  }
369  }
370  outputfM << " " << endl;
371  outputfM << "end Behavioral;" << endl;
372  outputfP << " " << endl;
373  outputfP << "end Behavioral;" << endl;
374  return;
375 }
376 
377 void
378 printMapping(ofstream& out, vector<vector<vector<int>>>& HPcell)
379 {
380  //...Date...
381  string ts = TRGUtil::dateString();
382  out << "// This file is generated by " << NAME << "(" << VERSION << ")" << endl;
383  out << "// " << ts << endl << endl;
384  //...Hough plane parameters...
385  out << nX << " " << nY << " " << minY << " " << maxY << endl;
386  //...TS patterns...
387  for (unsigned iy = 0; iy < nY; iy++) {
388  for (unsigned ix = 0; ix < nX; ix++) {
389  out << ix << " " << iy + 1 << " ";
390  for (unsigned its = 0; its < HPcell[iy][ix].size(); its++) {
391  out << HPcell[iy][ix][its] << " ";
392  }
393  out << " " << endl;
394  }
395  }
396 }
XY::y
double y
y of the point
Definition: TRGCDCHoughMapping.cc:21
XY::x
double x
x of the point
Definition: TRGCDCHoughMapping.cc:19
Plane::x
unsigned x
x of the plane
Definition: TRGCDCHoughMapping.cc:26
main
int main(int argc, char **argv)
Run all tests.
Definition: test_main.cc:77
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Plane::y
unsigned y
y of the point
Definition: TRGCDCHoughMapping.cc:28
Plane
Plane.
Definition: TRGCDCHoughMapping.cc:24
XY
Point.
Definition: TRGCDCHoughMapping.cc:17