14 #define TRG_SHORT_NAMES
15 #define TRGCDC_SHORT_NAMES
20 #include "framework/datastore/StoreArray.h"
21 #include "framework/datastore/RelationArray.h"
22 #include "rawdata/dataobjects/RawDataBlock.h"
23 #include "rawdata/dataobjects/RawCOPPER.h"
24 #include <rawdata/dataobjects/RawTRG.h>
25 #include "cdc/geometry/CDCGeometryPar.h"
26 #include "cdc/dataobjects/CDCHit.h"
27 #include "cdc/dataobjects/CDCSimHit.h"
28 #include "mdst/dataobjects/MCParticle.h"
29 #include "trg/trg/Debug.h"
30 #include "trg/trg/Time.h"
31 #include "trg/trg/State.h"
32 #include "trg/trg/Signal.h"
33 #include "trg/trg/Channel.h"
34 #include "trg/trg/Utilities.h"
35 #include "trg/cdc/dataobjects/CDCTriggerSegmentHit.h"
36 #include "trg/cdc/dataobjects/CDCTriggerTrack.h"
37 #include "trg/trg/dataobjects/TRGTiming.h"
38 #include "trg/cdc/TRGCDC.h"
39 #include "trg/cdc/Wire.h"
40 #include "trg/cdc/Layer.h"
41 #include "trg/cdc/WireHit.h"
42 #include "trg/cdc/WireHitMC.h"
43 #include "trg/cdc/TrackMC.h"
44 #include "trg/cdc/Relation.h"
45 #include "trg/cdc/Track.h"
46 #include "trg/cdc/Segment.h"
47 #include "trg/cdc/SegmentHit.h"
48 #include "trg/cdc/LUT.h"
49 #include "trg/cdc/FrontEnd.h"
50 #include "trg/cdc/Merger.h"
51 #include "trg/cdc/TrackSegmentFinder.h"
52 #include "trg/cdc/Tracker2D.h"
53 #include "trg/cdc/PerfectFinder.h"
54 #include "trg/cdc/HoughFinder.h"
55 #include "trg/cdc/Hough3DFinder.h"
56 #include "trg/cdc/Fitter3D.h"
57 #include "trg/cdc/Fitter3DUtility.h"
58 #include "trg/cdc/Link.h"
59 #include "trg/cdc/EventTime.h"
60 #include <framework/gearbox/Const.h>
62 #define NOT_USE_SOCKETLIB
68 #define SEND_BY_WRITEV
71 #include "trg/cdc/DisplayRphi.h"
72 #include "trg/cdc/DisplayHough.h"
73 namespace Belle2_TRGCDC {
74 Belle2::TRGCDCDisplayRphi* D = 0;
76 using namespace Belle2_TRGCDC;
79 #define P3D HepGeom::Point3D<double>
90 TRGCDC::name(
void)
const
96 TRGCDC::version(
void)
const
98 return string(
"TRGCDC 5.39");
105 TRGCDC::getTRGCDC(
const string& configFile,
106 unsigned simulationMode,
107 unsigned fastSimulationMode,
108 unsigned firmwareSimulationMode,
109 int firmwareSimulationStart,
110 int firmwareSimulationStop,
112 bool perfect2DFinder,
113 bool perfect3DFinder,
114 const string& innerTSLUTFile,
115 const string& outerTSLUTFile,
116 const string& rootTRGCDCFile,
117 const string& rootFitter3DFile,
118 unsigned houghFinderPeakMin,
119 const string& houghMappingFilePlus,
120 const string& houghMappingFileMinus,
125 bool fFitter3Ds2DFit,
126 bool fFitter3Ds2DFitDrift,
135 bool fXtSimpleFitter3D,
137 int trgCDCDataInputMode,
138 const string& cdchitCollectionName)
145 if (configFile !=
"good-bye") {
146 _cdc =
new TRGCDC(configFile,
149 firmwareSimulationMode,
150 firmwareSimulationStart,
151 firmwareSimulationStop,
160 houghMappingFilePlus,
161 houghMappingFileMinus,
167 fFitter3Ds2DFitDrift,
179 cdchitCollectionName);
181 cout <<
"TRGCDC::getTRGCDC ... good-bye" << endl;
190 TRGCDC::getTRGCDC(
void)
193 cout <<
"TRGCDC::getTRGCDC !!! TRGCDC is not created yet" << endl;
198 TRGCDC::getTrackList2D(
void)
204 TRGCDC::getTrackList2DFitted(
void)
206 return _trackList2DFitted;
210 TRGCDC::getTrackList3D(
void)
215 TRGCDC::TRGCDC(
const string& configFile,
216 unsigned simulationMode,
217 unsigned fastSimulationMode,
218 unsigned firmwareSimulationMode,
219 int firmwareSimulationStart,
220 int firmwareSimulationStop,
222 bool perfect2DFinder,
223 bool perfect3DFinder,
224 const string& innerTSLUTFile,
225 const string& outerTSLUTFile,
226 const string& rootTRGCDCFile,
227 const string& rootFitter3DFile,
228 unsigned houghFinderPeakMin,
229 const string& houghMappingFilePlus,
230 const string& houghMappingFileMinus,
235 bool fFitter3Ds2DFit,
236 bool fFitter3Ds2DFitDrift,
245 bool fXtSimpleFitter3D,
247 int trgCDCDataInputMode,
248 const string& cdchitCollectionName):
250 _configFilename(configFile),
251 _simulationMode(simulationMode),
252 _fastSimulationMode(fastSimulationMode),
253 _firmwareSimulationMode(firmwareSimulationMode),
254 _firmwareSimulationStart(firmwareSimulationStart),
255 _firmwareSimulationStop(firmwareSimulationStop),
257 _makeRootFile(makeRootFile),
258 _perfect2DFinder(perfect2DFinder),
259 _perfect3DFinder(perfect3DFinder),
260 _innerTSLUTFilename(innerTSLUTFile),
261 _outerTSLUTFilename(outerTSLUTFile),
262 _rootTRGCDCFilename(rootTRGCDCFile),
263 _rootFitter3DFilename(rootFitter3DFile),
264 _fLogicLUTTSF(fLogicLUTTSF),
266 _fFitter3Dsmclr(fFitter3Dsmclr),
267 _fFitter3Ds2DFit(fFitter3Ds2DFit),
268 _fFitter3Ds2DFitDrift(fFitter3Ds2DFitDrift),
269 _inefficiency(inefficiency),
273 _fprintFirmETF(fprintFirmETF),
274 _fileHough3D(fileHough3D),
275 _finder3DMode(finder3DMode),
276 _fileFitter3D(fileFitter3D),
277 _fXtSimpleFitter3D(fXtSimpleFitter3D),
282 _clock(
"CDCTrigger system clock", 0, 125. / TdcBinWidth),
283 _clockFE(
"CDCFE TDC clock", _clock, 8),
284 _clockTDC(
"CDCTrigger TDC clock (after mergers)", _clock, 4),
285 _clockD(
"CDCTrigger data clock", _clock, 1, 4),
286 _clockUser3125(
"CDCTrigger Aurora user clock (3.125Gbps)",
288 _clockUser6250(
"CDCTrigger Aurora user clock (6.250Gbps)",
296 _trgCDCDataInputMode(trgCDCDataInputMode),
297 _cdchitCollectionName(cdchitCollectionName)
302 #ifdef TRGCDC_DISPLAY
305 Gtk::Main main_instance(argc, argv);
307 D =
new TCDisplayRphi();
310 cout <<
"TRGCDC ... GTK initialized" << endl;
358 houghMappingFilePlus,
359 houghMappingFileMinus,
371 const string& houghMappingFilePlus,
372 const string& houghMappingFileMinus,
385 unsigned lastNWires = 0;
386 int lastShifts = -1000;
394 for (
unsigned i = 0; i <
nLayers; i++) {
403 unsigned axialStereoLayerId = 0;
406 axialStereoLayerId = ia;
409 axialStereoLayerId = is;
414 << nWiresInLayer <<
" axial: " << axial <<
" nShifts: "
419 if ((lastNWires != nWiresInLayer) || (lastShifts != nShifts)) {
432 lastNWires = nWiresInLayer;
433 lastShifts = nShifts;
438 << lastNWires <<
" " << nWiresInLayer << endl;
439 cout <<
TRGDebug::tab() <<
"(lastShifts,nShifts) " << lastShifts
440 <<
" " << nShifts << endl;
442 cout <<
"ia: " << ia <<
" is: " << is <<
" ias: " << ias
443 <<
" iss: " << iss << endl;
452 const float innerRadius = swr - (fwr - swr);
453 const float outerRadius = swr + (fwr - swr);
456 cout <<
TRGDebug::tab() <<
"lyr " << i <<
", in=" << innerRadius
457 <<
", out=" << outerRadius <<
", swr=" << swr <<
", fwr"
470 /
double(nWiresInLayer),
482 for (
unsigned j = 0; j < nWiresInLayer; j++) {
493 layer->push_back(tw);
502 const unsigned nWiresInTS[2] = {15, 11};
503 const int shape[2][30] = {
554 const int layerOffset[2] = {5, 2};
565 cout <<
"TRGCDC !!! can not create TS because "
566 <<
"#layers is less than 5 in super layer " << i
572 const TCCell& ww = *(*
_superLayers[i])[layerOffset[tsType]]->front();
577 const unsigned nWiresInLayer = ww.layer().nCells();
578 for (
unsigned j = 0; j < nWiresInLayer; j++) {
580 *(TCWire*)(*(*
_superLayers[i])[layerOffset[tsType]])[j];
582 const unsigned localId = w.localId();
583 const unsigned layerId = w.layerId();
584 vector<const TCWire*> cells;
586 for (
unsigned k = 0; k < nWiresInTS[tsType]; k++) {
587 const unsigned laid =
layerId + shape[tsType][k * 2];
588 const unsigned loid =
localId + shape[tsType][k * 2 + 1];
590 const TCWire* c =
wire(laid, loid);
592 cout <<
"TRGCDC !!! no such a wire for TS : "
593 <<
"layer id=" << laid <<
", local id=" << loid
600 if (w.superLayerId()) {
619 _tsSL[i].push_back(ts);
620 layer->push_back(ts);
626 if (
_r)
delete []
_r;
633 _width[i] = M_PI * 2 / float(slayer.back()->nCells());
634 _r[i] = slayer[0]->innerRadius();
637 _r[i + 1] = slayer.back()->outerRadius();
638 _r2[i + 1] =
_r[i + 1] *
_r[i + 1];
642 const TCCell& wi = *slayer[0]->front();
643 const unsigned layerId = wi.layerId();
646 cout <<
" super layer " << i <<
" radius=" <<
_r[i]
647 <<
"(r^2=" <<
_r2[i] <<
")" << endl;
655 _pFinder =
new TCPFinder(
"Perfect2DFinder", *
this);
661 _hFinder =
new TCHFinder(
"HoughFinder",
664 houghMappingFilePlus,
665 houghMappingFileMinus,
672 map<string, bool> flags = {
693 m_tree =
new TTree(
"m_tree",
"tree");
709 m_evtTime =
new TClonesArray(
"TVectorD");
712 _tree2D =
new TTree(
"tree2D",
"2D Tracks");
713 _tracks2D =
new TClonesArray(
"TVectorD");
769 if (msg.find(
"name") != string::npos ||
770 msg.find(
"version") != string::npos ||
771 msg.find(
"detail") != string::npos ||
776 if (msg.find(
"detail") != string::npos ||
777 msg.find(
"state") != string::npos) {
784 if (msg ==
"" || msg.find(
"geometry") != string::npos) {
786 unsigned nLayer =
_layers.size();
787 cout <<
" version : " <<
version() << endl;
788 cout <<
" cdc version: " <<
versionCDC() << endl;
789 cout <<
" # of wires : " <<
_wires.size() << endl;
790 cout <<
" # of layers: " << nLayer << endl;
791 cout <<
" super layer information" << endl;
792 cout <<
" # of super layers() = "
794 cout <<
" # of Axial super layers = "
796 cout <<
" # of Stereo super layers = "
799 if (msg.find(
"superLayers") != string::npos) {
800 cout <<
" super layer detail" << endl;
801 cout <<
" id #layers (stereo type)" << endl;
804 cout <<
" " << i <<
" " << n <<
" (";
805 for (
unsigned j = 0; j < n; j++) {
812 cout <<
" layer information" << endl;
813 cout <<
" # of Axial layers = "
815 cout <<
" # of Stereo layers = "
818 if (msg.find(
"layers") != string::npos) {
819 cout <<
" layer detail" << endl;
820 cout <<
" id type sId #wires lId asId assId"
822 for (
unsigned int i = 0; i <
nLayers(); ++i) {
825 <<
" " << l.stereoType()
826 <<
" " << l.superLayerId()
828 <<
" " << l.localLayerId()
829 <<
" " << l.axialStereoLayerId()
830 <<
" " << l.axialStereoSuperLayerId()
835 if (msg.find(
"wires") != string::npos) {
836 cout <<
" wire information" << endl;
837 for (
unsigned i = 0; i <
nWires(); i++)
838 (
_wires[i])->dump(
"neighbor", tab);
843 if (msg.find(
"hits") != string::npos) {
844 cout <<
" hits : " <<
_hits.size() << endl;
845 for (
unsigned i = 0; i < (unsigned)
_hits.size(); i++)
848 if (msg.find(
"axialHits") != string::npos) {
849 cout <<
" hits : " <<
_axialHits.size() << endl;
850 for (
unsigned i = 0; i < (unsigned)
_axialHits.size(); i++)
853 if (msg.find(
"stereoHits") != string::npos) {
855 for (
unsigned i = 0; i < (unsigned)
_stereoHits.size(); i++)
858 if (msg.find(
"trgWireHits") != string::npos) {
859 const string dumpOption =
"trigger detail";
860 cout <<
" wire hits" << endl;
861 for (
unsigned i = 0; i <
nWires(); i++) {
862 const TCWire& w = *
wire(i);
863 if (w.signal().active())
867 if (msg.find(
"trgWireCentralHits") != string::npos) {
868 const string dumpOption =
"trigger detail";
869 cout <<
" wire hits" << endl;
870 for (
unsigned i = 0; i <
nSegments(); i++) {
871 const TCSegment& s =
segment(i);
872 if (s.wires()[5]->signal().active())
876 if (msg.find(
"trgTSHits") != string::npos) {
877 const string dumpOption =
"trigger detail";
878 cout <<
" TS hits" << endl;
879 for (
unsigned i = 0; i <
nSegments(); i++) {
880 const TCSegment& s =
segment(i);
881 if (s.signal().active())
950 for (
unsigned i = 0; i <
_hitsMC.size(); i++)
957 for (
unsigned i = 0; i <
_wires.size(); i++) {
961 for (
unsigned i = 0; i <
_tss.size(); i++) {
962 TCSegment* s =
_tss[i];
972 for (
unsigned i = 0; i < 9; i++)
1002 cout <<
"TRGCDC !!! can not access to CDCSimHits" << endl;
1012 cout <<
"TRGCDC !!! can not access to CDCHits" << endl;
1020 if (! mcParticles) {
1022 cout <<
"TRGCDC !!! can not access to MCParticles" << endl;
1031 const unsigned nRelsMC = relsMC.
getEntries();
1034 for (
unsigned i = 0; i < nHits; i++) {
1035 const CDCHit& h = *CDCHits[i];
1036 double tmp = rand() / (double(RAND_MAX));
1044 unsigned iSimHit = 0;
1045 for (
unsigned j = 0; j < nRels; j++) {
1046 const unsigned k = rels[j].getToIndices().size();
1047 for (
unsigned l = 0; l < k; l++) {
1048 if (rels[j].getToIndex(l) == i)
1049 iSimHit = rels[j].getFromIndex();
1054 cout <<
"TRGCDC::update !!! CDCSimHit[" << iSimHit
1055 <<
"] has multiple CDCHit(" << k <<
" hits)" << endl;
1060 unsigned iMCPart = 0;
1061 for (
unsigned j = 0; j < nRelsMC; j++) {
1062 const unsigned k = relsMC[j].getToIndices().size();
1063 for (
unsigned l = 0; l < k; l++) {
1064 if (relsMC[j].getToIndex(l) == i) {
1065 iMCPart = relsMC[j].getFromIndex();
1072 cout <<
"TRGCDC::update !!! MCParticle[" << iMCPart
1073 <<
"] has multiple CDCHit(" << k <<
" hits)" << endl;
1078 if (h.getISuperLayer() == 0) t_layerId = h.getILayer();
1079 else t_layerId = h.getILayer() + 6 * h.getISuperLayer() + 2;
1080 const unsigned layerId = t_layerId;
1081 const unsigned wireId = h.getIWire();
1087 "tdc count:" << h.getTDCCount());
1089 - h.getTDCCount() + 0.5);
1093 const float driftLengthError = 0.013;
1105 w._signal.name(w.name());
1108 const int LRflag = SimHits[iSimHit]->getPosFlag();
1111 TCWHit* hit =
new TCWHit(w,
1121 hit->state(CellHitFindingValid | CellHitFittingValid);
1126 _hits.push_back(hit);
1134 << h.getTDCCount() << std::endl;
1145 cout <<
TRGDebug::tab() <<
"#CDCSimHit=" << n <<
",#CDCHit="
1154 for (
unsigned i = 0; i <
_hits.size(); i++) {
1155 const TCWHit& h = *
_hits[i];
1162 <<
" hits of a wire" << endl;
1163 for (
unsigned i = 0; i < n; i++) {
1164 const TCWHit& h = *
_hits[i];
1187 int nCDCWindows = 28 + 1;
1196 unsigned nTRGWindows = 48;
1197 int nTrgBitsInWindow = 352;
1199 vector<string> trgInformations;
1203 vector<vector<unsigned> > inBinaryData;
1205 if (inputMode == 1) {
1207 int nCdcBitsInWindow = 1536;
1220 int* temp_buf = raw_datablkarray[ iBlock ]->GetBuffer(iCopper);
1221 int nwords = raw_datablkarray[ iBlock ]->GetBlockNwords(iCopper);
1222 int malloc_flag = 0;
1227 raw_copper_buf.
SetBuffer(temp_buf, nwords, malloc_flag, num_nodes, num_events);
1228 raw_copper = &raw_copper_buf;
1254 int widthCDCFEData = 3 - 1 + nCdcBitsInWindow / 32 * nCDCWindows;
1255 int widthTRGData = nTrgBitsInWindow / 32 * nTRGWindows;
1260 vector<vector<unsigned> > rawBinaryData;
1261 for (
int iFe = 0; iFe < 2; iFe++) {
1262 int startCDCData = daqHeader + cdcHeader + (iFe * (widthCDCFEData + cdcTrailer + cdcHeader)) + 1;
1264 vector<unsigned> t_feData(widthCDCFEData + 1);
1266 for (
int iWord = 0; iWord < widthCDCFEData; iWord++) {
1267 t_feData[iWord] = (raw_copper->
GetBuffer(iCopper))[iWord + startCDCData];
1270 t_feData[widthCDCFEData] = 0;
1271 rawBinaryData.push_back(t_feData);
1273 vector<unsigned> t_trgData(widthTRGData);
1275 int startTRGData = daqHeader + 2 * (cdcHeader + widthCDCFEData + cdcTrailer) + trgHeader + 1;
1276 for (
int iWord = 0; iWord < widthTRGData; iWord++) {
1277 t_trgData[iWord] = (raw_copper->
GetBuffer(iCopper))[iWord + startTRGData];
1279 rawBinaryData.push_back(t_trgData);
1290 for (
unsigned iBoard = 0; iBoard < 2; iBoard++) {
1293 vector<unsigned> t_feData(widthCDCFEData + 1);
1294 for (
unsigned iWord = 0; iWord < rawBinaryData[iBoard].size(); iWord++) {
1295 unsigned t_value = t_buf + (rawBinaryData[iBoard][iWord] >> 8);
1296 t_buf = rawBinaryData[iBoard][iWord] << 24;
1297 t_feData[iWord] = t_value;
1299 inBinaryData.push_back(t_feData);
1301 inBinaryData.push_back(rawBinaryData[2]);
1310 }
else if (inputMode == 2) {
1313 for (
int i = 0; i < raw_trgarray.
getEntries(); i++) {
1314 cout <<
"\n===== DataBlock(RawTRG) : Block # " << i << endl;
1315 for (
int j = 0; j < raw_trgarray[ i ]->GetNumEntries(); j++) {
1316 RawCOPPER* raw_copper = raw_trgarray[ i ];
1317 vector<vector<unsigned> > rawBinaryData;
1319 for (
int iFinesse = 0; iFinesse < 4; iFinesse++) {
1321 if (bufferNwords > 0) {
1322 printf(
"===== Detector Buffer(FINESSE # %i) 0x%x words \n", iFinesse, bufferNwords);
1325 for (
int iWord = 0; iWord < bufferNwords; iWord++) {
1328 rawBinaryData.push_back(t_copperData);
1329 inBinaryData.push_back(t_copperData);
1345 cout <<
"[ERROR} TRGCDCDataInputMode is incorrect! No simulation will be done." << endl;
1352 vector<unsigned> cdcTrgTiming(2);
1353 for (
unsigned iFe = 0; iFe < 2; iFe++) {
1354 cdcTrgTiming[iFe] = (inBinaryData[iFe][1]) >> 16;
1355 if (cdcTrgTiming[iFe] >= 32768) cout <<
"CDC trigger timing error. TDC is larger than 0x8000" << endl;
1362 vector<vector<unsigned>> hitCdcData;
1364 for (
int iWindow = 0; iWindow < nCDCWindows; iWindow++) {
1365 for (
int iFE = 0; iFE < 2; iFE++) {
1366 for (
int iWire = 0; iWire < 48; iWire++) {
1369 unsigned t_layerId = 3 * iFE + ((iWire / 8) % 3);
1370 unsigned t_wireId = 8 * (iWire / 24) + iWire % 8;
1372 t_adc = ((inBinaryData[iFE][iWire / 2 + iWindow * 48 + 3] << (16 * (iWire % 2))) >> 16);
1373 t_tdc = ((inBinaryData[iFE][iWire / 2 + iWindow * 48 + 24 + 3] << (16 * (iWire % 2))) >> 16);
1375 if (t_tdc >= 32768) {
1379 unsigned startTiming;
1380 if (
unsigned(cdcTrgTiming[iFE] / 32) * 32 > (unsigned)(cdcDelay + nCDCWindows) * 32) {
1381 startTiming = unsigned(cdcTrgTiming[iFE] / 32) * 32 - (cdcDelay + nCDCWindows) * 32;
1383 startTiming = unsigned(cdcTrgTiming[iFE] / 32) * 32 + 32768 - (cdcDelay + nCDCWindows) * 32;
1385 if (t_tdc > startTiming) t_tdc -= startTiming;
1386 else t_tdc += 32768 - startTiming;
1388 vector<unsigned> t_hitCdcData(5);
1389 t_hitCdcData[0] = t_layerId;
1390 t_hitCdcData[1] = t_wireId;
1391 t_hitCdcData[2] = t_adc;
1392 t_hitCdcData[3] = t_tdc;
1393 t_hitCdcData[4] = cdcTrgTiming[iFE];
1394 hitCdcData.push_back(t_hitCdcData);
1396 cout <<
"CDC TDC data error. TDC is smaller than 0x8000. " << hex << t_tdc << dec << endl;
1409 stringstream t_trgWindow;
1410 for (
unsigned iWord = 0; iWord < inBinaryData[2].size(); iWord++) {
1411 t_trgWindow << setw(8) << setfill(
'0') << hex << inBinaryData[2][iWord];
1412 if (iWord % 11 == 10) {
1413 trgInformations.push_back(t_trgWindow.str());
1414 t_trgWindow.str(
"");
1420 vector<vector<int> > hitTrgData;
1421 for (
unsigned iWindow = 0; iWindow < nTRGWindows; iWindow++) {
1423 vector<unsigned> t_trgData(inBinaryData[2].begin() + 11 * iWindow, inBinaryData[2].begin() + 11 * iWindow + 11);
1438 for (
unsigned iWire = 0; iWire < t_hitMap.
size(); iWire++) {
1439 unsigned t_layerId = iWire / 16;
1440 unsigned t_wireId = iWire % 16;
1441 if (t_hitMap[iWire] == 1) {
1442 vector<int> t_hitTrgData(4);
1443 t_hitTrgData[0] = t_layerId;
1444 t_hitTrgData[1] = t_wireId;
1445 t_hitTrgData[2] = iWindow;
1447 if (t_layerId == 2) {
1448 t_hitTrgData[3] = (unsigned)t_priorityTiming.
subset(4 * t_wireId, 4);
1450 }
else if (t_layerId == 3) {
1456 if (t_wireId != 15) {
1457 t_priorityHits = t_priorityMap.
subset(t_wireId, 2);
1458 t_secondaryLR = t_secondHitFlag.
subset(t_wireId, 2);
1459 t_leftInformation.
set(1, t_priorityMap.
subset(t_wireId + 1, 1));
1460 t_leftInformation.
set(0, t_secondHitFlag.
subset(t_wireId + 1, 1));
1461 t_rightInformation.
set(1, t_priorityMap.
subset(t_wireId, 1));
1462 t_rightInformation.
set(0, t_secondHitFlag.
subset(t_wireId, 1));
1464 t_priorityHits.
set(0, t_priorityMap.
subset(t_wireId, 1));
1465 t_priorityHits.
set(1, 0);
1466 t_secondaryLR.
set(0, t_secondHitFlag.
subset(t_wireId, 1));
1467 t_secondaryLR.
set(1, 0);
1468 t_leftInformation.
set(1, 1);
1469 t_leftInformation.
set(0, 0);
1470 t_rightInformation.
set(1, t_priorityMap.
subset(t_wireId, 1));
1471 t_rightInformation.
set(0, t_secondHitFlag.
subset(t_wireId, 1));
1477 int secondaryFlag = -1;
1478 if ((
unsigned)t_leftInformation == (
unsigned)
TRGState(
"00", 0)) secondaryFlag = 0;
1479 if ((
unsigned)t_leftInformation == (
unsigned)
TRGState(
"01", 0)) secondaryFlag = 1;
1481 if (secondaryFlag != -1) {
1482 if (t_wireId + (1 - secondaryFlag) < 16) {
1484 unsigned t_secondaryTiming = (unsigned)t_priorityTiming.
subset(4 * (t_wireId + (1 - secondaryFlag)), 4);
1485 t_hitTrgData[3] = t_secondaryTiming;
1488 cout <<
"Error in TRGDATA for secondary priority hits" << endl;
1492 t_hitTrgData[3] = -1;
1496 t_hitTrgData[3] = -1;
1498 hitTrgData.push_back(t_hitTrgData);
1538 trgData =
new TRGSignalVector(
string(
"TRGData"), dClock, nTrgBitsInWindow);
1540 for (
unsigned iWindow = 0; iWindow < nTRGWindows; iWindow++) {
1541 vector<unsigned> t_trgData(inBinaryData[2].begin() + 11 * iWindow, inBinaryData[2].begin() + 11 * iWindow + 11);
1543 trgData->
set(t_trgState, iWindow);
1545 allTrgData->push_back(trgData);
1548 for (
unsigned i = 0; i < allTrgData->size(); i++)
delete(*allTrgData)[i];
1557 if (trgOutput == 2) {
1564 for (
unsigned iHit = 0; iHit < hitCdcData.size(); iHit++) {
1565 unsigned t_layerId = hitCdcData[iHit][0] + 50;
1566 unsigned t_wireId = hitCdcData[iHit][1];
1567 int t_rawTdc = hitCdcData[iHit][2];
1568 TCWire& currentWire = *(TCWire*)
wire(t_layerId, t_wireId);
1574 if (currentWire._signal.active()) {
1575 currentWire._signal = currentWire.signal() |
TRGSignal(rise & fall);
1576 currentWire._signal.name(currentWire.name());
1578 currentWire._signal =
TRGSignal(rise & fall);
1579 currentWire._signal.name(currentWire.name());
1583 TCWHit* hit =
new TCWHit(currentWire, 0, 0, 0, t_rawTdc, 0, t_rawTdc, 0, 0, 1);
1584 hit->state(CellHitFindingValid | CellHitFittingValid);
1585 ((TCWire*)(*
_layers[t_layerId])[t_wireId])->hit(hit);
1586 _hits.push_back(hit);
1587 if (currentWire.axial())
_axialHits.push_back(hit);
1606 unsigned n =
_hits.size();
1608 for (
unsigned i = 0; i < n; i++) {
1609 TCWHit* h =
_hits[i];
1610 const TCWire& w = h->wire();
1611 unsigned state = h->state();
1615 for (
unsigned j = 0; j < 7; j++)
neighbor[j] = w.neighbor(j);
1618 unsigned pattern = 0;
1619 for (
unsigned j = 0; j < 7; j++) {
1622 pattern += (1 << j);
1624 state |= (pattern << CellHitNeighborHit);
1627 const TCWHit* hr1 =
neighbor[2]->hit();
1628 const TCWHit* hl1 =
neighbor[3]->hit();
1629 if ((hr1 == 0) && (hl1 == 0)) {
1630 state |= CellHitIsolated;
1632 const TCWHit* hr2 =
neighbor[2]->neighbor(2)->hit();
1633 const TCWHit* hl2 =
neighbor[3]->neighbor(3)->hit();
1634 if (((hr2 == 0) && (hr1 != 0) && (hl1 == 0)) ||
1635 ((hl2 == 0) && (hl1 != 0) && (hr1 == 0)))
1636 state |= CellHitIsolated;
1641 bool previous =
false;
1643 if (
neighbor[0] == 0) previous =
true;
1657 if (previous || next) state |= CellHitContinuous;
1660 if ((pattern == 34) || (pattern == 42) ||
1661 (pattern == 40) || (pattern == 10) ||
1662 (pattern == 35) || (pattern == 50))
1663 state |= CellHitPatternRight;
1664 else if ((pattern == 17) || (pattern == 21) ||
1665 (pattern == 20) || (pattern == 5) ||
1666 (pattern == 19) || (pattern == 49))
1667 state |= CellHitPatternLeft;
1676 vector<const TCWHit*>
1679 vector<const TCWHit*> t;
1689 vector<const TCWHit*>
1692 vector<const TCWHit*> t;
1702 vector<const TCWHit*>
1705 vector<const TCWHit*> t;
1748 vector<const TCWHitMC*>
1751 vector<const TCWHitMC*> t;
1760 const TCWire*
const w =
wire(wireId);
1765 return "invalid_wire(" + TRGUtil::itostring(wireId) +
")";
1767 return TRGUtil::itostring(
layerId(wireId)) + as + TRGUtil::itostring(
localId(wireId));
1773 cout <<
"TRGCDC::localId !!! this function is not tested yet"
1775 unsigned iLayer = 0;
1777 bool nextLayer =
true;
1779 unsigned nWLast = nW;
1786 cout <<
"TRGCDC::localId !!! no such a wire (id=" <<
id << endl;
1787 return TRGCDC_UNDEFINED;
1793 cout <<
"TRGCDC::layerId !!! this function is not tested yet"
1795 unsigned iLayer = 0;
1797 bool nextLayer =
true;
1805 cout <<
"TRGCDC::layerId !!! no such a wire (id=" <<
id << endl;
1806 return TRGCDC_UNDEFINED;
1812 cout <<
"TRGCDC::layerId !!! this function is not implemented yet"
1814 return TRGCDC_UNDEFINED;
1820 unsigned iLayer = 0;
1822 bool nextLayer =
true;
1824 const vector<TRGCDCLayer*>& sl = *
superLayer(iLayer);
1825 const unsigned nLayers = sl.size();
1826 for (
unsigned i = 0; i <
nLayers; i++)
1827 nW += sl[i]->nCells();
1834 cout <<
"TRGCDC::superLayerId !!! no such a wire (id=" <<
id
1836 return TRGCDC_UNDEFINED;
1842 unsigned iLayer = 0;
1844 bool nextLayer =
true;
1846 const vector<TRGCDCLayer*>& sl = *
superLayer(iLayer);
1847 const unsigned nLayers = sl.size();
1848 for (
unsigned i = 0; i <
nLayers; i++) {
1849 nW += sl[i]->nCells();
1857 cout <<
"TRGCDC::localLayerId !!! no such a wire (id=" <<
id
1859 return TRGCDC_UNDEFINED;
1870 }
else if (i <= 13) {
1872 }
else if (i <= 19) {
1874 }
else if (i <= 25) {
1876 }
else if (i <= 31) {
1879 }
else if (aors == 1) {
1882 }
else if (i <= 11) {
1884 }
else if (i <= 17) {
1886 }
else if (i <= 23) {
2018 #ifdef TRGCDC_DISPLAY
2021 cout <<
"TRGCDC ... rphi displays deleted" << endl;
2030 const int lyr0 = w0.layerId();
2031 const int lyr1 = w1.layerId();
2032 const int lyr = lyr0 - lyr1;
2034 if (abs(lyr) > 1)
return false;
2035 if (w0.superLayerId() != w1.superLayerId())
return false;
2037 for (
unsigned i = 0; i < 7; i++) {
2038 if (w0.neighbor(i)) {
2039 if (w0.neighbor(i)->id() == w1.id())
2066 #ifdef TRGCDC_DISPLAY
2067 D->beginningOfEvent();
2076 trackSegmentClockSimulation,
2083 if (trackSegmentSimulationOnly) {
2088 #ifdef TRGCDC_DISPLAY
2090 string stg =
"fast simulation results";
2094 D->information(inf);
2095 D->area().append(
hits());
2117 cout <<
"> links=" << t.links().size() << endl;
2142 int t_debugValue = 0;
2143 for (
unsigned iTrack = 0; iTrack <
_trackList3D.size(); iTrack++) {
2144 t_debugValue |=
_trackList3D[iTrack]->getDebugValue(TRGCDCTrack::EDebugValueType::any);
2155 for (
unsigned iTrack = 0; iTrack <
_trackList3D.size(); iTrack++) {
2157 if (aTrack.fitted()) {
2158 double fitPt = aTrack.pt();
2159 double fitPhi0 = aTrack.helix().phi0();
2160 int fitCharge = aTrack.charge();
2161 if (fitCharge < 0) {
2163 if (fitPhi0 < 0) fitPhi0 += 2 * M_PI;
2165 double fitZ0 = aTrack.helix().dz();
2166 double fitCot = aTrack.helix().tanl();
2167 cout <<
TRGDebug::tab() <<
"Track[" << iTrack <<
"]: charge=" << fitCharge
2168 <<
" pt=" << fitPt <<
" phi_c=" << fitPhi0
2169 <<
" z0=" << fitZ0 <<
" cot=" << fitCot << endl;
2170 const TCRelation& trackRelation = aTrack.relation();
2171 const MCParticle& trackMCParticle = trackRelation.mcParticle(0);
2172 int mcCharge = trackMCParticle.
getCharge();
2174 double mcPhi0 = 0.0;
2175 if (mcCharge > 0) mcPhi0 = trackMCParticle.
getMomentum().Phi() - M_PI / 2;
2176 if (mcCharge < 0) mcPhi0 = trackMCParticle.
getMomentum().Phi() + M_PI / 2;
2178 if (mcPhi0 < 0) mcPhi0 += 2 * M_PI;
2180 TVector3 vertex = trackMCParticle.
getVertex();
2181 TLorentzVector vector4 = trackMCParticle.
get4Vector();
2182 TVector2 helixCenter;
2183 TVector3 impactPosition;
2185 double mcZ0 = impactPosition.Z();
2187 cout <<
TRGDebug::tab() <<
"Track[" << iTrack <<
"]: mcCharge=" << mcCharge
2188 <<
" mcPt=" << mcPt <<
" mcPhi_c=" << mcPhi0
2189 <<
" mcZ0=" << mcZ0 <<
" mcCot=" << mcCot << endl;
2203 #ifdef TRGCDC_DISPLAY
2205 vector<const TCTrack*> t2;
2207 vector<const TCTrack*> t2f;
2209 vector<const TCTrack*> t3;
2212 stg =
"fast simulation results";
2213 inf =
"red:2D, blue:2DF, green:3D";
2216 D->information(inf);
2217 D->area().append(
hits());
2219 D->area().append(t2, Gdk::Color(
"red"));
2220 D->area().append(t2f, Gdk::Color(
"blue"));
2221 D->area().append(t3, Gdk::Color(
"green"));
2247 string collection2Dfitter,
2248 string collection3Dfitter)
2256 for (
unsigned its = 0; its <
_segmentHits.size(); ++its) {
2258 const CDCHit* priorityHit = cdcHits[segmentHit->iCDCHit()];
2259 const TCSegment*
segment =
static_cast<const TCSegment*
>(&segmentHit->cell());
2261 storeSegmentHits.
appendNew(*priorityHit,
2270 for (
unsigned iw = 0; iw <
segment->
wires().size(); ++iw) {
2280 for (
unsigned imc = 0; imc < mcrel.
size(); ++imc) {
2281 mcrel[imc]->addRelationTo(storeHit, mcrel.
weight(imc));
2290 for (
unsigned itr = 0; itr <
_trackList2D.size(); ++itr) {
2292 double phi0 = remainder(track2D->helix().phi0() + M_PI_2, 2. * M_PI);
2293 double omega = track2D->charge() / track2D->helix().radius();
2295 storeTracks2Dfinder.
appendNew(phi0, omega, 0.);
2297 vector<TRGCDCLink*> links = track2D->links();
2298 for (
unsigned its = 0; its < links.size(); ++its) {
2300 const vector<const CDCTriggerSegmentHit*> storeHits =
segment->
storeHits();
2301 for (
unsigned ihit = 0; ihit < storeHits.size(); ++ihit) {
2302 track->addRelationTo(storeHits[ihit]);
2310 if (!track2D->fitted())
continue;
2311 double phi0 = remainder(track2D->helix().phi0() + M_PI_2, 2. * M_PI);
2312 double omega = track2D->charge() / track2D->helix().radius();
2313 double chi2 = track2D->get2DFitChi2();
2315 storeTracks2Dfitter.
appendNew(phi0, omega, chi2);
2317 vector<TRGCDCLink*> links = track2D->links();
2318 for (
unsigned its = 0; its < links.size(); ++its) {
2320 const vector<const CDCTriggerSegmentHit*> storeHits =
segment->
storeHits();
2321 for (
unsigned ihit = 0; ihit < storeHits.size(); ++ihit) {
2322 track->addRelationTo(storeHits[ihit]);
2327 storeTracks2Dfinder[itr]->addRelationTo(track);
2331 for (
unsigned itr = 0; itr <
_trackList3D.size(); ++itr) {
2333 if (!track3D->fitted())
continue;
2334 double phi0 = remainder(track3D->helix().phi0() + M_PI_2, 2. * M_PI);
2335 double omega = 1. / track3D->helix().radius();
2336 double chi2 = track3D->get2DFitChi2();
2337 double z = track3D->helix().dz();
2338 double cot = track3D->helix().tanl();
2339 double chi3 = track3D->get3DFitChi2();
2341 storeTracks3Dfitter.
appendNew(phi0, omega, chi2, z, cot, chi3);
2343 vector<TRGCDCLink*> links = track3D->links();
2344 for (
unsigned its = 0; its < links.size(); ++its) {
2346 const vector<const CDCTriggerSegmentHit*> storeHits =
segment->
storeHits();
2347 for (
unsigned ihit = 0; ihit < storeHits.size(); ++ihit) {
2348 track->addRelationTo(storeHits[ihit]);
2353 storeTracks2Dfinder[itr]->addRelationTo(track);
2363 #ifdef TRGCDC_DISPLAY
2364 D->beginningOfEvent();
2372 const unsigned nFronts =
_fronts.size();
2373 for (
unsigned i = 0; i < nFronts; i++) {
2383 const unsigned nMergers =
_mergers.size();
2384 for (
unsigned i = 0; i < nMergers; i++) {
2389 const unsigned nTSFBoards =
_tsfboards.size();
2390 for (
unsigned i = 0; i < nTSFBoards; i++) {
2398 int t_eventTime =
_eventTime.back()->getT0();
2402 for (
unsigned iHits = 0; iHits <
_hits.size(); iHits++) {
2404 _hits[iHits]->setDriftTime((
_hits[iHits]->drift(0) - t_eventTime) / 250, 0);
2405 _hits[iHits]->setDriftTime((
_hits[iHits]->drift(1) - t_eventTime) / 250, 1);
2412 for (
unsigned i = 0; i < nTracker2Ds; i++) {
2416 #ifdef TRGCDC_DISPLAY
2419 string stg =
"Firmware simulation";
2423 D->information(inf);
2424 D->area().append(
hits());
2454 if (infile.fail()) {
2455 cout <<
"TRGCDC !!! can not open file" << endl
2463 unsigned lastSl = 0;
2464 unsigned lastMergerLocalId = 0;
2465 while (! infile.eof()) {
2466 infile.getline(b, 800);
2476 if (lid != 0) mid = lid + tid;
2477 for (
unsigned i = 0; i < 5; i++) {
2478 const string car = TRGUtil::carstring(cdr);
2479 cdr = TRGUtil::cdrstring(cdr);
2484 }
else if (car ==
"CDC") {
2490 wid = atoi(car.c_str());
2491 }
else if (i == 1) {
2492 lid = atoi(car.c_str());
2493 }
else if (i == 2) {
2494 fid = atoi(car.c_str());
2495 }
else if (i == 3) {
2496 mid = atoi(car.c_str());
2497 }
else if (i == 4) {
2498 tid = atoi(car.c_str());
2508 const unsigned sl =
_wires[wid]->superLayerId();
2510 lastMergerLocalId = 0;
2513 bool newFrontEnd =
false;
2519 const string name =
"CDCFrontEnd" + TRGUtil::itostring(fid);
2520 TCFrontEnd::boardType t = TCFrontEnd::unknown;
2523 t = TCFrontEnd::innerInside;
2525 t = TCFrontEnd::innerOutside;
2528 t = TCFrontEnd::outerInside;
2530 t = TCFrontEnd::outerOutside;
2536 f->push_back(
_wires[wid]);
2544 bool newMerger =
false;
2552 const string name =
"CDCMerger" + TRGUtil::itostring(sl) +
2553 "-" + TRGUtil::itostring(lastMergerLocalId);
2554 TCMerger::unitType mt = TCMerger::unknown;
2556 mt = TCMerger::innerType;
2558 mt = TCMerger::outerType;
2559 m =
new TCMerger(
name,
2566 ++lastMergerLocalId;
2575 const string n = f->name() + string(
"-") + m->name();
2577 f->appendOutput(ch);
2589 const string name =
"CDCTSFBoard" + TRGUtil::itostring(tid);
2590 TSFinder::boardType tt = TSFinder::unknown;
2592 tt = TSFinder::innerType;
2594 tt = TSFinder::outerType;
2595 t =
new TSFinder(*
this,
2609 const string n = m->name() + string(
"-") + t->name();
2611 m->appendOutput(chmt);
2612 t->appendInput(chmt);
2621 for (
unsigned i = 0; i < 4; i++) {
2622 const string name =
"CDC2DBoard" + TRGUtil::itostring(i);
2623 TCTracker2D* t =
new TCTracker2D(
name,
2629 for (
unsigned j = 0; j < 9; j++) {
2631 const string n = tsf->name() + string(
"-") + t->name();
2633 tsf->appendOutput(ch);
2636 const string n = t->name() + string(
"-") +
"CDC3DBoard" + TRGUtil::itostring(i);
2639 t->appendOutput(ch);
2645 for (
unsigned i = 0; i <
_tsfboards.size(); i++) {
2649 cout <<
TRGDebug::tab() <<
"Tracker2D configuration" << endl;
2650 for (
unsigned i = 0; i <
_tracker2Ds.size(); i++) {
2760 rootCDCHitInformation.Clear();
2761 for (
unsigned iHit = 0; iHit < hitWiresFromCDC.size(); iHit++) {
2762 TVectorD t_rootCDCHitInformation(5);
2763 t_rootCDCHitInformation[0] = hitWiresFromCDC[iHit][0];
2764 t_rootCDCHitInformation[1] = hitWiresFromCDC[iHit][1];
2765 t_rootCDCHitInformation[2] = hitWiresFromCDC[iHit][2];
2766 t_rootCDCHitInformation[3] = hitWiresFromCDC[iHit][3];
2767 t_rootCDCHitInformation[4] = hitWiresFromCDC[iHit][4];
2768 new(rootCDCHitInformation[iHit]) TVectorD(t_rootCDCHitInformation);
2775 rootTRGHitInformation.Clear();
2776 for (
unsigned iHit = 0; iHit < hitWiresFromTRG.size(); iHit++) {
2777 TVectorD t_rootTRGHitInformation(4);
2778 t_rootTRGHitInformation[0] = hitWiresFromTRG[iHit][0];
2779 t_rootTRGHitInformation[1] = hitWiresFromTRG[iHit][1];
2780 t_rootTRGHitInformation[2] = hitWiresFromTRG[iHit][2];
2781 t_rootTRGHitInformation[3] = hitWiresFromTRG[iHit][3];
2782 new(rootTRGHitInformation[iHit]) TVectorD(t_rootTRGHitInformation);
2789 rootTRGRawInformation.Clear();
2790 for (
unsigned iWindow = 0; iWindow < trgInformations.size(); iWindow++) {
2791 TObjString t_rootTRGRawInformation;
2792 t_rootTRGRawInformation.SetString(trgInformations[iWindow].c_str());
2793 new(rootTRGRawInformation[iWindow]) TObjString(t_rootTRGRawInformation);