Belle II Software development
HoughFinder.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
9//-----------------------------------------------------------------------------
10// Description : A class to find tracks usning Hough algorithm
11//-----------------------------------------------------------------------------
12
13#define TRG_SHORT_NAMES
14#define TRGCDC_SHORT_NAMES
15
16#define ENV_PATH "BELLE2_LOCAL_DIR"
17
18#include <stdlib.h>
19#include <map>
20#include <iostream>
21#include <fstream>
22#include "cdc/geometry/CDCGeometryPar.h"
23#include "trg/trg/Debug.h"
24#include "trg/trg/Utilities.h"
25#include "trg/trg/Time.h"
26#include "trg/trg/Signal.h"
27#include "trg/trg/Constants.h"
28#include "trg/cdc/TRGCDC.h"
29#include "trg/cdc/Layer.h"
30#include "trg/cdc/Cell.h"
31#include "trg/cdc/Wire.h"
32#include "trg/cdc/WireHit.h"
33#include "trg/cdc/HoughFinder.h"
34#include "trg/cdc/Segment.h"
35#include "trg/cdc/SegmentHit.h"
36#include "trg/cdc/Circle.h"
37#include "trg/cdc/TRGCDCTrack.h"
38#include "trg/cdc/Link.h"
39#include "trg/cdc/Relation.h"
40#include "trg/cdc/Fitter3DUtility.h"
41#include "trg/cdc/Fitter3D.h"
42#include "trg/cdc/Helix.h"
43#include "trg/cdc/JLUT.h"
44#include "trg/cdc/JSignal.h"
45#include "trg/cdc/JSignalData.h"
46#include "trg/cdc/WireHitMC.h"
47#include "trg/cdc/TrackMC.h"
48#include "trg/cdc/FrontEnd.h"
49#include "trg/cdc/Merger.h"
50#include "trg/cdc/LUT.h"
51#include "trg/cdc/EventTime.h"
52
53#ifdef TRGCDC_DISPLAY_HOUGH
54#include "trg/cdc/DisplayRphi.h"
55#include "trg/cdc/DisplayHough.h"
56namespace Belle2_TRGCDC {
57 Belle2::TRGCDCDisplayHough* H0 = 0;
58 Belle2::TRGCDCDisplayHough* H1 = 0;
59}
60#endif
61
62using namespace std;
63#ifdef TRGCDC_DISPLAY_HOUGH
64using namespace Belle2_TRGCDC;
65#endif
66
67namespace Belle2 {
73 string
74 TRGCDCHoughFinder::version(void) const
75 {
76 return string("TRGCDCHoughFinder 5.24");
77 }
78
79 TRGCDCHoughFinder::TRGCDCHoughFinder(const string& name,
80 const TRGCDC& TRGCDC,
81 unsigned peakMin,
82 const string& mappingFilePlus,
83 const string& mappingFileMinus,
84 unsigned doit)
85 : _name(name),
86 _cdc(TRGCDC),
87 _circleH("CircleHough"),
88 _doit(doit),
89 _peakFinder("PeakFinder"),
90 _peakMin(peakMin)
91 {
92
93 _commonData = 0;
94
95 //...Read Hough plane parameters from file
96 const string fMap[2] = {mappingFilePlus, mappingFileMinus};
97 const string fName[2] = {string("circle hough plus"), string("circle hough minus")};
98
99 for (unsigned f = 0; f < 2; f++) {
100 const string fn = fMap[f];
101 ifstream infile(fn.c_str(), ios::in);
102 if (infile.fail()) {
103 B2FATAL("Cannot open Hough mapping file " << fn);
104 return;
105 }
106
107 //...Ignore lines not starting with a digit...
108 string ignores;
109 while (!isdigit(infile.peek())) {
110 getline(infile, ignores);
111 }
112
113 //...Read Hough plane cell number and limits...
114 unsigned nX = 0;
115 unsigned nY = 0;
116 double Ymin = 0.;
117 double Ymax = 0.;
118
119 infile >> nX >> nY >> Ymin >> Ymax;
120 infile.close();
121
122 _plane[f] = new TCHPlaneMulti2(fName[f], _circleH,
123 nX, 0, 2 * M_PI,
124 nY, Ymin, Ymax,
125 5);
126 }
127
128 //...Set charge...
129 _plane[0]->charge(1);
130 _plane[1]->charge(-1);
131
132#ifdef TRGCDC_DISPLAY_HOUGH
133 if (! H0)
134 H0 = new TCDisplayHough("Plus");
135 H0->link(* D);
136 H0->clear();
137 H0->show();
138 H0->move(630, 0);
139 if (! H1)
140 H1 = new TCDisplayHough("Minus");
141 H1->link(* D);
142 H1->clear();
143 H1->show();
144 H0->move(1260, 0);
145#endif
146
147 //...Old mapping (trasan methode)...
148 if (_doit == 0 || _doit == 1) {
149
150 //...Create patterns...
151 unsigned axialSuperLayerId = 0;
152 for (unsigned i = 0; i < _cdc.nSegmentLayers(); i++) {
153 const Belle2::TRGCDCLayer* l = _cdc.segmentLayer(i);
154 const unsigned nWires = l->nCells();
155
156 if (! nWires) continue;
157 if ((* l)[0]->stereo()) continue;
158
159 _plane[0]->preparePatterns(axialSuperLayerId, nWires);
160 _plane[1]->preparePatterns(axialSuperLayerId, nWires);
161 for (unsigned j = 0; j < nWires; j++) {
162 const TCCell& w = * (* l)[j];
163 const float x = w.xyPosition().x();
164 const float y = w.xyPosition().y();
165
166 _plane[0]->clear();
167 _plane[0]->vote(x, y, +1, axialSuperLayerId, 1);
168 _plane[0]->registerPattern(axialSuperLayerId, j);
169
170 _plane[1]->clear();
171 _plane[1]->vote(x, y, -1, axialSuperLayerId, 1);
172 _plane[1]->registerPattern(axialSuperLayerId, j);
173 }
174 ++axialSuperLayerId;
175 }
176 }
177
178 //...Kaiyu's methode...
179 else {
180 mappingByFile(mappingFilePlus, mappingFileMinus);
181 }
182 }
183
184 TRGCDCHoughFinder::~TRGCDCHoughFinder()
185 {
186#ifdef TRGCDC_DISPLAY_HOUGH
187 if (H0)
188 delete H0;
189 if (H1)
190 delete H1;
191 cout << "TRGCDCHoughFinder ... Hough displays deleted" << endl;
192#endif
193 }
194
195 int
196 TRGCDCHoughFinder::doFindingTrasan(vector<unsigned> peaks[],
197 vector<TCTrack*>& trackList2D) const
198 {
199
200 const string sn = "Hough Finder Finding (trasan version)";
201 TRGDebug::enterStage(sn);
202
203 //...Initialization...
204 _plane[0]->clear();
205 _plane[1]->clear();
206
207 //...Voting...
208 unsigned nLayers = _cdc.nAxialSuperLayers();
209 for (unsigned i = 0; i < nLayers; i++) {
210 const vector<const TCSHit*> hits = _cdc.axialSegmentHits(i);
211 for (unsigned j = 0; j < hits.size(); j++) {
212 _plane[0]->vote(i, hits[j]->cell().localId());
213 _plane[1]->vote(i, hits[j]->cell().localId());
214 }
215 }
216 _plane[0]->merge();
217 _plane[1]->merge();
218
219#ifdef TRGCDC_DISPLAY_HOUGH
220 string stg = "2D : Hough : Results of Peak Finding";
221 string inf = " ";
222 H0->stage(stg);
223 H0->information(inf);
224 H0->clear();
225 H0->area().append(_plane[0]);
226 H0->show();
227 H1->stage(stg);
228 H1->information(inf);
229 H1->clear();
230 H1->area().append(_plane[1]);
231 H1->show();
232#endif
233
234 //...Look for peaks which have 5 hits...
235 _peakFinder.findPeaksTrasan(* _plane[0], _peakMin, false, peaks[0]);
236 _peakFinder.findPeaksTrasan(* _plane[1], _peakMin, false, peaks[1]);
237
238 //...Peak loop to pick up segment hits...
239 // (no fit, using peak position only)
240 for (unsigned pm = 0; pm < 2; pm++) {
241 for (unsigned i = 0; i < peaks[pm].size(); i++) {
242 const unsigned peakId = peaks[pm][i];
243
244 //...Make a track...
245 TCTrack* t = makeTrack(peakId, pm);
246 trackList2D.push_back(t);
247
248#ifdef TRGCDC_DISPLAY_HOUGH
249 vector<const TCTrack*> cc;
250 cc.push_back(t);
251 const string stgNow = "doFinding : track made";
252 const string infNow = " ";
253 D->clear();
254 D->stage(stgNow);
255 D->information(infNow);
256 D->area().append(cc, Gdk::Color("#FF0066009900"));
257 D->area().append(_cdc.hits());
258 D->area().append(_cdc.segmentHits());
259 D->show();
260 D->run();
261#endif
262 }
263 }
264
265 TRGDebug::leaveStage(sn);
266 return 0;
267 }
268
269 vector<TCLink*>
270 TRGCDCHoughFinder::selectBestHits(const vector<TCLink*>& links) const
271 {
272 vector<TCLink*> bests;
273 vector<TCLink*> layers[9];
274 TCLink::separate(links, 9, layers);
275
276 if (TRGDebug::level()) {
277 for (unsigned i = 0; i < 9; i++) {
278 cout << TRGDebug::tab() << "layer " << i << endl;
279 TCLink::dump(layers[i], "", TRGDebug::tab(4));
280 }
281 }
282
283 //...Select links to be removed...
284 for (unsigned i = 0; i < 9; i++) {
285 if (layers[i].size() == 0) continue;
286 if (layers[i].size() == 1) {
287 bests.push_back(layers[i][0]);
288 continue;
289 }
290
291 TCLink* best = layers[i][0];
292 int timeMin = (layers[i][0]->cell()->signal())[0]->time();
293 for (unsigned j = 1; j < layers[i].size(); j++) {
294 const TRGTime& t = * (layers[i][j]->cell()->signal())[0];
295 if (t.time() < timeMin) {
296 timeMin = t.time();
297 best = layers[i][j];
298 }
299 }
300
301 bests.push_back(best);
302 }
303 return bests;
304 }
305
306 int
307 TRGCDCHoughFinder::doFittingTrasan(vector<unsigned> peaks[],
308 vector<TRGCDCTrack*>& trackList2DFitted) const
309 {
310
311 const string sn = "Hough Finder Fitting (trasan version)";
312 TRGDebug::enterStage(sn);
313
314 //...Peak loop...
315 unsigned nCircles = 0;
316 for (unsigned pm = 0; pm < 2; pm++) {
317 for (unsigned i = 0; i < peaks[pm].size(); i++) {
318 const unsigned peakId = peaks[pm][i];
319
320 //...Get segment hits...
321 vector<TCLink*> links;
322 vector<const TCSegment*> segments;
323 const unsigned nLayers = _plane[pm]->nLayers();
324 for (unsigned j = 0; j < nLayers; j++) {
325 const vector<unsigned>& ptn =
326 _plane[pm]->patternId(j, peakId);
327 for (unsigned k = 0; k < ptn.size(); k++) {
328 const TCSegment& s = _cdc.axialSegment(j, ptn[k]);
329 segments.push_back(& s);
330 if (s.hit()) {
331 TCLink* l = new TCLink(0, s.hit());
332 links.push_back(l);
333 }
334 }
335 }
336
337 //...Select best hits in each layer...
338 const vector<TCLink*> bests = selectBestHits(links);
339
340 //...Make a circle...
341 TCCircle c(bests);
342 c.fit();
343 c.name("CircleFitted_" + TRGUtil::itostring(nCircles));
344 ++nCircles;
345
346 if (TRGDebug::level()) {
347 cout << TRGDebug::tab() << "peak#" << nCircles << ":"
348 << "plane" << pm << ",serialId=" << peakId << endl;
349 cout << TRGDebug::tab() << "segments below" << endl;
350 cout << TRGDebug::tab(4);
351 for (unsigned j = 0; j < segments.size(); j++) {
352 cout << segments[j]->name();
353 if (j != (segments.size() - 1))
354 cout << ",";
355 }
356 cout << endl;
357 cout << TRGDebug::tab() << "best links below" << endl;
358 TCLink::dump(bests, "", TRGDebug::tab(1));
359 c.dump("detail", TRGDebug::tab() + "Circle> ");
360 }
361
362 //...Make a track...
363 TCTrack& t = * new TCTrack(c);
364 t.name("Track_" + TRGUtil::itostring(i));
365 trackList2DFitted.push_back(& t);
366
367 if (TRGDebug::level()) {
368 t.relation().dump("", TRGDebug::tab());
369 t.dump("detail", TRGDebug::tab(1));
370 }
371
372#ifdef TRGCDC_DISPLAY_HOUGH
373 vector<const TCCircle*> cc;
374 cc.push_back(& c);
375 const string stg = "doFitting : trasan method";
376 const string inf = " ";
377 D->clear();
378 D->stage(stg);
379 D->information(inf);
380 D->area().append(_cdc.hits());
381 D->area().append(_cdc.segmentHits());
382 D->area().append(cc, Gdk::Color("blue"));
383 D->show();
384 D->run();
385#endif
386 }
387 }
388
389 TRGDebug::leaveStage(sn);
390
391 return 0;
392 }
393
394 void
395 TRGCDCHoughFinder::terminate()
396 {
397 if (_commonData) delete _commonData;
398 }
399
400 double TRGCDCHoughFinder::calPhi(TRGCDCSegmentHit const* segmentHit, double eventTime)
401 {
402 const CDC::CDCGeometryPar& cdcp = CDC::CDCGeometryPar::Instance();
403 unsigned localId = segmentHit->segment().center().localId();
404 unsigned layerId = segmentHit->segment().center().layerId();
405 int nWires = cdcp.nWiresInLayer(layerId) * 2;
406 double rr = cdcp.senseWireR(layerId);
407 double driftTime = segmentHit->segment().priorityTime();
408 int lr = segmentHit->segment().LUT()->getValue(segmentHit->segment().lutPattern());
409 return Fitter3DUtility::calPhi(localId, nWires, driftTime, eventTime, rr, lr);
410 }
411
412 void TRGCDCHoughFinder::calCosPhi(std::map<std::string, Belle2::TRGCDCJSignal>& mSignalStorage,
413 std::map<std::string, Belle2::TRGCDCJLUT* >& mLutStorage)
414 {
415
417 Belle2::TRGCDCJSignalData* commonData = mSignalStorage.begin()->second.getCommonData();
418 {
419 mSignalStorage["invPhiAxMin"] = Belle2::TRGCDCJSignal(M_PI, mSignalStorage["phi_4"].getToReal(), commonData);
420 mSignalStorage["invPhiAxMax"] = Belle2::TRGCDCJSignal(2 * M_PI, mSignalStorage["phi_4"].getToReal(), commonData);
421 }
422
423 for (unsigned iSt = 0; iSt < 5; iSt++) {
424 string t_inputName = "phi_" + to_string(iSt);
425 string t_outputName = "cosPhi_" + to_string(iSt);
426
427 if (!mLutStorage.count(t_outputName)) {
428 mLutStorage[t_outputName] = new Belle2::TRGCDCJLUT(t_outputName);
429 mLutStorage[t_outputName]->setFloatFunction(
430 [ = ](double aValue) -> double{return cos(aValue);},
431 mSignalStorage[t_inputName],
432 mSignalStorage["invPhiAxMin"], mSignalStorage["invPhiAxMax"], mSignalStorage[t_inputName].getToReal(),
433 12, 12);
434 };//if name
435
436 mLutStorage[t_outputName]->operate(mSignalStorage[t_inputName], mSignalStorage[t_outputName]);
437 } // end for
438 }//calCosPhi
439
440 void TRGCDCHoughFinder::calSinPhi(std::map<std::string, Belle2::TRGCDCJSignal>& mSignalStorage,
441 std::map<std::string, Belle2::TRGCDCJLUT* >& mLutStorage)
442 {
443 Belle2::TRGCDCJSignalData* commonData = mSignalStorage.begin()->second.getCommonData();
444 {
445 mSignalStorage["invPhiAxMin"] = Belle2::TRGCDCJSignal(-M_PI / 2, mSignalStorage["phi_4"].getToReal(), commonData);
446 mSignalStorage["invPhiAxMax"] = Belle2::TRGCDCJSignal(M_PI / 2, mSignalStorage["phi_4"].getToReal(), commonData);
447 }
448
449 for (unsigned iSt = 0; iSt < 5; iSt++) {
450 string t_sininputName = "phi_" + to_string(iSt);
451 string t_sinoutputName = "sinPhi_" + to_string(iSt);
452
453 if (!mLutStorage.count(t_sinoutputName)) {
454 mLutStorage[t_sinoutputName] = new Belle2::TRGCDCJLUT(t_sinoutputName);
455
456 mLutStorage[t_sinoutputName]->setFloatFunction(
457 [ = ](double aValue) -> double{return sin(aValue);},
458 mSignalStorage[t_sininputName],
459 mSignalStorage["invPhiAxMin"], mSignalStorage["invPhiAxMax"], mSignalStorage[t_sininputName].getToReal(),
460 12, 12);
461 };//if name
462
463 mLutStorage[t_sinoutputName]->operate(mSignalStorage[t_sininputName], mSignalStorage[t_sinoutputName]);
464 } // end for
465 }//calSinPhi
466
467 void TRGCDCHoughFinder::rPhi(std::map<std::string, Belle2::TRGCDCJSignal>& mSignalStorage,
468 std::map<std::string, Belle2::TRGCDCJLUT* >& mLutStorage)
469 {
470 Belle2::TRGCDCJSignalData* commonData = mSignalStorage.begin()->second.getCommonData();
471 for (unsigned iSt = 0; iSt < 5; iSt++) {
472 mSignalStorage["cosPhiMul_" + to_string(iSt)] <= mSignalStorage["cosPhi_" + to_string(iSt)] * mSignalStorage["cosPhi_" + to_string(
473 iSt)];
474 }
475
476 for (unsigned iSt = 0; iSt < 5; iSt++) {
477 mSignalStorage["sinPhiMul_" + to_string(iSt)] <= mSignalStorage["sinPhi_" + to_string(iSt)] * mSignalStorage["sinPhi_" + to_string(
478 iSt)];
479 }
480
481 for (unsigned iSt = 0; iSt < 5; iSt++) {
482 mSignalStorage["cossinPhiMul_" + to_string(iSt)] <= mSignalStorage["cosPhi_" + to_string(iSt)] * mSignalStorage["sinPhi_" +
483 to_string(iSt)];
484 }
485
486 for (unsigned iSt = 0; iSt < 5; iSt++) {
487 mSignalStorage["rcosPhiMul_" + to_string(iSt)] <= mSignalStorage["cosPhi_" + to_string(iSt)] * mSignalStorage["2Drr_" + to_string(
488 iSt)];
489 }
490
491 for (unsigned iSt = 0; iSt < 5; iSt++) {
492 mSignalStorage["rsinPhiMul_" + to_string(iSt)] <= mSignalStorage["sinPhi_" + to_string(iSt)] * mSignalStorage["2Drr_" + to_string(
493 iSt)];
494 }
495 mSignalStorage["cossum_p1_0"] <= mSignalStorage["cosPhiMul_0"] + mSignalStorage["cosPhiMul_1"];
496 mSignalStorage["cossum_p1_1"] <= mSignalStorage["cosPhiMul_2"] + mSignalStorage["cosPhiMul_3"];
497 mSignalStorage["cossum_p1_2"] <= mSignalStorage["cosPhiMul_4"] ;
498 mSignalStorage["cossum_p2_0"] <= mSignalStorage["cossum_p1_0"] + mSignalStorage["cossum_p1_1"];
499 mSignalStorage["cossum_p2_1"] <= mSignalStorage["cossum_p1_2"] ;
500 mSignalStorage["cossum"] <= mSignalStorage["cossum_p2_0"] + mSignalStorage["cossum_p2_1"];
501 mSignalStorage["sinsum_p1_0"] <= mSignalStorage["sinPhiMul_0"] + mSignalStorage["sinPhiMul_1"];
502 mSignalStorage["sinsum_p1_1"] <= mSignalStorage["sinPhiMul_2"] + mSignalStorage["sinPhiMul_3"];
503 mSignalStorage["sinsum_p1_2"] <= mSignalStorage["sinPhiMul_4"] ;
504 mSignalStorage["sinsum_p2_0"] <= mSignalStorage["sinsum_p1_0"] + mSignalStorage["sinsum_p1_1"];
505 mSignalStorage["sinsum_p2_1"] <= mSignalStorage["sinsum_p1_2"] ;
506 mSignalStorage["sinsum"] <= mSignalStorage["sinsum_p2_0"] + mSignalStorage["sinsum_p2_1"];
507 mSignalStorage["cossinsum_p1_0"] <= mSignalStorage["cossinPhiMul_0"] + mSignalStorage["cossinPhiMul_1"];
508 mSignalStorage["cossinsum_p1_1"] <= mSignalStorage["cossinPhiMul_2"] + mSignalStorage["cossinPhiMul_3"];
509 mSignalStorage["cossinsum_p1_2"] <= mSignalStorage["cossinPhiMul_4"] ;
510 mSignalStorage["cossinsum_p2_0"] <= mSignalStorage["cossinsum_p1_0"] + mSignalStorage["cossinsum_p1_1"];
511 mSignalStorage["cossinsum_p2_1"] <= mSignalStorage["cossinsum_p1_2"] ;
512 mSignalStorage["cossinsum"] <= mSignalStorage["cossinsum_p2_0"] + mSignalStorage["cossinsum_p2_1"];
513 mSignalStorage["rcossum_p1_0"] <= mSignalStorage["rcosPhiMul_0"] + mSignalStorage["rcosPhiMul_1"];
514 mSignalStorage["rcossum_p1_1"] <= mSignalStorage["rcosPhiMul_2"] + mSignalStorage["rcosPhiMul_3"];
515 mSignalStorage["rcossum_p1_2"] <= mSignalStorage["rcosPhiMul_4"] ;
516 mSignalStorage["rcossum_p2_0"] <= mSignalStorage["rcossum_p1_0"] + mSignalStorage["rcossum_p1_1"];
517 mSignalStorage["rcossum_p2_1"] <= mSignalStorage["rcossum_p1_2"] ;
518 mSignalStorage["rcossum"] <= mSignalStorage["rcossum_p2_0"] + mSignalStorage["rcossum_p2_1"];
519 mSignalStorage["rsinsum_p1_0"] <= mSignalStorage["rsinPhiMul_0"] + mSignalStorage["rsinPhiMul_1"];
520 mSignalStorage["rsinsum_p1_1"] <= mSignalStorage["rsinPhiMul_2"] + mSignalStorage["rsinPhiMul_3"];
521 mSignalStorage["rsinsum_p1_2"] <= mSignalStorage["rsinPhiMul_4"] ;
522 mSignalStorage["rsinsum_p2_0"] <= mSignalStorage["rsinsum_p1_0"] + mSignalStorage["rsinsum_p1_1"];
523 mSignalStorage["rsinsum_p2_1"] <= mSignalStorage["rsinsum_p1_2"] ;
524 mSignalStorage["rsinsum"] <= mSignalStorage["rsinsum_p2_0"] + mSignalStorage["rsinsum_p2_1"];
525 mSignalStorage["hcx_p1_0"] <= mSignalStorage["rcossum"] * mSignalStorage["sinsum"];
526 mSignalStorage["hcx_p1_1"] <= mSignalStorage["rsinsum"] * mSignalStorage["cossinsum"];
527 mSignalStorage["hcx_p2"] <= mSignalStorage["hcx_p1_0"] - mSignalStorage["hcx_p1_1"];
528 mSignalStorage["hcy_p1_0"] <= mSignalStorage["rsinsum"] * mSignalStorage["cossum"];
529 mSignalStorage["hcy_p1_1"] <= mSignalStorage["rcossum"] * mSignalStorage["cossinsum"];
530 mSignalStorage["hcy_p2"] <= mSignalStorage["hcy_p1_0"] - mSignalStorage["hcy_p1_1"];
531 mSignalStorage["den_p1"] <= mSignalStorage["cossum"] * mSignalStorage["sinsum"];
532 mSignalStorage["den_p2"] <= mSignalStorage["cossinsum"] * mSignalStorage["cossinsum"];
533 mSignalStorage["den"] <= mSignalStorage["den_p1"] - mSignalStorage["den_p2"];
534 {
535 mSignalStorage["denMin"] = Belle2::TRGCDCJSignal(0.00315354, mSignalStorage["den"].getToReal(),
536 commonData);//den Min = 0.00630708*0.5
537 mSignalStorage["denMax"] = Belle2::TRGCDCJSignal(1.332705, mSignalStorage["den"].getToReal(), commonData);//den Max = 0.88847*1.5
538 }
539
540 // Constrain den.
541 {
542 // Make ifElse data.
543 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
544 // Compare
545 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["den"], ">", mSignalStorage["denMax"]);
546 // Assignments
547 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
548 make_pair(&mSignalStorage["den_c"], mSignalStorage["denMax"])
549 };
550 // Push to data.
551 t_data.push_back(make_pair(t_compare, t_assigns));
552 // Compare
553 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["den"], ">", mSignalStorage["denMin"]);
554 // Assignments
555 t_assigns = {
556 make_pair(&mSignalStorage["den_c"], mSignalStorage["den"].limit(mSignalStorage["denMin"], mSignalStorage["denMax"]))
557 };
558 // Push to data.
559 t_data.push_back(make_pair(t_compare, t_assigns));
560 // Compare
561 t_compare = Belle2::TRGCDCJSignal();
562 // Assignments
563 t_assigns = {
564 make_pair(&mSignalStorage["den_c"], mSignalStorage["denMin"])
565 };
566 // Push to data.
567 t_data.push_back(make_pair(t_compare, t_assigns));
568 Belle2::TRGCDCJSignal::ifElse(t_data);
569 }
570
571 {
572 unsigned long long t_int = mSignalStorage["den_c"].getMaxInt();
573 double t_toReal = mSignalStorage["den_c"].getToReal();
574 double t_actual = mSignalStorage["den_c"].getMaxActual();
575 mSignalStorage["iDenMin"] = Belle2::TRGCDCJSignal(t_int, t_toReal, t_int, t_int, t_actual, t_actual, t_actual, -1, commonData);
576 t_int = mSignalStorage["den_c"].getMinInt();
577 t_actual = mSignalStorage["den_c"].getMinActual();
578 mSignalStorage["iDenMax"] = Belle2::TRGCDCJSignal(t_int, t_toReal, t_int, t_int, t_actual, t_actual, t_actual, -1, commonData);
579 }
580
581 if (!mLutStorage.count("iDen")) {
582 mLutStorage["iDen"] = new Belle2::TRGCDCJLUT("iDen");
583 mLutStorage["iDen"]->setFloatFunction(
584 [ = ](double aValue) -> double{return 0.5 / aValue;},
585 mSignalStorage["den_c"],
586 mSignalStorage["iDenMin"], mSignalStorage["iDenMax"], pow(1 / mSignalStorage["cossinsum"].getToReal(), 2),
587 12, 12);
588 }
589
590 //Operate using LUT(iDen = 1/2/den)
591 mLutStorage["iDen"]->operate(mSignalStorage["den_c"], mSignalStorage["iDen"]);
592
593 mSignalStorage["hcx"] <= mSignalStorage["hcx_p2"] * mSignalStorage["iDen"];
594 mSignalStorage["hcy"] <= mSignalStorage["hcy_p2"] * mSignalStorage["iDen"];
595
596 mSignalStorage["hcx2"] <= mSignalStorage["hcx"] * mSignalStorage["hcx"];
597 mSignalStorage["hcy2"] <= mSignalStorage["hcy"] * mSignalStorage["hcy"];
598 mSignalStorage["sumhc"] <= mSignalStorage["hcx2"] + mSignalStorage["hcy2"];
599 {
600 mSignalStorage["sumhcMin"] = Belle2::TRGCDCJSignal(4489, mSignalStorage["sumhc"].getToReal(), commonData);
601 mSignalStorage["sumhcMax"] = Belle2::TRGCDCJSignal(2560000, mSignalStorage["sumhc"].getToReal(), commonData);
602 }
603
604 // Constrain den.
605 {
606 // Make ifElse data.
607 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
608 // Compare
609 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["sumhc"], ">", mSignalStorage["sumhcMax"]);
610 // Assignments
611 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
612 make_pair(&mSignalStorage["sumhc_c"], mSignalStorage["sumhcMax"])
613 };
614 // Push to data.
615 t_data.push_back(make_pair(t_compare, t_assigns));
616 // Compare
617 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["sumhc"], ">", mSignalStorage["sumhcMin"]);
618 // Assignments
619 t_assigns = {
620 make_pair(&mSignalStorage["sumhc_c"], mSignalStorage["sumhc"].limit(mSignalStorage["sumhcMin"], mSignalStorage["sumhcMax"]))
621 };
622 // Push to data.
623 t_data.push_back(make_pair(t_compare, t_assigns));
624 // Compare
625 t_compare = Belle2::TRGCDCJSignal();
626 // Assignments
627 t_assigns = {
628 make_pair(&mSignalStorage["sumhc_c"], mSignalStorage["sumhcMin"])
629 };
630 // Push to data.
631 t_data.push_back(make_pair(t_compare, t_assigns));
632 Belle2::TRGCDCJSignal::ifElse(t_data);
633 }
634 {
635 unsigned long long t_int = mSignalStorage["sumhc_c"].getMaxInt();
636 double t_toReal = mSignalStorage["sumhc_c"].getToReal();
637 double t_actual = mSignalStorage["sumhc_c"].getMaxActual();
638 mSignalStorage["iRhoMax"] = Belle2::TRGCDCJSignal(t_int, t_toReal, t_int, t_int, t_actual, t_actual, t_actual, -1, commonData);
639 t_int = mSignalStorage["sumhc_c"].getMinInt();
640 t_actual = mSignalStorage["sumhc_c"].getMinActual();
641 mSignalStorage["iRhoMin"] = Belle2::TRGCDCJSignal(t_int, t_toReal, t_int, t_int, t_actual, t_actual, t_actual, -1, commonData);
642 }
643 // Generate sqrt LUT sqrt(hcx*hcx + hcy*hcy)
644 if (!mLutStorage.count("iRho")) {
645 mLutStorage["iRho"] = new Belle2::TRGCDCJLUT("iRho");
646 mLutStorage["iRho"]->setFloatFunction(
647 [ = ](double aValue) -> double{return abs(sqrt(aValue));},
648 mSignalStorage["sumhc_c"],
649 mSignalStorage["iRhoMin"], mSignalStorage["iRhoMax"], mSignalStorage["sumhc"].getToReal(),
650 24, 12);
651 }
652 // Operate using LUT(iRho = sqrt(sumhc))
653 mLutStorage["iRho"]->operate(mSignalStorage["sumhc_c"], mSignalStorage["iRho"]);
654 if (TRGDebug::level()) {
655 cout << "<<<iRho>>>" << endl;
656 mSignalStorage["iRho"].dump();
657 }
658 }//rphi
659
660 void
661 TRGCDCHoughFinder::mappingByFile(const string& mappingFilePlus,
662 const string& mappingFileMinus)
663 {
664 const string sn = "mappingByFile";
665 TRGDebug::enterStage(sn);
666
667 const string* fMap[2] = {& mappingFilePlus, & mappingFileMinus};
668
669 for (unsigned f = 0; f < 2; f++) {
670
671 const string& fn = * fMap[f];
672 ifstream infile(fn.c_str(), ios::in);
673 if (infile.fail()) {
674 cout << " !!! can not open file" << endl
675 << " " << fn << endl
676 << " Mapping aborted" << endl;
677 return;
678 }
679
680 //...Map data storage...
681 vector<unsigned> x;
682 vector<unsigned> y;
683 vector<vector<unsigned>> tsf;
684
685 //...Ignore lines not starting with a digit...
686 string ignores;
687 while (!isdigit(infile.peek())) {
688 getline(infile, ignores);
689 }
690 //...Skip the first line (read in constructor)...
691 getline(infile, ignores);
692
693 //...Read map file...
694 string car;
695 while (getline(infile, car)) {
696 unsigned i = 0;
697 unsigned b;
698 istringstream in(car);
699 vector<unsigned> slts;
700 while (in >> b) {
701 ++i;
702 unsigned j = i % 2;
703 if (i == 1) // cell position x
704 x.push_back(b);
705 if (i == 2) // cell position y (stored with offset 1)
706 y.push_back(b - 1);
707 if (j != 0 && i != 1) // TSF SL
708 slts.push_back(b);
709 if (j != 1 && i != 2) // TSF local ID
710 slts.push_back(b);
711 }
712 tsf.push_back(slts);
713 }
714 infile.close();
715
716 unsigned axialSuperLayerId = 0;
717 for (unsigned i = 0; i < _cdc.nSegmentLayers(); i++) {
718 const Belle2::TRGCDCLayer* l = _cdc.segmentLayer(i);
719 const unsigned nWires = l->nCells();
720
721 if (! nWires) continue;
722 if ((* l)[0]->stereo()) continue;
723
724 _plane[f]->preparePatterns(axialSuperLayerId, nWires);
725 for (unsigned j = 0; j < nWires; j++) {
726
727 _plane[f]->clear();
728
729 //...loop over all hp cells...
730 const unsigned n = x.size();
731 for (unsigned k = 0; k < n; k++) {
732 const int ix = x[k];
733 const int iy = y[k];
734 const unsigned sid = _plane[f]->serialId(ix, iy);
735 for (unsigned itsf = 0; itsf < tsf[k].size();) {
736 const unsigned sl = tsf[k][itsf];
737 const unsigned id = tsf[k][itsf + 1];
738 if ((i == sl) && (j == id)) {
739 _plane[f]->setEntry(sid, axialSuperLayerId, 1);
740 }
741 itsf += 2;
742 }
743 }
744
745 _plane[f]->registerPattern(axialSuperLayerId, j);
746 }
747 ++axialSuperLayerId;
748 }
749 }
750
751 TRGDebug::leaveStage(sn);
752 }
753
754 int
755 TRGCDCHoughFinder::FindAndFit(std::vector<TRGCDCTrack*>& trackList2D,
756 std::vector<TRGCDCTrack*>& trackList2DFitted)
757 {
758
759 if (_doit == 0 || _doit == 1)
760 return doFindingAndFittingTrasan(trackList2D, trackList2DFitted);
761 else
762 return doFindingAndFitting(trackList2D, trackList2DFitted);
763 }
764
765 int
766 TRGCDCHoughFinder::doFindingAndFittingTrasan(std::vector<TRGCDCTrack*>& trackList2D,
767 std::vector<TRGCDCTrack*>& trackList2DFitted)
768 {
769
770 const string sn = "HoughFinder::doFindingAndFitting (trasan version)";
771 TRGDebug::enterStage(sn);
772
773 vector<unsigned> peaks[2];
774 doFindingTrasan(peaks, trackList2D);
775 doFittingTrasan(peaks, trackList2DFitted);
776
777 TRGDebug::leaveStage(sn);
778
779 return 0;
780 }
781
782 int
783 TRGCDCHoughFinder::doFindingAndFitting(std::vector<TRGCDCTrack*>& trackList2D,
784 std::vector<TRGCDCTrack*>& trackList2DFitted)
785 {
786
787 const string sn = "HoughFinder::doFindingAndFitting (kaiyu version)";
788 TRGDebug::enterStage(sn);
789
790 vector<vector<unsigned>> peaks[2];
791 doFinding(peaks, trackList2D);
792 doFitting(trackList2D, trackList2DFitted);
793
794 TRGDebug::leaveStage(sn);
795
796 return 0;
797 }
798
799 int
800 TRGCDCHoughFinder::doFinding(std::vector<std::vector<unsigned>> peaks[],
801 std::vector<TRGCDCTrack*>& trackList2D)
802 {
803
804 const string sn = "HoughFinder::doFinding";
805 TRGDebug::enterStage(sn);
806
807 //...Initialization...
808 _plane[0]->clear();
809 _plane[1]->clear();
810
811 //...Voting...
812 unsigned nLayers = _cdc.nAxialSuperLayers();
813 for (unsigned i = 0; i < nLayers; i++) {
814 const vector<const TCSHit*> hits = _cdc.axialSegmentHits(i);
815 for (unsigned j = 0; j < hits.size(); j++) {
816 _plane[0]->vote(i, hits[j]->cell().localId());
817 _plane[1]->vote(i, hits[j]->cell().localId());
818 }
819 }
820 _plane[0]->merge();
821 _plane[1]->merge();
822
823 //...Look for peaks which have 5 hits...
824 _peakFinder.findPeaks(* _plane[0], _peakMin, peaks[0]);
825 _peakFinder.findPeaks(* _plane[1], _peakMin, peaks[1]);
826
827 //...Peak loop to make tracks (no fit, using hough position only)...
828 for (unsigned pm = 0; pm < 2; pm++) {
829 for (unsigned i = 0; i < peaks[pm].size(); i++) {
830 const unsigned peakId = _plane[pm]->serialId(peaks[pm][i][0],
831 peaks[pm][i][1]);
832 TCTrack* t = makeTrack(peakId, pm);
833 trackList2D.push_back(t);
834
835#ifdef TRGCDC_DISPLAY_HOUGH
836 vector<const TCTrack*> cc;
837 cc.push_back(t);
838 const string stg = "doFinding : track made";
839 const string inf = " ";
840 D->clear();
841 D->stage(stg);
842 D->information(inf);
843 D->area().append(cc, Gdk::Color("#FF0066009900"));
844 D->area().append(_cdc.hits());
845 D->area().append(_cdc.segmentHits());
846 D->show();
847 D->run();
848#endif
849 }
850 }
851
852 TRGDebug::leaveStage(sn);
853 return 0;
854 }
855
856 int
857 TRGCDCHoughFinder::doFitting(std::vector<TRGCDCTrack*>& trackList2D,
858 std::vector<TRGCDCTrack*>& trackList2DFitted)
859 {
860
861 const string sn = "Hough Finder Fitting";
862 TRGDebug::enterStage(sn);
863
864 if (_commonData == 0) {
865 _commonData = new Belle2::TRGCDCJSignalData();
866 }
867
868 //...Event time...
869 bool fEvtTime = true;
870 // cppcheck-suppress knownConditionTrueFalse
871 double eventTime = fEvtTime ? _cdc.getEventTime() : 0;
872
873 // Loop over all tracks
874 for (unsigned int iTrack = 0; iTrack < trackList2D.size(); iTrack++) {
875 TCTrack& aTrack = * new TCTrack(* trackList2D[iTrack]);
876 trackList2DFitted.push_back(& aTrack);
877
878 CDC::CDCGeometryPar& cdcp = CDC::CDCGeometryPar::Instance();
879 double phi02D;
880 vector<double> nWires(9);
881 vector<double> rr(9);
882 vector<double> rr2D;
883 vector<double> wirePhi2DError(5);
884 wirePhi2DError[0] = 0.0085106;
885 wirePhi2DError[1] = 0.0039841;
886 wirePhi2DError[2] = 0.0025806;
887 wirePhi2DError[3] = 0.0019084;
888 wirePhi2DError[4] = 0.001514;
889 vector<double> driftPhi2DError(5);
890 driftPhi2DError[0] = 0.0085106;
891 driftPhi2DError[1] = 0.0039841;
892 driftPhi2DError[2] = 0.0025806;
893 driftPhi2DError[3] = 0.0019084;
894 driftPhi2DError[4] = 0.001514;
895 // Check if all superlayers have one TS
896 bool trackFull = 1;
897 for (unsigned iSL = 0; iSL < _cdc.nSuperLayers(); iSL = iSL + 2) {
898 // Check if all superlayers have one TS
899 const vector<TCLink*>& links = aTrack.links(iSL);
900 const unsigned nSegments = links.size();
901
902 if (TRGDebug::level())
903 cout << TRGDebug::tab() << "#segments in SL" << iSL << " : "
904 << nSegments << endl;
905
906 // Find if there is a TS with a priority hit.
907 // Loop over all TS in same superlayer.
908 bool priorityHitTS = 0;
909 for (unsigned iTS = 0; iTS < nSegments; iTS++) {
910 const TCSegment* _segment = dynamic_cast<const TCSegment*>(& links[iTS]->hit()->cell());
911 if (_segment->center().hit() != 0) priorityHitTS = 1;
912 }
913 if (nSegments != 1) {
914 if (nSegments == 0) {
915 trackFull = 0;
916 if (TRGDebug::level())
917 cout << TRGDebug::tab()
918 << "=> Not enough TS." << endl;
919 } else {
920 if (TRGDebug::level())
921 cout << TRGDebug::tab()
922 << "=> multiple TS are assigned." << endl;
923 }
924 } else {
925 if (priorityHitTS == 0) {
926 trackFull = 0;
927 if (TRGDebug::level())
928 cout << TRGDebug::tab()
929 << "=> There are no priority hit TS" << endl;
930 }
931 }
932 } // End superlayer loop
933 if (trackFull == 0) {
934 TRGCDCHelix helix(ORIGIN, CLHEP::HepVector(5, 0), CLHEP::HepSymMatrix(5, 0));
935 CLHEP::HepVector helixParameters(5);
936 helixParameters = aTrack.helix().a();
937 aTrack.setFitted(0);
938 aTrack.setHelix(helix);
939 continue;
940 }
941
942 for (unsigned iSL = 0; iSL < 9; iSL++) {
943 unsigned _layerId = _cdc.segment(iSL, 0).center().layerId();
944 rr[iSL] = cdcp.senseWireR(_layerId);
945 nWires[iSL] = cdcp.nWiresInLayer(_layerId) * 2;
946 }
947 rr2D = vector<double>({rr[0], rr[2], rr[4], rr[6], rr[8]});
948
949 // Get wirePhi, Drift information
950 vector<double>wirePhi(9);
951 vector<double>driftLength(9);
952 vector<double>LR(9);
953 vector<double>lutLR(9);
954 vector<double>mcLR(9);
955 bool fmcLR = false, fLRLUT = true;
956 for (unsigned iSL = 0; iSL < 9; iSL = iSL + 2) {
957 const vector<TCLink*>& links = aTrack.links(iSL);
958 const TCSegment* _segment = dynamic_cast<const TCSegment*>(& links[0]->hit()->cell());
959 wirePhi[iSL] = _segment->localId() / nWires[iSL] * 4 * M_PI;
960 lutLR[iSL] = _segment->LUT()->getValue(_segment->lutPattern());
961 mcLR[iSL] = _segment->hit()->mcLR();
962 driftLength[iSL] = _segment->hit()->drift();
963 // cppcheck-suppress knownConditionTrueFalse
964 if (fmcLR) LR[iSL] = mcLR[iSL];
965 // cppcheck-suppress knownConditionTrueFalse
966 else if (fLRLUT) LR[iSL] = lutLR[iSL];
967 else LR[iSL] = 3;
968 }//End superlayer loop
969 // 2D fit values from IW 2D fitter
970 phi02D = aTrack.helix().phi0();
971 if (aTrack.charge() < 0) {
972 phi02D = phi02D - M_PI;
973 if (phi02D < 0) phi02D = phi02D + 2 * M_PI;
974 }
975 // Get 2D fit values from JB 2D fitter
976 // Set phi2DError for 2D fit
977 vector<double>phi2DError(5);
978 for (unsigned iAx = 0; iAx < 5; iAx++) {
979 if (LR[2 * iAx] != 2) phi2DError[iAx] = driftPhi2DError[iAx];
980 else phi2DError[iAx] = wirePhi2DError[iAx];
981 }
982 // Calculate phi2D using driftTime.
983 vector<double>phi2D(5);
984 for (unsigned iAx = 0; iAx < 5; iAx++) {
985 phi2D[iAx] = Fitter3DUtility::calPhi(wirePhi[iAx * 2],
986 driftLength[iAx * 2],
987 eventTime,
988 rr[iAx * 2],
989 LR[iAx * 2]);
990 if (TRGDebug::level()) {
991 cout << TRGDebug::tab() << "eventTime: " << eventTime << endl;
992 for (int i = 0 ; i < 5 ; i++) {
993 cout << TRGDebug::tab() << "phi2D: : " << phi2D[i] << endl;
994 cout << TRGDebug::tab() << "wirePhi: " << wirePhi[i * 2]
995 << endl;
996 cout << TRGDebug::tab() << "driftLength: " <<
997 driftLength[i * 2] << endl;
998 cout << TRGDebug::tab() << "LR: " << LR[i * 2] << endl;
999 }
1000 }
1001 }
1002 // Fit2D
1003 double rho, phi0, pt, chi2;
1004 rho = 0;
1005 phi0 = 0;
1006 chi2 = 0;
1007 vector<double>phi2DInvError(5);
1008 for (unsigned iAx = 0; iAx < 5; iAx++) {
1009 phi2DInvError[iAx] = 1 / phi2DError[iAx];
1010 }
1011
1012 Fitter3DUtility::rPhiFitter(&rr2D[0], &phi2D[0], &phi2DInvError[0], rho, phi0, chi2); // By JB
1013
1014 pt = 0.3 * 1.5 * rho / 100;
1015
1016 if (TRGDebug::level()) {
1017 cout << TRGDebug::tab() << "rho: " << rho << endl;
1018 cout << TRGDebug::tab() << "pt: " << pt << endl;
1019 }
1020
1021 // update track parameters
1022 TRGCDCHelix helix(ORIGIN, CLHEP::HepVector(5, 0), CLHEP::HepSymMatrix(5, 0));
1023 CLHEP::HepVector a(5);
1024 a = aTrack.helix().a();
1025 a[1] = (aTrack.charge() < 0) ? fmod(phi0 + M_PI, 2 * M_PI) : phi0;
1026 a[2] = 1 / pt * aTrack.charge();
1027 helix.a(a);
1028 aTrack.setFitted(1);
1029 aTrack.setHelix(helix);
1030 aTrack.set2DFitChi2(chi2);
1031
1032 // Change to Signals
1033 {
1035 double rhoMin = 67;
1036 double rhoMax = 1600;
1037 int phiBitSize = 12;
1038 int rhoBitSize = 12;
1039 vector<tuple<string, double, int, double, double, int> > t_values = {
1040 make_tuple("phi_0", phi2D[0], phiBitSize, -2 * M_PI, 2 * M_PI, 0),
1041 make_tuple("phi_1", phi2D[1], phiBitSize, -2 * M_PI, 2 * M_PI, 0),
1042 make_tuple("phi_2", phi2D[2], phiBitSize, -2 * M_PI, 2 * M_PI, 0),
1043 make_tuple("phi_3", phi2D[3], phiBitSize, -2 * M_PI, 2 * M_PI, 0),
1044 make_tuple("phi_4", phi2D[4], phiBitSize, -2 * M_PI, 2 * M_PI, 0),
1045 make_tuple("rho", rho, rhoBitSize, rhoMin, rhoMax, 0),
1046 make_tuple("2Drr_0", rr[0], 12, 18.8, 105.68, 0),
1047 make_tuple("2Drr_1", rr[2], 12, 18.8, 105.68, 0),
1048 make_tuple("2Drr_2", rr[4], 12, 18.8, 105.68, 0),
1049 make_tuple("2Drr_3", rr[6], 12, 18.8, 105.68, 0),
1050 make_tuple("2Drr_4", rr[8], 12, 18.8, 105.68, 0),
1051 };
1052 TRGCDCJSignal::valuesToMapSignals(t_values, _commonData, m_mSignalStorage);
1053 }
1054 TRGCDCHoughFinder::calCosPhi(m_mSignalStorage, m_mLutStorage);
1055 TRGCDCHoughFinder::calSinPhi(m_mSignalStorage, m_mLutStorage);
1056 TRGCDCHoughFinder::rPhi(m_mSignalStorage, m_mLutStorage);
1057 // Name signals only once.
1058 bool updateSecondName = false;
1059 if ((*m_mSignalStorage.begin()).second.getName() == "") {
1060 for (auto it = m_mSignalStorage.begin(); it != m_mSignalStorage.end(); it++) {
1061 (*it).second.setName((*it).first);
1062 }
1063 updateSecondName = true;
1064 }
1065
1066 if (updateSecondName) {
1067 // Print Vhdl
1068 if ((*m_mSignalStorage.begin()).second.getName() != "") {
1069 if (_commonData->getPrintedToFile() == 0) {
1070 if (_commonData->getPrintVhdl() == 0) {
1071 _commonData->setVhdlOutputFile("Fitter2D.vhd");
1072 _commonData->setPrintVhdl(1);
1073 } else {
1074 _commonData->setPrintVhdl(0);
1075 _commonData->entryVhdlCode();
1076 _commonData->signalsVhdlCode();
1077 _commonData->buffersVhdlCode();
1078 _commonData->printToFile();
1079 // Print LUTs.
1080 for (map<string, TRGCDCJLUT*>::iterator it = m_mLutStorage.begin(); it != m_mLutStorage.end(); ++it) {
1081 it->second->makeCOE(it->first + ".coe");
1082 }
1083 }
1084 }
1085 }
1086 }
1087
1088#ifdef TRGCDC_DISPLAY_HOUGH
1089 vector<const TCTrack*> cc;
1090 cc.push_back(& aTrack);
1091 const string stg = "doFitting : Kaiyu method";
1092 const string inf = " ";
1093 D->clear();
1094 D->stage(stg);
1095 D->information(inf);
1096 D->area().append(_cdc.hits());
1097 D->area().append(_cdc.segmentHits());
1098 D->area().append(cc, Gdk::Color("blue"));
1099 D->show();
1100 D->run();
1101#endif
1102 }
1103
1104 TRGDebug::leaveStage(sn);
1105
1106 return 0;
1107 }
1108
1109 TCTrack*
1110 TRGCDCHoughFinder::makeTrack(const unsigned peakId, const unsigned pm) const
1111 {
1112
1113 //...Cal. pt and phi from the cell number...
1114 const TCHTransformationCircle& tc =
1115 (TCHTransformationCircle&) _plane[pm]->transformation();
1116 unsigned x, y;
1117 _plane[pm]->id(peakId, x, y);
1118 const TRGPoint2D hp = _plane[pm]->position(x, y);
1119 const TRGPoint2D cc = tc.circleCenter(hp);
1120 const double r = cc.y();
1121 const double phi = cc.x();
1122
1123 //...Make a circle...
1124 TCCircle c(r, phi, pm ? -1. : 1., * _plane[pm]);
1125 c.name("circle_" + TRGUtil::itostring(int(peakId) * (pm ? -1 : 1)));
1126
1127 if (TRGDebug::level()) {
1128 cout << TRGDebug::tab() << "plane" << pm << ",serialId=" << peakId
1129 << endl;
1130 c.dump("detail", TRGDebug::tab() + "Circle> ");
1131 }
1132
1133 //...Get segment hits...
1134 vector<TCLink*> links;
1135 vector<const TCSegment*> segments;
1136 const unsigned nLayers = _plane[pm]->nLayers();
1137 for (unsigned j = 0; j < nLayers; j++) {
1138 const vector<unsigned>& ptn =
1139 _plane[pm]->patternId(j, peakId);
1140 for (unsigned k = 0; k < ptn.size(); k++) {
1141 const TCSegment& s = _cdc.axialSegment(j, ptn[k]);
1142 segments.push_back(& s);
1143 if (s.hit()) {
1144 TCLink* l = new TCLink(0, s.hit());
1145 links.push_back(l);
1146 }
1147 }
1148 }
1149 c.append(links);
1150
1151 if (TRGDebug::level()) {
1152 cout << TRGDebug::tab() << "attched segments below" << endl;
1153 cout << TRGDebug::tab(4);
1154 for (unsigned j = 0; j < segments.size(); j++) {
1155 cout << segments[j]->name();
1156 if (j != (segments.size() - 1))
1157 cout << ",";
1158 }
1159 cout << endl;
1160 }
1161
1162 //...Make a track...
1163 TCTrack* t = new TCTrack(c);
1164 t->name("track_" + TRGUtil::itostring(int(peakId) * (pm ? -1 : 1)));
1165
1166 if (TRGDebug::level()) {
1167 t->relation().dump("", TRGDebug::tab());
1168 t->dump("detail");
1169 }
1170
1171 return t;
1172 }
1173
1175} // namespace Belle2
Belle2::CDC::CDCGeometryPar
The Class for CDC Geometry Parameters.
Definition: CDCGeometryPar.h:51
Belle2::CDC::CDCGeometryPar::nWiresInLayer
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
Definition: CDCGeometryPar.h:1254
Belle2::CDC::CDCGeometryPar::Instance
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
Definition: CDCGeometryPar.cc:39
Belle2::CDC::CDCGeometryPar::senseWireR
double senseWireR(int layerId) const
Returns radius of sense wire in each layer.
Definition: CDCGeometryPar.h:1274
Belle2::TRGCDCHelix
TRGCDCHelix parameter class.
Definition: Helix.h:34
Belle2::TRGCDCHoughFinder::_commonData
TRGCDCJSignalData * _commonData
For VHDL code.
Definition: HoughFinder.h:184
Belle2::TRGCDCHoughFinder::_plane
TRGCDCHoughPlaneMulti2 * _plane[2]
Hough planes, for + and - charges.
Definition: HoughFinder.h:145
Belle2::TRGCDCHoughFinder::_cdc
const TRGCDC & _cdc
CDCTRG.
Definition: HoughFinder.h:142
Belle2::TRGCDCHoughFinder::_circleH
TRGCDCHoughTransformationCircle _circleH
Circle Hough transformtion.
Definition: HoughFinder.h:148
Belle2::TRGCDCHoughFinder::m_mLutStorage
std::map< std::string, TRGCDCJLUT * > m_mLutStorage
Map to hold JLuts.
Definition: HoughFinder.h:181
Belle2::TRGCDCHoughFinder::_doit
const unsigned _doit
Doit version.
Definition: HoughFinder.h:151
Belle2::TRGCDCHoughFinder::m_mSignalStorage
std::map< std::string, TRGCDCJSignal > m_mSignalStorage
Map to hold JSignals.
Definition: HoughFinder.h:178
Belle2::TRGCDCHoughFinder::_peakFinder
TRGCDCPeakFinder _peakFinder
Peak finder.
Definition: HoughFinder.h:154
Belle2::TRGCDCHoughFinder::_peakMin
const unsigned _peakMin
Min. peak height for the peak finder.
Definition: HoughFinder.h:157
Belle2::TRGCDCHoughPlaneMulti2::patternId
const std::vector< unsigned > & patternId(unsigned cellId) const
returns pattern ID which activates specified cell.
Belle2::TRGCDCJLUT
A class to use LUTs for TRGCDC.
Definition: JLUT.h:35
Belle2::TRGCDCJSignalData
A class to hold common data for JSignals.
Definition: JSignalData.h:33
Belle2::TRGCDCJSignal
A class to use Signals for TRGCDC 3D tracker.
Definition: JSignal.h:35
Belle2::TRGCDCLayer
A class to represent a cell layer.
Definition: Layer.h:33
Belle2::TRGCDCSegmentHit
A class to represent a track segment hit in CDC.
Definition: SegmentHit.h:31
Belle2::TRGCDC
The instance of TRGCDC is a singleton.
Definition: TRGCDC.h:69
Belle2::TRGPoint2D
A class to represent a point in 2D.
Definition: Point2D.h:27
Belle2::TRGTime
A class to represent a signal timing in the trigger system.
Definition: Time.h:25
Fitter3DUtility::rPhiFitter
static void rPhiFitter(double *rr, double *phi2, double *invphierror, double &rho, double &myphi0)
A circle fitter with invPhiError without fit chi2 output.
Definition: Fitter3DUtility.cc:72
Fitter3DUtility::calPhi
static double calPhi(double wirePhi, double driftLength, double rr, int lr)
Pre 3D fitter functions. rr is in cm scale. driftLength is in cm scale.
Definition: Fitter3DUtility.cc:239
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2::TRGCDCHoughFinder::doFindingTrasan
int doFindingTrasan(std::vector< unsigned > peaks[], std::vector< TRGCDCTrack * > &trackList2D) const
do track finding. (trasan version)
Definition: HoughFinder.cc:196
Belle2::TRGCDC::segmentHits
std::vector< const TRGCDCSegmentHit * > segmentHits(void) const
returns a list of TRGCDCSegmentHit.
Definition: TRGCDC.h:996
Belle2::TRGCDCHoughFinder::calCosPhi
static void calCosPhi(std::map< std::string, Belle2::TRGCDCJSignal > &mSignalStorage, std::map< std::string, Belle2::TRGCDCJLUT * > &mLutStorage)
Calculate Cos Sin ?
Definition: HoughFinder.cc:412
Belle2::TRGCDCHoughPlaneMulti2::preparePatterns
void preparePatterns(unsigned layerId, unsigned nPatterns)
allocate memory for patterns.
Definition: HoughPlaneMulti2.h:246
Belle2::TRGDebug::tab
static std::string tab(void)
returns tab spaces.
Definition: Debug.cc:47
Belle2::TRGCDC::axialSegmentHits
std::vector< const TRGCDCSegmentHit * > axialSegmentHits(unsigned) const
returns a list of TRGCDCSegmentHit in a axial super layer N.
Definition: TRGCDC.h:1014
Belle2::TRGCDC::nSegmentLayers
unsigned nSegmentLayers(void) const
returns # of track segment layers.
Definition: TRGCDC.h:1070
Belle2::TRGCDCSegment::center
const TRGCDCWire & center(void) const
returns a center wire.
Definition: Segment.h:265
Belle2::TRGCDCJSignalData::printToFile
void printToFile()
Utilities Function to print VHDL code.
Definition: JSignalData.cc:106
Belle2::TRGCDCJSignal::valuesToMapSignals
static void valuesToMapSignals(std::vector< std::tuple< std::string, double, int, double, double, int > > const &inValues, Belle2::TRGCDCJSignalData *inCommonData, std::map< std::string, Belle2::TRGCDCJSignal > &outMap)
Values => [name, value, bitwidth, min, max, clock] Changes values to signals.
Definition: JSignal.cc:2208
Belle2::TRGCDCHoughFinder::calSinPhi
static void calSinPhi(std::map< std::string, Belle2::TRGCDCJSignal > &mSignalStorage, std::map< std::string, Belle2::TRGCDCJLUT * > &mLutStorage)
Calculate Cos Sin ?
Definition: HoughFinder.cc:440
Belle2::TRGCDCHoughFinder::calPhi
static double calPhi(TRGCDCSegmentHit const *segmentHit, double eventTime)
Utility functions.
Definition: HoughFinder.cc:400
Belle2::TRGCDCHoughFinder::makeTrack
TRGCDCTrack * makeTrack(const unsigned serialID, const unsigned pm) const
Make a track from serial ID in Hough plane.
Definition: HoughFinder.cc:1110
Belle2::TRGCDCHoughFinder::FindAndFit
int FindAndFit(std::vector< TRGCDCTrack * > &trackList2D, std::vector< TRGCDCTrack * > &trackList2DFitted)
do track finding and fitting (wrapper that can choose between different versions).
Definition: HoughFinder.cc:755
Belle2::TRGCDCHoughFinder::mappingByFile
void mappingByFile(const std::string &mappingFilePlus, const std::string &mappingFileMinus)
creates mappings by a file.
Definition: HoughFinder.cc:661
Belle2::TRGCDCHoughPlaneBase::charge
float charge(void) const
returns charge for this plane.
Definition: HoughPlaneBase.h:229
Belle2::TRGCDCJSignalData::getPrintVhdl
bool getPrintVhdl() const
Gets the status of m_printVhdl.
Definition: JSignalData.cc:75
Belle2::TRGCDCJSignal::ifElse
static void ifElse(std::vector< std::pair< TRGCDCJSignal, std::vector< std::pair< TRGCDCJSignal *, TRGCDCJSignal > > > > &data, int targetClock)
If else implementation with target clock.
Definition: JSignal.cc:800
Belle2::TRGCDCCell::layerId
unsigned layerId(void) const
returns layer id.
Definition: Cell.h:214
Belle2::TRGCDCHoughFinder::selectBestHits
std::vector< TRGCDCLink * > selectBestHits(const std::vector< TRGCDCLink * > &links) const
selects the best(fastest) hits in each super layer.
Definition: HoughFinder.cc:270
Belle2::TRGCDCHoughPlaneBase::id
void id(unsigned serialId, unsigned &x, unsigned &y) const
returns x and y for serialID.
Definition: HoughPlaneBase.h:483
Belle2::TRGCDCHoughFinder::~TRGCDCHoughFinder
virtual ~TRGCDCHoughFinder()
Destructor.
Definition: HoughFinder.cc:184
Belle2::TRGCDCLUT::getValue
int getValue(unsigned) const
get LUT Values
Definition: LUT.cc:47
Belle2::TRGCDC::segmentLayer
const TRGCDCLayer * segmentLayer(unsigned id) const
returns a pointer to a track segment layer.
Definition: TRGCDC.h:1061
Belle2::TRGCDCJSignalData::getPrintedToFile
bool getPrintedToFile() const
Gets the status of m_printedToFile.
Definition: JSignalData.cc:80
Belle2::TRGCDCHoughPlaneBase::serialId
unsigned serialId(unsigned x, unsigned y) const
returns serial ID for position (x, y).
Definition: HoughPlaneBase.h:342
Belle2::TRGDebug::enterStage
static void enterStage(const std::string &stageName)
Declare that you enter new stage.
Definition: Debug.cc:24
Belle2::TRGCDCHoughPlaneMulti2::setEntry
unsigned setEntry(unsigned serialId, unsigned layerId, unsigned n)
Sets entry.
Definition: HoughPlaneMulti2.h:260
Belle2::TRGCDCHoughFinder::terminate
void terminate()
termination.
Definition: HoughFinder.cc:395
Belle2::TRGCDCPeakFinder::findPeaks
void findPeaks(const TRGCDCHoughPlaneMulti2 &hp, const unsigned threshold, std::vector< std::vector< unsigned > > &peaks) const
do peak finding. Kaiyu's version using p1p2Methode.
Definition: PeakFinder.cc:477
Belle2::TRGCDC::nAxialSuperLayers
unsigned nAxialSuperLayers(void) const
return # of axial super layers.
Definition: TRGCDC.h:891
Belle2::TRGCDCHelix::a
const CLHEP::HepVector & a(void) const
returns helix parameters.
Definition: Helix.h:313
Belle2::TRGCDCHoughFinder::doFindingAndFitting
int doFindingAndFitting(std::vector< TRGCDCTrack * > &trackList2D, std::vector< TRGCDCTrack * > &trackList2DFitted)
do track finding and fitting (Kaiyu version).
Definition: HoughFinder.cc:783
Belle2::TRGCDC::nSuperLayers
unsigned nSuperLayers(void) const
returns # of super layers.
Definition: TRGCDC.h:870
Belle2::TRGCDCSegment::lutPattern
unsigned lutPattern(void) const
hit pattern containing bit for priority position
Definition: Segment.cc:571
Belle2::TRGCDC::segment
const TRGCDCSegment & segment(unsigned id) const
returns a track segment.
Definition: TRGCDC.h:933
Belle2::TRGCDCHoughFinder::doFindingAndFittingTrasan
int doFindingAndFittingTrasan(std::vector< TRGCDCTrack * > &trackList2D, std::vector< TRGCDCTrack * > &trackList2DFitted)
do track finding and fitting (Trasan version).
Definition: HoughFinder.cc:766
Belle2::TRGCDCHoughPlaneMulti2::merge
void merge(void)
Merge layers into one.
Definition: HoughPlaneMulti2.cc:65
Belle2::TRGCDCPeakFinder::findPeaksTrasan
void findPeaksTrasan(TRGCDCHoughPlaneMulti2 &hp, const unsigned threshold, const bool centerIsPeak, std::vector< unsigned > &peakSerialIds) const
do peak finding. This is a copy from "trasan".
Definition: PeakFinder.cc:359
Belle2::TRGCDCSegment::LUT
const TRGCDCLUT * LUT(void) const
returns LUT
Definition: Segment.h:258
Belle2::TRGCDCHoughPlaneMulti2::nLayers
unsigned nLayers(void) const
returns # of Hough Boolean layers.
Definition: HoughPlaneMulti2.h:253
Belle2::TRGCDC::getEventTime
int getEventTime(void) const
returns bad hits(finding invalid hits).
Definition: TRGCDC.cc:2743
Belle2::TRGCDCHoughFinder::doFittingTrasan
int doFittingTrasan(std::vector< unsigned > peaks[], std::vector< TRGCDCTrack * > &trackList2DFitted) const
do track fitting. (old trasan version)
Definition: HoughFinder.cc:307
Belle2::TRGCDCHoughPlaneMulti2::clear
void clear(void) override
Clears all entries and regions.
Definition: HoughPlaneMulti2.h:143
Belle2::TRGCDCHoughFinder::doFitting
int doFitting(std::vector< TRGCDCTrack * > &trackList2D, std::vector< TRGCDCTrack * > &trackList2DFitted)
do track fitting. (kaiyu original)
Definition: HoughFinder.cc:857
Belle2::TRGCDC::axialSegment
const TRGCDCSegment & axialSegment(unsigned lyrId, unsigned id) const
returns a track segment in axial layers. (lyrId is axial #)
Definition: TRGCDC.h:940
Belle2::TRGCDCHoughFinder::rPhi
static void rPhi(std::map< std::string, Belle2::TRGCDCJSignal > &mSignalStorage, std::map< std::string, Belle2::TRGCDCJLUT * > &mLutStorage)
Calculate r * phi ?
Definition: HoughFinder.cc:467
Belle2::TRGDebug::level
static int level(void)
returns the debug level.
Definition: Debug.cc:67
Belle2::TRGCDC::hits
std::vector< const TRGCDCWireHit * > hits(void) const
returns a list of TRGCDCWireHit.
Definition: TRGCDC.cc:1705
Belle2::TRGCDCSegment::priorityTime
float priorityTime(void) const
return priority time in TSHit.
Definition: Segment.cc:463
Belle2::TRGDebug::leaveStage
static void leaveStage(const std::string &stageName)
Declare that you leave a stage.
Definition: Debug.cc:34
Belle2::TRGCDCHoughFinder::doFinding
int doFinding(std::vector< std::vector< unsigned > > peaks[], std::vector< TRGCDCTrack * > &trackList2D)
do track finding. (kaiyu version)
Definition: HoughFinder.cc:800
Belle2::TRGCDCHoughPlaneMulti2::registerPattern
void registerPattern(unsigned layerId, unsigned id)
registers a pattern..
Definition: HoughPlaneMulti2.h:231
Belle2::TRGCDCHoughFinder::TRGCDCHoughFinder
TRGCDCHoughFinder(const std::string &name, const TRGCDC &, unsigned peakMin, const std::string &mappingFilePlus, const std::string &mappingFileMinus, unsigned doit)
Contructor.
Definition: HoughFinder.cc:79
Belle2::TRGCDCLayer::nCells
unsigned nCells(void) const
returns # of cells.
Definition: Layer.h:194
Belle2::TRGCDCJSignalData::entryVhdlCode
void entryVhdlCode()
Function to print entry VHDL code.
Definition: JSignalData.cc:195
Belle2::TRGCDCHoughFinder::version
std::string version(void) const
returns version.
Definition: HoughFinder.cc:74
Belle2::TRGCDCSegmentHit::segment
const TRGCDCSegment & segment(void) const
returns a pointer to a track segment.
Definition: SegmentHit.cc:78
Belle2::TRGCDCJSignalData::signalsVhdlCode
void signalsVhdlCode()
Function to print definition of signal VHDL code.
Definition: JSignalData.cc:178
Belle2::TRGCDCHoughPlaneBase::position
TRGPoint2D position(unsigned x, unsigned y) const
returns position in Hough plain for a cell (x, y)..
Definition: HoughPlaneBase.h:363
Belle2::TRGCDCJSignalData::setVhdlOutputFile
void setVhdlOutputFile(const std::string &)
Sets the filename for VHDL output.
Definition: JSignalData.cc:45
Belle2::ORIGIN
const HepGeom::Point3D< double > ORIGIN
Origin 3D point.
Belle2::TRGCDCHoughPlaneMulti2::vote
void vote(float rx, float ry, int charge, unsigned layerId, int weight=1)
Voting.
Definition: HoughPlaneMulti2.h:175
Belle2::TRGCDCJSignalData::buffersVhdlCode
void buffersVhdlCode()
Function to print buffer VHDL code.
Definition: JSignalData.cc:150
Belle2::TRGCDCJSignalData::setPrintVhdl
void setPrintVhdl(bool)
Sets if to print VHDL output.
Definition: JSignalData.cc:50
Belle2::TRGCDCCell::localId
unsigned localId(void) const
returns local id in a layer.
Definition: Cell.h:207
Belle2::TRGCDCJSignal::comp
static TRGCDCJSignal comp(TRGCDCJSignal const &lhs, const std::string &operate, TRGCDCJSignal const &rhs)
Compare two signals.
Definition: JSignal.cc:1169
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17
std
STL namespace.