Belle II Software  release-05-01-25
Hough3DFinder.cc
1 //-----------------------------------------------------------------------------
2 // $Id$
3 //-----------------------------------------------------------------------------
4 // Filename : Hough3DFinder.cc
5 // Section : TRG CDC
6 // Owner : Jaebak Kim
7 // Email : jbkim@hep.korea.ac.kr
8 //-----------------------------------------------------------------------------
9 // Description : A class to find 3D tracks using Hough algorithm
10 //-----------------------------------------------------------------------------
11 // $Log$
12 //-----------------------------------------------------------------------------
13 
14 #define TRG_SHORT_NAMES
15 #define TRGCDC_SHORT_NAMES
16 
17 #include <map>
18 #include "TFile.h"
19 #include "TTree.h"
20 #include <framework/dataobjects/EventMetaData.h>
21 #include "framework/datastore/StoreArray.h"
22 #include "framework/datastore/RelationArray.h"
23 #include "cdc/dataobjects/CDCHit.h"
24 #include "cdc/dataobjects/CDCSimHit.h"
25 #include "cdc/geometry/CDCGeometryPar.h"
26 #include "trg/trg/Debug.h"
27 #include "trg/trg/Utilities.h"
28 #include "trg/cdc/TRGCDC.h"
29 #include "trg/cdc/Layer.h"
30 #include "trg/cdc/Cell.h"
31 #include "trg/cdc/Wire.h"
32 #include "trg/cdc/WireHit.h"
33 #include "trg/cdc/Hough3DFinder.h"
34 #include "trg/cdc/Segment.h"
35 #include "trg/cdc/SegmentHit.h"
36 #include "trg/cdc/Circle.h"
37 #include "trg/cdc/Track.h"
38 #include "trg/cdc/Link.h"
39 #include "trg/cdc/Relation.h"
40 #include "trg/cdc/Fitter3DUtility.h"
41 #include "trg/cdc/Fitter3D.h"
42 #include "trg/cdc/HandleRoot.h"
43 
44 using namespace std;
45 
46 namespace Belle2 {
52  TRGCDCHough3DFinder::TRGCDCHough3DFinder(const TRGCDC& TRGCDC, bool makeRootFile, int finderMode)
53  : _cdc(TRGCDC), m_makeRootFile(makeRootFile) , m_finderMode(finderMode)
54  {
55 
56  m_fileFinder3D = 0;
57  m_treeTrackFinder3D = 0;
58  m_treeConstantsFinder3D = 0;
59 
60  m_mBool["fMc"] = 1;
61  m_mBool["fVerbose"] = 0;
62  m_mBool["fIsPrintError"] = 0;
63  m_mBool["fIsIntegerEffect"] = 1;
64  m_mBool["fmcLR"] = 0;
65  m_mBool["fLRLUT"] = 1;
66  m_mBool["f2DFitDrift"] = 1;
67  m_mBool["f2DFit"] = 1;
68 
69  // For finder3D.
70  m_mBool["debugEfficiency"] = 1;
71  m_mBool["debugNTs"] = 1;
72 
73  m_mConstD["Trg_PI"] = 3.141592653589793;
74  // Get rr,zToStraw,angleSt,nWire
75  CDC::CDCGeometryPar& cdcp = CDC::CDCGeometryPar::Instance();
76  m_mConstV["rr"] = vector<double> (9);
77  m_mConstV["nWires"] = vector<double> (9);
78  m_mConstV["nTSs"] = vector<double> (9);
79  for (unsigned iSL = 0; iSL < 9; iSL++) {
80  unsigned t_layerId = _cdc.segment(iSL, 0).center().layerId();
81  m_mConstV["rr"][iSL] = cdcp.senseWireR(t_layerId);
82  m_mConstV["nWires"][iSL] = cdcp.nWiresInLayer(t_layerId) * 2;
83  m_mConstV["nTSs"][iSL] = cdcp.nWiresInLayer(t_layerId);
84  }
85  m_mConstV["nTSs2D"] = vector<double> (5);
86  for (unsigned iAx = 0; iAx < 5; iAx++) {
87  m_mConstV["nTSs2D"][iAx] = m_mConstV["nTSs"][2 * iAx];
88  }
89 
90  m_mConstV["zToStraw"] = vector<double> (4);
91  m_mConstV["zToOppositeStraw"] = vector<double> (4);
92  m_mConstV["angleSt"] = vector<double> (4);
93  m_mConstV["nShift"] = vector<double> (4);
94  for (int iSt = 0; iSt < 4; iSt++) {
95  unsigned t_layerId = _cdc.stereoSegment(iSt, 0).center().layerId();
96  m_mConstV["zToStraw"][iSt] = cdcp.senseWireBZ(t_layerId);
97  m_mConstV["zToOppositeStraw"][iSt] = cdcp.senseWireFZ(t_layerId);
98  m_mConstV["angleSt"][iSt] = 2 * m_mConstV["rr"][2 * iSt + 1] * sin(m_mConstD["Trg_PI"] * cdcp.nShifts(t_layerId) /
99  (2 * cdcp.nWiresInLayer(t_layerId))) / (cdcp.senseWireFZ(t_layerId) - cdcp.senseWireBZ(t_layerId));
100  m_mConstV["nShift"][iSt] = cdcp.nShifts(t_layerId);
101  }
102  m_mConstV["rr2D"] = vector<double> (5);
103  m_mConstV["rr3D"] = vector<double> (4);
104  for (int iAx = 0; iAx < 5; iAx++) m_mConstV["rr2D"][iAx] = m_mConstV["rr"][iAx * 2];
105  for (int iSt = 0; iSt < 4; iSt++) m_mConstV["rr3D"][iSt] = m_mConstV["rr"][iSt * 2 + 1];
106  m_mConstV["wirePhi2DError"] = vector<double> (5);
107  m_mConstV["driftPhi2DError"] = vector<double> (5);
108 
110  //m_mConstD["driftResolution"] = 2 * 40 * 0.0001; // (cm)
111  //m_mConstV["cellResolution"] = vector<double> (5);
112  //for(unsigned iAx=0; iAx<5; iAx++){
113  // m_mConstV["cellResolution"][iAx] = m_mConstV["rr2D"][iAx] * 2 * m_mConstD["Trg_PI"] / (m_mConstV["nWires"][iAx*2]/2); // (cm)
114  //}
115  //for(unsigned iAx=0; iAx<5; iAx++){
116  // m_mConstV["wirePhi2DError"][iAx] = m_mConstV["cellResolution"][iAx]/m_mConstV["rr2D"][iAx] / sqrt(12);
117  //}
118  //for(unsigned iAx=0; iAx<5; iAx++){
119  // m_mConstV["driftPhi2DError"][iAx] = m_mConstD["driftResolution"]/m_mConstV["rr2D"][iAx] / sqrt(12);
120  //}
121  m_mConstV["wirePhi2DError"][0] = 0.00085106;
122  m_mConstV["wirePhi2DError"][1] = 0.00039841;
123  m_mConstV["wirePhi2DError"][2] = 0.00025806;
124  m_mConstV["wirePhi2DError"][3] = 0.00019084;
125  m_mConstV["wirePhi2DError"][4] = 0.0001514;
126  m_mConstV["driftPhi2DError"][0] = 0.00085106;
127  m_mConstV["driftPhi2DError"][1] = 0.00039841;
128  m_mConstV["driftPhi2DError"][2] = 0.00025806;
129  m_mConstV["driftPhi2DError"][3] = 0.00019084;
130  m_mConstV["driftPhi2DError"][4] = 0.0001514;
135  //m_mConstV["wireZError"] = vector<double> (4);
136  //m_mConstV["wireZError"][0] = 3.19263;
137  //m_mConstV["wireZError"][1] = 2.8765;
138  //m_mConstV["wireZError"][2] = 2.90057;
139  //m_mConstV["wireZError"][3] = 3.96206;
140  //m_mConstV["driftZError"] = vector<double> (4);
141  //m_mConstV["driftZError"][0] = 3.19263;
142  //m_mConstV["driftZError"][1] = 2.8765;
143  //m_mConstV["driftZError"][2] = 2.90057;
144  //m_mConstV["driftZError"][3] = 3.96206;
145 
146  // (2016.06.07) study
147  //m_mConstV["wireZError"] = vector<double> ({4.752, 6.393, 6.578, 6.418});
148  //m_mConstV["driftZError"] = vector<double> ({0.4701, 0.7203, 0.8058, 0.9382});
149  m_mConstV["driftZError"] = vector<double> ({0.7676, 0.9753, 1.029, 1.372});
150  m_mConstV["wireZError"] = vector<double> ({0.7676, 0.9753, 1.029, 1.372});
151 
152  // Make driftLength table for each superlayer. Up to 511 clock ticks.
153  // driftLengthTableSLX[ tdcCount (~2ns unit) ] = drift length (cm)
154  for (unsigned iSl = 0; iSl < 9; iSl++) {
155  string tableName = "driftLengthTableSL" + to_string(iSl);
156  unsigned tableSize = 512;
157  m_mConstV[tableName] = vector<double> (tableSize);
158  unsigned t_layer = _cdc.segment(iSl, 0).center().layerId();
159  for (unsigned iTick = 0; iTick < tableSize; iTick++) {
160  double t_driftTime = iTick * 2 * cdcp.getTdcBinWidth();
161  double avgDriftLength = 0;
162  if (m_mBool["fXtSimple"] == 1) {
163  avgDriftLength = cdcp.getNominalDriftV() * t_driftTime;
164  } else {
165  double driftLength_0 = cdcp.getDriftLength(t_driftTime, t_layer, 0);
166  double driftLength_1 = cdcp.getDriftLength(t_driftTime, t_layer, 1);
167  avgDriftLength = (driftLength_0 + driftLength_1) / 2;
168  }
169  m_mConstV[tableName][iTick] = avgDriftLength;
170  }
171  }
172 
173  // Save geometry to root file.
174  TVectorD geometryHough3D(16);
175  for (int i = 0; i < 4; i++) {
176  geometryHough3D[i] = m_mConstV["rr"][2 * i + 1] / 100;
177  geometryHough3D[i + 4] = m_mConstV["angleSt"][i];
178  geometryHough3D[i + 8] = m_mConstV["zToStraw"][i] / 100;
179  geometryHough3D[i + 12] = m_mConstV["nWires"][2 * i + 1];
180  }
181 
182  m_Hough3DFinder = new Hough3DFinder();
183  // 1: Hough3DFinder 2: GeoFinder 3: VHDL GeoFinder
184  m_Hough3DFinder->setMode(m_finderMode);
185  // Set input file name for VHDL GeoFinder.
186  m_Hough3DFinder->setInputFileName(string(std::getenv("BELLE2_LOCAL_DIR")) + "/data/trg/cdc/GeoFinder.input");
187 
188  // For VHDL GEoFinder
189  //m_Hough3DFinder->setInputFileName("GeoFinder.input");
190  // cotStart, cotEnd, z0Start, z0End, cotSteps, z0Steps
191  float tempInitVariables[] = { -3, 3, -2, 2, 1001, 1001};
192  vector<float > initVariables(tempInitVariables, tempInitVariables + sizeof(tempInitVariables) / sizeof(tempInitVariables[0]));
193  // Save the init variables
194  m_Hough3DFinder->initialize(geometryHough3D, initVariables);
195 
196  m_mConstD["modeHough3D"] = m_Hough3DFinder->getMode();
197  m_mConstV["initVariablesHough3D"] = vector<double> (6);
198  m_mConstV["initVariablesHough3D"][0] = tempInitVariables[0];
199  m_mConstV["initVariablesHough3D"][1] = tempInitVariables[1];
200  m_mConstV["initVariablesHough3D"][2] = tempInitVariables[2];
201  m_mConstV["initVariablesHough3D"][3] = tempInitVariables[3];
202  m_mConstV["initVariablesHough3D"][4] = tempInitVariables[4];
203  m_mConstV["initVariablesHough3D"][5] = tempInitVariables[5];
204 
205  }
206 
207  TRGCDCHough3DFinder::~TRGCDCHough3DFinder()
208  {
209  }
210 
211  void TRGCDCHough3DFinder::terminate(void)
212  {
213  if (m_makeRootFile) {
214  HandleRoot::writeRoot(m_fileFinder3D);
215  HandleRoot::terminateRoot(m_mRunTVectorD, m_mEventTVectorD, m_mTClonesArray);
216  delete m_fileFinder3D;
217  }
218  m_Hough3DFinder->destruct();
219  }
220 
221  void TRGCDCHough3DFinder::doit(vector<TCTrack*> const& trackList2D, vector<TCTrack*>& trackList3D)
222  {
223  // Loop over trackList2D and copy to make a new trackList3D. Will delete it at TRGCDC.cc.
224  for (unsigned int iTrack = 0; iTrack < trackList2D.size(); iTrack++) {
225  TCTrack& aTrack = * new TCTrack(* trackList2D[iTrack]);
226  trackList3D.push_back(& aTrack);
227  }
228  doit(trackList3D);
229  }
230 
231  void TRGCDCHough3DFinder::doit(vector<TCTrack*>& trackList)
232  {
233 
234  // Assign track ID's.
235  for (unsigned iTrack = 0; iTrack < trackList.size(); iTrack++) {
236  TCTrack& aTrack = *trackList[iTrack];
237  aTrack.setTrackID(iTrack + 1);
238  }
239 
240  if (m_finderMode == 0) doitPerfectly(trackList);
241  if (m_finderMode != 0) doitFind(trackList);
242 
243  }
244 
245  void TRGCDCHough3DFinder::doitFind(vector<TCTrack*>& trackList)
246  {
247  TRGDebug::enterStage("3D finder");
248 
249  // For saving to root file.
250  if (m_makeRootFile) m_mDouble["iSave"] = 0;
251  if (m_makeRootFile) HandleRoot::initializeEvent(m_mEventTVectorD, m_mTClonesArray);
252 
253  // Get event number.
254  StoreObjPtr<EventMetaData> eventMetaDataPtr;
255  // Event starts from 0.
256  m_mDouble["eventNumber"] = eventMetaDataPtr->getEvent();;
257 
258  // Generate arrays for TS candidates.
259  vector<vector<double> > stTSs(4);
260  vector<vector<int> > stTSDrift(4);
261  vector<vector<const TCSHit*> > p_stTSs(4);
262  for (unsigned iSL = 0; iSL < 4; iSL++) {
263  vector<const TCSHit*> hits = _cdc.stereoSegmentHits(iSL);
264  // Initialize vectors.
265  string slName = "st" + to_string(iSL);
266  m_mEventV[slName + "_hit"] = vector<double> ();
267  m_mEventV[slName + "_driftHit"] = vector<double> ();
268  stTSs[iSL] = vector<double> ();
269  stTSDrift[iSL] = vector<int> ();
270  p_stTSs[iSL] = vector<const TCSHit*> ();
271  // Fill vectors
272  for (unsigned iHit = 0; iHit < hits.size(); iHit++) {
273  if (hits[iHit] == 0) continue;
276  //if(hits[iHit]->segment().priorityPosition() != 3) {
277  // continue;
278  //}
279  double t_wirePhi = ((double)hits[iHit]->cell().localId()) / m_mConstV["nWires"][2 * iSL + 1] * 4 * m_mConstD["Trg_PI"];
280  m_mEventV[slName + "_hit"].push_back(t_wirePhi);
281  m_mEventV[slName + "_driftHit"].push_back(TRGCDCFitter3D::calPhi(hits[iHit], _cdc.getEventTime()));
282  int t_tdc = hits[iHit]->segment().priorityTime();
283  int t_lr = hits[iHit]->segment().LUT()->getValue(hits[iHit]->segment().lutPattern());;
284  int t_priorityPosition = hits[iHit]->segment().priorityPosition();
285  int t_driftInfo = (t_tdc << 4) + (t_lr << 2) + t_priorityPosition;
286  stTSs[iSL].push_back(t_wirePhi);
287  p_stTSs[iSL].push_back(hits[iHit]);
288  stTSDrift[iSL].push_back(t_driftInfo);
289  //cout<<"iSL:"<<iSL<<" iHit:"<<iHit<<" t_tdc:"<<t_tdc<<" t_lr:"<<t_lr<<" t_priorityPosition:"<<t_priorityPosition<<" t_driftInfo:"<<t_driftInfo<<endl;
290  }
291  }
292 
293  // Get MC values related with finding
294  // numberTSsForParticle[mcId] = # superlayer hits.
295  map<unsigned, unsigned> numberTSsForParticle;
296  if (m_mBool["fMc"]) findNumberOfHitSuperlayersForMcParticles(p_stTSs, numberTSsForParticle);
297 
298  // Loop over all the tracks.
299  m_mEventD["nTracks"] = trackList.size();
300  for (unsigned iTrack = 0; iTrack < m_mEventD["nTracks"]; iTrack++) {
301 
302  TCTrack& aTrack = * trackList[iTrack];
303 
304  // Get MC values related with fitting
305  if (m_mBool["fMc"]) TRGCDCFitter3D::getMCValues(_cdc, &aTrack, m_mConstD, m_mDouble, m_mVector);
306 
307  // Get track ID
308  m_mDouble["trackId"] = aTrack.getTrackID();
309 
310  // 2D Fitter
311  int fit2DResult = TRGCDCFitter3D::do2DFit(aTrack, m_mBool, m_mConstD, m_mConstV, m_mDouble, m_mVector);
312  if (fit2DResult != 0) continue;
313 
314  // Set input of finder
315  //vector<double > trackVariables = { m_mDouble["charge"], m_mDouble["rho"]/100, m_mDouble["phi0"] } ;
316  vector<double > trackVariables = { m_mDouble["charge2D"], m_mDouble["rho"] / 100, m_mDouble["phi0"] } ;
317 
318  // Run finder
319  m_Hough3DFinder->runFinder(trackVariables, stTSs, stTSDrift);
320 
321  // Get results of finder
322  m_Hough3DFinder->getValues("bestTSIndex", m_mVector["bestTSIndex"]);
323  const TCSHit* p_bestTS[4] = {0, 0, 0, 0};
324  for (int iSt = 0; iSt < 4; iSt++) {
325  if (m_mVector["bestTSIndex"][iSt] == 999) p_bestTS[iSt] = 0;
326  else p_bestTS[iSt] = p_stTSs[iSt][(int)m_mVector["bestTSIndex"][iSt]];
327  }
328  m_Hough3DFinder->getValues("bestTS", m_mVector["bestTS"]);
329  // Find and append TS to track.
330  for (unsigned iSt = 0; iSt < 4; iSt++) {
331  if (m_mVector["bestTS"][iSt] != 999) aTrack.append(new TCLink(0, p_bestTS[iSt], p_bestTS[iSt]->cell().xyPosition()));
332  }
333 
334  // For saving values from finder.
335  if (m_Hough3DFinder->getMode() == 1) {
336  m_Hough3DFinder->getValues("bestZ0", m_mVector["bestZ0"]);
337  m_Hough3DFinder->getValues("bestCot", m_mVector["bestCot"]);
338  m_Hough3DFinder->getValues("houghMax", m_mVector["houghMax"]);
339  m_Hough3DFinder->getValues("minDiffHough", m_mVector["minDiffHough"]);
340  }
341  if (m_Hough3DFinder->getMode() == 2) {
342  m_Hough3DFinder->getValues("st0GeoCandidatesPhi", m_mVector["st0GeoCandidatesPhi"]);
343  m_Hough3DFinder->getValues("st1GeoCandidatesPhi", m_mVector["st1GeoCandidatesPhi"]);
344  m_Hough3DFinder->getValues("st2GeoCandidatesPhi", m_mVector["st2GeoCandidatesPhi"]);
345  m_Hough3DFinder->getValues("st3GeoCandidatesPhi", m_mVector["st3GeoCandidatesPhi"]);
346  m_Hough3DFinder->getValues("st0GeoCandidatesDiffStWires", m_mVector["st0GeoCandidatesDiffStWires"]);
347  m_Hough3DFinder->getValues("st1GeoCandidatesDiffStWires", m_mVector["st1GeoCandidatesDiffStWires"]);
348  m_Hough3DFinder->getValues("st2GeoCandidatesDiffStWires", m_mVector["st2GeoCandidatesDiffStWires"]);
349  m_Hough3DFinder->getValues("st3GeoCandidatesDiffStWires", m_mVector["st3GeoCandidatesDiffStWires"]);
350  m_Hough3DFinder->getValues("stAxPhi", m_mVector["stAxPhi"]);
351  }
352 
353  // Get MC values.
354  if (m_mBool["fMc"]) {
355  // Call storage array.
356  StoreArray<CDCSimHit> SimHits("CDCSimHits");
357  StoreArray<CDCHit> CDCHits("CDCHits");
358  RelationArray relationCDCHits(SimHits, CDCHits);
359 
360  // save performance values (purity, efficiency)
361  const TCRelation& trackRelation3D = aTrack.relation3D();
362  m_mDouble["purity"] = trackRelation3D.purity3D(aTrack.relation2D().contributor(0));
363  m_mDouble["efficiency"] = trackRelation3D.efficiency3D(aTrack.relation2D().contributor(0), numberTSsForParticle);
364  // Find one mc stereo track segment per layer.
365  m_mVector["mcTSs"] = vector<double> (4, 999);
366  vector<const TCSHit*> mcTSList;
367  perfectFinder(trackList, iTrack, mcTSList);
368  for (unsigned iTS = 0; iTS < mcTSList.size(); iTS++) {
369  int iSuperLayer = (int)(double(mcTSList[iTS]->cell().superLayerId()) - 1) / 2;
370  m_mVector["mcTSs"][iSuperLayer] = (double)mcTSList[iTS]->cell().localId() / m_mConstV["nWires"][2 * iSuperLayer + 1] * 4 *
371  m_mConstD["Trg_PI"];
372  }
373  // Save MC ture CDC's hit position
374  m_mVector["mcTSsX"] = vector<double> (4);
375  m_mVector["mcTSsY"] = vector<double> (4);
376  for (unsigned iTS = 0; iTS < mcTSList.size(); iTS++) {
377  unsigned iCDCSimHit = mcTSList[iTS]->iCDCSimHit();
378  CDCSimHit* aCDCSimHit = SimHits[iCDCSimHit];
379  TVector3 posWire = aCDCSimHit->getPosWire();
380  m_mVector["mcTSsX"][iTS] = posWire.X();
381  m_mVector["mcTSsY"][iTS] = posWire.Y();
382  }
383  // Calculate diff from perfect
384  m_mVector["perfectWireDiff"] = vector<double> (4);
385  m_mVector["perfectCalZ"] = vector<double> (4);
386  for (unsigned iSt = 0; iSt < 4; iSt++) {
387  if (m_mVector["mcTSs"][iSt] == 999) {
388  m_mVector["perfectWireDiff"][iSt] = 999;
389  m_mVector["perfectCalZ"][iSt] = 999;
390  } else {
391  m_mVector["perfectWireDiff"][iSt] = Fitter3DUtility::calDeltaPhi(m_mDouble["mcCharge"], m_mConstV["angleSt"][iSt],
392  m_mConstV["zToStraw"][iSt] / 100, m_mConstV["rr"][2 * iSt + 1] / 100, m_mVector["mcTSs"][iSt], m_mDouble["rho"] / 100,
393  m_mDouble["phi0"]) / 4 / m_mConstD["Trg_PI"] * m_mConstV["nWires"][2 * iSt + 1];
394  m_mVector["perfectCalZ"][iSt] = Fitter3DUtility::calZ(m_mDouble["mcCharge"], m_mConstV["angleSt"][iSt],
395  m_mConstV["zToStraw"][iSt] / 100, m_mConstV["rr"][2 * iSt + 1] / 100, m_mVector["mcTSs"][iSt], m_mDouble["rho"] / 100,
396  m_mDouble["phi0"]);
397  }
398  }
399  // Find multiple mc stereo track segments per layer.
400  unsigned int mcParticleId = aTrack.relation2D().contributor(0);
401  for (unsigned iSt = 0; iSt < 4; iSt++) {
402  string mcTsName = "mcTsSt" + to_string(iSt);
403  m_mVector[mcTsName] = vector<double> ();
404  vector<const TCSHit*> hits = _cdc.stereoSegmentHits(iSt);
405  // Find TS which match mc index
406  for (unsigned iTS = 0; iTS < hits.size(); iTS++) {
407  if (hits[iTS]->iMCParticle() == mcParticleId) m_mVector[mcTsName].push_back((double)hits[iTS]->cell().localId() /
408  m_mConstV["nWires"][2 * iSt + 1] * 4 * m_mConstD["Trg_PI"]);
409  }
410  }
411 
412  //cout<<"----newStart----"<<endl;
413  //cout<<"mcCharge:"<<m_mDouble["mcCharge"]<<endl;
414  //cout<<"rho:"<<m_mDouble["rho"]<<" phi0:"<<m_mDouble["phi0"]<<endl;
415  //cout<<"purity:"<<m_mDouble["purity"]<<" efficiency:"<<m_mDouble["efficiency"]<<endl;
416  //cout<<"mcTSs[0]:"<<m_mVector["mcTSs"][0]<<" mcTSs[1]:"<<m_mVector["mcTSs"][1]<<" mcTSs[2]:"<<m_mVector["mcTSs"][2]<<" mcTSs[3]:"<<m_mVector["mcTSs"][3]<<endl;
417  //cout<<"mcTSsX[0]:"<<m_mVector["mcTSsX"][0]<<" mcTSsX[1]:"<<m_mVector["mcTSsX"][1]<<" mcTSsX[2]:"<<m_mVector["mcTSsX"][2]<<" mcTSsX[3]:"<<m_mVector["mcTSsX"][3]<<endl;
418  //cout<<"mcTSsY[0]:"<<m_mVector["mcTSsY"][0]<<" mcTSsY[1]:"<<m_mVector["mcTSsY"][1]<<" mcTSsY[2]:"<<m_mVector["mcTSsY"][2]<<" mcTSsY[3]:"<<m_mVector["mcTSsY"][3]<<endl;
419  //cout<<"perfectWireDiff[0]:"<<m_mVector["perfectWireDiff"][0]<<" perfectWireDiff[1]:"<<m_mVector["perfectWireDiff"][1]<<" perfectWireDiff[2]:"<<m_mVector["perfectWireDiff"][2]<<" perfectWireDiff[3]:"<<m_mVector["perfectWireDiff"][3]<<endl;
420  //cout<<"perfectCalZ[0]:"<<m_mVector["perfectCalZ"][0]<<" perfectCalZ[1]:"<<m_mVector["perfectCalZ"][1]<<" perfectCalZ[2]:"<<m_mVector["perfectCalZ"][2]<<" perfectCalZ[3]:"<<m_mVector["perfectCalZ"][3]<<endl;
421  //cout<<"----newEnd----"<<endl;
422  }
423 
424  if (m_makeRootFile) {
425  if (m_fileFinder3D == 0) {
426  m_fileFinder3D = new TFile("Finder3D.root", "RECREATE");
427  HandleRoot::initializeRoot("hough3D", &m_treeConstantsFinder3D, &m_treeTrackFinder3D,
428  m_mRunTVectorD, m_mEventTVectorD, m_mTClonesArray,
429  m_mConstD, m_mConstV,
430  m_mEventD, m_mEventV,
431  m_mDouble, m_mVector
432  );
433  }
434  HandleRoot::saveTrackValues("hough3D",
435  m_mTClonesArray, m_mDouble, m_mVector
436  );
437  }
438 
439  } // End of loop over all the tracks.
440 
441  // Will ignore events until event has a track.
442  if (m_makeRootFile && m_fileFinder3D) {
443  HandleRoot::saveEventValues("hough3D",
444  m_mEventTVectorD, m_mEventD, m_mEventV
445  );
446  m_treeTrackFinder3D->Fill();
447  }
448 
449 
450  // Set debug values.
451  if (m_mBool["debugEfficiency"]) {
452  // Notify when efficiency is not 1.
453  for (unsigned iTrack = 0; iTrack < trackList.size(); iTrack++) {
454  /* cppcheck-suppress variableScope */
455  TCTrack& aTrack = * trackList[iTrack];
456  // Find number of super layers that have priority layer hit.
457  int nPriorityHitSL = 0;
458  for (unsigned iSt = 0; iSt < 4; iSt++) {
459  int priorityHitSL = 0;
460  for (unsigned iTS = 0; iTS < stTSDrift[iSt].size(); iTS++) {
461  int t_priorityPosition = (stTSDrift[iSt][iTS] & 3);
462  if (t_priorityPosition == 3) priorityHitSL = 1;
463  //cout<<"iSt:"<<iSt<<" iTS:"<<iTS<<" priorityPosition:"<<t_priorityPosition<<" priorityHitSL:"<<priorityHitSL<<endl;
464  }
465  if (priorityHitSL) nPriorityHitSL++;
466  }
467  //cout<<"nPriorityHitSL:"<<nPriorityHitSL<<endl;
468  if (m_mDouble["efficiency"] != 1) {
469  // Remove case when all hits are secondary priority.
470  if (m_mDouble["efficiency"] != nPriorityHitSL * 1. / 4) {
471  aTrack.setDebugValue(TRGCDCTrack::EDebugValueType::find3D, 1);
472  }
473  }
474  }
475  }
476  if (m_mBool["debugNTs"]) {
477  // Notify when not enough TS are found
478  for (unsigned iTrack = 0; iTrack < trackList.size(); iTrack++) {
479  unsigned nHitStSl = 0;
480  for (unsigned iSt = 0; iSt < 4; iSt++) {
481  if (trackList[iTrack]->links(2 * iSt + 1).size() != 0) nHitStSl++;
482  }
483  if (nHitStSl < 2) trackList[iTrack]->setDebugValue(TRGCDCTrack::EDebugValueType::find3D, 1);
484  }
485  }
486 
487  TRGDebug::leaveStage("3D finder");
488 
489  }
490 
491  void TRGCDCHough3DFinder::perfectFinder(vector<TCTrack*>& trackList, unsigned j, vector<const TCSHit*>& mcTSList)
492  {
493 
494  //Just a test
495  StoreArray<CDCHit> CDCHits("CDCHits");
496  StoreArray<CDCSimHit> SimHits("CDCSimHits");
497  RelationArray rels(SimHits, CDCHits);
498 
499  //...G4 trackID...
500  unsigned id = trackList[j]->relation().contributor(0);
501  vector<const TCSHit*> tsList[9];
502  //cout<<"[JB] id: "<<id<<endl;
503 
504  //...Segment loop...
505  const vector<const TCSHit*> hits = _cdc.segmentHits();
506  for (unsigned i = 0; i < hits.size(); i++) {
507  const TCSHit& ts = * hits[i];
508  if (ts.segment().axial()) continue;
509  if (! ts.signal().active()) continue;
510  const TCWHit* wh = ts.segment().center().hit();
511  if (! wh) continue;
512  const unsigned trackId = wh->iMCParticle();
513  // Try to track down the mcParticle another way.
514  //cout<<"[CDCTRG] trackId:"<<trackId<<" "<<ts.cell().name()<<endl;
515 //iw commented out because simind is not used
516 // int ind=wh->iCDCHit();
517 // int simind=rels[ind].getFromIndex();
518  //CDCSimHit &h=*SimHits[simind];
519  //cout<<"[CDCTRG] simTrackId: "<<wh->simHit()->getTrackId()<<endl;;
520  //cout<<"[CDC] simTrackId: "<<h.getTrackId()<<" from CDCHit: "<<CDCHits[ind]->getIWire()<<endl;;
521 
522  if (id == trackId)
523  tsList[wh->wire().superLayerId()].push_back(& ts);
524  }
525 
526  if (TRGDebug::level()) {
527  for (unsigned k = 0; k < 9; k++) {
528  if (k % 2) {
529  cout << TRGDebug::tab(4) << "superlayer " << k << ":";
530  for (unsigned l = 0; l < tsList[k].size(); l++) {
531  if (l)
532  cout << ",";
533  cout << tsList[k][l]->cell().name();
534  }
535  cout << endl;
536  }
537  }
538  }
539 
540  //...Select best one in each super layer...
541  for (unsigned i = 0; i < 9; i++) {
542  const TCSHit* best = 0;
543  if (tsList[i].size() == 0) {
544  continue;
545  } else if (tsList[i].size() == 1) {
546  best = tsList[i][0];
547  } else {
548  int timeMin = 99999;
549  for (unsigned k = 0; k < tsList[i].size(); k++) {
550  const TRGSignal& timing = tsList[i][k]->signal();
551  const TRGTime& t = * timing[0];
552  if (t.time() < timeMin) {
553  timeMin = t.time();
554  best = tsList[i][k];
555  }
556  }
557  }
558  mcTSList.push_back(best);
559  }
560 
561  }
562 
563 
564  void TRGCDCHough3DFinder::doitPerfectly(vector<TCTrack*>& trackList)
565  {
566 
567  TRGDebug::enterStage("Perfect 3D Finder");
568  if (TRGDebug::level())
569  cout << TRGDebug::tab() << "givenTrk#=" << trackList.size() << endl;
570 
571 
572  //...Track loop....
573  for (unsigned j = 0; j < trackList.size(); j++) {
574  //...G4 trackID...
575  TCTrack* trk = trackList[j];
576 
577  vector<const TCSHit*> mcTSList;
578  perfectFinder(trackList, j, mcTSList);
579  for (unsigned iTS = 0; iTS < mcTSList.size(); iTS++) {
580  trk->append(new TCLink(0, mcTSList[iTS], mcTSList[iTS]->cell().xyPosition()));
581  }
582 
583  if (TRGDebug::level())
584  trk->dump("", "> ");
585 
586  }
587 
588  TRGDebug::leaveStage("Perfect 3D Finder");
589 
590  }
591 
592  void TRGCDCHough3DFinder::findNumberOfHitSuperlayersForMcParticles(vector<vector<const TCSHit*> >& p_stTSs,
593  map<unsigned, unsigned>& numberTSsForParticle)
594  {
595  vector<unsigned> mcParticleList;
596  for (unsigned iLayer = 0; iLayer < 4; iLayer++) {
597  // Find what mc particles there are in a layer
598  mcParticleList.clear();
599  for (unsigned iTS = 0; iTS < p_stTSs[iLayer].size(); iTS++) {
600  unsigned iMCParticle = p_stTSs[iLayer][iTS]->iMCParticle();
601  if (find(mcParticleList.begin(), mcParticleList.end(), iMCParticle) == mcParticleList.end()) {
602  mcParticleList.push_back(iMCParticle);
603  }
604  }
605  // Loop over mcParticleList and add to numberTSsForParticle
606  for (unsigned iMCPart = 0; iMCPart < mcParticleList.size(); iMCPart++) {
607  map<unsigned, unsigned>::iterator it = numberTSsForParticle.find(mcParticleList[iMCPart]);
608  if (it != numberTSsForParticle.end()) ++it->second;
609  else numberTSsForParticle[mcParticleList[iMCPart]] = 1;
610  }
611  }
612  }
613 
614 
616 }
617 
Belle2::TRGSignal
A class to represent a digitized signal. Unit is nano second.
Definition: Signal.h:28
Belle2::RelationArray
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:72
Belle2::TRGCDCFitter3D::getMCValues
static void getMCValues(const TRGCDC &m_cdc_in, TRGCDCTrack *aTrack, std::map< std::string, double > &m_mConstD_in, std::map< std::string, double > &m_mDouble_in, std::map< std::string, std::vector< double > > &m_mVector_in)
Function for mc debugging.
Definition: Fitter3D.cc:961
Belle2::TRGCDCHough3DFinder::m_mEventD
std::map< std::string, double > m_mEventD
Map to hold event values for Fitter3D.
Definition: Hough3DFinder.h:89
Hough3DFinder::runFinder
void runFinder(std::vector< double > &trackVariables, std::vector< std::vector< double > > &stTSs, std::vector< std::vector< int > > &stTSDrift)
Track variables.
Definition: Hough3DUtility.cc:128
Hough3DFinder::setInputFileName
void setInputFileName(std::string inputFileName)
Sets the config file for the GeoFinder.
Definition: Hough3DUtility.cc:255
Belle2::TRGCDCHough3DFinder::~TRGCDCHough3DFinder
~TRGCDCHough3DFinder()
Destructor.
Definition: Hough3DFinder.cc:207
Belle2::TRGCDCHough3DFinder::m_makeRootFile
bool m_makeRootFile
Choose whether to save root file.
Definition: Hough3DFinder.h:65
Belle2::TRGCDCHough3DFinder::m_mTClonesArray
std::map< std::string, TClonesArray * > m_mTClonesArray
TClonesArray map for saving track values to root file.
Definition: Hough3DFinder.h:99
Belle2::TRGCDC::stereoSegmentHits
std::vector< const TRGCDCSegmentHit * > stereoSegmentHits(unsigned) const
returns a list of TRGCDCSegmentHit in a stereo super layer N.
Definition: TRGCDC.h:1024
Belle2::TRGCDC::segmentHits
std::vector< const TRGCDCSegmentHit * > segmentHits(void) const
returns a list of TRGCDCSegmentHit.
Definition: TRGCDC.h:997
Belle2::TRGDebug::level
static int level(void)
returns the debug level.
Definition: Debug.cc:72
Belle2::TRGCDC::stereoSegment
const TRGCDCSegment & stereoSegment(unsigned lyrId, unsigned id) const
returns a track segment in stereo layers. (lyrId is stereo #)
Definition: TRGCDC.h:948
Belle2::CDC::CDCGeometryPar::senseWireBZ
double senseWireBZ(int layerId) const
Returns backward z position of sense wire in each layer.
Definition: CDCGeometryPar.h:1241
Fitter3DUtility::calDeltaPhi
static double calDeltaPhi(int charge, double anglest, double ztostraw, double rr, double phi, double rho, double phi0)
Calculates the phi difference between fitted axial phi and stereo phi.
Definition: Fitter3DUtility.cc:808
Belle2::TRGCDCHough3DFinder::m_mBool
std::map< std::string, bool > m_mBool
Map to hold input options.
Definition: Hough3DFinder.h:93
Belle2::TRGCDCHough3DFinder::perfectFinder
void perfectFinder(std::vector< TRGCDCTrack * > &trackList, unsigned j, std::vector< const TRGCDCSegmentHit * > &mcTSList)
Perfect 3D finder for a track.
Definition: Hough3DFinder.cc:491
Belle2::CDC::CDCGeometryPar::getTdcBinWidth
double getTdcBinWidth() const
Return TDC bin width (nsec).
Definition: CDCGeometryPar.h:732
Belle2::TRGCDCHough3DFinder::m_mConstV
std::map< std::string, std::vector< double > > m_mConstV
Map to hold run vectcors for Fitter3D.
Definition: Hough3DFinder.h:87
Belle2::TRGCDCHough3DFinder::m_Hough3DFinder
Hough3DFinder * m_Hough3DFinder
Hough Variables.
Definition: Hough3DFinder.h:71
Belle2::CDC::CDCGeometryPar::senseWireFZ
double senseWireFZ(int layerId) const
Returns forward z position of sense wire in each layer.
Definition: CDCGeometryPar.h:1236
Belle2::CDCSimHit
Example Detector.
Definition: CDCSimHit.h:33
Hough3DFinder::initialize
void initialize(TVectorD &geometryVariables, std::vector< float > &initVariables)
Geometry variables.
Definition: Hough3DUtility.cc:63
Belle2::TRGCDCHough3DFinder::doitPerfectly
void doitPerfectly(std::vector< TRGCDCTrack * > &trackList)
Perfect 3D finder for a tracklist.
Definition: Hough3DFinder.cc:564
Belle2::TRGCDCHough3DFinder::m_mRunTVectorD
std::map< std::string, TVectorD * > m_mRunTVectorD
TVectorD map for saving run values to root file.
Definition: Hough3DFinder.h:95
Fitter3DUtility::calZ
static double calZ(int charge, double anglest, double ztostraw, double rr, double phi, double rho, double phi0)
Calculates z.
Definition: Fitter3DUtility.cc:896
Belle2::TRGCDCFitter3D::do2DFit
static int do2DFit(TRGCDCTrack &aTrack, std::map< std::string, bool > &m_mBool_in, std::map< std::string, double > &m_mConstD_in, std::map< std::string, std::vector< double > > &m_mConstV_in, std::map< std::string, double > &m_mDouble_in, std::map< std::string, std::vector< double > > &m_mVector_in)
Does 2D fit. Returns 0 if fit is done successfully. m_mBool should have fIsPrintError,...
Definition: Fitter3D.cc:1257
Belle2::CDC::CDCGeometryPar::getDriftLength
double getDriftLength(double dt, unsigned short layer, unsigned short lr, double alpha=0., double theta=0.5 *M_PI, bool calculateMinTime=true, double minTime=0.) const
Return the drift dength to the sense wire.
Definition: CDCGeometryPar.cc:2168
Belle2::CDC::CDCGeometryPar
The Class for CDC Geometry Parameters.
Definition: CDCGeometryPar.h:75
Belle2::TRGCDCHough3DFinder::doitFind
void doitFind(std::vector< TRGCDCTrack * > &trackList)
Finds tracks using tracklist.
Definition: Hough3DFinder.cc:245
Belle2::TRGCDCHough3DFinder::m_mEventTVectorD
std::map< std::string, TVectorD * > m_mEventTVectorD
TVectorD map for saving event values to root file.
Definition: Hough3DFinder.h:97
Hough3DFinder::getMode
int getMode(void)
Gets which 3D finder is used.
Definition: Hough3DUtility.cc:58
Belle2::CDC::CDCGeometryPar::senseWireR
double senseWireR(int layerId) const
Returns radius of sense wire in each layer.
Definition: CDCGeometryPar.h:1231
Belle2::TRGCDCHough3DFinder::m_finderMode
int m_finderMode
0: perfect finder, 1: Hough3DFinder, 2: GeoFinder, 3: VHDL GeoFinder Choose what finder to use.
Definition: Hough3DFinder.h:68
Belle2::CDC::CDCGeometryPar::Instance
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
Definition: CDCGeometryPar.cc:41
Hough3DFinder
A class to finded stereo TS hits related to 2D tracks.
Definition: Hough3DUtility.h:18
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TRGCDCHough3DFinder::_cdc
const TRGCDC & _cdc
Members.
Definition: Hough3DFinder.h:63
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
Hough3DFinder::destruct
void destruct(void)
Destructs the 3D finder.
Definition: Hough3DUtility.cc:111
Belle2::TRGCDCHough3DFinder::m_mDouble
std::map< std::string, double > m_mDouble
Map to hold track double values for Fitter3D.
Definition: Hough3DFinder.h:81
Belle2::TRGCDCFitter3D::calPhi
static double calPhi(TRGCDCSegmentHit const *segmentHit, double eventTime)
Utility functions.
Definition: Fitter3D.cc:893
Belle2::TRGCDCHough3DFinder::m_mEventV
std::map< std::string, std::vector< double > > m_mEventV
Map to hold event vectcors for Fitter3D.
Definition: Hough3DFinder.h:91
Belle2::CDC::CDCGeometryPar::nShifts
int nShifts(int layerId) const
Returns number shift.
Definition: CDCGeometryPar.h:1201
Belle2::CDC::CDCGeometryPar::getNominalDriftV
double getNominalDriftV() const
Return the nominal drift velocity of He-ethane gas (default: 4.0x10^-3 cm/nsec).
Definition: CDCGeometryPar.h:742
Belle2::TRGCDC::getEventTime
int getEventTime(void) const
returns bad hits(finding invalid hits).
Definition: TRGCDC.cc:2740
Belle2::TRGCDC::segment
const TRGCDCSegment & segment(unsigned id) const
returns a track segment.
Definition: TRGCDC.h:934
Belle2::TRGCDC
The instance of TRGCDC is a singleton.
Definition: TRGCDC.h:70
Belle2::CDCSimHit::getPosWire
TVector3 getPosWire() const
The method to get position on wire.
Definition: CDCSimHit.h:216
Belle2::TRGCDCHough3DFinder::m_treeConstantsFinder3D
TTree * m_treeConstantsFinder3D
TTree for constants of Hough3D.
Definition: Hough3DFinder.h:79
Belle2::TRGCDCSegment::center
const TRGCDCWire & center(void) const
returns a center wire.
Definition: Segment.h:250
Belle2::TRGCDCHough3DFinder::m_mVector
std::map< std::string, std::vector< double > > m_mVector
Map to hold track vector values for Fitter3D.
Definition: Hough3DFinder.h:83
Belle2::TRGDebug::tab
static std::string tab(void)
returns tab spaces.
Definition: Debug.cc:52
Belle2::CDC::CDCGeometryPar::nWiresInLayer
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
Definition: CDCGeometryPar.h:1211
Belle2::TRGCDCHough3DFinder::m_fileFinder3D
TFile * m_fileFinder3D
Tfile for Hough3D root file.
Definition: Hough3DFinder.h:75
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::TRGCDCHough3DFinder::terminate
void terminate(void)
Termination method.
Definition: Hough3DFinder.cc:211
Hough3DFinder::setMode
void setMode(int mode)
Sets which 3D finder to use.
Definition: Hough3DUtility.cc:53
Belle2::TRGCDCCell::layerId
unsigned layerId(void) const
returns layer id.
Definition: Cell.h:212
Belle2::TRGDebug::leaveStage
static void leaveStage(const std::string &stageName)
Declare that you leave a stage.
Definition: Debug.cc:39
Belle2::TRGCDCHough3DFinder::m_mConstD
std::map< std::string, double > m_mConstD
Map to hold run values for Fitter3D.
Definition: Hough3DFinder.h:85
Belle2::TRGDebug::enterStage
static void enterStage(const std::string &stageName)
Declare that you enter new stage.
Definition: Debug.cc:29
Belle2::TRGTime
A class to represent a signal timing in the trigger system.
Definition: Time.h:30
Hough3DFinder::getValues
void getValues(const std::string &input, std::vector< double > &result)
Gets results from the 3D finder.
Definition: Hough3DUtility.cc:1132
Belle2::TRGCDCHough3DFinder::doit
void doit(std::vector< TRGCDCTrack * > const &trackList2D, std::vector< TRGCDCTrack * > &trackList3D)
Member functions.
Belle2::TRGCDCHough3DFinder::findNumberOfHitSuperlayersForMcParticles
void findNumberOfHitSuperlayersForMcParticles(std::vector< std::vector< const TRGCDCSegmentHit * > > &p_stTSs, std::map< unsigned, unsigned > &numberTSsForParticle)
Finds number of hit superlayers for each mc particle.
Definition: Hough3DFinder.cc:592
Belle2::TRGCDCHough3DFinder::m_treeTrackFinder3D
TTree * m_treeTrackFinder3D
TTree for tracks of Hough3D.
Definition: Hough3DFinder.h:77