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