14 #define TRG_SHORT_NAMES
15 #define TRGCDC_SHORT_NAMES
17 #define ENV_PATH "BELLE2_LOCAL_DIR"
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"
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;
64 #ifdef TRGCDC_DISPLAY_HOUGH
65 using namespace Belle2_TRGCDC;
75 TRGCDCHoughFinder::version(
void)
const
77 return string(
"TRGCDCHoughFinder 5.24");
80 TRGCDCHoughFinder::TRGCDCHoughFinder(
const string& name,
83 const string& mappingFilePlus,
84 const string& mappingFileMinus,
88 _circleH(
"CircleHough"),
90 _peakFinder(
"PeakFinder"),
97 const string fMap[2] = {mappingFilePlus, mappingFileMinus};
98 const string fName[2] = {string(
"circle hough plus"), string(
"circle hough minus")};
100 for (
unsigned f = 0; f < 2; f++) {
101 const string fn = fMap[f];
102 ifstream infile(fn.c_str(), ios::in);
104 B2FATAL(
"Cannot open Hough mapping file " << fn);
110 while (!isdigit(infile.peek())) {
111 getline(infile, ignores);
120 infile >> nX >> nY >> Ymin >> Ymax;
133 #ifdef TRGCDC_DISPLAY_HOUGH
135 H0 =
new TCDisplayHough(
"Plus");
141 H1 =
new TCDisplayHough(
"Minus");
152 unsigned axialSuperLayerId = 0;
155 const unsigned nWires = l->nCells();
157 if (! nWires)
continue;
158 if ((* l)[0]->stereo())
continue;
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();
168 _plane[0]->
vote(x, y, +1, axialSuperLayerId, 1);
172 _plane[1]->
vote(x, y, -1, axialSuperLayerId, 1);
187 #ifdef TRGCDC_DISPLAY_HOUGH
192 cout <<
"TRGCDCHoughFinder ... Hough displays deleted" << endl;
198 vector<TCTrack*>& trackList2D)
const
201 const string sn =
"Hough Finder Finding (trasan version)";
210 for (
unsigned i = 0; i < nLayers; 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());
220 #ifdef TRGCDC_DISPLAY_HOUGH
221 string stg =
"2D : Hough : Results of Peak Finding";
224 H0->information(inf);
226 H0->area().append(
_plane[0]);
229 H1->information(inf);
231 H1->area().append(
_plane[1]);
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];
247 trackList2D.push_back(t);
249 #ifdef TRGCDC_DISPLAY_HOUGH
250 vector<const TCTrack*> cc;
252 const string stg =
"doFinding : track made";
253 const string inf =
" ";
257 D->area().append(cc, Gdk::Color(
"#FF0066009900"));
273 vector<TCLink*> bests;
274 vector<TCLink*> layers[9];
275 TCLink::separate(links, 9, layers);
278 for (
unsigned i = 0; i < 9; i++) {
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]);
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) {
302 bests.push_back(best);
309 vector<TRGCDCTrack*>& trackList2DFitted)
const
312 const string sn =
"Hough Finder Fitting (trasan version)";
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];
322 vector<TCLink*> links;
323 vector<const TCSegment*> segments;
325 for (
unsigned j = 0; j < nLayers; j++) {
326 const vector<unsigned>& ptn =
328 for (
unsigned k = 0; k < ptn.size(); k++) {
330 segments.push_back(& s);
332 TCLink* l =
new TCLink(0, s.hit());
344 c.name(
"CircleFitted_" + TRGUtil::itostring(nCircles));
349 <<
"plane" << pm <<
",serialId=" << peakId << endl;
352 for (
unsigned j = 0; j < segments.size(); j++) {
353 cout << segments[j]->name();
354 if (j != (segments.size() - 1))
364 TCTrack& t = *
new TCTrack(c);
365 t.name(
"Track_" + TRGUtil::itostring(i));
366 trackList2DFitted.push_back(& t);
373 #ifdef TRGCDC_DISPLAY_HOUGH
374 vector<const TCCircle*> cc;
376 const string stg =
"doFitting : trasan method";
377 const string inf =
" ";
383 D->area().append(cc, Gdk::Color(
"blue"));
414 std::map<std::string, Belle2::TRGCDCJLUT* >& mLutStorage)
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);
424 for (
unsigned iSt = 0; iSt < 5; iSt++) {
425 string t_inputName =
"phi_" + to_string(iSt);
426 string t_outputName =
"cosPhi_" + to_string(iSt);
428 if (mLutStorage.find(t_outputName) == mLutStorage.end()) {
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(),
437 mLutStorage[t_outputName]->operate(mSignalStorage[t_inputName], mSignalStorage[t_outputName]);
442 std::map<std::string, Belle2::TRGCDCJLUT* >& mLutStorage)
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);
450 for (
unsigned iSt = 0; iSt < 5; iSt++) {
451 string t_sininputName =
"phi_" + to_string(iSt);
452 string t_sinoutputName =
"sinPhi_" + to_string(iSt);
454 if (mLutStorage.find(t_sinoutputName) == mLutStorage.end()) {
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(),
464 mLutStorage[t_sinoutputName]->operate(mSignalStorage[t_sininputName], mSignalStorage[t_sinoutputName]);
469 std::map<std::string, Belle2::TRGCDCJLUT* >& mLutStorage)
472 for (
unsigned iSt = 0; iSt < 5; iSt++) {
473 mSignalStorage[
"cosPhiMul_" + to_string(iSt)] <= mSignalStorage[
"cosPhi_" + to_string(iSt)] * mSignalStorage[
"cosPhi_" + to_string(
477 for (
unsigned iSt = 0; iSt < 5; iSt++) {
478 mSignalStorage[
"sinPhiMul_" + to_string(iSt)] <= mSignalStorage[
"sinPhi_" + to_string(iSt)] * mSignalStorage[
"sinPhi_" + to_string(
482 for (
unsigned iSt = 0; iSt < 5; iSt++) {
483 mSignalStorage[
"cossinPhiMul_" + to_string(iSt)] <= mSignalStorage[
"cosPhi_" + to_string(iSt)] * mSignalStorage[
"sinPhi_" +
487 for (
unsigned iSt = 0; iSt < 5; iSt++) {
488 mSignalStorage[
"rcosPhiMul_" + to_string(iSt)] <= mSignalStorage[
"cosPhi_" + to_string(iSt)] * mSignalStorage[
"2Drr_" + to_string(
492 for (
unsigned iSt = 0; iSt < 5; iSt++) {
493 mSignalStorage[
"rsinPhiMul_" + to_string(iSt)] <= mSignalStorage[
"sinPhi_" + to_string(iSt)] * mSignalStorage[
"2Drr_" + to_string(
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"];
538 mSignalStorage[
"denMax"] =
Belle2::TRGCDCJSignal(1.332705, mSignalStorage[
"den"].getToReal(), commonData);
544 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
548 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
549 make_pair(&mSignalStorage[
"den_c"], mSignalStorage[
"denMax"])
552 t_data.push_back(make_pair(t_compare, t_assigns));
557 make_pair(&mSignalStorage[
"den_c"], mSignalStorage[
"den"].limit(mSignalStorage[
"denMin"], mSignalStorage[
"denMax"]))
560 t_data.push_back(make_pair(t_compare, t_assigns));
565 make_pair(&mSignalStorage[
"den_c"], mSignalStorage[
"denMin"])
568 t_data.push_back(make_pair(t_compare, t_assigns));
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);
582 if (mLutStorage.find(
"iDen") == mLutStorage.end()) {
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),
592 mLutStorage[
"iDen"]->operate(mSignalStorage[
"den_c"], mSignalStorage[
"iDen"]);
594 mSignalStorage[
"hcx"] <= mSignalStorage[
"hcx_p2"] * mSignalStorage[
"iDen"];
595 mSignalStorage[
"hcy"] <= mSignalStorage[
"hcy_p2"] * mSignalStorage[
"iDen"];
597 mSignalStorage[
"hcx2"] <= mSignalStorage[
"hcx"] * mSignalStorage[
"hcx"];
598 mSignalStorage[
"hcy2"] <= mSignalStorage[
"hcy"] * mSignalStorage[
"hcy"];
599 mSignalStorage[
"sumhc"] <= mSignalStorage[
"hcx2"] + mSignalStorage[
"hcy2"];
601 mSignalStorage[
"sumhcMin"] =
Belle2::TRGCDCJSignal(4489, mSignalStorage[
"sumhc"].getToReal(), commonData);
602 mSignalStorage[
"sumhcMax"] =
Belle2::TRGCDCJSignal(2560000, mSignalStorage[
"sumhc"].getToReal(), commonData);
608 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
612 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
613 make_pair(&mSignalStorage[
"sumhc_c"], mSignalStorage[
"sumhcMax"])
616 t_data.push_back(make_pair(t_compare, t_assigns));
621 make_pair(&mSignalStorage[
"sumhc_c"], mSignalStorage[
"sumhc"].limit(mSignalStorage[
"sumhcMin"], mSignalStorage[
"sumhcMax"]))
624 t_data.push_back(make_pair(t_compare, t_assigns));
629 make_pair(&mSignalStorage[
"sumhc_c"], mSignalStorage[
"sumhcMin"])
632 t_data.push_back(make_pair(t_compare, t_assigns));
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);
645 if (mLutStorage.find(
"iRho") == mLutStorage.end()) {
647 mLutStorage[
"iRho"]->setFloatFunction(
648 [ = ](
double aValue) ->
double{
return abs(sqrt(aValue));},
649 mSignalStorage[
"sumhc_c"],
650 mSignalStorage[
"iRhoMin"], mSignalStorage[
"iRhoMax"], mSignalStorage[
"sumhc"].getToReal(),
654 mLutStorage[
"iRho"]->operate(mSignalStorage[
"sumhc_c"], mSignalStorage[
"iRho"]);
656 cout <<
"<<<iRho>>>" << endl;
657 mSignalStorage[
"iRho"].dump();
663 const string& mappingFileMinus)
665 const string sn =
"mappingByFile";
668 const string* fMap[2] = {& mappingFilePlus, & mappingFileMinus};
670 for (
unsigned f = 0; f < 2; f++) {
672 const string& fn = * fMap[f];
673 ifstream infile(fn.c_str(), ios::in);
675 cout <<
" !!! can not open file" << endl
677 <<
" Mapping aborted" << endl;
684 vector<vector<unsigned>> tsf;
688 while (!isdigit(infile.peek())) {
689 getline(infile, ignores);
692 getline(infile, ignores);
696 while (getline(infile, car)) {
699 istringstream in(car);
700 vector<unsigned> slts;
708 if (j != 0 && i != 1)
710 if (j != 1 && i != 2)
717 unsigned axialSuperLayerId = 0;
720 const unsigned nWires = l->nCells();
722 if (! nWires)
continue;
723 if ((* l)[0]->stereo())
continue;
726 for (
unsigned j = 0; j < nWires; j++) {
731 const unsigned n = x.size();
732 for (
unsigned k = 0; k < n; k++) {
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)) {
757 std::vector<TRGCDCTrack*>& trackList2DFitted)
768 std::vector<TRGCDCTrack*>& trackList2DFitted)
771 const string sn =
"HoughFinder::doFindingAndFitting (trasan version)";
774 vector<unsigned> peaks[2];
785 std::vector<TRGCDCTrack*>& trackList2DFitted)
788 const string sn =
"HoughFinder::doFindingAndFitting (kaiyu version)";
791 vector<vector<unsigned>> peaks[2];
793 doFitting(trackList2D, trackList2DFitted);
802 std::vector<TRGCDCTrack*>& trackList2D)
805 const string sn =
"HoughFinder::doFinding";
814 for (
unsigned i = 0; i < nLayers; 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());
829 for (
unsigned pm = 0; pm < 2; pm++) {
830 for (
unsigned i = 0; i < peaks[pm].size(); i++) {
834 trackList2D.push_back(t);
836 #ifdef TRGCDC_DISPLAY_HOUGH
837 vector<const TCTrack*> cc;
839 const string stg =
"doFinding : track made";
840 const string inf =
" ";
844 D->area().append(cc, Gdk::Color(
"#FF0066009900"));
859 std::vector<TRGCDCTrack*>& trackList2DFitted)
862 const string sn =
"Hough Finder Fitting";
870 bool fEvtTime =
true;
874 for (
unsigned int iTrack = 0; iTrack < trackList2D.size(); iTrack++) {
875 TCTrack& aTrack = *
new TCTrack(* trackList2D[iTrack]);
876 trackList2DFitted.push_back(& aTrack);
879 double phi02D, Trg_PI;
880 vector<double> nWires(9);
881 vector<double> rr(9);
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;
900 const vector<TCLink*>& links = aTrack.links(iSL);
901 const unsigned nSegments = links.size();
905 << nSegments << endl;
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;
914 if (nSegments != 1) {
915 if (nSegments == 0) {
919 <<
"=> Not enough TS." << endl;
923 <<
"=> multiple TS are assigned." << endl;
926 if (priorityHitTS == 0) {
930 <<
"=> There are no priority hit TS" << endl;
934 if (trackFull == 0) {
936 CLHEP::HepVector helixParameters(5);
937 helixParameters = aTrack.helix().a();
939 aTrack.setHelix(helix);
943 for (
unsigned iSL = 0; iSL < 9; iSL++) {
948 rr2D = vector<double>({rr[0], rr[2], rr[4], rr[6], rr[8]});
951 vector<double>wirePhi(9);
952 vector<double>driftLength(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];
969 phi02D = aTrack.helix().phi0();
970 if (aTrack.charge() < 0) {
971 phi02D = phi02D - Trg_PI;
972 if (phi02D < 0) phi02D = phi02D + 2 * Trg_PI;
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];
982 vector<double>phi2D(5);
983 for (
unsigned iAx = 0; iAx < 5; iAx++) {
985 driftLength[iAx * 2],
990 cout <<
TRGDebug::tab() <<
"eventTime: " << eventTime << endl;
991 for (
int i = 0 ; i < 5 ; i++) {
996 driftLength[i * 2] << endl;
1002 double rho, phi0, pt, chi2;
1007 vector<double>phi2DInvError(5);
1008 for (
unsigned iAx = 0; iAx < 5; iAx++) {
1009 phi2DInvError[iAx] = 1 / phi2DError[iAx];
1014 pt = 0.3 * 1.5 * rho / 100;
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();
1028 aTrack.setFitted(1);
1029 aTrack.setHelix(helix);
1030 aTrack.set2DFitChi2(chi2);
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),
1058 bool updateSecondName =
false;
1061 (*it).second.setName((*it).first);
1063 updateSecondName =
true;
1066 if (updateSecondName) {
1081 it->second->makeCOE(it->first +
".coe");
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 =
" ";
1095 D->information(inf);
1098 D->area().append(cc, Gdk::Color(
"blue"));
1114 const TCHTransformationCircle& tc =
1115 (TCHTransformationCircle&)
_plane[pm]->transformation();
1120 const double r = cc.y();
1121 const double phi = cc.x();
1124 TCCircle c(r, phi, pm ? -1. : 1., *
_plane[pm]);
1125 c.name(
"circle_" + TRGUtil::itostring(
int(peakId) * (pm ? -1 : 1)));
1128 cout <<
TRGDebug::tab() <<
"plane" << pm <<
",serialId=" << peakId
1134 vector<TCLink*> links;
1135 vector<const TCSegment*> segments;
1137 for (
unsigned j = 0; j < nLayers; j++) {
1138 const vector<unsigned>& ptn =
1140 for (
unsigned k = 0; k < ptn.size(); k++) {
1142 segments.push_back(& s);
1144 TCLink* l =
new TCLink(0, s.hit());
1152 cout <<
TRGDebug::tab() <<
"attched segments below" << endl;
1154 for (
unsigned j = 0; j < segments.size(); j++) {
1155 cout << segments[j]->name();
1156 if (j != (segments.size() - 1))
1163 TCTrack* t =
new TCTrack(c);
1164 t->name(
"track_" + TRGUtil::itostring(
int(peakId) * (pm ? -1 : 1)));