Belle II Software  release-08-01-10
TRGCDC.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 //-----------------------------------------------------------------------------
10 // Description : A class to represent CDC.
11 //-----------------------------------------------------------------------------
12 
13 #define TRG_SHORT_NAMES
14 #define TRGCDC_SHORT_NAMES
15 
16 #include <math.h>
17 #include <TRandom.h>
18 #include <iostream>
19 #include <fstream>
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/TRGCDCTrack.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>
61 #include <framework/geometry/B2Vector3.h>
62 #include <TObjString.h>
63 
64 #define NOT_USE_SOCKETLIB
65 //#define NOT_SEND
66 //#define DUMMY_DATA
67 #define TIME_MONITOR
68 //#define MULTIPLE_SEND
69 //#define MEMCPY_TO_ONE_BUFFER
70 #define SEND_BY_WRITEV
71 
72 #ifdef TRGCDC_DISPLAY
73 #include "trg/cdc/DisplayRphi.h"
74 #include "trg/cdc/DisplayHough.h"
75 namespace Belle2_TRGCDC {
76  Belle2::TRGCDCDisplayRphi* D = 0;
77 }
78 using namespace Belle2_TRGCDC;
79 #endif
80 
81 #define P3D HepGeom::Point3D<double>
82 
83 using namespace std;
84 
85 namespace Belle2 {
91  string
92  TRGCDC::name(void) const
93  {
94  return "TRGCDC";
95  }
96 
97  string
98  TRGCDC::version(void) const
99  {
100  return string("TRGCDC 5.39");
101  }
102 
103  TRGCDC*
104  TRGCDC::_cdc = 0;
105 
106  TRGCDC*
107  TRGCDC::getTRGCDC(const string& configFile,
108  unsigned simulationMode,
109  unsigned fastSimulationMode,
110  unsigned firmwareSimulationMode,
111  int firmwareSimulationStart,
112  int firmwareSimulationStop,
113  bool makeRootFile,
114  bool perfect2DFinder,
115  bool perfect3DFinder,
116  const string& innerTSLUTFile,
117  const string& outerTSLUTFile,
118  const string& rootTRGCDCFile,
119  const string& rootFitter3DFile,
120  unsigned houghFinderPeakMin,
121  const string& houghMappingFilePlus,
122  const string& houghMappingFileMinus,
123  unsigned houghDoit,
124  bool fLogicLUTTSF,
125  bool fLRLUT,
126  bool fFitter3Dsmclr,
127  bool fFitter3Ds2DFit,
128  bool fFitter3Ds2DFitDrift,
129  double inefficiency,
130  bool fileTSF,
131  bool fileETF,
132  int fverETF,
133  bool fprintFirmETF,
134  bool fileHough3D,
135  int finder3DMode,
136  bool fileFitter3D,
137  bool fXtSimpleFitter3D,
138  double TdcBinWidth,
139  int trgCDCDataInputMode,
140  const string& cdchitCollectionName)
141  {
142  if (_cdc) {
143  //delete _cdc;
144  _cdc = 0;
145  }
146 
147  if (configFile != "good-bye") {
148  _cdc = new TRGCDC(configFile,
149  simulationMode,
150  fastSimulationMode,
151  firmwareSimulationMode,
152  firmwareSimulationStart,
153  firmwareSimulationStop,
154  makeRootFile,
155  perfect2DFinder,
156  perfect3DFinder,
157  innerTSLUTFile,
158  outerTSLUTFile,
159  rootTRGCDCFile,
160  rootFitter3DFile,
161  houghFinderPeakMin,
162  houghMappingFilePlus,
163  houghMappingFileMinus,
164  houghDoit,
165  fLogicLUTTSF,
166  fLRLUT,
167  fFitter3Dsmclr,
168  fFitter3Ds2DFit,
169  fFitter3Ds2DFitDrift,
170  inefficiency,
171  fileTSF,
172  fileETF,
173  fverETF,
174  fprintFirmETF,
175  fileHough3D,
176  finder3DMode,
177  fileFitter3D,
178  fXtSimpleFitter3D,
179  TdcBinWidth,
180  trgCDCDataInputMode,
181  cdchitCollectionName);
182  } else {
183  cout << "TRGCDC::getTRGCDC ... good-bye" << endl;
184  // delete _cdc;
185  _cdc = 0;
186  }
187 
188  return _cdc;
189  }
190 
191  TRGCDC*
192  TRGCDC::getTRGCDC(void)
193  {
194  if (! _cdc)
195  cout << "TRGCDC::getTRGCDC !!! TRGCDC is not created yet" << endl;
196  return _cdc;
197  }
198 
199  vector<TCTrack*>
200  TRGCDC::getTrackList2D(void)
201  {
202  return _trackList2D;
203  }
204 
205  vector<TCTrack*>
206  TRGCDC::getTrackList2DFitted(void)
207  {
208  return _trackList2DFitted;
209  }
210 
211  vector<TCTrack*>
212  TRGCDC::getTrackList3D(void)
213  {
214  return _trackList3D;
215  }
216 
217  TRGCDC::TRGCDC(const string& configFile,
218  unsigned simulationMode,
219  unsigned fastSimulationMode,
220  unsigned firmwareSimulationMode,
221  int firmwareSimulationStart,
222  int firmwareSimulationStop,
223  bool makeRootFile,
224  bool perfect2DFinder,
225  bool perfect3DFinder,
226  const string& innerTSLUTFile,
227  const string& outerTSLUTFile,
228  const string& rootTRGCDCFile,
229  const string& rootFitter3DFile,
230  unsigned houghFinderPeakMin,
231  const string& houghMappingFilePlus,
232  const string& houghMappingFileMinus,
233  unsigned houghDoit,
234  bool fLogicLUTTSF,
235  bool fLRLUT,
236  bool fFitter3Dsmclr,
237  bool fFitter3Ds2DFit,
238  bool fFitter3Ds2DFitDrift,
239  double inefficiency,
240  bool fileTSF,
241  bool fileETF,
242  int fverETF,
243  bool fprintFirmETF,
244  bool fileHough3D,
245  int finder3DMode,
246  bool fileFitter3D,
247  bool fXtSimpleFitter3D,
248  double TdcBinWidth,
249  int trgCDCDataInputMode,
250  const string& cdchitCollectionName):
251  _debugLevel(0),
252  _configFilename(configFile),
253  _simulationMode(simulationMode),
254  _fastSimulationMode(fastSimulationMode),
255  _firmwareSimulationMode(firmwareSimulationMode),
256  _firmwareSimulationStart(firmwareSimulationStart),
257  _firmwareSimulationStop(firmwareSimulationStop),
258  _returnValue(0),
259  _makeRootFile(makeRootFile),
260  _perfect2DFinder(perfect2DFinder),
261  _perfect3DFinder(perfect3DFinder),
262  _innerTSLUTFilename(innerTSLUTFile),
263  _outerTSLUTFilename(outerTSLUTFile),
264  _rootTRGCDCFilename(rootTRGCDCFile),
265  _rootFitter3DFilename(rootFitter3DFile),
266  _fLogicLUTTSF(fLogicLUTTSF),
267  _fLRLUT(fLRLUT),
268  _fFitter3Dsmclr(fFitter3Dsmclr),
269  _fFitter3Ds2DFit(fFitter3Ds2DFit),
270  _fFitter3Ds2DFitDrift(fFitter3Ds2DFitDrift),
271  _inefficiency(inefficiency),
272  _fileTSF(fileTSF),
273  _fileETF(fileETF),
274  _fverETF(fverETF),
275  _fprintFirmETF(fprintFirmETF),
276  _fileHough3D(fileHough3D),
277  _finder3DMode(finder3DMode),
278  _fileFitter3D(fileFitter3D),
279  _fXtSimpleFitter3D(fXtSimpleFitter3D),
280  _fudgeFactor(1.),
281  _width(0),
282  _r(0),
283  _r2(0),
284  _clock("CDCTrigger system clock", 0, 125. / TdcBinWidth),
285  _clockFE("CDCFE TDC clock", _clock, 8),
286  _clockTDC("CDCTrigger TDC clock (after mergers)", _clock, 4),
287  _clockD("CDCTrigger data clock", _clock, 1, 4),
288  _clockUser3125("CDCTrigger Aurora user clock (3.125Gbps)",
289  _clock, 25, 20),
290  _clockUser6250("CDCTrigger Aurora user clock (6.250Gbps)",
291  _clock, 50, 20),
292  _offset(5.3),
293  _pFinder(0),
294  _hFinder(0),
295  _h3DFinder(0),
296  _fitter3D(0),
297  _eventTime(0),
298  _trgCDCDataInputMode(trgCDCDataInputMode),
299  _cdchitCollectionName(cdchitCollectionName)
300  {
301 
302  TRGDebug::enterStage("TRGCDC constructor");
303 
304 #ifdef TRGCDC_DISPLAY
305  int argc = 0;
306  char** argv = 0;
307  Gtk::Main main_instance(argc, argv);
308  if (! D)
309  D = new TCDisplayRphi();
310  D->clear();
311  D->show();
312  cout << "TRGCDC ... GTK initialized" << endl;
313 #endif
314 
315  //...Clock summary...
316  _clock.dump();
317  _clockFE.dump();
318  _clockTDC.dump();
319  _clockD.dump();
320 
321  //...Firmware simulation time window...
323  TRGTime fall = rise;
325  _firmwareSimulationWindow = TRGSignal(rise & fall);
326  _firmwareSimulationWindow.name("Firmware simulation window in FE clock");
328 
329  //...Time window in other clocks...
334 
335  cout << "In " << _clock.name() << endl;
336  cout << " Start : "
338  << endl;
339  cout << " Stop : "
341  << endl;
342  cout << "In " << _clockTDC.name() << endl;
343  cout << " Start : "
345  << endl;
346  cout << " Stop : "
348  << endl;
349  cout << "In " << _clockD.name() << endl;
350  cout << " Start : " << _firmwareSimulationStartDataClock << endl;
351  cout << " Stop : " << _firmwareSimulationStopDataClock << endl;
352 
353  if (TRGDebug::level()) {
354  cout << "TRGCDC ... TRGCDC initializing with " << _configFilename
355  << endl
356  << "mode=0x" << hex << _simulationMode << dec << endl;
357  }
358 
359  initialize(houghFinderPeakMin,
360  houghMappingFilePlus,
361  houghMappingFileMinus,
362  houghDoit);
363 
364  if (TRGDebug::level()) {
365  cout << "TRGCDC ... TRGCDC created with " << _configFilename << endl;
366  }
367 
368  TRGDebug::leaveStage("TRGCDC constructor");
369  }
370 
371  void
372  TRGCDC::initialize(unsigned houghFinderPeakMin,
373  const string& houghMappingFilePlus,
374  const string& houghMappingFileMinus,
375  unsigned houghDoit)
376  {
377 
378  TRGDebug::enterStage("TRGCDC initialize");
379 
380  //...CDC...
382  const unsigned nLayers = m_cdcp->nWireLayers();
383 
384  //...Loop over layers...
385  int superLayerId = -1;
386  vector<TRGCDCLayer*>* superLayer;
387  unsigned lastNWires = 0;
388  int lastShifts = -1000;
389  int ia = -1;
390  int is = -1;
391  int ias = -1;
392  int iss = -1;
393  unsigned nWires = 0;
394  float fwr = 0;
395  unsigned axialStereoSuperLayerId = 0;
396  for (unsigned i = 0; i < nLayers; i++) {
397  const unsigned nWiresInLayer = m_cdcp->nWiresInLayer(i);
398 
399  //...Axial or stereo?...
400  int nShifts = m_cdcp->nShifts(i);
401  bool axial = true;
402  if (nShifts != 0)
403  axial = false;
404 
405  unsigned axialStereoLayerId = 0;
406  if (axial) {
407  ++ia;
408  axialStereoLayerId = ia;
409  } else {
410  ++is;
411  axialStereoLayerId = is;
412  }
413 
414  if (TRGDebug::level() > 1) {
415  cout << TRGDebug::tab() << "No " << i << " nWiresInLayer: "
416  << nWiresInLayer << " axial: " << axial << " nShifts: "
417  << nShifts << endl;
418  }
419 
420  //...Is this in a new super layer?...
421  if ((lastNWires != nWiresInLayer) || (lastShifts != nShifts)) {
422  ++superLayerId;
423  superLayer = new vector<TRGCDCLayer*>;
424  _superLayers.push_back(superLayer);
425  if (axial) {
426  ++ias;
428  _axialSuperLayers.push_back(superLayer);
429  } else {
430  ++iss;
432  _stereoSuperLayers.push_back(superLayer);
433  }
434  lastNWires = nWiresInLayer;
435  lastShifts = nShifts;
436  }
437 
438  if (TRGDebug::level() > 1) {
439  cout << TRGDebug::tab() << "(lastNWires,nWiresInLayer) "
440  << lastNWires << " " << nWiresInLayer << endl;
441  cout << TRGDebug::tab() << "(lastShifts,nShifts) " << lastShifts
442  << " " << nShifts << endl;
443  cout << "superLayerId: " << superLayerId << endl;
444  cout << "ia: " << ia << " is: " << is << " ias: " << ias
445  << " iss: " << iss << endl;
446  }
447 
448  //...Calculate radius...
449  const float swr = m_cdcp->senseWireR(i);
450  if (i < nLayers - 1)
451  fwr = m_cdcp->fieldWireR(i);
452  else
453  fwr = swr + (swr - m_cdcp->fieldWireR(i - 1));
454  const float innerRadius = swr - (fwr - swr);
455  const float outerRadius = swr + (fwr - swr);
456 
457  if (TRGDebug::level() > 1)
458  cout << TRGDebug::tab() << "lyr " << i << ", in=" << innerRadius
459  << ", out=" << outerRadius << ", swr=" << swr << ", fwr"
460  << fwr << endl;
461 
462  //...New layer...
463  TRGCDCLayer* layer = new TRGCDCLayer(i,
464  superLayerId,
465  _superLayers[superLayerId]->size(),
466  axialStereoLayerId,
469  nShifts,
470  M_PI * m_cdcp->senseWireR(i)
471  * m_cdcp->senseWireR(i)
472  / double(nWiresInLayer),
473  nWiresInLayer,
474  innerRadius,
475  outerRadius);
476  _layers.push_back(layer);
477  superLayer->push_back(layer);
478  if (axial)
479  _axialLayers.push_back(layer);
480  else
481  _stereoLayers.push_back(layer);
482 
483  //...Loop over all wires in a layer...
484  for (unsigned j = 0; j < nWiresInLayer; j++) {
485  const P3D fp = P3D(m_cdcp->wireForwardPosition(i, j).x(),
486  m_cdcp->wireForwardPosition(i, j).y(),
487  m_cdcp->wireForwardPosition(i, j).z());
488  const P3D bp = P3D(m_cdcp->wireBackwardPosition(i, j).x(),
489  m_cdcp->wireBackwardPosition(i, j).y(),
490  m_cdcp->wireBackwardPosition(i, j).z());
491  TCWire* tw = new TCWire(nWires++, j, *layer, fp, bp, _clockFE);
492  if (_simulationMode & 1)
493  tw->_signal.clock(_clockTDC);
494  _wires.push_back(tw);
495  layer->push_back(tw);
496  }
497  }
498 
499  //...event Time...
500  _eventTime.push_back(new TCEventTime(*this, _fileETF));
501  _eventTime.back()->initialize();
502 
503  //...Make TSF's...
504  const unsigned nWiresInTS[2] = {15, 11};
505  const int shape[2][30] = {
506  {
507  -2, 0, // relative layer id, relative wire id
508  -1, -1, // assuming layer offset 0.0, not 0.5
509  -1, 0,
510  0, -1,
511  0, 0,
512  0, 1,
513  1, -2,
514  1, -1,
515  1, 0,
516  1, 1,
517  2, -2,
518  2, -1,
519  2, 0,
520  2, 1,
521  2, 2
522  //-2, 0, // relative layer id, relative wire id
523  //-1, 0, // assuming layer offset 0.5, not 0.0
524  //-1, 1,
525  //0, -1,
526  //0, 0,
527  //0, 1,
528  //1, -1,
529  //1, -0,
530  //1, 1,
531  //1, 2,
532  //2, -2,
533  //2, -1,
534  //2, 0,
535  //2, 1,
536  //2, 2
537  },
538  {
539  -2, -1,
540  -2, 0,
541  -2, 1,
542  -1, -1,
543  -1, 0,
544  0, 0,
545  1, -1,
546  1, 0,
547  2, -1,
548  2, 0,
549  2, 1,
550  0, 0,
551  0, 0,
552  0, 0,
553  0, 0
554  }
555  };
556  const int layerOffset[2] = {5, 2};
557  unsigned id = 0;
558  unsigned idTS = 0;
559  for (unsigned i = 0; i < nSuperLayers(); i++) {
560  unsigned tsType = 0;
561  if (i)
562  tsType = 1;
563 
564 // const unsigned nLayers = _superLayers[i]->size(); //'nLayers' shadows a previous local
565 // if (nLayers < 5) {
566  if (_superLayers[i]->size() < 5) {
567  cout << "TRGCDC !!! can not create TS because "
568  << "#layers is less than 5 in super layer " << i
569  << endl;
570  continue;
571  }
572 
573  //...TS layer... w is a central wire
574  const TCCell& ww = *(*_superLayers[i])[layerOffset[tsType]]->front();
575  TRGCDCLayer* layer = new TRGCDCLayer(id++, ww);
576  _tsLayers.push_back(layer);
577 
578  //...Loop over all wires in a central wire layer...
579  const unsigned nWiresInLayer = ww.layer().nCells();
580  for (unsigned j = 0; j < nWiresInLayer; j++) {
581  const TCWire& w =
582  *(TCWire*)(*(*_superLayers[i])[layerOffset[tsType]])[j];
583 
584  const unsigned localId = w.localId();
585  const unsigned layerId = w.layerId();
586  vector<const TCWire*> cells;
587 
588  for (unsigned k = 0; k < nWiresInTS[tsType]; k++) {
589  const unsigned laid = layerId + shape[tsType][k * 2];
590  const unsigned loid = localId + shape[tsType][k * 2 + 1];
591 
592  const TCWire* c = wire(laid, loid);
593  if (! c)
594  cout << "TRGCDC !!! no such a wire for TS : "
595  << "layer id=" << laid << ", local id=" << loid
596  << endl;
597 
598  cells.push_back(c);
599  }
600 
601  TRGCDCSegment* ts;
602  if (w.superLayerId()) {
603  ts = new TRGCDCSegment(idTS++,
604  *layer,
605  w,
606  _clockD,
608  cells);
609  } else {
610  ts = new TRGCDCSegment(idTS++,
611  *layer,
612  w,
613  _clockD,
615  cells);
616  }
617  ts->initialize();
618 
619  //...Store it...
620  _tss.push_back(ts);
621  _tsSL[i].push_back(ts);
622  layer->push_back(ts);
623  }
624  }
625 
626  //...Fill caches...
627  if (_width) delete [] _width;
628  if (_r) delete [] _r;
629  if (_r2) delete [] _r2;
630  _width = new float[nSuperLayers()];
631  _r = new float[nSuperLayers() + 1];
632  _r2 = new float[nSuperLayers() + 1];
633  for (unsigned i = 0; i < nSuperLayers(); i++) {
634  const vector<TRGCDCLayer*>& slayer = *_superLayers[i];
635  _width[i] = M_PI * 2 / float(slayer.back()->nCells());
636  _r[i] = slayer[0]->innerRadius();
637  _r2[i] = _r[i] * _r[i];
638  if (i == (nSuperLayers() - 1)) {
639  _r[i + 1] = slayer.back()->outerRadius();
640  _r2[i + 1] = _r[i + 1] * _r[i + 1];
641  }
642 
643  if (TRGDebug::level() > 9) {
644  const TCCell& wi = *slayer[0]->front();
645  const unsigned layerId = wi.layerId();
646  cout << layerId << "," << m_cdcp->senseWireR(layerId) << ","
647  << m_cdcp->fieldWireR(layerId) << endl;
648  cout << " super layer " << i << " radius=" << _r[i]
649  << "(r^2=" << _r2[i] << ")" << endl;
650  }
651  }
652 
653  //...Track Segment Finder...
654  _tsFinder = new TSFinder(*this, _fileTSF, _fLogicLUTTSF);
655 
656  //...Perfect 2D Finder...
657  _pFinder = new TCPFinder("Perfect2DFinder", *this);
658 
659  //...Perfect 3D Finder...
660  _p3DFinder = 0;
661 
662  //...Hough 2D Finder...
663  _hFinder = new TCHFinder("HoughFinder",
664  *this,
665  houghFinderPeakMin,
666  houghMappingFilePlus,
667  houghMappingFileMinus,
668  houghDoit);
669 
670  //...Hough 3D Finder...
671  _h3DFinder = new TCH3DFinder(*this, _fileHough3D, _finder3DMode);
672 
673  //...3D fitter...
674  map<string, bool> flags = {
675  {"fLRLUT", _fLRLUT},
676  {"fmcLR", _fFitter3Dsmclr},
677  {"f2DFit", _fFitter3Ds2DFit},
678  {"f2DFitDrift", _fFitter3Ds2DFitDrift},
679  {"fRootFile", _fileFitter3D},
680  {"fXtSimple", _fXtSimpleFitter3D}
681  };
682  _fitter3D = new TCFitter3D("Fitter3D",
684  *this,
685  flags);
687 
688  //...For module simulation (Front-end)...
689  configure();
690 
691  //...Initialize root file...
692  if (_makeRootFile) {
693  m_file = new TFile((char*)_rootTRGCDCFilename.c_str(), "RECREATE");
694  //m_file = new TFile("TRGCDC.root", "RECREATE");
695  m_tree = new TTree("m_tree", "tree");
696  m_treeAllTracks = new TTree("m_treeAllTracks", "treeAllTracks");
697 
698  m_fitParameters = new TClonesArray("TVectorD");
699  m_mcParameters = new TClonesArray("TVectorD");
700  m_mcTrack4Vector = new TClonesArray("TLorentzVector");
701  m_mcTrackVertexVector = new TClonesArray("TVector3");
702  m_mcTrackStatus = new TClonesArray("TVectorD");
703  // m_parameters2D = new TClonesArray("TVectorD");
704 
705  m_tree->Branch("fitParameters", &m_fitParameters, 32000, 0);
706  m_tree->Branch("mcParameters", &m_mcParameters, 32000, 0);
707  m_treeAllTracks->Branch("mcTrack4Vector", &m_mcTrack4Vector, 32000, 0);
708  m_treeAllTracks->Branch("mcTrackVertexVector", &m_mcTrackVertexVector, 32000, 0);
709  m_treeAllTracks->Branch("mcTrackStatus", &m_mcTrackStatus, 32000, 0);
710 
711  m_evtTime = new TClonesArray("TVectorD");
712  m_treeAllTracks->Branch("evtTime", &m_evtTime, 32000, 0);
713 
714  _tree2D = new TTree("tree2D", "2D Tracks");
715  _tracks2D = new TClonesArray("TVectorD");
716  _tree2D->Branch("track parameters", & _tracks2D, 32000, 0);
717 
718  //...Initialize firmware ROOT input
719  //m_minCDCTdc = 9999;
720  //m_maxCDCTdc = 0;
721  //m_minTRGTdc = 9999;
722  //m_maxTRGTdc = 0;
723 
724  m_treeROOTInput = new TTree("m_treeROOTInput", "treeRootInput");
725  //m_CDCTRGTimeMatch = new TClonesArray("TVectorD");
726  m_rootCDCHitInformation = new TClonesArray("TVectorD");
727  m_rootTRGHitInformation = new TClonesArray("TVectorD");
728  m_rootTRGRawInformation = new TClonesArray("TObjString");
729  //m_treeROOTInput->Branch("CDCTRGTimeMatch", &m_CDCTRGTimeMatch,32000,0);
730  m_treeROOTInput->Branch("rootCDCHitInformation", &m_rootCDCHitInformation, 32000, 0);
731  m_treeROOTInput->Branch("rootTRGHitInformation", &m_rootTRGHitInformation, 32000, 0);
732  m_treeROOTInput->Branch("rootTRGRawInformation", &m_rootTRGRawInformation, 32000, 0);
733  }
734 
735  TRGDebug::leaveStage("TRGCDC initialize");
736  }
737 
738  void
740  {
741  TRGDebug::enterStage("TRGCDC terminate");
742 
743  if (_tsFinder) {
744  _tsFinder->terminate();
745  }
746  if (_eventTime.back()) {
747  _eventTime.back()->terminate();
748  }
749  if (_fitter3D) {
750  _fitter3D->terminate();
751  }
752  if (_hFinder) {
753  _hFinder->terminate();
754  }
755  if (_h3DFinder) {
757  }
758  if (_makeRootFile) {
759  m_file->Write();
760  m_file->Close();
761  }
762 
763  TRGDebug::leaveStage("TRGCDC terminate");
764  }
765 
766  void
767  TRGCDC::dump(const string& msg) const
768  {
769  TRGDebug::enterStage("TRGCDC dump");
770 
771  if (msg.find("name") != string::npos ||
772  msg.find("version") != string::npos ||
773  msg.find("detail") != string::npos ||
774  msg == "") {
775  cout << name() << "(CDC version=" << versionCDC() << ", "
776  << version() << ") ";
777  }
778  if (msg.find("detail") != string::npos ||
779  msg.find("state") != string::npos) {
780  cout << "Debug Level=" << _debugLevel;
781  }
782  cout << endl;
783 
784  string tab(" ");
785 
786  if (msg == "" || msg.find("geometry") != string::npos) {
787  //...Get information..."
788  unsigned nLayer = _layers.size();
789  cout << " version : " << version() << endl;
790  cout << " cdc version: " << versionCDC() << endl;
791  cout << " # of wires : " << _wires.size() << endl;
792  cout << " # of layers: " << nLayer << endl;
793  cout << " super layer information" << endl;
794  cout << " # of super layers() = "
795  << nSuperLayers() << endl;
796  cout << " # of Axial super layers = "
797  << nAxialSuperLayers() << endl;
798  cout << " # of Stereo super layers = "
799  << nStereoSuperLayers() << endl;
800 
801  if (msg.find("superLayers") != string::npos) {
802  cout << " super layer detail" << endl;
803  cout << " id #layers (stereo type)" << endl;
804  for (unsigned i = 0; i < nSuperLayers(); ++i) {
805  const unsigned n = _superLayers[i]->size();
806  cout << " " << i << " " << n << " (";
807  for (unsigned j = 0; j < n; j++) {
808  cout << (* _superLayers[i])[0]->stereoType();
809  }
810  cout << ")" << endl;
811  }
812  }
813 
814  cout << " layer information" << endl;
815  cout << " # of Axial layers = "
816  << nAxialLayers() << endl;
817  cout << " # of Stereo layers = "
818  << nStereoLayers() << endl;
819 
820  if (msg.find("layers") != string::npos) {
821  cout << " layer detail" << endl;
822  cout << " id type sId #wires lId asId assId"
823  << endl;
824  for (unsigned int i = 0; i < nLayers(); ++i) {
825  const TRGCDCLayer& l = *_layers[i];
826  cout << " " << i
827  << " " << l.stereoType()
828  << " " << l.superLayerId()
829  << " " << l.nCells()
830  << " " << l.localLayerId()
831  << " " << l.axialStereoLayerId()
832  << " " << l.axialStereoSuperLayerId()
833  << endl;
834  }
835  }
836 
837  if (msg.find("wires") != string::npos) {
838  cout << " wire information" << endl;
839  for (unsigned i = 0; i < nWires(); i++)
840  (_wires[i])->dump("neighbor", tab);
841  }
842 
843  return;
844  }
845  if (msg.find("hits") != string::npos) {
846  cout << " hits : " << _hits.size() << endl;
847  for (unsigned i = 0; i < (unsigned) _hits.size(); i++)
848  _hits[i]->dump("mc drift", tab);
849  }
850  if (msg.find("axialHits") != string::npos) {
851  cout << " hits : " << _axialHits.size() << endl;
852  for (unsigned i = 0; i < (unsigned) _axialHits.size(); i++)
853  _axialHits[i]->dump("mc drift", tab);
854  }
855  if (msg.find("stereoHits") != string::npos) {
856  cout << " hits : " << _stereoHits.size() << endl;
857  for (unsigned i = 0; i < (unsigned) _stereoHits.size(); i++)
858  _stereoHits[i]->dump("mc drift", tab);
859  }
860  if (msg.find("trgWireHits") != string::npos) {
861  const string dumpOption = "trigger detail";
862  cout << " wire hits" << endl;
863  for (unsigned i = 0; i < nWires(); i++) {
864  const TCWire& w = *wire(i);
865  if (w.signal().active())
866  w.dump(dumpOption, TRGDebug::tab(4));
867  }
868  }
869  if (msg.find("trgWireCentralHits") != string::npos) {
870  const string dumpOption = "trigger detail";
871  cout << " wire hits" << endl;
872  for (unsigned i = 0; i < nSegments(); i++) {
873  const TCSegment& s = segment(i);
874  if (s.wires()[5]->signal().active())
875  s.wires()[5]->dump(dumpOption, TRGDebug::tab(4));
876  }
877  }
878  if (msg.find("trgTSHits") != string::npos) {
879  const string dumpOption = "trigger detail";
880  cout << " TS hits" << endl;
881  for (unsigned i = 0; i < nSegments(); i++) {
882  const TCSegment& s = segment(i);
883  if (s.signal().active())
884  s.dump(dumpOption, TRGDebug::tab(4));
885  }
886  }
887 
888  TRGDebug::leaveStage("TRGCDC dump");
889  }
890 
891  const TCWire*
892  TRGCDC::wire(unsigned id) const
893  {
894  if (id < nWires())
895  return _wires[id];
896  return 0;
897  }
898 
899  const TCWire*
900  TRGCDC::wire(unsigned layerId, int localId) const
901  {
902  if (layerId < nLayers())
903  return (TCWire*) & _layers[layerId]->cell(localId);
904  return 0;
905  }
906 
907 // const TCWire *
908 // TRGCDC::wire(const HepGeom::Point3D<double> & p) const {
909 // float r = p.mag();
910 // float phi = p.phi();
911 // return wire(r, phi);
912 // }
913 
914  const TCWire*
915  TRGCDC::wire(float, float) const
916  {
917  //...Not implemented yet...
918  return _wires[0];
919 
920  // // cout << "r,phi = " << r << "," << p << endl;
921 
922  // // unsigned id = 25;
923  // // bool ok = false;
924  // // const TRGCDCLayer * l;
925  // // while (! ok) {
926  // // l = layer(id);
927  // // if (! l) return 0;
928 
929  // // const geocdc_layer * geo = l->geocdc();
930  // // if (geo->m_r + geo->m_rcsiz2 < r) ++id;
931  // // else if (geo->m_r - geo->m_rcsiz1 > r) --id;
932  // // else ok = true;
933  // // }
934  // // float dPhi = 2. * M_PI / float(l->nCells());
935  // // if (l->geocdc()->m_offset > 0.) p -= dPhi / 2.;
936  // // unsigned localId = unsigned(phi(p) / dPhi);
937  // // return l->wire(localId);
938  // }
939  }
940 
941  void
943  {
944  TRGDebug::enterStage("TRGCDC clear");
945 
946  TCWHit::removeAll();
947  TCSHit::removeAll();
948  TCLink::removeAll();
949 
950  // for (unsigned i = 0; i < _hits.size(); i++)
951  // delete _hits[i];
952  for (unsigned i = 0; i < _hitsMC.size(); i++)
953  delete _hitsMC[i];
954  // for (unsigned i = 0; i < _badHits.size(); i++)
955  // delete _badHits[i];
956  // for (unsigned i = 0; i < _segmentHits.size(); i++)
957  // delete _segmentHits[i];
958 
959  for (unsigned i = 0; i < _wires.size(); i++) {
960  TCWire* w = _wires[i];
961  w->clear();
962  }
963  for (unsigned i = 0; i < _tss.size(); i++) {
964  TCSegment* s = _tss[i];
965  s->clear();
966  }
967  _hitWires.clear();
968  _hits.clear();
969  _axialHits.clear();
970  _stereoHits.clear();
971  _badHits.clear();
972  _hitsMC.clear();
973  _segmentHits.clear();
974  for (unsigned i = 0; i < 9; i++)
975  _segmentHitsSL[i].clear();
976 
977  TRGDebug::leaveStage("TRGCDC clear");
978  }
979 
980  void
982 
983  void
985  {
986  _trackList2D.clear();
987  _trackList2DFitted.clear();
988  _trackList3D.clear();
989 
990  TRGDebug::enterStage("TRGCDC update");
991 
992  if (_trgCDCDataInputMode != 0) {
994  } else {
995 
996  //...Clear old information...
997  // fastClear();
998  clear();
999 
1000  //...CDCSimHit...
1001  StoreArray<CDCSimHit> SimHits;
1002  if (! SimHits) {
1003  if (TRGDebug::level())
1004  cout << "TRGCDC !!! can not access to CDCSimHits" << endl;
1005  TRGDebug::leaveStage("TRGCDC update");
1006  return;
1007  }
1008  const unsigned n = SimHits.getEntries();
1009 
1010  //...CDCHit...
1012  if (! CDCHits) {
1013  if (TRGDebug::level())
1014  cout << "TRGCDC !!! can not access to CDCHits" << endl;
1015  TRGDebug::leaveStage("TRGCDC update");
1016  return;
1017  }
1018  const unsigned nHits = CDCHits.getEntries();
1019 
1020  //...MCParticle...
1021  StoreArray<MCParticle> mcParticles;
1022  if (! mcParticles) {
1023  if (TRGDebug::level())
1024  cout << "TRGCDC !!! can not access to MCParticles" << endl;
1025  TRGDebug::leaveStage("TRGCDC update");
1026  return;
1027  }
1028 
1029  //...Relations...
1030  RelationArray rels(SimHits, CDCHits);
1031  const unsigned nRels = rels.getEntries();
1032  RelationArray relsMC(mcParticles, CDCHits);
1033  const unsigned nRelsMC = relsMC.getEntries();
1034 
1035  //...Loop over CDCHits...
1036  for (unsigned i = 0; i < nHits; i++) {
1037  const CDCHit& h = *CDCHits[i];
1038  double tmp = gRandom->Uniform(0.0, 1.0) / (double(RAND_MAX));
1039  if (tmp < _inefficiency)
1040  continue;
1041 
1042  // //...Check validity (skip broken channel)...
1043  // if (! (h->m_stat & CellHitFindingValid)) continue;
1044 
1045  //...Get CDCSimHit... This is expensive. Should be moved outside.
1046  unsigned iSimHit = 0;
1047  for (unsigned j = 0; j < nRels; j++) {
1048  const unsigned k = rels[j].getToIndices().size();
1049  for (unsigned l = 0; l < k; l++) {
1050  if (rels[j].getToIndex(l) == i)
1051  iSimHit = rels[j].getFromIndex();
1052  }
1053 
1054  if (TRGDebug::level())
1055  if (k > 1)
1056  cout << "TRGCDC::update !!! CDCSimHit[" << iSimHit
1057  << "] has multiple CDCHit(" << k << " hits)" << endl;
1058  }
1059 
1060  //...Get MCParticle... This is expensive, again.
1061  // (Getting the first MCParticle only)
1062  unsigned iMCPart = 0;
1063  for (unsigned j = 0; j < nRelsMC; j++) {
1064  const unsigned k = relsMC[j].getToIndices().size();
1065  for (unsigned l = 0; l < k; l++) {
1066  if (relsMC[j].getToIndex(l) == i) {
1067  iMCPart = relsMC[j].getFromIndex();
1068  break;
1069  }
1070  }
1071 
1072  if (TRGDebug::level())
1073  if (k > 1)
1074  cout << "TRGCDC::update !!! MCParticle[" << iMCPart
1075  << "] has multiple CDCHit(" << k << " hits)" << endl;
1076  }
1077 
1078  //...Wire...
1079  int t_layerId;
1080  if (h.getISuperLayer() == 0) t_layerId = h.getILayer();
1081  else t_layerId = h.getILayer() + 6 * h.getISuperLayer() + 2;
1082  const unsigned layerId = t_layerId;
1083  const unsigned wireId = h.getIWire();
1084  const TCWire& w = * static_cast<const TCWire*>(wire(layerId, wireId));
1085 
1086  //...TDC count...
1087  B2INFO("t0:" << m_cdcp->getT0(WireID(h.getID())) <<
1088  "binwidth:" << m_cdcp->getTdcBinWidth() <<
1089  "tdc count:" << h.getTDCCount());
1090  const int tdcCount = floor(m_cdcp->getT0(WireID(h.getID())) / m_cdcp->getTdcBinWidth()
1091  - h.getTDCCount() + 0.5);
1092 
1093  //...Drift length from TDC...
1094  const float driftLength = tdcCount * m_cdcp->getTdcBinWidth() * m_cdcp->getNominalDriftV();
1095  const float driftLengthError = 0.013;
1096 
1097  //...Trigger timing...
1098  TRGTime rise = TRGTime(tdcCount, true, _clockFE, w.name());
1099  if (_simulationMode & 1)
1100  rise.clock(_clockTDC);
1101  TRGTime fall = rise;
1102  fall.shift(1).reverse();
1103  if (_simulationMode & 1)
1104  w._signal |= TRGSignal(rise & fall);
1105  else
1106  w._signal |= (TRGSignal(rise & fall) & _firmwareSimulationWindow);
1107  w._signal.name(w.name());
1108 
1109  //...Left/right...
1110  const int LRflag = SimHits[iSimHit]->getPosFlag();
1111 
1112  //...TCWireHit...
1113  TCWHit* hit = new TCWHit(w,
1114  i,
1115  iSimHit,
1116  iMCPart,
1117  driftLength,
1118  driftLengthError,
1119  driftLength,
1120  driftLengthError,
1121  LRflag,
1122  1);
1123  hit->state(CellHitFindingValid | CellHitFittingValid);
1124 
1125  //...Store a hit...
1126  if (!(*_layers[layerId])[wireId]->hit())
1127  ((TCWire*)(*_layers[layerId])[wireId])->hit(hit);
1128  _hits.push_back(hit);
1129  if (w.axial()) _axialHits.push_back(hit);
1130  else _stereoHits.push_back(hit);
1131 
1132  //...Debug...
1133  if (TRGDebug::level() > 2) {
1134  w._signal.dump("", TRGDebug::tab());
1135  cout << TRGDebug::tab(4) << "CDCHit TDC count="
1136  << h.getTDCCount() << std::endl;
1137  }
1138  }
1139 
1140  //...Track segment... This part is moved to ::simulate().
1141 
1142  //...Hit classification...
1143  // _hits.sort(TCWHit::sortByWireId);
1144  classification();
1145 
1146  if (TRGDebug::level()) {
1147  cout << TRGDebug::tab() << "#CDCSimHit=" << n << ",#CDCHit="
1148  << nHits << endl;
1149  }
1150 
1151  if (TRGDebug::level() > 2) {
1152  // StoreArray<CDCSimHit> simHits("CDCSimHits"); //TODO variable is not used, is it ok?
1153  _clock.dump("detail", TRGDebug::tab());
1154  _clockFE.dump("detail", TRGDebug::tab());
1155  if (TRGDebug::level() > 10) {
1156  for (unsigned i = 0; i < _hits.size(); i++) {
1157  const TCWHit& h = *_hits[i];
1158  h.dump("detail", TRGDebug::tab(4));
1159  }
1160  } else {
1161  //unsigned nh = 10; //TODO these two lines have no impact
1162  //if (nh > _hits.size()) nh = _hits.size();
1163  cout << TRGDebug::tab() << "Dump of the first " << n
1164  << " hits of a wire" << endl;
1165  for (unsigned i = 0; i < n; i++) {
1166  const TCWHit& h = *_hits[i];
1167  h.dump("detail", TRGDebug::tab(4));
1168  }
1169  }
1170  }
1171 
1172  }
1173 
1174  TRGDebug::leaveStage("TRGCDC update");
1175  }
1176 
1177 
1178  void
1179  TRGCDC::updateByData(int inputMode)
1180  {
1181  TRGDebug::enterStage("TRGCDC updateByData");
1182 
1183  //...Clear old information...
1184  clear();
1185 
1186  // There is a 1 window difference between slow control parameter and actual data.
1187  // For SPrint-8 data
1188  int cdcDelay = 101;
1189  int nCDCWindows = 28 + 1;
1190  // For Fuji hall data
1191  //int cdcDelay = 87;
1192  //int nCDCWindows = 32+1;
1193  // For Fuji hall matching data
1194  //int cdcDelay = 79;
1195  //int nCDCWindows = 47+1;
1196 
1197  // For TRG parameters
1198  unsigned nTRGWindows = 48;
1199  int nTrgBitsInWindow = 352;
1200 
1201  vector<string> trgInformations;
1202 
1203  // ##### Store data into a common binary format #####
1204  // inBinaryData[iBoard][iWord]
1205  vector<vector<unsigned> > inBinaryData;
1206  // Get data using Belle2Link
1207  if (inputMode == 1) {
1208 
1209  int nCdcBitsInWindow = 1536;
1210  //[FIXME] Getting data from ROOT file is fixed to 2 CDC FE, 1 TRG board.
1211  // Data should be in CDCFE0, CDCFE1, TRG sequence.
1212 
1213  // For Fuji data
1214  // Get data from root file.
1215  StoreArray<RawDataBlock> raw_datablkarray;
1216  // One block is one event. Take only first block = event.
1217  // Copper board 0 is for two CDC FE and one TRG.
1218  int iBlock = 0;
1219  int iCopper = 0;
1220  //cout<<"Number of blocks: "<<raw_cdcarray.getEntries()<<endl;
1221  //cout<<"Number of Copper: "<<raw_cdcarray[iBlock]->GetNumEntries()<<endl;
1222  int* temp_buf = raw_datablkarray[ iBlock ]->GetBuffer(iCopper);
1223  int nwords = raw_datablkarray[ iBlock ]->GetBlockNwords(iCopper);
1224  int malloc_flag = 0;
1225  int num_nodes = 1;
1226  int num_events = 1;
1227  RawCOPPER* raw_copper;
1228  RawCOPPER raw_copper_buf;
1229  raw_copper_buf.SetBuffer(temp_buf, nwords, malloc_flag, num_nodes, num_events);
1230  raw_copper = &raw_copper_buf;
1231 
1233  //StoreArray<RawCDC> raw_cdcarray;
1234  //int iBlock = 0;
1235  //int iCopper = 0;
1236  //RawCOPPER* raw_copper = raw_cdcarray[iBlock];
1237 
1239  //printf("*******Start of BODY**********\n");
1240  //printf("\n%.8d : ", 0);
1241  //for (int j = 0; j < raw_copper->GetBlockNwords(iCopper); j++) {
1242  // printf("0x%.8x ", (raw_copper->GetBuffer(iCopper))[ j ]);
1243  // if ((j + 1) % 10 == 0) {
1244  // printf("\n%.8d : ", j + 1);
1245  // }
1246  //}
1247  //printf("*******End of BODY**********\n");
1248 
1249  // Number of words for data
1250  int daqHeader = 32;
1251  int cdcHeader = 5;
1252  int cdcTrailer = 2;
1253  int trgHeader = 1;
1254  //int trgTrailer = 1;
1255  // 3 for event data. -1 for lost word in DAQ. Each window is 48 words.
1256  int widthCDCFEData = 3 - 1 + nCdcBitsInWindow / 32 * nCDCWindows;
1257  int widthTRGData = nTrgBitsInWindow / 32 * nTRGWindows;
1258 
1259  // Remove headers and trailer from data.
1260  // rawBinaryData[iBoard][iWord]
1261  // CDCFE0, CDCFE1, TSF sequence.
1262  vector<vector<unsigned> > rawBinaryData;
1263  for (int iFe = 0; iFe < 2; iFe++) {
1264  int startCDCData = daqHeader + cdcHeader + (iFe * (widthCDCFEData + cdcTrailer + cdcHeader)) + 1;
1265  // 1 is due to missing data.
1266  vector<unsigned> t_feData(widthCDCFEData + 1);
1267  // Loop over all the words for CDC FE data
1268  for (int iWord = 0; iWord < widthCDCFEData; iWord++) {
1269  t_feData[iWord] = (raw_copper->GetBuffer(iCopper))[iWord + startCDCData];
1270  }
1271  // Fill with 0 for last missing data
1272  t_feData[widthCDCFEData] = 0;
1273  rawBinaryData.push_back(t_feData);
1274  }
1275  vector<unsigned> t_trgData(widthTRGData);
1276  // Loop over all the words for TRG data
1277  int startTRGData = daqHeader + 2 * (cdcHeader + widthCDCFEData + cdcTrailer) + trgHeader + 1;
1278  for (int iWord = 0; iWord < widthTRGData; iWord++) {
1279  t_trgData[iWord] = (raw_copper->GetBuffer(iCopper))[iWord + startTRGData];
1280  }
1281  rawBinaryData.push_back(t_trgData);
1283  //for(unsigned iBoard=0; iBoard<rawBinaryData.size(); iBoard++){
1284  // cout<<"Board"<<iBoard<<endl;
1285  // for(unsigned iWord=0; iWord<rawBinaryData[iBoard].size(); iWord++){
1286  // cout<<setfill('0')<<setw(8)<<hex<<rawBinaryData[iBoard][iWord];
1287  // }
1288  // cout<<dec<<endl;
1289  //}
1290 
1291  // Shift all binary data by 8 bits to the right for two CDC FE data
1292  for (unsigned iBoard = 0; iBoard < 2; iBoard++) {
1293  unsigned t_buf = 0;
1294  // 1 is due to missing data
1295  vector<unsigned> t_feData(widthCDCFEData + 1);
1296  for (unsigned iWord = 0; iWord < rawBinaryData[iBoard].size(); iWord++) {
1297  unsigned t_value = t_buf + (rawBinaryData[iBoard][iWord] >> 8);
1298  t_buf = rawBinaryData[iBoard][iWord] << 24;
1299  t_feData[iWord] = t_value;
1300  }
1301  inBinaryData.push_back(t_feData);
1302  }
1303  inBinaryData.push_back(rawBinaryData[2]);
1305  //for(unsigned iBoard=0; iBoard<inBinaryData.size(); iBoard++){
1306  // cout<<"Board"<<iBoard<<endl;
1307  // for(unsigned iWord=0; iWord<inBinaryData[iBoard].size(); iWord++){
1308  // cout<<setfill('0')<<setw(8)<<hex<<inBinaryData[iBoard][iWord];
1309  // }
1310  // cout<<dec<<endl;
1311  //}
1312  } else if (inputMode == 2) {
1313  // get data from Belle2Link
1314  StoreArray<RawTRG> raw_trgarray;
1315  for (int i = 0; i < raw_trgarray.getEntries(); i++) {
1316  cout << "\n===== DataBlock(RawTRG) : Block # " << i << endl;
1317  for (int j = 0; j < raw_trgarray[ i ]->GetNumEntries(); j++) {
1318  RawCOPPER* raw_copper = raw_trgarray[ i ];
1319  vector<vector<unsigned> > rawBinaryData;
1320  // Loop over 4 FINESSEs
1321  for (int iFinesse = 0; iFinesse < 4; iFinesse++) {
1322  int bufferNwords = raw_copper->GetDetectorNwords(j, iFinesse);
1323  if (bufferNwords > 0) {
1324  printf("===== Detector Buffer(FINESSE # %i) 0x%x words \n", iFinesse, bufferNwords);
1325  vector<unsigned> t_copperData(raw_copper->GetDetectorNwords(j, iFinesse));
1326  // Loop over all the words for board data
1327  for (int iWord = 0; iWord < bufferNwords; iWord++) {
1328  t_copperData[iWord] = (raw_copper->GetDetectorBuffer(j, iFinesse))[iWord];
1329  }
1330  rawBinaryData.push_back(t_copperData);
1331  inBinaryData.push_back(t_copperData);
1332 
1333  }
1334  }
1335  }
1336  }
1337 
1339  //for(unsigned iBoard=0; iBoard<inBinaryData.size(); iBoard++){
1340  // cout<<"Board"<<iBoard<<endl;
1341  // for(unsigned iWord=0; iWord<inBinaryData[iBoard].size(); iWord++){
1342  // cout<<setfill('0')<<setw(8)<<hex<<inBinaryData[iBoard][iWord];
1343  // }
1344  // cout<<dec<<endl;
1345  //}
1346  } else {
1347  cout << "[ERROR} TRGCDCDataInputMode is incorrect! No simulation will be done." << endl;
1348  return;
1349  }
1350 
1351  // #####Separate binary into meaningful data#####
1352  // CDC
1353  // cdcTrgTiming[iFE]
1354  vector<unsigned> cdcTrgTiming(2);
1355  for (unsigned iFe = 0; iFe < 2; iFe++) {
1356  cdcTrgTiming[iFe] = (inBinaryData[iFe][1]) >> 16;
1357  if (cdcTrgTiming[iFe] >= 32768) cout << "CDC trigger timing error. TDC is larger than 0x8000" << endl;
1358  }
1360  //for(unsigned iFe=0; ieE<2; iFe++){
1361  // cout<<"FE"<<iFe<<" TRG Timing: "<<setfill('0')<<setw(8)<<hex<<cdcTrgTiming[0]<<dec<<endl;
1362  //}
1363  // [iHit] [0]: layerId, [1]: wireId, [2]: adc, [3]: tdc, [4]: FE trg timing
1364  vector<vector<unsigned>> hitCdcData;
1365  // Store information according to window,FE,Wire
1366  for (int iWindow = 0; iWindow < nCDCWindows; iWindow++) {
1367  for (int iFE = 0; iFE < 2; iFE++) {
1368  for (int iWire = 0; iWire < 48; iWire++) {
1369  unsigned t_adc;
1370  unsigned t_tdc;
1371  unsigned t_layerId = 3 * iFE + ((iWire / 8) % 3);
1372  unsigned t_wireId = 8 * (iWire / 24) + iWire % 8;
1373  // Take 16 bits according to wire.
1374  t_adc = ((inBinaryData[iFE][iWire / 2 + iWindow * 48 + 3] << (16 * (iWire % 2))) >> 16);
1375  t_tdc = ((inBinaryData[iFE][iWire / 2 + iWindow * 48 + 24 + 3] << (16 * (iWire % 2))) >> 16);
1376  if (t_tdc != 0) {
1377  if (t_tdc >= 32768) {
1378  // Calculate TDC.
1379  t_tdc -= 32768;
1380  // Round down and go back cdcDelay+nCDCWindows windows
1381  unsigned startTiming;
1382  if (unsigned(cdcTrgTiming[iFE] / 32) * 32 > (unsigned)(cdcDelay + nCDCWindows) * 32) {
1383  startTiming = unsigned(cdcTrgTiming[iFE] / 32) * 32 - (cdcDelay + nCDCWindows) * 32;
1384  } else {
1385  startTiming = unsigned(cdcTrgTiming[iFE] / 32) * 32 + 32768 - (cdcDelay + nCDCWindows) * 32;
1386  }
1387  if (t_tdc > startTiming) t_tdc -= startTiming;
1388  else t_tdc += 32768 - startTiming;
1389  // Fill data
1390  vector<unsigned> t_hitCdcData(5);
1391  t_hitCdcData[0] = t_layerId;
1392  t_hitCdcData[1] = t_wireId;
1393  t_hitCdcData[2] = t_adc;
1394  t_hitCdcData[3] = t_tdc;
1395  t_hitCdcData[4] = cdcTrgTiming[iFE];
1396  hitCdcData.push_back(t_hitCdcData);
1397  } else {
1398  cout << "CDC TDC data error. TDC is smaller than 0x8000. " << hex << t_tdc << dec << endl;
1399  continue;
1400  }
1401  } // end tdc hit
1402  } // wire loop
1403  } // fe loop
1404  } // window loop
1406  //for(unsigned iHit=0; iHit<hitCdcData.size(); iHit++){
1407  // cout<<"["<<hitCdcData[iHit][0]<<"]["<<hitCdcData[iHit][1]<<"]: "<<hitCdcData[iHit][2]<<", "<<hitCdcData[iHit][3]<<endl;
1408  //}
1409  // TRG
1410  // Save data into string
1411  stringstream t_trgWindow;
1412  for (unsigned iWord = 0; iWord < inBinaryData[2].size(); iWord++) {
1413  t_trgWindow << setw(8) << setfill('0') << hex << inBinaryData[2][iWord];
1414  if (iWord % 11 == 10) {
1415  trgInformations.push_back(t_trgWindow.str());
1416  t_trgWindow.str("");
1417  }
1418  }
1419  // Save data into hits
1420  // [iHit] [0]: layerId, [1]: wireId, [2]: window timing, [3]: priority timing
1421  // priority timing = -1 : No timing
1422  vector<vector<int> > hitTrgData;
1423  for (unsigned iWindow = 0; iWindow < nTRGWindows; iWindow++) {
1424  // Change binary into TRGState
1425  vector<unsigned> t_trgData(inBinaryData[2].begin() + 11 * iWindow, inBinaryData[2].begin() + 11 * iWindow + 11);
1426  TRGState t_trgState(t_trgData, 1);
1427  TRGState t_mergerState = t_trgState.subset(96, 256);
1428  TRGState t_hitMap = t_mergerState.subset(0, 80);
1429  TRGState t_priorityMap = t_mergerState.subset(32, 16);
1430  TRGState t_secondPriorityMap = t_mergerState.subset(48, 16);
1431  TRGState t_secondHitFlag = t_mergerState.subset(208, 16);
1432  TRGState t_priorityTiming = t_mergerState.subset(80, 64);
1433  TRGState t_fastestTiming = t_mergerState.subset(144, 64);
1434  //cout<<"["<<iWindow<<"] HM: "<<t_hitMap<<endl;
1435  //cout<<"["<<iWindow<<"] PHM: "<<t_priorityMap<<endl;
1436  //cout<<"["<<iWindow<<"] SHM: "<<t_secondPriorityMap<<endl;
1437  //cout<<"["<<iWindow<<"] SPF: "<<t_secondHitFlag<<endl;
1438  //cout<<"["<<iWindow<<"] PT: "<<setfill('0')<<setw(16)<<hex<<(unsigned long long) (t_priorityTiming)<<dec<<endl;
1439  //cout<<"["<<iWindow<<"] FT: "<<setfill('0')<<setw(16)<<hex<<(unsigned long long) (t_fastestTiming)<<dec<<endl;
1440  for (unsigned iWire = 0; iWire < t_hitMap.size(); iWire++) {
1441  unsigned t_layerId = iWire / 16;
1442  unsigned t_wireId = iWire % 16;
1443  if (t_hitMap[iWire] == 1) {
1444  vector<int> t_hitTrgData(4);
1445  t_hitTrgData[0] = t_layerId;
1446  t_hitTrgData[1] = t_wireId;
1447  t_hitTrgData[2] = iWindow;
1448  // Calculate priority timing
1449  if (t_layerId == 2) {
1450  t_hitTrgData[3] = (unsigned)t_priorityTiming.subset(4 * t_wireId, 4);
1451  //cout<<"["<<t_layerId<<"]["<<t_wireId<<"]: "<<iWindow<<", "<<t_hitTrgData[3]<<endl;
1452  } else if (t_layerId == 3) {
1453  // [1] Priority hit, [0] secondaryLR
1454  TRGState t_leftInformation(2);
1455  TRGState t_rightInformation(2);
1456  TRGState t_priorityHits(2);
1457  TRGState t_secondaryLR(2);
1458  if (t_wireId != 15) {
1459  t_priorityHits = t_priorityMap.subset(t_wireId, 2);
1460  t_secondaryLR = t_secondHitFlag.subset(t_wireId, 2);
1461  t_leftInformation.set(1, t_priorityMap.subset(t_wireId + 1, 1));
1462  t_leftInformation.set(0, t_secondHitFlag.subset(t_wireId + 1, 1));
1463  t_rightInformation.set(1, t_priorityMap.subset(t_wireId, 1));
1464  t_rightInformation.set(0, t_secondHitFlag.subset(t_wireId, 1));
1465  } else {
1466  t_priorityHits.set(0, t_priorityMap.subset(t_wireId, 1));
1467  t_priorityHits.set(1, 0);
1468  t_secondaryLR.set(0, t_secondHitFlag.subset(t_wireId, 1));
1469  t_secondaryLR.set(1, 0);
1470  t_leftInformation.set(1, 1);
1471  t_leftInformation.set(0, 0);
1472  t_rightInformation.set(1, t_priorityMap.subset(t_wireId, 1));
1473  t_rightInformation.set(0, t_secondHitFlag.subset(t_wireId, 1));
1474  }
1475  //cout<<"["<<iWindow<<"]["<<t_layerId<<"]["<<t_wireId<<"] PriorityHits: "<<t_priorityHits<<endl;
1476  //cout<<"["<<iWindow<<"]["<<t_layerId<<"]["<<t_wireId<<"] SecondaryLR: "<<t_secondaryLR<<endl;
1477  //cout<<"["<<iWindow<<"]["<<t_layerId<<"]["<<t_wireId<<"] Combine0: "<<t_leftInformation<<endl;
1478  //cout<<"["<<iWindow<<"]["<<t_layerId<<"]["<<t_wireId<<"] Combine1: "<<t_rightInformation<<endl;
1479  int secondaryFlag = -1;
1480  if ((unsigned)t_leftInformation == (unsigned)TRGState("00", 0)) secondaryFlag = 0;
1481  if ((unsigned)t_leftInformation == (unsigned)TRGState("01", 0)) secondaryFlag = 1;
1482  //cout<<"["<<iWindow<<"]["<<t_layerId<<"]["<<t_wireId<<"] SecondaryFlag: "<<secondaryFlag<<" timing: ";
1483  if (secondaryFlag != -1) {
1484  if (t_wireId + (1 - secondaryFlag) < 16) {
1485  //int t_priorityTiming = strtol(priorityTiming.substr(15-t_wireId-(1-secondaryFlag),1).c_str(),0,16);
1486  unsigned t_secondaryTiming = (unsigned)t_priorityTiming.subset(4 * (t_wireId + (1 - secondaryFlag)), 4);
1487  t_hitTrgData[3] = t_secondaryTiming;
1488  //cout<<t_secondaryTiming<<endl;
1489  } else {
1490  cout << "Error in TRGDATA for secondary priority hits" << endl;
1491  }
1492  //cout<<"Layer3end"<<endl;
1493  } else {
1494  t_hitTrgData[3] = -1;
1495  //cout<<"-1"<<endl;
1496  }
1497  } else {
1498  t_hitTrgData[3] = -1;
1499  }
1500  hitTrgData.push_back(t_hitTrgData);
1501  } // wire is hit
1502  } // wire loop
1503  } // window loop
1505  //if(hitTrgData.size()) cout<<"TRG: "<<endl;
1506  //for(unsigned iHit=0; iHit<hitTrgData.size(); iHit++){
1507  // cout<<"["<<hitTrgData[iHit][0]<<"]["<<hitTrgData[iHit][1]<<"]: "<<hitTrgData[iHit][2]*32+hitTrgData[iHit][3]*2<<"("<<hitTrgData[iHit][2]<<","<<hitTrgData[iHit][3]<<")"<<endl;;
1508  //}
1509  //if(hitTrgData.size()) {
1510  // for(unsigned iWindow=0; iWindow<nTRGWindows; iWindow++){
1511  // vector<unsigned> t_trgData(inBinaryData[2].begin()+11*iWindow, inBinaryData[2].begin()+11*iWindow+11);
1512  // TRGState t_trgState(t_trgData,1);
1513  // TRGState t_mergerState = t_trgState.subset(96, 256);
1514  // TRGState t_hitMap = t_mergerState.subset(0,80);
1515  // TRGState t_priorityMap = t_mergerState.subset(32,16);
1516  // TRGState t_secondPriorityMap = t_mergerState.subset(48,16);
1517  // TRGState t_secondHitFlag= t_mergerState.subset(208,16);
1518  // TRGState t_priorityTiming = t_mergerState.subset(80,64);
1519  // TRGState t_fastestTiming = t_mergerState.subset(144,64);
1520  // cout<<"["<<iWindow<<"] HM: "<<t_hitMap<<endl;
1521  // cout<<"["<<iWindow<<"] PHM: "<<t_priorityMap<<endl;
1522  // cout<<"["<<iWindow<<"] SHM: "<<t_secondPriorityMap<<endl;
1523  // cout<<"["<<iWindow<<"] SPF: "<<t_secondHitFlag<<endl;
1524  // cout<<"["<<iWindow<<"] PT: "<<setfill('0')<<setw(16)<<hex<<(unsigned long long) (t_priorityTiming)<<dec<<endl;
1525  // cout<<"["<<iWindow<<"] FT: "<<setfill('0')<<setw(16)<<hex<<(unsigned long long) (t_fastestTiming)<<dec<<endl;
1526  // }
1527  //}
1528 
1529  const int trgOutput = 2;
1530  // Store TRG data to TRGSignalBundle.
1531  // 1 is FE output. 2 is Merger output. 3 is TRG output.
1532  //...Clock...
1533  const TRGClock& dClock = TRGCDC::getTRGCDC()->dataClock();
1534  // Create signal bundle
1535  TRGSignalBundle* allTrgData;
1536  allTrgData = new TRGSignalBundle(string("TRGData"), dClock);
1537  // Create signal vector. Need to create signals for of signal vector.
1538  TRGSignalVector* trgData;
1539  trgData = new TRGSignalVector(string("TRGData"), dClock, nTrgBitsInWindow);
1540  // Loop over data and fill TRGSignalVector
1541  for (unsigned iWindow = 0; iWindow < nTRGWindows; iWindow++) {
1542  vector<unsigned> t_trgData(inBinaryData[2].begin() + 11 * iWindow, inBinaryData[2].begin() + 11 * iWindow + 11);
1543  TRGState t_trgState(t_trgData, 1);
1544  trgData->set(t_trgState, iWindow);
1545  }
1546  if (allTrgData) {
1547  allTrgData->push_back(trgData);
1548  // Clean up memory
1549  for (unsigned i = 0; i < allTrgData->size(); i++) delete (*allTrgData)[i];
1550  delete allTrgData;
1551  }
1552 
1553 
1554  // Save CDC hit data and TRG hit data to ROOT file.
1555  if (_makeRootFile) {
1556  saveTRGRawInformation(trgInformations);
1557  saveCDCHitInformation(hitCdcData);
1558  // cppcheck-suppress knownConditionTrueFalse
1559  if (trgOutput == 2) {
1560  saveTRGHitInformation(hitTrgData);
1561  }
1562  m_treeROOTInput->Fill();
1563  }
1564 
1565  // Update CDC wires using data. Start from layer 50.
1566  for (unsigned iHit = 0; iHit < hitCdcData.size(); iHit++) {
1567  unsigned t_layerId = hitCdcData[iHit][0] + 50;
1568  unsigned t_wireId = hitCdcData[iHit][1];
1569  int t_rawTdc = hitCdcData[iHit][2];
1570  const TCWire& currentWire = *static_cast<const TCWire*>(wire(t_layerId, t_wireId));
1571 
1572  // Update currentWire._signal
1573  TRGTime rise = TRGTime(t_rawTdc, true, _clockFE, currentWire.name());
1574  TRGTime fall = rise;
1575  fall.shift(1).reverse();
1576  if (currentWire._signal.active()) {
1577  currentWire._signal = currentWire.signal() | TRGSignal(rise & fall);
1578  currentWire._signal.name(currentWire.name());
1579  } else {
1580  currentWire._signal = TRGSignal(rise & fall);
1581  currentWire._signal.name(currentWire.name());
1582 
1583  // Make TCWireHit and store in _layers.
1584  // [FIXME] TCWireHit can not accept wire being hit two times.
1585  TCWHit* hit = new TCWHit(currentWire, 0, 0, 0, t_rawTdc, 0, t_rawTdc, 0, 0, 1);
1586  hit->state(CellHitFindingValid | CellHitFittingValid);
1587  ((TCWire*)(*_layers[t_layerId])[t_wireId])->hit(hit);
1588  _hits.push_back(hit);
1589  if (currentWire.axial()) _axialHits.push_back(hit);
1590  else _stereoHits.push_back(hit);
1591  }
1592  //cout<<"iLayer, iWire, iHit: "<<t_layerId<<" "<<t_wireId<<" "<<iHit<<endl;
1593  //currentWire.signal().dump();
1594 
1595  } // End updating CDC wires
1596 
1597  //...Hit classification...
1598  classification();
1599 
1600  TRGDebug::leaveStage("TRGCDC updateByData");
1601  }
1602 
1603  void
1605  {
1606  TRGDebug::enterStage("TRGCDC classification");
1607 
1608  unsigned n = _hits.size();
1609 
1610  for (unsigned i = 0; i < n; i++) {
1611  TCWHit* h = _hits[i];
1612  const TCWire& w = h->wire();
1613  unsigned state = h->state();
1614 
1615  //...Cache pointers to a neighbor...
1616  const TCWire* neighbor[7];
1617  for (unsigned j = 0; j < 7; j++) neighbor[j] = w.neighbor(j);
1618 
1619  //...Decide hit pattern...
1620  unsigned pattern = 0;
1621  for (unsigned j = 0; j < 7; j++) {
1622  if (neighbor[j])
1623  if (neighbor[j]->hit())
1624  pattern += (1 << j);
1625  }
1626  state |= (pattern << CellHitNeighborHit);
1627 
1628  //...Check isolation...
1629  const TCWHit* hr1 = neighbor[2]->hit();
1630  const TCWHit* hl1 = neighbor[3]->hit();
1631  if ((hr1 == 0) && (hl1 == 0)) {
1632  state |= CellHitIsolated;
1633  } else {
1634  const TCWHit* hr2 = neighbor[2]->neighbor(2)->hit();
1635  const TCWHit* hl2 = neighbor[3]->neighbor(3)->hit();
1636  if (((hr2 == 0) && (hr1 != 0) && (hl1 == 0)) ||
1637  ((hl2 == 0) && (hl1 != 0) && (hr1 == 0)))
1638  state |= CellHitIsolated;
1639  }
1640 
1641  //...Check continuation...
1642  // unsigned superLayer = w.superLayerId();
1643  bool previous = false;
1644  bool next = false;
1645  if (neighbor[0] == 0) previous = true;
1646  else {
1647  if ((neighbor[0]->hit()) || neighbor[1]->hit())
1648  previous = true;
1649  // if (m_smallcell && w.layerId() == 3)
1650  // if (neighbor[6]->hit())
1651  // previous = true;
1652  }
1653  if (neighbor[5] == 0) next = true;
1654  else {
1655  if ((neighbor[4]->hit()) || neighbor[5]->hit())
1656  next = true;
1657  }
1658  // if (previous && next) state |= CellHitContinuous;
1659  if (previous || next) state |= CellHitContinuous;
1660 
1661  //...Solve LR locally...
1662  if ((pattern == 34) || (pattern == 42) ||
1663  (pattern == 40) || (pattern == 10) ||
1664  (pattern == 35) || (pattern == 50))
1665  state |= CellHitPatternRight;
1666  else if ((pattern == 17) || (pattern == 21) ||
1667  (pattern == 20) || (pattern == 5) ||
1668  (pattern == 19) || (pattern == 49))
1669  state |= CellHitPatternLeft;
1670 
1671  //...Store it...
1672  h->state(state);
1673  }
1674 
1675  TRGDebug::leaveStage("TRGCDC classification");
1676  }
1677 
1678  vector<const TCWHit*>
1679  TRGCDC::axialHits(void) const
1680  {
1681  vector<const TCWHit*> t;
1682  t.assign(_axialHits.begin(), _axialHits.end());
1683  return t;
1684 
1685  // if (! mask) return _axialHits;
1686  // else if (mask == CellHitFindingValid) return _axialHits;
1687  // cout << "TRGCDC::axialHits !!! unsupported mask given" << endl;
1688  // return _axialHits;
1689  }
1690 
1691  vector<const TCWHit*>
1693  {
1694  vector<const TCWHit*> t;
1695  t.assign(_stereoHits.begin(), _stereoHits.end());
1696  return t;
1697 
1698  // if (! mask) return _stereoHits;
1699  // else if (mask == CellHitFindingValid) return _stereoHits;
1700  // cout << "TRGCDC::stereoHits !!! unsupported mask given" << endl;
1701  // return _stereoHits;
1702  }
1703 
1704  vector<const TCWHit*>
1705  TRGCDC::hits(void) const
1706  {
1707  vector<const TCWHit*> t;
1708  t.assign(_hits.begin(), _hits.end());
1709  return t;
1710 
1711  // if (! mask) return _hits;
1712  // else if (mask == CellHitFindingValid) return _hits;
1713  // cout << "TRGCDC::hits !!! unsupported mask given" << endl;
1714  // return _hits;
1715  }
1716 
1717 // vector<const TCWHit *>
1718 // TRGCDC::badHits(void) const {
1719 // vector<const TCWHit *> t;
1720 // t.assign(_badHits.begin(), _badHits.end());
1721 // return t;
1722 
1723 // //cnv if (! updated()) update();
1724 // // if (_badHits.length()) return _badHits;
1725 
1726 // // //...Loop over RECCDC_WIRHIT bank...
1727 // // x unsigned nReccdc = BsCouTab(RECCDC_WIRHIT);
1728 // // for (unsigned i = 0; i < nReccdc; i++) {
1729 // // x struct reccdc_wirhit * h =
1730 // // (struct reccdc_wirhit *)
1731 // // BsGetEnt(RECCDC_WIRHIT, i + 1, BBS_No_Index);
1732 
1733 // // //...Check validity...
1734 // // if (h->m_stat & CellHitFindingValid) continue;
1735 
1736 // // //...Obtain a pointer to GEOCDC...
1737 // // x geocdc_wire * g =
1738 // // (geocdc_wire *) BsGetEnt(GEOCDC_WIRE, h->m_geo, BBS_No_Index);
1739 
1740 // // //...Get a pointer to a TRGCDCWire...
1741 // // TCWire * w = _wires[g->m_ID - 1];
1742 
1743 // // //...Create TCWHit...
1744 // // _badHits.append(new TCWHit(w, h, _fudgeFactor));
1745 // // }
1746 
1747 // // return _badHits;
1748 // }
1749 
1750  vector<const TCWHitMC*>
1751  TRGCDC::hitsMC(void) const
1752  {
1753  vector<const TCWHitMC*> t;
1754  t.assign(_hitsMC.begin(), _hitsMC.end());
1755  return t;
1756  }
1757 
1758  string
1759  TRGCDC::wireName(unsigned wireId) const
1760  {
1761  string as = "-";
1762  const TCWire* const w = wire(wireId);
1763  if (w) {
1764  if (w->stereo())
1765  as = "=";
1766  } else {
1767  return "invalid_wire(" + TRGUtil::itostring(wireId) + ")";
1768  }
1769  return TRGUtil::itostring(layerId(wireId)) + as + TRGUtil::itostring(localId(wireId));
1770  }
1771 
1772  unsigned
1773  TRGCDC::localId(unsigned id) const
1774  {
1775  cout << "TRGCDC::localId !!! this function is not tested yet"
1776  << endl;
1777  unsigned iLayer = 0;
1778  unsigned nW = 0;
1779  bool nextLayer = true;
1780  while (nextLayer) {
1781  unsigned nWLast = nW;
1782  nW += layer(iLayer++)->nCells();
1783  if (id < (nW - 1))
1784  return id - nWLast;
1785  if (nW >= nWires())
1786  nextLayer = false;
1787  }
1788  cout << "TRGCDC::localId !!! no such a wire (id=" << id << endl;
1789  return TRGCDC_UNDEFINED;
1790  }
1791 
1792  unsigned
1793  TRGCDC::layerId(unsigned id) const
1794  {
1795  cout << "TRGCDC::layerId !!! this function is not tested yet"
1796  << endl;
1797  unsigned iLayer = 0;
1798  unsigned nW = 0;
1799  bool nextLayer = true;
1800  while (nextLayer) {
1801  nW += layer(iLayer++)->nCells();
1802  if (id < (nW - 1))
1803  return iLayer - 1;
1804  if (nW >= nWires())
1805  nextLayer = false;
1806  }
1807  cout << "TRGCDC::layerId !!! no such a wire (id=" << id << endl;
1808  return TRGCDC_UNDEFINED;
1809  }
1810 
1811  unsigned
1812  TRGCDC::layerId(unsigned, unsigned) const
1813  {
1814  cout << "TRGCDC::layerId !!! this function is not implemented yet"
1815  << endl;
1816  return TRGCDC_UNDEFINED;
1817  }
1818 
1819  unsigned
1820  TRGCDC::superLayerId(unsigned id) const
1821  {
1822  unsigned iLayer = 0;
1823  unsigned nW = 0;
1824  bool nextLayer = true;
1825  while (nextLayer) {
1826  const vector<TRGCDCLayer*>& sl = *superLayer(iLayer);
1827  const unsigned nLayers = sl.size();
1828  for (unsigned i = 0; i < nLayers; i++)
1829  nW += sl[i]->nCells();
1830 
1831  if (id < (nW - 1))
1832  return iLayer;
1833  if (nW >= nWires())
1834  nextLayer = false;
1835  }
1836  cout << "TRGCDC::superLayerId !!! no such a wire (id=" << id
1837  << endl;
1838  return TRGCDC_UNDEFINED;
1839  }
1840 
1841  unsigned
1842  TRGCDC::localLayerId(unsigned id) const
1843  {
1844  unsigned iLayer = 0;
1845  unsigned nW = 0;
1846  bool nextLayer = true;
1847  while (nextLayer) {
1848  const vector<TRGCDCLayer*>& sl = *superLayer(iLayer);
1849  const unsigned nLayers = sl.size();
1850  for (unsigned i = 0; i < nLayers; i++) {
1851  nW += sl[i]->nCells();
1852  if (id < (nW - 1))
1853  return i;
1854  }
1855 
1856  if (nW >= nWires())
1857  nextLayer = false;
1858  }
1859  cout << "TRGCDC::localLayerId !!! no such a wire (id=" << id
1860  << endl;
1861  return TRGCDC_UNDEFINED;
1862  }
1863 
1864  unsigned
1865  TRGCDC::axialStereoSuperLayerId(unsigned aors, unsigned i) const
1866  {
1867  unsigned is = 99;
1868  // cout << "aors,i= " << aors <<" "<< i << std::endl;
1869  if (aors == 0) { //axial
1870  if (i <= 7) {
1871  is = 0;
1872  } else if (i <= 13) {
1873  is = 1;
1874  } else if (i <= 19) {
1875  is = 2;
1876  } else if (i <= 25) {
1877  is = 3;
1878  } else if (i <= 31) {
1879  is = 4;
1880  }
1881  } else if (aors == 1) { //stereo
1882  if (i <= 5) {
1883  is = 0;
1884  } else if (i <= 11) {
1885  is = 1;
1886  } else if (i <= 17) {
1887  is = 2;
1888  } else if (i <= 23) {
1889  is = 3;
1890  }
1891  }
1892 
1893  assert(is != 99);
1894  return is;
1895  }
1896 
1897 // void
1898 // TRGCDC::driftDistance(TLink & l,
1899 // const TTrack & t,
1900 // unsigned flag,
1901 // float t0Offset) {
1902 
1903 // //...No correction...
1904 // if (flag == 0) {
1905 // if (l.hit()) {
1906 // // l.drift(0, l.hit()->drift(0));
1907 // // l.drift(1, l.hit()->drift(1));
1908 // // l.dDrift(0, l.hit()->dDrift(0));
1909 // // l.dDrift(1, l.hit()->dDrift(1));
1910 // l.drift(l.hit()->drift(0), 0);
1911 // l.drift(l.hit()->drift(1), 1);
1912 // l.dDrift(l.hit()->dDrift(0), 0);
1913 // l.dDrift(l.hit()->dDrift(1), 1);
1914 // }
1915 // else {
1916 // // l.drift(0, 0.);
1917 // // l.drift(1, 0.);
1918 // // l.dDrift(0, 0.);
1919 // // l.dDrift(1, 0.);
1920 // l.drift(0., 0);
1921 // l.drift(0., 1);
1922 // l.dDrift(0., 0);
1923 // l.dDrift(0., 1);
1924 // }
1925 
1926 // return;
1927 // }
1928 
1929 // //...TOF correction...
1930 // float tof = 0.;
1931 // if (flag && 1) {
1932 // int imass = 3;
1933 // float tl = t.helix().a()[4];
1934 // float f = sqrt(1. + tl * tl);
1935 // float s = fabs(t.helix().curv()) * fabs(l.dPhi()) * f;
1936 // float p = f / fabs(t.helix().a()[2]);
1937 // calcdc_tof2_(& imass, & p, & s, & tof);
1938 // }
1939 
1940 // //...T0 correction....
1941 // if (! (flag && 2)) t0Offset = 0.;
1942 
1943 // //...Propagation corrections...
1944 // const TCWHit & h = * l.hit();
1945 // int wire = h.wire()->id();
1946 // HepGeom::Vector3D<double> tp = t.helix().momentum(l.dPhi());
1947 // float p[3] = {tp.x(), tp.y(), tp.z()};
1948 // const HepGeom::Point3D<double> & onWire = l.positionOnWire();
1949 // float x[3] = {onWire.x(), onWire.y(), onWire.z()};
1950 // //cnv float time = h.reccdc()->m_tdc + t0Offset - tof;
1951 // float time = 0;
1952 // float dist;
1953 // float edist;
1954 // int prop = (flag & 4);
1955 
1956 // //...Calculation with left side...
1957 // int side = -1;
1958 // if (side == 0) side = -1;
1959 // calcdc_driftdist_(& prop,
1960 // & wire,
1961 // & side,
1962 // p,
1963 // x,
1964 // & time,
1965 // & dist,
1966 // & edist);
1967 // // l.drift(0, dist);
1968 // // l.dDrift(0, edist);
1969 // l.drift(dist, 0);
1970 // l.dDrift(edist, 0);
1971 
1972 // //...Calculation with left side...
1973 // side = 1;
1974 // calcdc_driftdist_(& prop,
1975 // & wire,
1976 // & side,
1977 // p,
1978 // x,
1979 // & time,
1980 // & dist,
1981 // & edist);
1982 // // l.drift(1, dist);
1983 // // l.dDrift(1, edist);
1984 // l.drift(dist, 1);
1985 // l.dDrift(edist, 1);
1986 
1987 // //...tan(lambda) correction...
1988 // if (flag && 8) {
1989 // float tanl = abs(p[2]) / tp.perp();
1990 // float c;
1991 // if ((tanl >= 0.0) && (tanl < 0.5)) c = -0.48 * tanl + 1.3;
1992 // else if ((tanl >= 0.5) && (tanl < 1.0)) c = -0.28 * tanl + 1.2;
1993 // else if ((tanl >= 1.0) && (tanl < 1.5)) c = -0.16 * tanl + 1.08;
1994 // else c = 0.84;
1995 
1996 // // l.dDrift(0, l.dDrift(0) * c);
1997 // // l.dDrift(1, l.dDrift(1) * c);
1998 // l.dDrift(l.dDrift(0) * c, 0);
1999 // l.dDrift(l.dDrift(1) * c, 1);
2000 // }
2001 // }
2002 
2004  {
2005  TRGDebug::enterStage("TRGCDC destructor");
2006  clear();
2007 
2008  delete [] _width;
2009  delete [] _r;
2010  delete [] _r2;
2011 
2012  if (_tsFinder) delete _tsFinder;
2013  if (_hFinder)
2014  delete _hFinder;
2015  if (_h3DFinder)
2016  delete _h3DFinder;
2017  if (_fitter3D)
2018  delete _fitter3D;
2019 
2020 #ifdef TRGCDC_DISPLAY
2021  if (D)
2022  delete D;
2023  cout << "TRGCDC ... rphi displays deleted" << endl;
2024 #endif
2025 
2026  TRGDebug::leaveStage("TRGCDC destructor");
2027  }
2028 
2029  bool
2030  TRGCDC::neighbor(const TCWire& w0, const TCWire& w1) const
2031  {
2032  const int lyr0 = w0.layerId();
2033  const int lyr1 = w1.layerId();
2034  const int lyr = lyr0 - lyr1;
2035 
2036  if (abs(lyr) > 1) return false;
2037  if (w0.superLayerId() != w1.superLayerId()) return false;
2038 
2039  for (unsigned i = 0; i < 7; i++) {
2040  if (w0.neighbor(i)) {
2041  if (w0.neighbor(i)->id() == w1.id())
2042  return true;
2043  }
2044  }
2045  return false;
2046  }
2047 
2048  void
2050  {
2051  TRGDebug::enterStage("TRGCDC simulate");
2052 
2053  const bool fast = (_simulationMode & 1);
2054  const bool firm = (_simulationMode & 2);
2055  if (fast)
2056  fastSimulation();
2057  if (firm)
2059 
2060  TRGDebug::leaveStage("TRGCDC simulate");
2061  }
2062 
2063  void
2065  {
2066  TRGDebug::enterStage("TRGCDC fastSimulation");
2067 
2068 #ifdef TRGCDC_DISPLAY
2069  D->beginningOfEvent();
2070 #endif
2071 
2072  //...Options...
2073  const bool trackSegmentSimulationOnly = _fastSimulationMode & 1;
2074  const bool trackSegmentClockSimulation = _fastSimulationMode & 2;
2075 
2076  //...TSF simulation...
2077  _tsFinder->doit(_tss,
2078  trackSegmentClockSimulation,
2079  _segmentHits,
2080  _segmentHitsSL);
2081  _eventTime.back()->doit(_fverETF, _fprintFirmETF);
2082 
2083  if (_segmentHits.size() == 0) setReturnValue(TSF, 1);
2084 
2085  if (trackSegmentSimulationOnly) {
2086  TRGDebug::leaveStage("TRGCDC fastSimulation");
2087  return;
2088  }
2089 
2090 #ifdef TRGCDC_DISPLAY
2091  // dump("hits");
2092  string stg = "fast simulation results";
2093  string inf = "#segments=" + TRGUtilities::itostring(segmentHits().size());
2094  D->clear();
2095  D->stage(stg);
2096  D->information(inf);
2097  D->area().append(hits());
2098  D->area().append(segmentHits());
2099  D->show();
2100  D->run();
2101 #endif
2102 
2103  //...2D finder and fitter...
2104  if (_perfect2DFinder)
2106  else
2108  if (_trackList2D.size() == 0) setReturnValue(find2D, 1);
2109  if (_trackList2DFitted.size() == 0) setReturnValue(fit2D, 1);
2110 
2111  //cout<<"s3D"<<endl;
2112  //...Stereo finder...
2114 
2115  //...Check tracks...
2116  if (TRGDebug::level()) {
2117  for (unsigned i = 0; i < _trackList3D.size(); i++) {
2118  const TCTrack& t = *_trackList3D[i];
2119  cout << "> links=" << t.links().size() << endl;
2120  t.dump();
2121  }
2122  }
2123 
2124  //...Check relations...
2125  if (TRGDebug::level()) {
2126  for (unsigned i = 0; i < _trackList3D.size(); i++) {
2127  // trackList[i]->relation().dump();
2128  _trackList3D[i]->dump();
2129  }
2130  }
2131 
2132  //...Event Time... In ns scale.
2133  // [FIXME] This does nothing.
2134  _eventTime.back()->getT0();
2135 
2136  //...3D tracker...
2138  //_fitter3D->doitComplex(_trackList3D);
2139  //cout<<"e3D"<<endl;
2140 
2141  // Notify if track's debugValue is not 0
2142  {
2143  // For any track
2144  int t_debugValue = 0;
2145  for (unsigned iTrack = 0; iTrack < _trackList3D.size(); iTrack++) {
2146  t_debugValue |= _trackList3D[iTrack]->getDebugValue(TRGCDCTrack::EDebugValueType::any);
2147  }
2149  //int t_debugValue = TRGCDCTrack::EDebugValueType::any;
2150  //for(unsigned iTrack = 0; iTrack<_trackList3D.size(); iTrack++) {
2151  // t_debugValue &= _trackList3D[iTrack]->getDebugValue(TRGCDCTrack::EDebugValueType::any);
2152  //}
2153  setReturnValue(t_debugValue | getReturnValue());
2154  }
2155 
2156  if (TRGDebug::level() > 1) {
2157  for (unsigned iTrack = 0; iTrack < _trackList3D.size(); iTrack++) {
2158  const TCTrack& aTrack = *_trackList3D[iTrack];
2159  if (aTrack.fitted()) {
2160  double fitPt = aTrack.pt();
2161  double fitPhi0 = aTrack.helix().phi0();
2162  int fitCharge = aTrack.charge();
2163  if (fitCharge < 0) {
2164  fitPhi0 -= M_PI;
2165  if (fitPhi0 < 0) fitPhi0 += 2 * M_PI;
2166  }
2167  double fitZ0 = aTrack.helix().dz();
2168  double fitCot = aTrack.helix().tanl();
2169  cout << TRGDebug::tab() << "Track[" << iTrack << "]: charge=" << fitCharge
2170  << " pt=" << fitPt << " phi_c=" << fitPhi0
2171  << " z0=" << fitZ0 << " cot=" << fitCot << endl;
2172  const TCRelation& trackRelation = aTrack.relation();
2173  const MCParticle& trackMCParticle = trackRelation.mcParticle(0);
2174  int mcCharge = trackMCParticle.getCharge();
2175  double mcPt = trackMCParticle.getMomentum().Rho();
2176  double mcPhi0 = 0.0;
2177  if (mcCharge > 0) mcPhi0 = trackMCParticle.getMomentum().Phi() - M_PI / 2;
2178  if (mcCharge < 0) mcPhi0 = trackMCParticle.getMomentum().Phi() + M_PI / 2;
2179  // Change range to [0,2pi]
2180  if (mcPhi0 < 0) mcPhi0 += 2 * M_PI;
2181  // Calculated impact position
2182  TVector3 vertex = B2Vector3D(trackMCParticle.getVertex());
2183  ROOT::Math::PxPyPzEVector vector4 = trackMCParticle.get4Vector();
2184  TVector2 helixCenter;
2185  TVector3 impactPosition;
2186  Fitter3DUtility::findImpactPosition(&vertex, &vector4, mcCharge, helixCenter, impactPosition);
2187  double mcZ0 = impactPosition.Z();
2188  double mcCot = trackMCParticle.getMomentum().z() / trackMCParticle.getMomentum().Rho();
2189  cout << TRGDebug::tab() << "Track[" << iTrack << "]: mcCharge=" << mcCharge
2190  << " mcPt=" << mcPt << " mcPhi_c=" << mcPhi0
2191  << " mcZ0=" << mcZ0 << " mcCot=" << mcCot << endl;
2192  }
2193  }
2194  }
2195 
2196  if (TRGDebug::level()) {
2197  cout << TRGDebug::tab() << "Number of 2D tracks : "
2198  << _trackList2D.size() << endl;
2199  cout << TRGDebug::tab() << "Number of 2D fitted tracks : "
2200  << _trackList2DFitted.size() << endl;
2201  cout << TRGDebug::tab() << "Number of 3D tracks : "
2202  << _trackList3D.size() << endl;
2203  }
2204 
2205 #ifdef TRGCDC_DISPLAY
2206  // dump("hits");
2207  vector<const TCTrack*> t2;
2208  t2.assign(_trackList2D.begin(), _trackList2D.end());
2209  vector<const TCTrack*> t2f;
2210  t2f.assign(_trackList2DFitted.begin(), _trackList2DFitted.end());
2211  vector<const TCTrack*> t3;
2212  t3.assign(_trackList3D.begin(), _trackList3D.end());
2213  D->endOfEvent();
2214  stg = "fast simulation results";
2215  inf = "red:2D, blue:2DF, green:3D";
2216  D->clear();
2217  D->stage(stg);
2218  D->information(inf);
2219  D->area().append(hits());
2220  D->area().append(segmentHits());
2221  D->area().append(t2, Gdk::Color("red"));
2222  D->area().append(t2f, Gdk::Color("blue"));
2223  D->area().append(t3, Gdk::Color("green"));
2224  D->show();
2225  D->run();
2226 
2227  // unsigned iFront = 0;
2228  // while (const TCFrontEnd * f = _cdc.frontEnd(iFront++)) {
2229  // D->clear();
2230  // D->beginEvent();
2231  // D->area().append(* f);
2232  // D->run();
2233  // }
2234  // unsigned iMerger = 0;
2235  // while (const TCMerger * f = _cdc.merger(iMerger++)) {
2236  // D->clear();
2237  // D->beginEvent();
2238  // D->area().append(* f);
2239  // D->run();
2240  // }
2241 #endif
2242 
2243  TRGDebug::leaveStage("TRGCDC fastSimulation");
2244  return;
2245  }
2246 
2247  void
2248  TRGCDC::storeSimulationResults(string collection2Dfinder,
2249  string collection2Dfitter,
2250  string collection3Dfitter)
2251  {
2252  TRGDebug::enterStage("TRGCDC storeSimulationResults");
2253 
2254  // hits
2255  if (!(_fastSimulationMode & 2)) {
2256  StoreArray<CDCTriggerSegmentHit> storeSegmentHits;
2258  for (unsigned its = 0; its < _segmentHits.size(); ++its) {
2259  const TCSHit* segmentHit = _segmentHits[its];
2260  const CDCHit* priorityHit = cdcHits[segmentHit->iCDCHit()];
2261  const TCSegment* segment = static_cast<const TCSegment*>(&segmentHit->cell());
2262  const CDCTriggerSegmentHit* storeHit =
2263  storeSegmentHits.appendNew(*priorityHit,
2264  segment->id(),
2267  segment->priorityTime(),
2268  segment->fastestTime(),
2269  segment->foundTime());
2270  _tss[segment->id()]->addStoreHit(storeHit);
2271  // relation to all CDCHits in segment
2272  for (unsigned iw = 0; iw < segment->wires().size(); ++iw) {
2273  const TRGCDCWire* wire = (TRGCDCWire*)(*segment)[iw];
2274  if (wire->signal().active()) {
2275  // priority wire has relation weight 2
2276  double weight = (wire == &(segment->priority())) ? 2. : 1.;
2277  storeHit->addRelationTo(cdcHits[wire->hit()->iCDCHit()], weight);
2278  }
2279  }
2280  // relation to MCParticles (same as priority hit)
2281  RelationVector<MCParticle> mcrel = priorityHit->getRelationsFrom<MCParticle>();
2282  for (unsigned imc = 0; imc < mcrel.size(); ++imc) {
2283  mcrel[imc]->addRelationTo(storeHit, mcrel.weight(imc));
2284  }
2285  }
2286  }
2287  // event time
2288  StoreObjPtr<TRGTiming> eventTime("CDCTriggerEventTime");
2289  eventTime.construct(Const::CDC, getEventTime());
2290  // 2D finder tracks
2291  StoreArray<CDCTriggerTrack> storeTracks2Dfinder(collection2Dfinder);
2292  for (unsigned itr = 0; itr < _trackList2D.size(); ++itr) {
2293  const TCTrack* track2D = _trackList2D[itr];
2294  double phi0 = remainder(track2D->helix().phi0() + M_PI_2, 2. * M_PI);
2295  double omega = track2D->charge() / track2D->helix().radius();
2296  const CDCTriggerTrack* track =
2297  storeTracks2Dfinder.appendNew(phi0, omega, 0.);
2298  // relation to SegmentHits
2299  vector<TRGCDCLink*> links = track2D->links();
2300  for (unsigned its = 0; its < links.size(); ++its) {
2301  const TRGCDCSegment* segment = static_cast<const TRGCDCSegment*>(links[its]->cell());
2302  const vector<const CDCTriggerSegmentHit*> storeHits = segment->storeHits();
2303  for (unsigned ihit = 0; ihit < storeHits.size(); ++ihit) {
2304  track->addRelationTo(storeHits[ihit]);
2305  }
2306  }
2307  }
2308  // 2D fitter tracks
2309  StoreArray<CDCTriggerTrack> storeTracks2Dfitter(collection2Dfitter);
2310  for (unsigned itr = 0; itr < _trackList2DFitted.size(); ++itr) {
2311  const TCTrack* track2D = _trackList2DFitted[itr];
2312  if (!track2D->fitted()) continue;
2313  double phi0 = remainder(track2D->helix().phi0() + M_PI_2, 2. * M_PI);
2314  double omega = track2D->charge() / track2D->helix().radius();
2315  double chi2 = track2D->get2DFitChi2();
2316  const CDCTriggerTrack* track =
2317  storeTracks2Dfitter.appendNew(phi0, omega, chi2);
2318  // relation to SegmentHits
2319  vector<TRGCDCLink*> links = track2D->links();
2320  for (unsigned its = 0; its < links.size(); ++its) {
2321  const TRGCDCSegment* segment = static_cast<const TRGCDCSegment*>(links[its]->cell());
2322  const vector<const CDCTriggerSegmentHit*> storeHits = segment->storeHits();
2323  for (unsigned ihit = 0; ihit < storeHits.size(); ++ihit) {
2324  track->addRelationTo(storeHits[ihit]);
2325  }
2326  }
2327  // relation to 2D finder, assuming same order in tracklist
2328  if (storeTracks2Dfinder.getEntries() > 0)
2329  storeTracks2Dfinder[itr]->addRelationTo(track);
2330  }
2331  // 3D tracks
2332  StoreArray<CDCTriggerTrack> storeTracks3Dfitter(collection3Dfitter);
2333  for (unsigned itr = 0; itr < _trackList3D.size(); ++itr) {
2334  const TCTrack* track3D = _trackList3D[itr];
2335  if (!track3D->fitted()) continue;
2336  double phi0 = remainder(track3D->helix().phi0() + M_PI_2, 2. * M_PI);
2337  double omega = 1. / track3D->helix().radius();
2338  double chi2 = track3D->get2DFitChi2();
2339  double z = track3D->helix().dz();
2340  double cot = track3D->helix().tanl();
2341  double chi3 = track3D->get3DFitChi2();
2342  const CDCTriggerTrack* track =
2343  storeTracks3Dfitter.appendNew(phi0, omega, chi2, z, cot, chi3);
2344  // relation to SegmentHits
2345  vector<TRGCDCLink*> links = track3D->links();
2346  for (unsigned its = 0; its < links.size(); ++its) {
2347  const TRGCDCSegment* segment = static_cast<const TRGCDCSegment*>(links[its]->cell());
2348  const vector<const CDCTriggerSegmentHit*> storeHits = segment->storeHits();
2349  for (unsigned ihit = 0; ihit < storeHits.size(); ++ihit) {
2350  track->addRelationTo(storeHits[ihit]);
2351  }
2352  }
2353  // relation to 2D finder, assuming same order in tracklist
2354  if (storeTracks2Dfinder.getEntries() > 0)
2355  storeTracks2Dfinder[itr]->addRelationTo(track);
2356  }
2357 
2358  TRGDebug::leaveStage("TRGCDC storeSimulationResults");
2359  return;
2360  }
2361 
2362  void
2364  {
2365 #ifdef TRGCDC_DISPLAY
2366  D->beginningOfEvent();
2367 #endif
2368 
2369  TRGDebug::enterStage("TRGCDC firmwareSimulation");
2370 
2371  //...Wire level simulation is in update().
2372 
2373  //...Front ends...
2374  const unsigned nFronts = _fronts.size();
2375  for (unsigned i = 0; i < nFronts; i++) {
2376  //...Skip inner-most FE because they have no contribution to TRG
2377  if (i < 10)
2378  continue;
2379 
2380  //...FE simulation...
2381  _fronts[i]->simulate();
2382  }
2383 
2384  //...Mergers...
2385  const unsigned nMergers = _mergers.size();
2386  for (unsigned i = 0; i < nMergers; i++) {
2387  _mergers[i]->simulate();
2388  }
2389 
2390  //...TSFs...
2391  const unsigned nTSFBoards = _tsfboards.size();
2392  for (unsigned i = 0; i < nTSFBoards; i++) {
2393  //iw _tsfboards[i]->simulateBoard();
2394  // _tsfboards[i]->simulate();
2395  _tsfboards[i]->simulate2();
2396  }
2397 
2398  //...Event Time... In ns scale.
2399  // [FIXME] Not really firmware mode. Just for display.
2400  int t_eventTime = _eventTime.back()->getT0();
2401  // Update WireHit driftTime using eventTime
2402  // [FIXME] This method is a temorary method. Next time change everything inside EventTime class.
2403  //cout<<"Hit timing ";
2404  for (unsigned iHits = 0; iHits < _hits.size(); iHits++) {
2405  //cout<<"["<<_hits[iHits]->wire().layerId()<<"]["<<_hits[iHits]->wire().localId()<<"]: "<<(_hits[iHits]->drift(0))<<" ";
2406  _hits[iHits]->setDriftTime((_hits[iHits]->drift(0) - t_eventTime) / 250, 0);
2407  _hits[iHits]->setDriftTime((_hits[iHits]->drift(1) - t_eventTime) / 250, 1);
2408  }
2409  //cout<< endl;
2410  //cout<<"Event timing: "<<t_eventTime<<endl;
2411 
2412  //...Tracker2D...
2413  const unsigned nTracker2Ds = _tracker2Ds.size();
2414  for (unsigned i = 0; i < nTracker2Ds; i++) {
2415  _tracker2Ds[i]->simulate();
2416  }
2417 
2418 #ifdef TRGCDC_DISPLAY
2419  dump("hits");
2420  D->endOfEvent();
2421  string stg = "Firmware simulation";
2422  string inf = " ";
2423  D->clear();
2424  D->stage(stg);
2425  D->information(inf);
2426  D->area().append(hits());
2427  D->area().append(segmentHits());
2428  D->show();
2429  D->run();
2430  // unsigned iFront = 0;
2431  // while (const TCFrontEnd * f = _cdc.frontEnd(iFront++)) {
2432  // D->clear();
2433  // D->beginEvent();
2434  // D->area().append(* f);
2435  // D->run();
2436  // }
2437  // unsigned iMerger = 0;
2438  // while (const TCMerger * f = _cdc.merger(iMerger++)) {
2439  // D->clear();
2440  // D->beginEvent();
2441  // D->area().append(* f);
2442  // D->run();
2443  // }
2444 #endif
2445 
2446  TRGDebug::leaveStage("TRGCDC firmwareSimulation");
2447  }
2448 
2449  void
2451  {
2452  TRGDebug::enterStage("TRGCDC configure");
2453 
2454  //...Open configuration file...
2455  ifstream infile(_configFilename.c_str(), ios::in);
2456  if (infile.fail()) {
2457  cout << "TRGCDC !!! can not open file" << endl
2458  << " " << _configFilename << endl;
2459  return;
2460  }
2461 
2462  //...Read configuration data...
2463  char b[800];
2464  unsigned lines = 0;
2465  unsigned lastSl = 0;
2466  unsigned lastMergerLocalId = 0;
2467  while (! infile.eof()) {
2468  infile.getline(b, 800);
2469  const string l(b);
2470  string cdr = l;
2471 
2472  bool skip = false;
2473  unsigned wid = 0;
2474  unsigned lid = 0;
2475  unsigned fid = 0;
2476  unsigned mid = 0;
2477  unsigned tid = 0;
2478  // cppcheck-suppress knownConditionTrueFalse
2479  if (lid != 0) mid = lid + tid; //jb
2480  for (unsigned i = 0; i < 5; i++) {
2481  const string car = TRGUtil::carstring(cdr);
2482  cdr = TRGUtil::cdrstring(cdr);
2483 
2484  if (car == "#") {
2485  skip = true;
2486  break;
2487  } else if (car == "CDC") {
2488  skip = true;
2489  break;
2490  }
2491 
2492  if (i == 0) {
2493  wid = atoi(car.c_str());
2494  } else if (i == 1) {
2495  lid = atoi(car.c_str());
2496  } else if (i == 2) {
2497  fid = atoi(car.c_str());
2498  } else if (i == 3) {
2499  mid = atoi(car.c_str());
2500  } else if (i == 4) {
2501  tid = atoi(car.c_str());
2502  }
2503  }
2504 
2505  if (skip)
2506  continue;
2507  if (lines != wid)
2508  continue;
2509 
2510  //...Super layer ID...
2511  const unsigned sl = _wires[wid]->superLayerId();
2512  if (sl != lastSl)
2513  lastMergerLocalId = 0;
2514 
2515  //...Make a front-end board if necessary...
2516  bool newFrontEnd = false;
2517  TCFrontEnd* f = 0;
2518  if (fid < _fronts.size())
2519  f = _fronts[fid];
2520  if (! f) {
2521  newFrontEnd = true;
2522  const string name = "CDCFrontEnd" + TRGUtil::itostring(fid);
2523  TCFrontEnd::boardType t;
2524  if (sl == 0) {
2525  if (_wires[wid]->localLayerId() < 5)
2526  t = TCFrontEnd::innerInside;
2527  else
2528  t = TCFrontEnd::innerOutside;
2529  } else {
2530  if (_wires[wid]->localLayerId() < 3)
2531  t = TCFrontEnd::outerInside;
2532  else
2533  t = TCFrontEnd::outerOutside;
2534  }
2535  f = new TCFrontEnd(name, t, _clock, _clockD, _clockUser3125);
2536 
2537  _fronts.push_back(f);
2538  }
2539  f->push_back(_wires[wid]);
2540 
2541  //...Make a merger board if necessary... 2013,0908: physjg: I
2542  // think this should be done only when a new frontboard is
2543  // created. i.e., inlcuded in the if(!f) { .... } block.
2544  // Or I can put this part inside the new FrontEnd creation
2545  // part, well, seems not good in coding
2546  //
2547  bool newMerger = false;
2548  TCMerger* m = nullptr;
2549  if (newFrontEnd) {
2550  if (mid != 99999) {
2551  if (mid < _mergers.size())
2552  m = _mergers[mid];
2553  if (! m) {
2554  newMerger = true;
2555  const string name = "CDCMerger" + TRGUtil::itostring(sl) +
2556  "-" + TRGUtil::itostring(lastMergerLocalId);
2557  TCMerger::unitType mt;
2558  if (sl == 0)
2559  mt = TCMerger::innerType;
2560  else
2561  mt = TCMerger::outerType;
2562  m = new TCMerger(name,
2563  mt,
2564  _clock,
2565  _clockD,
2567  _clockUser6250);
2568  _mergers.push_back(m);
2569  ++lastMergerLocalId;
2570  lastSl = sl;
2571  //cout << "new merger : " << name << endl;
2572  }
2573  m->push_back(f);
2574  }
2575 
2576  //...Make Aurora channel...
2577  if (f && m) {
2578  const string n = f->name() + string("-") + m->name();
2579  TRGChannel* ch = new TRGChannel(n, *f, *m);
2580  f->appendOutput(ch);
2581  m->appendInput(ch);
2582  }
2583  }
2584 
2585  //...Make a TSF board if necessary...
2586  if (newMerger) {
2587  TSFinder* t = 0;
2588  if (tid != 99999) {
2589  if (tid < _tsfboards.size())
2590  t = _tsfboards[tid];
2591  if (!t) {
2592  const string name = "CDCTSFBoard" + TRGUtil::itostring(tid);
2593  TSFinder::boardType tt;
2594  if (_wires[wid]->superLayerId() == 0)
2595  tt = TSFinder::innerType;
2596  else
2597  tt = TSFinder::outerType;
2598  t = new TSFinder(*this,
2599  name,
2600  tt,
2601  _clock,
2602  _clockD,
2605  _tsSL[tid]);
2606  _tsfboards.push_back(t);
2607  }
2608  t->push_back(m);
2609  }
2610 
2611  if (m && t) {
2612  const string n = m->name() + string("-") + t->name();
2613  TRGChannel* chmt = new TRGChannel(n, *m, *t);
2614  m->appendOutput(chmt);
2615  t->appendInput(chmt);
2616  }
2617  }
2618 
2619  ++lines;
2620  }
2621  infile.close();
2622 
2623  //...Make a 2D finder...
2624  for (unsigned i = 0; i < 4; i++) {
2625  const string name = "CDC2DBoard" + TRGUtil::itostring(i);
2626  TCTracker2D* t = new TCTracker2D(name,
2627  _clock,
2628  _clockD,
2630  _clockUser6250);
2631  _tracker2Ds.push_back(t);
2632  for (unsigned j = 0; j < 9; j++) {
2633  TCTSFinder* tsf = _tsfboards[j];
2634  const string n = tsf->name() + string("-") + t->name();
2635  TRGChannel* ch = new TRGChannel(n, *tsf, *t);
2636  tsf->appendOutput(ch);
2637  t->appendInput(ch);
2638  }
2639  const string n = t->name() + string("-") + "CDC3DBoard" + TRGUtil::itostring(i);
2640  TRGChannel* ch = new TRGChannel(n, *t, *t);
2641  // last (* t) should be 3D tracker in future
2642  t->appendOutput(ch);
2643  }
2644 
2645  //...For debug...
2646  if (TRGDebug::level()) {
2647  cout << TRGDebug::tab() << "TSF configuration" << endl;
2648  for (unsigned i = 0; i < _tsfboards.size(); i++) {
2649  const TSFinder& t = *_tsfboards[i];
2650  t.dump("detail", TRGDebug::tab() + " ");
2651  }
2652  cout << TRGDebug::tab() << "Tracker2D configuration" << endl;
2653  for (unsigned i = 0; i < _tracker2Ds.size(); i++) {
2654  const TCTracker2D& t = *_tracker2Ds[i];
2655  t.dump("detail", TRGDebug::tab() + " ");
2656  }
2657  }
2658 
2659  TRGDebug::leaveStage("TRGCDC configure");
2660  }
2661 
2662  const TRGCDCSegment&
2663  TRGCDC::segment(unsigned lid, unsigned id) const
2664  {
2665  return *(const TRGCDCSegment*)(* _tsLayers[lid])[id];
2666  }
2667 
2668 //jb void
2669 //jb TRGCDC::perfect3DFinder(vector<TCTrack*> trackList) const {
2670 //jb
2671 //jb TRGDebug::enterStage("Perfect 3D Finder");
2672 //jb if (TRGDebug::level())
2673 //jb cout << TRGDebug::tab() << "givenTrk#=" << trackList.size() << endl;
2674 //jb
2675 //jb //...Track loop....
2676 //jb for (unsigned j = 0; j < trackList.size(); j++) {
2677 //jb
2678 //jb //...G4 trackID...
2679 //jb TCTrack * trk = trackList[j];
2680 //jb unsigned id = trackList[j]->relation().contributor(0);
2681 //jb vector<const TCSHit*> tsList[9];
2682 //jb
2683 //jb //...Segment loop...
2684 //jb const vector<const TCSHit*> hits = segmentHits();
2685 //jb for (unsigned i = 0; i < hits.size(); i++) {
2686 //jb const TCSHit & ts = * hits[i];
2687 //jb if (ts.segment().axial()) continue;
2688 //jb if (! ts.signal().active()) continue;
2689 //jb const TCWHit * wh = ts.segment().center().hit();
2690 //jb if (! wh) continue;
2691 //jb const unsigned trackId = wh->iMCParticle();
2692 //jb
2693 //jb if (id == trackId)
2694 //jb tsList[wh->wire().superLayerId()].push_back(& ts);
2695 //jb }
2696 //jb
2697 //jb if (TRGDebug::level()) {
2698 //jb cout << TRGDebug::tab() << "trk#" << j << endl;
2699 //jb for (unsigned k = 0; k < 9; k++) {
2700 //jb if (k % 2) {
2701 //jb cout << TRGDebug::tab(4) << "superlayer " << k << ":";
2702 //jb for (unsigned l = 0; l < tsList[k].size(); l++) {
2703 //jb if (l)
2704 //jb cout << ",";
2705 //jb cout << tsList[k][l]->cell().name();
2706 //jb }
2707 //jb cout << endl;
2708 //jb }
2709 //jb }
2710 //jb }
2711 //jb
2712 //jb //...Select best one in each super layer...
2713 //jb for (unsigned i = 0; i < 9; i++) {
2714 //jb const TCSHit * best = 0;
2715 //jb if (tsList[i].size() == 0) {
2716 //jb continue;
2717 //jb } else if (tsList[i].size() == 1) {
2718 //jb best = tsList[i][0];
2719 //jb } else {
2720 //jb int timeMin = 99999;
2721 //jb for (unsigned k = 0; k < tsList[i].size(); k++) {
2722 //jb const TRGSignal & signal = tsList[i][k]->signal();
2723 //jb const TRGTime & t = * signal[0];
2724 //jb if (t.time() < timeMin) {
2725 //jb timeMin = t.time();
2726 //jb best = tsList[i][k];
2727 //jb }
2728 //jb }
2729 //jb }
2730 //jb //trk->append(new TCLink(0,
2731 //jb // best,
2732 //jb // best->cell().xyPosition()));
2733 //jb }
2734 //jb
2735 //jb if (TRGDebug::level())
2736 //jb trk->dump("", "> ");
2737 //jb
2738 //jb }
2739 //jb
2740 //jb TRGDebug::leaveStage("Perfect 3D Finder");
2741 //jb }
2742 
2743  int TRGCDC::getEventTime(void) const
2744  {
2745  return _eventTime.back()->getT0();
2746  }
2747 
2748  void TRGCDC::setReturnValue(EReturnValueType const& moduleName, bool flag)
2749  {
2750  if (flag) _returnValue |= moduleName;
2751  else _returnValue &= ~moduleName;
2752  }
2753 
2754  int TRGCDC::getReturnValue(EReturnValueType const& moduleName) const
2755  {
2756  return _returnValue & moduleName;
2757  }
2758 
2759 
2760  void TRGCDC::saveCDCHitInformation(vector<vector<unsigned>>& hitWiresFromCDC)
2761  {
2762  TClonesArray& rootCDCHitInformation = *m_rootCDCHitInformation;
2763  rootCDCHitInformation.Clear();
2764  for (unsigned iHit = 0; iHit < hitWiresFromCDC.size(); iHit++) {
2765  TVectorD t_rootCDCHitInformation(5);
2766  t_rootCDCHitInformation[0] = hitWiresFromCDC[iHit][0];
2767  t_rootCDCHitInformation[1] = hitWiresFromCDC[iHit][1];
2768  t_rootCDCHitInformation[2] = hitWiresFromCDC[iHit][2];
2769  t_rootCDCHitInformation[3] = hitWiresFromCDC[iHit][3];
2770  t_rootCDCHitInformation[4] = hitWiresFromCDC[iHit][4];
2771  new (rootCDCHitInformation[iHit]) TVectorD(t_rootCDCHitInformation);
2772  } // End of hit loop
2773  }
2774 
2775  void TRGCDC::saveTRGHitInformation(vector<vector<int>>& hitWiresFromTRG)
2776  {
2777  TClonesArray& rootTRGHitInformation = *m_rootTRGHitInformation;
2778  rootTRGHitInformation.Clear();
2779  for (unsigned iHit = 0; iHit < hitWiresFromTRG.size(); iHit++) {
2780  TVectorD t_rootTRGHitInformation(4);
2781  t_rootTRGHitInformation[0] = hitWiresFromTRG[iHit][0];
2782  t_rootTRGHitInformation[1] = hitWiresFromTRG[iHit][1];
2783  t_rootTRGHitInformation[2] = hitWiresFromTRG[iHit][2];
2784  t_rootTRGHitInformation[3] = hitWiresFromTRG[iHit][3];
2785  new (rootTRGHitInformation[iHit]) TVectorD(t_rootTRGHitInformation);
2786  } // End of hit loop
2787  }
2788 
2789  void TRGCDC::saveTRGRawInformation(vector<string >& trgInformations)
2790  {
2791  TClonesArray& rootTRGRawInformation = *m_rootTRGRawInformation;
2792  rootTRGRawInformation.Clear();
2793  for (unsigned iWindow = 0; iWindow < trgInformations.size(); iWindow++) {
2794  TObjString t_rootTRGRawInformation;
2795  t_rootTRGRawInformation.SetString(trgInformations[iWindow].c_str());
2796  new (rootTRGRawInformation[iWindow]) TObjString(t_rootTRGRawInformation);
2797  } // End of hit loop
2798  }
2799 
2800 //void TRGCDC::saveCDCTRGTimeMatchInformation(vector<vector<map<int, int> > >& hitTimingComp) {
2801 //
2802 // TClonesArray& CDCTRGTimeMatch = *m_CDCTRGTimeMatch;
2803 // CDCTRGTimeMatch.Clear();
2804 // int nTotalHits = 0;
2805 // for(unsigned iLayer=0; iLayer<5; iLayer++){
2806 // for(unsigned iWire=0; iWire<16; iWire++){
2807 // for(map<int, int>::iterator it = hitTimingComp[iLayer][iWire].begin(); it != hitTimingComp[iLayer][iWire].end(); it++){
2808 // TVectorD t_CDCTRGTimeMatch(4);
2809 // t_CDCTRGTimeMatch[0] = iLayer;
2810 // t_CDCTRGTimeMatch[1] = iWire;
2811 // t_CDCTRGTimeMatch[2] = (*it).first;
2812 // t_CDCTRGTimeMatch[3] = (*it).second;
2813 // new(CDCTRGTimeMatch[nTotalHits++]) TVectorD(t_CDCTRGTimeMatch);
2814 // }
2815 // } // End of wire loop
2816 // } // End of layer loop
2817 //
2818 //}
2819 
2820  unsigned
2821  TRGCDC::nSegments(unsigned id) const
2822  {
2823  if (id < _tsLayers.size())
2824  return _tsLayers[id]->size();
2825  return 0;
2826  }
2827 
2829 } // namespace Belle2
DataType y() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:427
DataType z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:429
DataType x() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:425
Class containing the result of the unpacker in raw data and the result of the digitizer in simulation...
Definition: CDCHit.h:40
Combination of several CDCHits to a track segment hit for the trigger.
This class is the interface between TSim/basf2 TSF module and the firmware simulation core of XSim/IS...
Track created by the CDC trigger.
int nShifts(int layerId) const
Returns number shift.
double fieldWireR(int layerId) const
Returns radius of field wire in each layer.
const B2Vector3D wireForwardPosition(uint layerId, int cellId, EWirePosition set=c_Base) const
Returns the forward position of the input sense wire.
const B2Vector3D wireBackwardPosition(uint layerId, int cellId, EWirePosition set=c_Base) const
Returns the backward position of the input sense wire.
double getTdcBinWidth() const
Return TDC bin width (nsec).
float getT0(const WireID &wireID) const
Returns t0 parameter of the specified sense wire.
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
double getNominalDriftV() const
Return the nominal drift velocity of He-ethane gas (default: 4.0x10^-3 cm/nsec).
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
unsigned nWireLayers() const
Returns a number of wire layers.
double senseWireR(int layerId) const
Returns radius of sense wire in each layer.
double zOffsetWireLayer(unsigned i) const
Returns the offset of z of the wire layer i.
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
ROOT::Math::XYZVector getVertex() const
Return production vertex position, shorthand for getProductionVertex().
Definition: MCParticle.h:183
float getCharge() const
Return the particle charge defined in TDatabasePDG.
Definition: MCParticle.cc:36
ROOT::Math::PxPyPzEVector get4Vector() const
Return 4Vector of particle.
Definition: MCParticle.h:207
ROOT::Math::XYZVector getMomentum() const
Return momentum.
Definition: MCParticle.h:198
The Raw COPPER class This class stores data received by COPPER via belle2linkt Data from all detector...
Definition: RawCOPPER.h:52
void SetBuffer(int *bufin, int nwords, int delete_flag, int num_events, int num_nodes) OVERRIDE_CPP17
set buffer ( delete_flag : m_buffer is freeed( = 0 )/ not freeed( = 1 ) in Destructer )
Definition: RawCOPPER.cc:141
virtual int * GetBuffer(int n)
get nth buffer pointer
Definition: RawDataBlock.h:53
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:62
int getEntries() const
Get the number of elements.
Class for type safe access to objects that are referred to in relations.
size_t size() const
Get number of relations.
float weight(int index) const
Get weight with index.
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
RelationVector< FROM > getRelationsFrom(const std::string &name="", const std::string &namedRelation="") const
Get the relations that point from another store array to this object.
Accessor to arrays stored in the data store.
Definition: StoreArray.h:113
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:246
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:216
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:96
bool construct(Args &&... params)
Construct an object of type T in this StoreObjPtr, using the provided constructor arguments.
Definition: StoreObjPtr.h:119
A class to represent a cell layer.
Definition: Layer.h:33
A class to represent a wire in CDC.
Definition: Segment.h:39
A class to represent a wire in CDC.
Definition: Wire.h:56
The instance of TRGCDC is a singleton.
Definition: TRGCDC.h:69
TClonesArray * m_fitParameters
3D fit
Definition: TRGCDC.h:747
bool _fileFitter3D
Switch for Fitter3D.root file.
Definition: TRGCDC.h:602
std::vector< TRGCDCTrack * > _trackList3D
Track list for 3D fitted tracks.
Definition: TRGCDC.h:542
float * _r
R of cell.
Definition: TRGCDC.h:668
std::vector< TRGCDCTracker2D * > _tracker2Ds
CDC 2D finder boards.
Definition: TRGCDC.h:708
bool _fFitter3Ds2DFitDrift
Switch to us wire 2D fit or drift 2D fit.
Definition: TRGCDC.h:578
int _firmwareSimulationStopDataClock
Firmware simulation stop clock in CDCTRG data clock.
Definition: TRGCDC.h:530
std::vector< TRGCDCWireHit * > _axialHits
CDC hits on axial wires.
Definition: TRGCDC.h:635
bool _fFitter3Ds2DFit
Switch to us 2D fit or Hough finder results.
Definition: TRGCDC.h:575
unsigned _simulationMode
Simulation mode.
Definition: TRGCDC.h:502
TClonesArray * m_rootTRGRawInformation
[0]: iLayer, [1]: iWire, [2]: Timing(CDC), [3]: MatchNumber MatchNumber: 1: Only CDC,...
Definition: TRGCDC.h:788
std::string _configFilename
CDC trigger configuration filename.
Definition: TRGCDC.h:499
TTree * m_tree
root tree for reconstructed 3D tracks
Definition: TRGCDC.h:738
std::vector< std::vector< TRGCDCLayer * > * > _superLayers
Super layers.
Definition: TRGCDC.h:608
unsigned _fastSimulationMode
Fast simulation mode.
Definition: TRGCDC.h:505
TTree * m_treeAllTracks
root tree for MC tracks
Definition: TRGCDC.h:741
std::vector< TRGCDCWire * > _hitWires
Wires with a hit.
Definition: TRGCDC.h:629
std::vector< TRGCDCSegment * > _tsSL[9]
Track Segments.
Definition: TRGCDC.h:650
int _trgCDCDataInputMode
Switch for TRG CDC input mode.
Definition: TRGCDC.h:768
std::string _cdchitCollectionName
name of the CDCHit DataStore array used as input
Definition: TRGCDC.h:771
int _debugLevel
Debug level.
Definition: TRGCDC.h:493
std::vector< TRGCDCTrack * > _trackList2DFitted
Track list for 2D fitted tracks.
Definition: TRGCDC.h:539
int _finder3DMode
Sets mode of 3DFinder.
Definition: TRGCDC.h:599
std::vector< TRGCDCTrackSegmentFinder * > _tsfboards
CDC trigger tsf boards.
Definition: TRGCDC.h:705
bool _fXtSimpleFitter3D
Switch for using simple x-t curve or non-linear x-t curve. 0: non-linear 1: simple.
Definition: TRGCDC.h:605
TClonesArray * m_mcParameters
MC.
Definition: TRGCDC.h:750
float * _width
Cell width in radian.
Definition: TRGCDC.h:665
std::string _outerTSLUTFilename
The filename of LUT for outer track segments.
Definition: TRGCDC.h:557
std::string _innerTSLUTFilename
The filename of LUT for the inner-most track segments.
Definition: TRGCDC.h:554
bool _fLogicLUTTSF
Switch for logic or LUT TSF.
Definition: TRGCDC.h:566
int _fverETF
Switch for selecting ETF version.
Definition: TRGCDC.h:590
std::vector< TRGCDCWireHitMC * > _hitsMC
MC info. of CDC hits.
Definition: TRGCDC.h:644
TClonesArray * m_mcTrackVertexVector
MC vertex.
Definition: TRGCDC.h:756
std::vector< TRGCDCSegment * > _tss
Track Segments.
Definition: TRGCDC.h:647
const bool _perfect2DFinder
Switch to activate perfect 2D finder.
Definition: TRGCDC.h:548
bool _fileHough3D
Switch for Hough3D.root file.
Definition: TRGCDC.h:596
std::vector< TRGCDCLayer * > _stereoLayers
Stereo layers.
Definition: TRGCDC.h:623
const TRGClock _clockD
CDC trigger data clock.
Definition: TRGCDC.h:684
bool _fileTSF
Switch for TSF.root file.
Definition: TRGCDC.h:584
bool _fileETF
Switch for ETF.root file.
Definition: TRGCDC.h:587
std::vector< std::vector< TRGCDCLayer * > * > _stereoSuperLayers
Stereo super layers.
Definition: TRGCDC.h:614
std::string _rootTRGCDCFilename
The filename of root file for TRGCDC.
Definition: TRGCDC.h:560
std::vector< TRGCDCSegmentHit * > _segmentHits
Track Segments with hits.
Definition: TRGCDC.h:656
std::vector< TRGCDCLayer * > _tsLayers
Track Segment layers.
Definition: TRGCDC.h:653
TTree * m_treeROOTInput
Debugging members for firmware ROOT input.
Definition: TRGCDC.h:779
std::vector< TRGCDCWireHit * > _badHits
Bad CDC hits.(not used now)
Definition: TRGCDC.h:641
bool _makeRootFile
Switch for TRGCDC.root file.
Definition: TRGCDC.h:545
EReturnValueType
Enum for returnValue types.
Definition: TRGCDC.h:74
int _returnValue
Return value for trg cdc module;.
Definition: TRGCDC.h:533
std::vector< TRGCDCTrack * > _trackList2D
Track list by 2D finding.
Definition: TRGCDC.h:536
bool _fprintFirmETF
Switch for printing Firmware inputs for ETF.
Definition: TRGCDC.h:593
TRGCDCHoughFinder * _hFinder
Hough finder.
Definition: TRGCDC.h:720
TRGCDCTrackSegmentFinder * _tsFinder
Track Segement Finder.
Definition: TRGCDC.h:711
std::vector< TRGCDCEventTime * > _eventTime
EventTime.
Definition: TRGCDC.h:729
TClonesArray * m_mcTrack4Vector
MC track.
Definition: TRGCDC.h:753
TRGSignal _firmwareSimulationWindow
Firmware simulation time window in FE.
Definition: TRGCDC.h:524
std::vector< TRGCDCWire * > _wires
All wires.
Definition: TRGCDC.h:626
TClonesArray * m_rootTRGHitInformation
[0]: iLayer, [1]: iWire, [2]: window number, [3]: priority timing
Definition: TRGCDC.h:783
std::vector< TRGCDCLayer * > _layers
All layers.
Definition: TRGCDC.h:617
const TRGClock _clock
CDC trigger system clock.
Definition: TRGCDC.h:674
TClonesArray * m_mcTrackStatus
MC track status.
Definition: TRGCDC.h:759
double _inefficiency
Hit inefficiency parameter.
Definition: TRGCDC.h:581
TClonesArray * m_evtTime
Event time.
Definition: TRGCDC.h:762
TFile * m_file
root file
Definition: TRGCDC.h:735
TRGCDCFitter3D * _fitter3D
3D fitter.
Definition: TRGCDC.h:726
std::vector< TRGCDCFrontEnd * > _fronts
CDC front-end boards.
Definition: TRGCDC.h:699
float * _r2
R^2 of cell.
Definition: TRGCDC.h:671
std::vector< TRGCDCSegmentHit * > _segmentHitsSL[9]
Track Segments with hits in each super layer.
Definition: TRGCDC.h:659
std::vector< TRGCDCWireHit * > _stereoHits
CDC hits on stereo wires.
Definition: TRGCDC.h:638
int _firmwareSimulationStop
Fimrware simulation stop clock in FE.
Definition: TRGCDC.h:521
int _firmwareSimulationStart
Fimrware simulation start clock in FE.
Definition: TRGCDC.h:518
TRGCDCPerfectFinder * _pFinder
Perfect 2D finder.
Definition: TRGCDC.h:714
std::vector< std::vector< TRGCDCLayer * > * > _axialSuperLayers
Axial super layers.
Definition: TRGCDC.h:611
TClonesArray * _tracks2D
2D track information
Definition: TRGCDC.h:765
std::vector< TRGCDCLayer * > _axialLayers
Axial layers.
Definition: TRGCDC.h:620
TRGCDCPerfectFinder * _p3DFinder
Perfect 3D finder.
Definition: TRGCDC.h:717
CDC::CDCGeometryPar * m_cdcp
returns a pointer to CDCGeometryPar
Definition: TRGCDC.h:252
TRGCDCHough3DFinder * _h3DFinder
Hough 3D finder.
Definition: TRGCDC.h:723
const TRGClock _clockUser3125
CDC trigger user clock for Aurora 3.125 Gbps.
Definition: TRGCDC.h:687
std::vector< TRGCDCWireHit * > _hits
CDC hits.
Definition: TRGCDC.h:632
const TRGClock _clockFE
CDC front end clock. Resolution is CDC TdcBinWidth.
Definition: TRGCDC.h:677
std::string _rootFitter3DFilename
The filename of root file for Fitter3D.
Definition: TRGCDC.h:563
bool _fFitter3Dsmclr
Switch for MC L/R information in Fitter3D.
Definition: TRGCDC.h:572
std::vector< TRGCDCMerger * > _mergers
CDC trigger merger boards.
Definition: TRGCDC.h:702
bool _fLRLUT
Switch for the LR LUT in Fitter3D.
Definition: TRGCDC.h:569
const TRGClock _clockUser6250
CDC trigger user clock for Aurora 6.250 Gbps.
Definition: TRGCDC.h:690
TTree * _tree2D
root tree for 2D tracks
Definition: TRGCDC.h:744
const TRGClock _clockTDC
CDC trigger TDC clock.
Definition: TRGCDC.h:681
int _firmwareSimulationStartDataClock
Firmware simulation start clock in CDCTRG data clock.
Definition: TRGCDC.h:527
TClonesArray * m_rootCDCHitInformation
[0]: iLayer, [1]: iWire, [2]: CDCADC, [3]: CDCTDC, [4]: CDC FE TRG timing
Definition: TRGCDC.h:781
A class to represent a serial link between trigger hardware modules.
Definition: Channel.h:24
A class to represent a digitized signal. Unit is nano second.
Definition: Clock.h:38
A class to represent a bundle of SignalVectors.
Definition: SignalBundle.h:26
A class to represent a bundle of digitized signals.
Definition: SignalVector.h:26
A class to represent a digitized signal. Unit is nano second.
Definition: Signal.h:23
A class to represent a state of multi bits.
Definition: State.h:24
A class to represent a signal timing in the trigger system.
Definition: Time.h:25
Class to identify a wire inside the CDC.
Definition: WireID.h:34
static void findImpactPosition(TVector3 *mcPosition, ROOT::Math::PxPyPzEVector *mcMomentum, int charge, TVector2 &helixCenter, TVector3 &impactPosition)
MC calculation functions Calculates the impact position of track.
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
int GetDetectorNwords(int n, int finesse_num)
get Detector buffer length
Definition: RawCOPPER.h:657
int * GetDetectorBuffer(int n, int finesse_num)
get Detector buffer
Definition: RawCOPPER.h:681
std::vector< const TRGCDCSegmentHit * > segmentHits(void) const
returns a list of TRGCDCSegmentHit.
Definition: TRGCDC.h:996
unsigned nAxialLayers(void) const
return # of axial layers.
Definition: TRGCDC.h:884
unsigned size(void) const
returns bit size.
Definition: State.h:149
unsigned localId(unsigned wireId) const
returns local ID in a layer. This function is expensive.
Definition: TRGCDC.cc:1773
void configure(void)
configures trigger modules for firmware simulation.
Definition: TRGCDC.cc:2450
static std::string tab(void)
returns tab spaces.
Definition: Debug.cc:47
const TRGSignalVector & set(const TRGState &, int clockPosition)
sets state at given clock.
static TRGCDC * getTRGCDC(void)
returns TRGCDC object.
Definition: TRGCDC.cc:192
void setReturnValue(EReturnValueType const &moduleName, bool flag)
sets return value for trg cdc module.
Definition: TRGCDC.cc:2748
unsigned layerId(unsigned wireId) const
returns layer ID. This function is expensive.
Definition: TRGCDC.cc:1793
std::string name(void) const
simulates track segment decisions.
Definition: TRGCDC.cc:92
void initialize()
Initialization.
Definition: Fitter3D.cc:78
unsigned localLayerId(unsigned wireId) const
returns local layer ID in a super layer. This function is expensive.
Definition: TRGCDC.cc:1842
void fastClear(void)
clears TRGCDC information.
Definition: TRGCDC.cc:981
double absoluteTime(int clockPosition) const
returns absolute time of clock position
Definition: Clock.cc:128
int priorityPosition(void) const
return priority cell position in TSHit. 0: no hit, 3: 1st priority, 1: 2nd right, 2: 2nd left
Definition: Segment.cc:473
std::string wireName(unsigned wireId) const
returns wire name.
Definition: TRGCDC.cc:1759
int FindAndFit(std::vector< TRGCDCTrack * > &trackList2D, std::vector< TRGCDCTrack * > &trackList2DFitted)
do track finding and fitting (wrapper that can choose between different versions).
Definition: HoughFinder.cc:755
TRGState subset(unsigned i, unsigned n) const
returns subset from i with n bits.
Definition: State.cc:356
float fastestTime(void) const
return fastest time in TSHit.
Definition: Segment.cc:387
unsigned id(void) const
returns id.
Definition: Cell.h:197
int doit(std::vector< TRGCDCTrack * > &trackList)
Does track fitting.
Definition: Fitter3D.cc:213
const TRGClock & clock(void) const
returns clock.
Definition: Time.h:149
int doit(std::vector< TRGCDCTrack * > &trackListClone, std::vector< TRGCDCTrack * > &trackList)
do track finding.
unsigned localLayerId(void) const
returns local layer id in a super layer.
Definition: Layer.h:173
void terminate(void)
terminates when run is finished
Definition: TRGCDC.cc:739
TRGTime & reverse(void)
reverse edge.
Definition: Time.h:141
int getValue(unsigned) const
get LUT Values
Definition: LUT.cc:47
const std::string & name(void) const
returns name.
Definition: Signal.h:188
static void enterStage(const std::string &stageName)
Declare that you enter new stage.
Definition: Debug.cc:24
bool neighbor(const TRGCDCWire &w0, const TRGCDCWire &w1) const
returns true if w0 and w1 are neighbor.
Definition: TRGCDC.cc:2030
unsigned superLayerId(void) const
returns super layer id.
Definition: Layer.h:166
virtual ~TRGCDC()
Destructor.
Definition: TRGCDC.cc:2003
void saveTRGHitInformation(std::vector< std::vector< int > > &)
Save functions for ROOT.
Definition: TRGCDC.cc:2775
void terminate()
Termination.
Definition: Fitter3D.cc:1530
const std::string stereoType(void) const
returns "A" or "U" or "V" depending on stereo type.
Definition: Layer.h:238
const TRGSignal & signal(void) const override
returns an input to the trigger. This is sync'ed to 1GHz clock.
Definition: Wire.h:226
static std::string itostring(int i)
converts int to string. (Use boost::lexical_cast)
Definition: Utilities.cc:50
const TRGCDCWire & priority(void) const
returns priority wire.
Definition: Segment.cc:525
unsigned nAxialSuperLayers(void) const
return # of axial super layers.
Definition: TRGCDC.h:891
unsigned axialStereoSuperLayerId(void) const
returns id of axial or stereo super layer id.
Definition: Layer.h:224
const TRGState & set(unsigned position, bool state=true)
sets state at bit i.
Definition: State.h:305
std::vector< const TRGCDCWireHitMC * > hitsMC(void) const
returns a list of TRGCDCWireHitMC.
Definition: TRGCDC.cc:1751
const TRGClock & dataClock(void) const
returns the data clock.
Definition: TRGCDC.h:982
unsigned nSegments(void) const
returns # of track segments.
Definition: TRGCDC.h:954
TRGTime & shift(int unit)
delays by clock unit.
Definition: Time.h:163
const TRGCDCWireHit * hit(void) const
returns a pointer to a TRGCDCWireHit.
Definition: Wire.h:144
unsigned axialStereoSuperLayerId(unsigned axialStereo, unsigned axialStereoLayerId) const
returns axialStereo super layer ID. This function is expensive.
Definition: TRGCDC.cc:1865
unsigned nSuperLayers(void) const
returns # of super layers.
Definition: TRGCDC.h:870
unsigned nStereoLayers(void) const
returns # of stereo layers.
Definition: TRGCDC.h:877
void initialize(void)
initilize variables.
Definition: Segment.cc:65
unsigned lutPattern(void) const
hit pattern containing bit for priority position
Definition: Segment.cc:559
std::string versionCDC(void) const
returns CDC version.
Definition: TRGCDC.h:856
void initialize(unsigned houghFinderPeakMin, const std::string &houghMappingFilePlus, const std::string &houghMappingFileMinus, unsigned houghDoit)
initializes CDC geometry.
Definition: TRGCDC.cc:372
const TRGCDCSegment & segment(unsigned id) const
returns a track segment.
Definition: TRGCDC.h:933
void saveTRGRawInformation(std::vector< std::string > &)
Save functions for ROOT.
Definition: TRGCDC.cc:2789
const TRGCDCLUT * LUT(void) const
returns LUT
Definition: Segment.h:243
unsigned nLayers(void) const
return # of layers.
Definition: TRGCDC.h:905
int getEventTime(void) const
returns bad hits(finding invalid hits).
Definition: TRGCDC.cc:2743
const std::vector< TRGCDCLayer * > * superLayer(unsigned id) const
returns a pointer to a super-layer.
Definition: TRGCDC.h:835
void dump(const std::string &message) const
dumps debug information.
Definition: TRGCDC.cc:767
std::vector< const TRGCDCWireHit * > stereoHits(void) const
returns a list of stereo hits.
Definition: TRGCDC.cc:1692
const std::vector< const CDCTriggerSegmentHit * > storeHits(void) const
returns a pointer to a CDCTriggerSegmentHit.
Definition: Segment.h:236
unsigned iCDCHit(void) const
returns an index to CDCHit.
Definition: CellHit.h:360
std::vector< const TRGCDCWireHit * > axialHits(void) const
returns a list of axial hits.
Definition: TRGCDC.cc:1679
unsigned superLayerId(unsigned wireId) const
returns super layer ID. This function is expensive.
Definition: TRGCDC.cc:1820
void doit(std::vector< TRGCDCSegment * > &tss, const bool trackSegmentClockSimulation, std::vector< TRGCDCSegmentHit * > &segmentHits, std::vector< TRGCDCSegmentHit * > *segmentHitsSL)
Member functions of doing TSF.
static int level(void)
returns the debug level.
Definition: Debug.cc:67
unsigned axialStereoLayerId(void) const
returns id of axial or stereo id.
Definition: Layer.h:217
std::vector< const TRGCDCWireHit * > hits(void) const
returns a list of TRGCDCWireHit.
Definition: TRGCDC.cc:1705
int getReturnValue(void) const
gets return value for trg cdc module.
Definition: TRGCDC.h:1098
bool active(void) const
returns true if there is a signal.
Definition: Signal.h:277
const std::vector< const TRGCDCWire * > & wires(void) const
returns a vector containing pointers to a wire.
Definition: Segment.h:190
float priorityTime(void) const
return priority time in TSHit.
Definition: Segment.cc:463
const TRGCDCLayer * layer(unsigned id) const
returns a pointer to a layer. 0 will be returned if 'id' is invalid.
Definition: TRGCDC.h:828
void update()
updates TRGCDC wire information. clear() is called in this function.
Definition: TRGCDC.cc:984
int position(double timing) const
returns clock position.
Definition: Clock.cc:114
static void leaveStage(const std::string &stageName)
Declare that you leave a stage.
Definition: Debug.cc:34
unsigned nCells(void) const
returns # of cells.
Definition: Layer.h:194
unsigned nStereoSuperLayers(void) const
returns # of stereo super layers.
Definition: TRGCDC.h:898
void fastSimulation(void)
Fast simulation.
Definition: TRGCDC.cc:2064
float foundTime(void) const
return found time in TSHit.
Definition: Segment.cc:423
std::string version(void) const
returns version.
Definition: TRGCDC.cc:98
unsigned nWires(void) const
return # of wires.
Definition: TRGCDC.h:863
void simulate(void)
fast trigger simulation.
Definition: TRGCDC.cc:2049
void clear(void)
clears all TRGCDC hit information.
Definition: TRGCDC.cc:942
void saveCDCHitInformation(std::vector< std::vector< unsigned > > &)
Save functions for ROOT.
Definition: TRGCDC.cc:2760
void doit(std::vector< TRGCDCTrack * > const &trackList2D, std::vector< TRGCDCTrack * > &trackList3D)
Member functions.
void firmwareSimulation(void)
Firmware simulation.
Definition: TRGCDC.cc:2363
void classification(void)
classify hits.
Definition: TRGCDC.cc:1604
void updateByData(int inputMode)
updates TRGCDC wire information by Hardware data 0: From CDC FE ASCII file (Implementing) 1: From CDC...
Definition: TRGCDC.cc:1179
void storeSimulationResults(std::string collection2Dfinder, std::string collection2Dfitter, std::string collection3Dfitter)
Save results of fast simulation to data store (segment hits & tracks).
Definition: TRGCDC.cc:2248
void dump(const std::string &message="", const std::string &pre="") const
dumps contents.
Definition: Clock.cc:88
const TRGCDCWire * wire(unsigned wireId) const
returns a pointer to a wire.
Definition: TRGCDC.cc:892
Abstract base class for different kinds of events.