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