Belle II Software  release-06-02-00
Fitter3D.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 fit tracks in 3D
11 //-----------------------------------------------------------------------------
12 
13 #define TRG_SHORT_NAMES
14 #define TRGCDC_SHORT_NAMES
15 
16 #include <iostream>
17 #include "trg/trg/Debug.h"
18 #include "trg/cdc/Fitter3D.h"
19 #include "trg/cdc/Segment.h"
20 #include "trg/cdc/TRGCDCTrack.h"
21 #include "trg/cdc/Link.h"
22 #include <cmath>
23 
24 #include <framework/dataobjects/EventMetaData.h>
25 #include "cdc/dataobjects/CDCSimHit.h"
26 #include "cdc/geometry/CDCGeometryPar.h"
27 #include "trg/trg/Time.h"
28 #include "trg/trg/Signal.h"
29 #include "trg/trg/Utilities.h"
30 #include "trg/cdc/TRGCDC.h"
31 #include "trg/cdc/Wire.h"
32 #include "trg/cdc/Layer.h"
33 #include "trg/cdc/WireHit.h"
34 #include "trg/cdc/WireHitMC.h"
35 #include "trg/cdc/SegmentHit.h"
36 #include "trg/cdc/TrackMC.h"
37 #include "trg/cdc/Relation.h"
38 #include "mdst/dataobjects/MCParticle.h"
39 #include "trg/cdc/FrontEnd.h"
40 #include "trg/cdc/Merger.h"
41 #include "trg/cdc/LUT.h"
42 #include "trg/trg/Constants.h"
43 #include "trg/cdc/Helix.h"
44 #include "trg/cdc/Fitter3DUtility.h"
45 #include "trg/cdc/EventTime.h"
46 #include "trg/cdc/JSignal.h"
47 #include "trg/cdc/JLUT.h"
48 #include "trg/cdc/JSignalData.h"
49 #include "trg/cdc/FpgaUtility.h"
50 #include "trg/cdc/HandleRoot.h"
51 
52 using namespace std;
53 namespace Belle2 {
59  TRGCDCFitter3D::TRGCDCFitter3D(const std::string& name,
60  const std::string& rootFitter3DFile,
61  const TRGCDC& TRGCDC,
62  const std::map<std::string, bool>& flags)
63  : m_name(name), m_cdc(TRGCDC),
64  m_mBool(flags), m_rootFitter3DFileName(rootFitter3DFile),
65  m_treeTrackFitter3D(), m_treeConstantsFitter3D() // 2019/07/31 by ytlai
66  {
67  m_fileFitter3D = 0;
68  m_commonData = 0;
69  }
70 
71  TRGCDCFitter3D::~TRGCDCFitter3D()
72  {
73  }
74 
75  void TRGCDCFitter3D::initialize()
76  {
77 
78  StoreObjPtr<EventMetaData> evtMetaData;
79  evtMetaData.isRequired();
80 
81  // If we are using Monte Carlo information.
82  m_mBool["fMc"] = 1;
83  m_mBool["fVerbose"] = 0;
84  m_mBool["fIsPrintError"] = 0;
85  m_mBool["fIsIntegerEffect"] = 1;
86  m_mBool["debugFitted"] = 1;
87  m_mBool["debugLargeZ0"] = 0;
88 
89  // Init values
90  m_mConstD["Trg_PI"] = 3.141592653589793;
91 
92  // Get rr,zToStraw,angleSt,nWire
93  const CDC::CDCGeometryPar& cdcp = CDC::CDCGeometryPar::Instance();
94  m_mConstV["rr"] = vector<double> (9);
95  m_mConstV["nWires"] = vector<double> (9);
96  m_mConstV["nTSs"] = vector<double> (9);
97  for (unsigned iSL = 0; iSL < 9; iSL++) {
98  unsigned t_layerId = m_cdc.segment(iSL, 0).center().layerId();
99  m_mConstV["rr"][iSL] = cdcp.senseWireR(t_layerId);
100  m_mConstV["nWires"][iSL] = cdcp.nWiresInLayer(t_layerId) * 2;
101  m_mConstV["nTSs"][iSL] = cdcp.nWiresInLayer(t_layerId);
102  }
103  m_mConstV["nTSs2D"] = vector<double> (5);
104  for (unsigned iAx = 0; iAx < 5; iAx++) {
105  m_mConstV["nTSs2D"][iAx] = m_mConstV["nTSs"][2 * iAx];
106  }
107 
108  m_mConstV["zToStraw"] = vector<double> (4);
109  m_mConstV["zToOppositeStraw"] = vector<double> (4);
110  m_mConstV["angleSt"] = vector<double> (4);
111  m_mConstV["nShift"] = vector<double> (4);
112  for (int iSt = 0; iSt < 4; iSt++) {
113  unsigned t_layerId = m_cdc.stereoSegment(iSt, 0).center().layerId();
114  m_mConstV["zToStraw"][iSt] = cdcp.senseWireBZ(t_layerId);
115  m_mConstV["zToOppositeStraw"][iSt] = cdcp.senseWireFZ(t_layerId);
116  m_mConstV["angleSt"][iSt] = 2 * m_mConstV["rr"][2 * iSt + 1] * sin(m_mConstD["Trg_PI"] * cdcp.nShifts(t_layerId) /
117  (2 * cdcp.nWiresInLayer(t_layerId))) / (cdcp.senseWireFZ(t_layerId) - cdcp.senseWireBZ(t_layerId));
118  m_mConstV["nShift"][iSt] = cdcp.nShifts(t_layerId);
119  }
120 
121  m_mConstV["rr2D"] = vector<double> (5);
122  m_mConstV["rr3D"] = vector<double> (4);
123  for (int iAx = 0; iAx < 5; iAx++) m_mConstV["rr2D"][iAx] = m_mConstV["rr"][iAx * 2];
124  for (int iSt = 0; iSt < 4; iSt++) m_mConstV["rr3D"][iSt] = m_mConstV["rr"][iSt * 2 + 1];
125 
126  m_mConstV["wirePhi2DError"] = vector<double> (5);
127  m_mConstV["driftPhi2DError"] = vector<double> (5);
129  //m_mConstD["driftResolution"] = 2 * 40 * 0.0001; // (cm)
130  //m_mConstV["cellResolution"] = vector<double> (5);
131  //for(unsigned iAx=0; iAx<5; iAx++){
132  // m_mConstV["cellResolution"][iAx] = m_mConstV["rr2D"][iAx] * 2 * m_mConstD["Trg_PI"] / (m_mConstV["nWires"][iAx*2]/2); // (cm)
133  //}
134  //for(unsigned iAx=0; iAx<5; iAx++){
135  // m_mConstV["wirePhi2DError"][iAx] = m_mConstV["cellResolution"][iAx]/m_mConstV["rr2D"][iAx] / sqrt(12);
136  //}
137  //for(unsigned iAx=0; iAx<5; iAx++){
138  // m_mConstV["driftPhi2DError"][iAx] = m_mConstD["driftResolution"]/m_mConstV["rr2D"][iAx] / sqrt(12);
139  //}
140 
141  m_mConstV["wirePhi2DError"][0] = 0.00085106;
142  m_mConstV["wirePhi2DError"][1] = 0.00039841;
143  m_mConstV["wirePhi2DError"][2] = 0.00025806;
144  m_mConstV["wirePhi2DError"][3] = 0.00019084;
145  m_mConstV["wirePhi2DError"][4] = 0.0001514;
146  m_mConstV["driftPhi2DError"][0] = 0.00085106;
147  m_mConstV["driftPhi2DError"][1] = 0.00039841;
148  m_mConstV["driftPhi2DError"][2] = 0.00025806;
149  m_mConstV["driftPhi2DError"][3] = 0.00019084;
150  m_mConstV["driftPhi2DError"][4] = 0.0001514;
151 
152  // (2016.06.07) study
153  //m_mConstV["wireZError"] = vector<double> ({4.752, 6.393, 6.578, 6.418});
154  //m_mConstV["driftZError"] = vector<double> ({0.4701, 0.7203, 0.8058, 0.9382});
155  m_mConstV["driftZError"] = vector<double> ({0.7676, 0.9753, 1.029, 1.372});
156  m_mConstV["wireZError"] = vector<double> ({0.7676, 0.9753, 1.029, 1.372});
157 
158  // Make driftLength table for each superlayer. Up to 511 clock ticks.
159  // driftLengthTableSLX[ tdcCount (~2ns unit) ] = drift length (cm)
160  for (unsigned iSl = 0; iSl < 9; iSl++) {
161  string tableName = "driftLengthTableSL" + to_string(iSl);
162  unsigned tableSize = 512;
163  m_mConstV[tableName] = vector<double> (tableSize);
164  unsigned t_layer = m_cdc.segment(iSl, 0).center().layerId();
165  for (unsigned iTick = 0; iTick < tableSize; iTick++) {
166  double t_driftTime = iTick * 2 * cdcp.getTdcBinWidth();
167  double avgDriftLength = 0;
168  if (m_mBool["fXtSimple"] == 1) {
169  avgDriftLength = cdcp.getNominalDriftV() * t_driftTime;
170  } else {
171  double driftLength_0 = cdcp.getDriftLength(t_driftTime, t_layer, 0);
172  double driftLength_1 = cdcp.getDriftLength(t_driftTime, t_layer, 1);
173  avgDriftLength = (driftLength_0 + driftLength_1) / 2;
174  }
175  m_mConstV[tableName][iTick] = avgDriftLength;
176  }
177  }
178 
179  if (m_mBool["fVerbose"]) {
180  cout << "fLRLUT: " << m_mBool["fLRLUT"] << endl;
181  cout << "fMc: " << m_mBool["fMc"] << endl;
182  cout << "fVerbose: " << m_mBool["fVerbose"] << endl;
183  if (m_mBool["fMc"]) cout << "fmcLR: " << m_mBool["fmcLR"] << endl;
184  cout << "fRoot: " << m_mBool["fRootFile"] << endl;
185  cout << "PI: " << m_mConstD["Trg_PI"] << endl;
186  cout << "rr: " << m_mConstV["rr"][0] << " " << m_mConstV["rr"][1] << " " << m_mConstV["rr"][2] << " " <<
187  m_mConstV["rr"][3] << " " << m_mConstV["rr"][4] << " " << m_mConstV["rr"][5] << " " << m_mConstV["rr"][6] << " " <<
188  m_mConstV["rr"][7] << " " << m_mConstV["rr"][8] << endl;
189  cout << "nWires: " << int(m_mConstV["nWires"][0]) << " " << int(m_mConstV["nWires"][1]) << " " << int(
190  m_mConstV["nWires"][2]) << " " << int(m_mConstV["nWires"][3]) << " " << int(m_mConstV["nWires"][4]) << " " << int(
191  m_mConstV["nWires"][5]) << " " << int(m_mConstV["nWires"][6]) << " " << int(m_mConstV["nWires"][7]) << " " << int(
192  m_mConstV["nWires"][8]) << endl;
193  cout << "zToStraw: " << m_mConstV["zToStraw"][0] << " " << m_mConstV["zToStraw"][1] << " " << m_mConstV["zToStraw"][2] << " " <<
194  m_mConstV["zToStraw"][3] << endl;
195  cout << "angleSt: " << m_mConstV["angleSt"][0] << " " << m_mConstV["angleSt"][1] << " " << m_mConstV["angleSt"][2] << " " <<
196  m_mConstV["angleSt"][3] << endl;
197  cout << "wireZError: " << m_mConstV["wireZError"][0] << " " << m_mConstV["wireZError"][1] << " " << m_mConstV["wireZError"][2] <<
198  " " << m_mConstV["wireZError"][3] << endl;
199  cout << "driftZError: " << m_mConstV["driftZError"][0] << " " << m_mConstV["driftZError"][1] << " " << m_mConstV["driftZError"][2]
200  << " " << m_mConstV["driftZError"][3] << endl;
201  cout << "wirePhiError: " << m_mConstV["wirePhi2DError"][0] << " " << m_mConstV["wirePhi2DError"][1] << " " <<
202  m_mConstV["wirePhi2DError"][2] << " " << m_mConstV["wirePhi2DError"][3] << " " << m_mConstV["wirePhi2DError"][4] << endl;
203  cout << "driftPhiError: " << m_mConstV["driftPhi2DError"][0] << " " << m_mConstV["driftPhi2DError"][1] << " " <<
204  m_mConstV["driftPhi2DError"][2] << " " << m_mConstV["driftPhi2DError"][3] << " " << m_mConstV["driftPhi2DError"][4] << endl;
205  }
206 
207 
208  }
209 
210  int TRGCDCFitter3D::doit(std::vector<TRGCDCTrack*>& trackList)
211  {
212 
213  TRGDebug::enterStage("Fitter 3D");
214 
215  // Set values for saving data.
216  if (m_mBool["fRootFile"]) {
217  m_mDouble["iSave"] = 0;
218  HandleRoot::initializeEvent(m_mTClonesArray);
219  }
220 
221  // Get common values for event.
222  m_mDouble["eventTime"] = m_cdc.getEventTime();
223  StoreObjPtr<EventMetaData> eventMetaDataPtr;
224  // Event starts from 0.
225  m_mDouble["eventNumber"] = eventMetaDataPtr->getEvent();
226 
227  // Fitter3D
228  // Loop over all tracks
229  for (unsigned iTrack = 0; iTrack < trackList.size(); iTrack++) {
230 
231  TCTrack& aTrack = * trackList[iTrack];
232 
234  // Get MC values for fitter.
235  if (m_mBool["fMc"]) getMCValues(m_cdc, &aTrack, m_mConstD, m_mDouble, m_mVector);
236  // Get track ID
237  m_mDouble["trackId"] = aTrack.getTrackID();
238 
241  int fit2DResult = do2DFit(aTrack, m_mBool, m_mConstD, m_mConstV, m_mDouble, m_mVector);
242  if (fit2DResult != 0) {
243  aTrack.setFitted(0);
244  if (m_mBool["fVerbose"]) cout << "Exit due to 2D fit result:" << fit2DResult << endl;
245  continue;
246  }
247 
249  // 3D Fitter
250  // Print input TSs
251  if (m_mBool["fVerbose"]) {
252  for (unsigned iSt = 0; iSt < 4; iSt++) {
253  const vector<TCLink*>& links = aTrack.links(iSt * 2 + 1);
254  const unsigned nSegments = links.size();
255  cout << "iSt:" << iSt << " nSegments:" << nSegments << endl;
256  for (unsigned iTS = 0; iTS < nSegments; iTS++) {
257  const TCSegment* t_segment = dynamic_cast<const TCSegment*>(& links[iTS]->hit()->cell());
258  cout << " tsId:" << t_segment->localId()
259  << " tdc:" << t_segment->priorityTime() << " lr:" << t_segment->LUT()->getValue(t_segment->lutPattern())
260  << " priorityPosition:" << t_segment->priorityPosition() << endl;
261  }
262  }
263  }
264  // Check which stereo super layers should be used.
265  m_mVector["useStSl"] = vector<double> (4);
266  findHitStereoSuperlayers(aTrack, m_mVector["useStSl"], m_mBool["fIsPrintError"]);
267  removeImpossibleStereoSuperlayers(m_mVector["useStSl"]);
268  m_mDouble["nHitStSl"] = m_mVector["useStSl"][0] + m_mVector["useStSl"][1] + m_mVector["useStSl"][2] + m_mVector["useStSl"][3];
269 
270  m_mVector["useSl"] = vector<double> (9);
271  for (unsigned iAx = 0; iAx < 5; iAx++) m_mVector["useSl"][2 * iAx] = m_mVector["useAxSl"][iAx];
272  for (unsigned iSt = 0; iSt < 4; iSt++) m_mVector["useSl"][2 * iSt + 1] = m_mVector["useAxSl"][iSt];
273 
274  // Fill information for stereo layers
275  for (unsigned iSt = 0; iSt < 4; iSt++) {
276  if (m_mVector["useStSl"][iSt] == 1) {
277  const vector<TCLink*>& links = aTrack.links(iSt * 2 + 1);
278  const TCSegment* t_segment = dynamic_cast<const TCSegment*>(& links[0]->hit()->cell());
279  m_mVector["tsId"][iSt * 2 + 1] = t_segment->localId();
280  m_mVector["wirePhi"][iSt * 2 + 1] = (double) t_segment->localId() / m_mConstV["nWires"][iSt * 2 + 1] * 4 * m_mConstD["Trg_PI"];
281  m_mVector["lutLR"][iSt * 2 + 1] = t_segment->LUT()->getValue(t_segment->lutPattern());
282  if (m_mBool["fMc"]) m_mVector["mcLR"][iSt * 2 + 1] = t_segment->hit()->mcLR() + 1;
283  m_mVector["driftLength"][iSt * 2 + 1] = t_segment->hit()->drift();
284  m_mVector["tdc"][iSt * 2 + 1] = t_segment->priorityTime();
285  if (m_mBool["fmcLR"] == 1) m_mVector["LR"][iSt * 2 + 1] = m_mVector["mcLR"][iSt * 2 + 1];
286  else if (m_mBool["fLRLUT"] == 1) m_mVector["LR"][iSt * 2 + 1] = m_mVector["lutLR"][iSt * 2 + 1];
287  else m_mVector["LR"][iSt * 2 + 1] = 3;
288  } else {
289  m_mVector["tsId"][iSt * 2 + 1] = 9999;
290  m_mVector["wirePhi"][iSt * 2 + 1] = 9999;
291  m_mVector["lutLR"][iSt * 2 + 1] = 0;
292  if (m_mBool["fMc"]) m_mVector["mcLR"][iSt * 2 + 1] = 9999;
293  m_mVector["driftLength"][iSt * 2 + 1] = 9999;
294  m_mVector["tdc"][iSt * 2 + 1] = 9999;
295  if (m_mBool["fmcLR"] == 1) m_mVector["LR"][iSt * 2 + 1] = 9999;
296  else if (m_mBool["fLRLUT"] == 1) m_mVector["LR"][iSt * 2 + 1] = 9999;
297  else m_mVector["LR"][iSt * 2 + 1] = 9999;
298  }
299  } // End superlayer loop
300 
301  // Calculate phi3D.
302  m_mVector["phi3D"] = vector<double> (4);
303  if (m_mDouble["eventTime"] == 9999) {
304  for (unsigned iSt = 0; iSt < 4; iSt++) {
305  m_mVector["phi3D"][iSt] = m_mVector["wirePhi"][iSt * 2 + 1];
306  }
307  } else {
308  for (unsigned iSt = 0; iSt < 4; iSt++) {
309  if (m_mVector["useStSl"][iSt] == 1) {
310  // Get drift length from table.
311  string tableName = "driftLengthTableSL" + to_string(iSt * 2 + 1);
312  double t_driftTime = m_mVector["tdc"][iSt * 2 + 1] - m_mDouble["eventTime"];
313  if (t_driftTime < 0) t_driftTime = 0;
314  double t_driftLength = m_mConstV[tableName][(unsigned)t_driftTime];
315  m_mVector["phi3D"][iSt] = Fitter3DUtility::calPhi(m_mVector["wirePhi"][iSt * 2 + 1], t_driftLength, m_mConstV["rr3D"][iSt],
316  m_mVector["LR"][iSt * 2 + 1]);
317  } else {
318  m_mVector["phi3D"][iSt] = 9999;
319  }
320  }
321  }
322  // Get zerror for 3D fit
323  m_mVector["zError"] = vector<double> (4);
324  for (unsigned iSt = 0; iSt < 4; iSt++) {
325  if (m_mVector["useStSl"][iSt] == 1) {
326  // Check LR.
327  if (m_mVector["LR"][2 * iSt + 1] != 3) m_mVector["zError"][iSt] = m_mConstV["driftZError"][iSt];
328  else m_mVector["zError"][iSt] = m_mConstV["wireZError"][iSt];
329  // Check eventTime
330  if (m_mDouble["eventTime"] == 9999) m_mVector["zError"][iSt] = m_mConstV["wireZError"][iSt];
331  } else {
332  m_mVector["zError"][iSt] = 9999;
333  }
334  }
335  // Get inverse zerror ^ 2
336  m_mVector["iZError2"] = vector<double> (4);
337  for (unsigned iSt = 0; iSt < 4; iSt++) {
338  if (m_mVector["useStSl"][iSt] == 1) {
339  m_mVector["iZError2"][iSt] = 1 / pow(m_mVector["zError"][iSt], 2);
340  } else {
341  m_mVector["iZError2"][iSt] = 0;
342  }
343  }
344 
345  // Calculate zz
346  m_mVector["zz"] = vector<double> (4);
347  for (unsigned iSt = 0; iSt < 4; iSt++) {
348  if (m_mVector["useStSl"][iSt] == 1) {
349  //m_mVector["zz"][iSt] = Fitter3DUtility::calZ(m_mDouble["charge"], m_mConstV["angleSt"][iSt], m_mConstV["zToStraw"][iSt], m_mConstV["rr3D"][iSt], m_mVector["phi3D"][iSt], m_mDouble["rho"], m_mDouble["phi0"]);
350  m_mVector["zz"][iSt] = Fitter3DUtility::calZ(m_mDouble["charge2D"], m_mConstV["angleSt"][iSt], m_mConstV["zToStraw"][iSt],
351  m_mConstV["rr3D"][iSt], m_mVector["phi3D"][iSt], m_mDouble["rho"], m_mDouble["phi0"]);
352  } else {
353  m_mVector["zz"][iSt] = 0;
354  }
355  }
356  // Calculate arcS
357  m_mVector["arcS"] = vector<double> (4);
358  for (unsigned iSt = 0; iSt < 4; iSt++) {
359  if (m_mVector["useStSl"][iSt] == 1) {
360  m_mVector["arcS"][iSt] = Fitter3DUtility::calS(m_mDouble["rho"], m_mConstV["rr3D"][iSt]);
361  } else {
362  m_mVector["arcS"][iSt] = 0;
363  }
364  }
365  // Fit3D
366  m_mDouble["z0"] = 0;
367  m_mDouble["cot"] = 0;
368  m_mDouble["zChi2"] = 0;
369  Fitter3DUtility::rSFit(&m_mVector["iZError2"][0], &m_mVector["arcS"][0], &m_mVector["zz"][0], m_mDouble["z0"], m_mDouble["cot"],
370  m_mDouble["zChi2"]);
371  // Change to deg
372  m_mDouble["theta"] = m_mConstD["Trg_PI"] / 2 - atan(m_mDouble["cot"]);
373  m_mDouble["theta"] = 180 / m_mConstD["Trg_PI"];
374 
375  if (m_mBool["fVerbose"]) print3DInformation(iTrack);
376 
377  // For failed fits. When cot is 0 or nan.
378  if (m_mDouble["cot"] == 0 || std::isnan(m_mDouble["cot"])) {
379  aTrack.setFitted(0);
380  aTrack.setDebugValue(TRGCDCTrack::EDebugValueType::fit3D, 1);
381  if (m_mBool["fVerbose"]) cout << "Exit due to 3D fit cot result:" << m_mDouble["cot"] << endl;
382  continue;
383  }
384 
385  // Set track in trackList
386  if (m_mBool["fVerbose"]) cout << "Fit was done successfully." << endl;
387  // Set Helix parameters
388  TRGCDCHelix helix(ORIGIN, CLHEP::HepVector(5, 0), CLHEP::HepSymMatrix(5, 0));
389  CLHEP::HepVector a(5);
390  a = aTrack.helix().a();
391  aTrack.setFitted(1);
392  //if(m_mDouble["charge"]<0)
393  if (m_mDouble["charge2D"] < 0) {
394  a[1] = fmod(m_mDouble["phi0"] + m_mConstD["Trg_PI"], 2 * m_mConstD["Trg_PI"]);
395  } else {
396  a[1] = m_mDouble["phi0"];
397  }
398  //a[2] = 1/m_mDouble["pt"]*m_mDouble["charge"];
399  a[2] = 1 / m_mDouble["pt"] * m_mDouble["charge2D"];
400  a[3] = m_mDouble["z0"];
401  a[4] = m_mDouble["cot"];
402  helix.a(a);
403  aTrack.set2DFitChi2(m_mDouble["fit2DChi2"]);
404  aTrack.set3DFitChi2(m_mDouble["zChi2"]);
405  aTrack.setHelix(helix);
406 
408  // Save values
409  if (m_mBool["fRootFile"]) {
410  if (m_fileFitter3D == 0) {
411  m_fileFitter3D = new TFile(m_rootFitter3DFileName.c_str(), "RECREATE");
412  HandleRoot::initializeRoot("fitter3D", &m_treeConstantsFitter3D, &m_treeTrackFitter3D,
413  m_mTVectorD, m_mTClonesArray,
414  m_mConstD, m_mConstV,
415  m_mDouble, m_mVector
416  );
417  }
418  HandleRoot::saveTrackValues("fitter3D",
419  m_mTClonesArray, m_mDouble, m_mVector
420  );
421  }
422 
423  } // End track loop
424 
425  // Save values to file
426  // To prevent crash when there are no tracks in first event.
427  // If there is no track in first case then the ROOT saving functions might fail.
428  // This is due to bad software design.
429  // The first event that have tracks will be event 0 in ROOT file.
430  if (m_mBool["fRootFile"]) {
431  if (m_fileFitter3D != 0) m_treeTrackFitter3D->Fill();
432  }
433 
434  if (m_mBool["debugFitted"]) {
435  for (unsigned iTrack = 0; iTrack < trackList.size(); iTrack++) {
436  TCTrack& aTrack = * trackList[iTrack];
437  if (aTrack.fitted() == 0) aTrack.setDebugValue(TRGCDCTrack::EDebugValueType::fit3D, 1);
438  }
439  }
440  if (m_mBool["debugLargeZ0"]) {
441  // If 3D fit z0 is larger or smaller than expected.
442  for (unsigned iTrack = 0; iTrack < trackList.size(); iTrack++) {
443  TCTrack& aTrack = *trackList[iTrack];
444  if (aTrack.fitted()) {
445  double fitZ0 = aTrack.helix().dz();
446  if (fitZ0 > 20 || fitZ0 < -20) aTrack.setDebugValue(TRGCDCTrack::EDebugValueType::fit3D, 1);
447  }
448  }
449  }
450 
451 
452  TRGDebug::leaveStage("Fitter 3D");
453  return 1;
454 
455  }
456 
457  int TRGCDCFitter3D::doitComplex(std::vector<TRGCDCTrack*>& trackList)
458  {
459 
460  TRGDebug::enterStage("Fitter 3D");
461 
462  // Set values for saving data.
463  if (m_mBool["fRootFile"]) {
464  m_mDouble["iSave"] = 0;
465  HandleRoot::initializeEvent(m_mTClonesArray);
466  }
467  if (m_commonData == 0) {
468  m_commonData = new Belle2::TRGCDCJSignalData();
469  }
470 
471  // Get common values for event.
472  m_mDouble["eventTime"] = m_cdc.getEventTime();
473  StoreObjPtr<EventMetaData> eventMetaDataPtr;
474  // Event starts from 0.
475  m_mDouble["eventNumber"] = eventMetaDataPtr->getEvent();
476 
477  // Fitter3D
478  // Loop over all tracks
479  for (unsigned iTrack = 0; iTrack < trackList.size(); iTrack++) {
480 
481  TCTrack& aTrack = * trackList[iTrack];
482 
484  // Get MC values for 3D fitter
485  if (m_mBool["fMc"]) getMCValues(m_cdc, &aTrack, m_mConstD, m_mDouble, m_mVector);
486  // Get input values for 3D fitter
487  // Get event and track ID
488  m_mDouble["trackId"] = aTrack.getTrackID();
489 
492  int fit2DResult = do2DFit(aTrack, m_mBool, m_mConstD, m_mConstV, m_mDouble, m_mVector);
493  if (fit2DResult != 0) {
494  aTrack.setFitted(0);
495  continue;
496  }
497 
499  // Start of 3D fitter
500  // Print input TSs
501  if (m_mBool["fVerbose"]) {
502  for (unsigned iSt = 0; iSt < 4; iSt++) {
503  const vector<TCLink*>& links = aTrack.links(iSt * 2 + 1);
504  const unsigned nSegments = links.size();
505  cout << "iSt:" << iSt << " nSegments:" << nSegments << endl;
506  for (unsigned iTS = 0; iTS < nSegments; iTS++) {
507  const TCSegment* t_segment = dynamic_cast<const TCSegment*>(& links[iTS]->hit()->cell());
508  cout << " tsId:" << t_segment->localId()
509  << " tdc:" << t_segment->priorityTime() << " lr:" << t_segment->LUT()->getValue(t_segment->lutPattern())
510  << " priorityPosition:" << t_segment->priorityPosition() << endl;
511  }
512  }
513  }
514  m_mVector["useStSl"] = vector<double> (4);
515  findHitStereoSuperlayers(aTrack, m_mVector["useStSl"], m_mBool["fIsPrintError"]);
516  removeImpossibleStereoSuperlayers(m_mVector["useStSl"]);
517  m_mDouble["nHitStSl"] = m_mVector["useStSl"][0] + m_mVector["useStSl"][1] + m_mVector["useStSl"][2] + m_mVector["useStSl"][3];
518 
519  // Fill information for stereo layers
520  for (unsigned iSt = 0; iSt < 4; iSt++) {
521  if (m_mVector["useStSl"][iSt] == 1) {
522  const vector<TCLink*>& links = aTrack.links(iSt * 2 + 1);
523  const TCSegment* t_segment = dynamic_cast<const TCSegment*>(& links[0]->hit()->cell());
524  m_mVector["tsId"][iSt * 2 + 1] = t_segment->localId();
525  m_mVector["wirePhi"][iSt * 2 + 1] = (double) t_segment->localId() / m_mConstV["nWires"][iSt * 2 + 1] * 4 * m_mConstD["Trg_PI"];
526  m_mVector["lutLR"][iSt * 2 + 1] = t_segment->LUT()->getValue(t_segment->lutPattern());
527  if (m_mBool["fMc"]) m_mVector["mcLR"][iSt * 2 + 1] = t_segment->hit()->mcLR() + 1;
528  m_mVector["driftLength"][iSt * 2 + 1] = t_segment->hit()->drift();
529  m_mVector["tdc"][iSt * 2 + 1] = t_segment->priorityTime();
530  if (m_mBool["fmcLR"] == 1) m_mVector["LR"][iSt * 2 + 1] = m_mVector["mcLR"][iSt * 2 + 1];
531  else if (m_mBool["fLRLUT"] == 1) m_mVector["LR"][iSt * 2 + 1] = m_mVector["lutLR"][iSt * 2 + 1];
532  else m_mVector["LR"][iSt * 2 + 1] = 3;
533  } else {
534  m_mVector["tsId"][iSt * 2 + 1] = 0;
535  m_mVector["wirePhi"][iSt * 2 + 1] = 9999;
536  m_mVector["lutLR"][iSt * 2 + 1] = 0;
537  if (m_mBool["fMc"]) m_mVector["mcLR"][iSt * 2 + 1] = 9999;
538  m_mVector["driftLength"][iSt * 2 + 1] = 9999;
539  m_mVector["tdc"][iSt * 2 + 1] = 0;
540  if (m_mBool["fmcLR"] == 1) m_mVector["LR"][iSt * 2 + 1] = 9999;
541  else if (m_mBool["fLRLUT"] == 1) m_mVector["LR"][iSt * 2 + 1] = 9999;
542  else m_mVector["LR"][iSt * 2 + 1] = 9999;
543  }
544  } // End superlayer loop
546  //int minTSTdc = 9999;
547  //for(unsigned iSl=0; iSl<9; iSl++){
548  // if (minTSTdc > m_mVector["tdc"][iSl]) minTSTdc = m_mVector["tdc"][iSl];
549  //}
550  //m_mDouble["eventTime"] = minTSTdc;
551 
552  // Calculate phi3D.
553  m_mVector["phi3D"] = vector<double> (4);
554  if (m_mDouble["eventTime"] == 9999) {
555  for (unsigned iSt = 0; iSt < 4; iSt++) {
556  m_mVector["phi3D"][iSt] = m_mVector["wirePhi"][iSt * 2 + 1];
557  }
558  } else {
559  for (unsigned iSt = 0; iSt < 4; iSt++) {
560  // Get drift length from table.
561  string tableName = "driftLengthTableSL" + to_string(iSt * 2 + 1);
562  double t_driftTime = m_mVector["tdc"][iSt * 2 + 1] - m_mDouble["eventTime"];
563  if (t_driftTime < 0) t_driftTime = 0;
564  double t_driftLength = m_mConstV[tableName][(unsigned)t_driftTime];
565  m_mVector["phi3D"][iSt] = Fitter3DUtility::calPhi(m_mVector["wirePhi"][iSt * 2 + 1], t_driftLength, m_mConstV["rr3D"][iSt],
566  m_mVector["LR"][iSt * 2 + 1]);
567  }
568  }
569 
570  // Get zerror for 3D fit
571  m_mVector["zError"] = vector<double> (4);
572  for (unsigned iSt = 0; iSt < 4; iSt++) {
573  // Check LR.
574  if (m_mVector["LR"][2 * iSt + 1] != 3) m_mVector["zError"][iSt] = m_mConstV["driftZError"][iSt];
575  else m_mVector["zError"][iSt] = m_mConstV["wireZError"][iSt];
576  // Check eventTime
577  if (m_mDouble["eventTime"] == 9999) m_mVector["zError"][iSt] = m_mConstV["wireZError"][iSt];
578  }
579  // Get inverse zerror ^ 2
580  m_mVector["iZError2"] = vector<double> (4);
581  for (unsigned iSt = 0; iSt < 4; iSt++) {
582  if (m_mVector["useStSl"][iSt] == 1) {
583  m_mVector["iZError2"][iSt] = 1 / pow(m_mVector["zError"][iSt], 2);
584  } else {
585  m_mVector["iZError2"][iSt] = 0;
586  }
587  }
588 
589  // Simple version of Fitter3D
591  //m_mVector["zz"] = vector<double> (4);
592  //for (unsigned iSt = 0; iSt < 4; iSt++) m_mVector["zz"][iSt] = Fitter3DUtility::calZ(m_mDouble["charge"], m_mConstV["angleSt"][iSt], m_mConstV["zToStraw"][iSt], m_mConstV["rr3D"][iSt], m_mVector["phi3D"][iSt], m_mDouble["rho"], m_mDouble["phi0"]);
594  //m_mVector["arcS"] = vector<double> (4);
595  //for (unsigned iSt = 0; iSt < 4; iSt++) m_mVector["arcS"][iSt] = Fitter3DUtility::calS(m_mDouble["rho"], m_mConstV["rr3D"][iSt]);
597  //m_mDouble["z0"] = 0;
598  //m_mDouble["cot"] = 0;
599  //m_mDouble["zChi2"] = 0;
600  //Fitter3DUtility::rSFit(&m_mVector["iZError2"][0], &m_mVector["arcS"][0], &m_mVector["zz"][0], m_mDouble["z0"], m_mDouble["cot"], m_mDouble["zChi2"]);
602  //m_mDouble["theta"] = m_mConstD["Trg_PI"]/2 - atan(m_mDouble["cot"]);
603  //m_mDouble["theta"] = 180 / m_mConstD["Trg_PI"];
604 
605  double phiMax = m_mConstD["Trg_PI"];
606  double phiMin = -m_mConstD["Trg_PI"];
607  int phiBitSize = 13;
608  // pt = 0.3*1.5*rho*0.01;
609  //double rhoMin = 48;
610  /* cppcheck-suppress variableScope */
611  double rhoMin = 20;
612  double rhoMax = 2500;
613  //double rhoMax = 1600;
614  // 5bit (clock counter) + 4 bit (2ns resolution)
615  m_mConstD["tdcBitSize"] = 9;
616  m_mConstD["rhoBitSize"] = 11;
617  m_mConstD["iError2BitSize"] = 8;
618  m_mConstD["iError2Max"] = 1 / pow(m_mConstV["wireZError"][0], 2);
619  // LUT values
620  m_mConstD["JB"] = 0;
621  m_mConstD["driftPhiLUTOutBitSize"] = phiBitSize - 1;
622  m_mConstD["driftPhiLUTInBitSize"] = m_mConstD["tdcBitSize"];
623  m_mConstD["acosLUTOutBitSize"] = phiBitSize - 1;
624  m_mConstD["acosLUTInBitSize"] = m_mConstD["rhoBitSize"];
625  m_mConstD["zLUTInBitSize"] = phiBitSize;
626  m_mConstD["zLUTOutBitSize"] = 9;
627  m_mConstD["iDenLUTInBitSize"] = 13;
628  m_mConstD["iDenLUTOutBitSize"] = 11; // To increase wireZError = 2.5*driftZError
629  // Rotate by quadrants depending on charge and cc(circle center)
630  m_mDouble["relRefPhi"] = 0;
631  int t_quadrant = Fitter3DUtility::findQuadrant(m_mDouble["phi0"]);
632  //if (m_mDouble["charge"] == -1)
633  if (m_mDouble["charge2D"] == -1) {
634  if (t_quadrant == 1) m_mDouble["relRefPhi"] = 0;
635  else if (t_quadrant == 2) m_mDouble["relRefPhi"] = -m_mConstD["Trg_PI"] / 2;
636  else if (t_quadrant == 3) m_mDouble["relRefPhi"] = -m_mConstD["Trg_PI"];
637  else if (t_quadrant == 4) m_mDouble["relRefPhi"] = m_mConstD["Trg_PI"] / 2;
638  } else {
639  if (t_quadrant == 1) m_mDouble["relRefPhi"] = m_mConstD["Trg_PI"] / 2;
640  else if (t_quadrant == 2) m_mDouble["relRefPhi"] = 0;
641  else if (t_quadrant == 3) m_mDouble["relRefPhi"] = -m_mConstD["Trg_PI"] / 2;
642  else if (t_quadrant == 4) m_mDouble["relRefPhi"] = -m_mConstD["Trg_PI"];
643  }
644  // Rotate phi
645  m_mDouble["relPhi0"] = Fitter3DUtility::rotatePhi(m_mDouble["phi0"], m_mDouble["relRefPhi"]);
646  m_mVector["relPhi3D"] = vector<double> (4);
647  for (unsigned iSt = 0; iSt < 4;
648  iSt++) m_mVector["relPhi3D"][iSt] = Fitter3DUtility::rotatePhi(m_mVector["phi3D"][iSt], m_mDouble["relRefPhi"]);
649 
650  // Rotate wirePhi
651  m_mVector["relWirePhi3D"] = vector<double> (4);
652  for (unsigned iSt = 0; iSt < 4; iSt++) {
653  double t_relWirePhi = Fitter3DUtility::rotatePhi(m_mVector["wirePhi"][2 * iSt + 1], m_mDouble["relRefPhi"]);
654  bool rangeOk = 0;
655  while (rangeOk == 0) {
656  if (t_relWirePhi < 0) t_relWirePhi += 2 * m_mConstD["Trg_PI"];
657  else if (t_relWirePhi > 2 * m_mConstD["Trg_PI"]) t_relWirePhi -= 2 * m_mConstD["Trg_PI"];
658  else rangeOk = 1;
659  }
660  m_mVector["relWirePhi3D"][iSt] = t_relWirePhi;
661  }
662 
663  // Rotate tsId
664  m_mVector["relTsId3D"] = vector<double> (4);
665  for (unsigned iSt = 0; iSt < 4;
666  iSt++) m_mVector["relTsId3D"][iSt] = Fitter3DUtility::rotateTsId(m_mVector["tsId"][2 * iSt + 1],
667  m_mDouble["relRefPhi"] / m_mConstD["Trg_PI"] / 2 * m_mConstV["nTSs"][2 * iSt + 1], m_mConstV["nTSs"][2 * iSt + 1]);
668 
669  // Constrain rho to rhoMax. For removing warnings when changing to signals.
670  if (m_mDouble["rho"] > rhoMax) {
671  m_mDouble["rho"] = rhoMax;
672  m_mDouble["pt"] = rhoMax * 0.3 * 1.5 * 0.01;
673  }
674 
675  // Constrain event time
676  m_mDouble["eventTimeValid"] = 1;
677  if (m_mDouble["eventTime"] == 9999) m_mDouble["eventTimeValid"] = 0;
678  // Change tdc and eventTime to 9 bit unsigned value. (Ex: -1 => 511, -2 => 510)
679  m_mVector["unsignedTdc"] = vector<double> (9);
680  for (unsigned iSt = 0; iSt < 4; iSt++) {
681  //cout<<"iSt:"<<iSt<<" tdc:"<<m_mVector["tdc"][2*iSt+1]<<" unsignedTdc:"<<Fitter3DUtility::toUnsignedTdc(m_mVector["tdc"][2*iSt+1], 9)<<endl;
682  m_mVector["unsignedTdc"][2 * iSt + 1] = Fitter3DUtility::toUnsignedTdc(m_mVector["tdc"][2 * iSt + 1], m_mConstD["tdcBitSize"]);
683  }
684  m_mDouble["unsignedEventTime"] = Fitter3DUtility::toUnsignedTdc(m_mDouble["eventTime"], m_mConstD["tdcBitSize"]);
685  //cout<<"eventTime:"<<m_mDouble["eventTime"]<<" unsignedEventTime:"<<Fitter3DUtility::toUnsignedTdc(m_mDouble["eventTime"], 9)<<endl;
686 
687  // Change to Signals.
688  {
689  vector<tuple<string, double, int, double, double, int> > t_values = {
690  make_tuple("phi0", m_mDouble["relPhi0"], phiBitSize, phiMin, phiMax, 0),
691  make_tuple("rho", m_mDouble["rho"], m_mConstD["rhoBitSize"], rhoMin, rhoMax, 0),
692  //make_tuple("charge",(int) (m_mDouble["charge"]==1 ? 1 : 0), 1, 0, 1.5, 0),
693  make_tuple("charge", (int)(m_mDouble["charge2D"] == 1 ? 1 : 0), 1, 0, 1.5, 0),
694  make_tuple("lr_0", m_mVector["lutLR"][1], 2, 0, 3.5, 0),
695  make_tuple("lr_1", m_mVector["lutLR"][3], 2, 0, 3.5, 0),
696  make_tuple("lr_2", m_mVector["lutLR"][5], 2, 0, 3.5, 0),
697  make_tuple("lr_3", m_mVector["lutLR"][7], 2, 0, 3.5, 0),
698  make_tuple("tsId_0", m_mVector["relTsId3D"][0], ceil(log(m_mConstV["nTSs"][1]) / log(2)), 0, pow(2, ceil(log(m_mConstV["nTSs"][1]) / log(2))) - 0.5, 0),
699  make_tuple("tsId_1", m_mVector["relTsId3D"][1], ceil(log(m_mConstV["nTSs"][3]) / log(2)), 0, pow(2, ceil(log(m_mConstV["nTSs"][3]) / log(2))) - 0.5, 0),
700  make_tuple("tsId_2", m_mVector["relTsId3D"][2], ceil(log(m_mConstV["nTSs"][5]) / log(2)), 0, pow(2, ceil(log(m_mConstV["nTSs"][5]) / log(2))) - 0.5, 0),
701  make_tuple("tsId_3", m_mVector["relTsId3D"][3], ceil(log(m_mConstV["nTSs"][7]) / log(2)), 0, pow(2, ceil(log(m_mConstV["nTSs"][7]) / log(2))) - 0.5, 0),
702  make_tuple("tdc_0", m_mVector["unsignedTdc"][1], m_mConstD["tdcBitSize"], 0, pow(2, m_mConstD["tdcBitSize"]) - 0.5, 0),
703  make_tuple("tdc_1", m_mVector["unsignedTdc"][3], m_mConstD["tdcBitSize"], 0, pow(2, m_mConstD["tdcBitSize"]) - 0.5, 0),
704  make_tuple("tdc_2", m_mVector["unsignedTdc"][5], m_mConstD["tdcBitSize"], 0, pow(2, m_mConstD["tdcBitSize"]) - 0.5, 0),
705  make_tuple("tdc_3", m_mVector["unsignedTdc"][7], m_mConstD["tdcBitSize"], 0, pow(2, m_mConstD["tdcBitSize"]) - 0.5, 0),
706  make_tuple("eventTime", m_mDouble["unsignedEventTime"], m_mConstD["tdcBitSize"], 0, pow(2, m_mConstD["tdcBitSize"]) - 0.5, 0),
707  make_tuple("eventTimeValid", (int) m_mDouble["eventTimeValid"], 1, 0, 1.5, 0),
708  };
709  TRGCDCJSignal::valuesToMapSignals(t_values, m_commonData, m_mSignalStorage);
710  }
711 
712  Fitter3DUtility::setError(m_mConstD, m_mConstV, m_mSignalStorage);
713  Fitter3DUtility::calPhi(m_mConstD, m_mConstV, m_mSignalStorage, m_mLutStorage);
714  Fitter3DUtility::constrainRPerStSl(m_mConstV, m_mSignalStorage);
715  Fitter3DUtility::calZ(m_mConstD, m_mConstV, m_mSignalStorage, m_mLutStorage);
716  Fitter3DUtility::calS(m_mConstD, m_mConstV, m_mSignalStorage, m_mLutStorage);
717  Fitter3DUtility::rSFit(m_mConstD, m_mConstV, m_mSignalStorage, m_mLutStorage);
718  // Change to values.
719  vector<tuple<string, double, int, double, double, int> > resultValues;
720  {
721  vector<pair<string, int> > t_chooseSignals = {
722  make_pair("z0", 0), make_pair("cot", 0), make_pair("chi2Sum", 0)
723  };
724  TRGCDCJSignal::mapSignalsToValues(m_mSignalStorage, t_chooseSignals, resultValues);
725  }
726 
727  // Post handling of signals.
728  // Name all signals.
729  if ((*m_mSignalStorage.begin()).second.getName() == "") {
730  for (auto it = m_mSignalStorage.begin(); it != m_mSignalStorage.end(); ++it) {
731  (*it).second.setName((*it).first);
732  }
733  }
742 
743 
744  // Takes some time.
745  // Save values.
746  if (m_mBool["fRootFile"]) {
747  // Check if there is a name.
748  if ((*m_mSignalStorage.begin()).second.getName() != "") {
749 
750  // Save values to root file.
751  {
752  vector<pair<string, int> > chooseValues = {
753  make_pair("zz_0", m_mBool["fIsIntegerEffect"]),
754  make_pair("zz_1", m_mBool["fIsIntegerEffect"]),
755  make_pair("zz_2", m_mBool["fIsIntegerEffect"]),
756  make_pair("zz_3", m_mBool["fIsIntegerEffect"]),
757  make_pair("arcS_0", m_mBool["fIsIntegerEffect"]),
758  make_pair("arcS_1", m_mBool["fIsIntegerEffect"]),
759  make_pair("arcS_2", m_mBool["fIsIntegerEffect"]),
760  make_pair("arcS_3", m_mBool["fIsIntegerEffect"]),
761  make_pair("z0", m_mBool["fIsIntegerEffect"]),
762  make_pair("cot", m_mBool["fIsIntegerEffect"]),
763  make_pair("zChi2", m_mBool["fIsIntegerEffect"])
764  };
765  // outValues => [name, value, bitwidth, min, max, clock]
766  vector<tuple<string, double, int, double, double, int> > outValues;
767  TRGCDCJSignal::mapSignalsToValues(m_mSignalStorage, chooseValues, outValues);
768  // Changes names with "_" to m_mVector.
769  HandleRoot::convertSignalValuesToMaps(outValues, m_mDouble, m_mVector);
770  }
771  }
772  }
773 
774  m_mBool["fVHDLFile"] = 0;
775  if (m_mBool["fVHDLFile"]) {
776  // Check if there is a name.
777  if ((*m_mSignalStorage.begin()).second.getName() != "") {
778  // Saves to file only one time.
779  saveVhdlAndCoe();
780  // Saves values to memory. Wrote to file at terminate().
781  saveAllSignals();
782  saveIoSignals();
783  }
784  }
785 
787  // Calculate zz
788  m_mVector["float_zz"] = vector<double> (4);
789  //for (unsigned iSt = 0; iSt < 4; iSt++) m_mVector["float_zz"][iSt] = Fitter3DUtility::calZ(m_mDouble["charge"], m_mConstV["angleSt"][iSt], m_mConstV["zToStraw"][iSt], m_mConstV["rr3D"][iSt], m_mVector["phi3D"][iSt], m_mDouble["rho"], m_mDouble["phi0"]);
790  for (unsigned iSt = 0; iSt < 4;
791  iSt++) m_mVector["float_zz"][iSt] = Fitter3DUtility::calZ(m_mDouble["charge2D"], m_mConstV["angleSt"][iSt],
792  m_mConstV["zToStraw"][iSt], m_mConstV["rr3D"][iSt], m_mVector["phi3D"][iSt], m_mDouble["rho"], m_mDouble["phi0"]);
793  // Calculate arcS
794  m_mVector["float_arcS"] = vector<double> (4);
795  for (unsigned iSt = 0; iSt < 4;
796  iSt++) m_mVector["float_arcS"][iSt] = Fitter3DUtility::calS(m_mDouble["rho"], m_mConstV["rr3D"][iSt]);
797  // Fit3D
798  m_mDouble["float_z0"] = 0;
799  m_mDouble["float_cot"] = 0;
800  m_mDouble["float_zChi2"] = 0;
801  Fitter3DUtility::rSFit(&m_mVector["iZError2"][0], &m_mVector["float_arcS"][0], &m_mVector["float_zz"][0], m_mDouble["float_z0"],
802  m_mDouble["float_cot"], m_mDouble["float_zChi2"]);
803 
804  if (m_mBool["fVerbose"]) {
805  for (unsigned iSt = 0; iSt < 4; iSt++) cout << "float_zz[" << iSt << "] : " << m_mVector["float_zz"][iSt] << " ";
806  cout << endl;
807  for (unsigned iSt = 0; iSt < 4; iSt++) cout << "float_arcS[" << iSt << "] : " << m_mVector["float_arcS"][iSt] << " ";
808  cout << endl;
809  cout << "float_z0: " << m_mDouble["float_z0"] << endl;
810  cout << "float_zChi2: " << m_mDouble["float_zChi2"] << endl;
811 
812  print3DInformation(iTrack);
813  }
814 
815  // For failed fits.
816  if (m_mDouble["cot"] == 0) {
817  aTrack.setFitted(0);
818  continue;
819  }
820 
821  // Set track in trackListOut.
822  // Set Helix parameters
823  TRGCDCHelix helix(ORIGIN, CLHEP::HepVector(5, 0), CLHEP::HepSymMatrix(5, 0));
824  CLHEP::HepVector a(5);
825  a = aTrack.helix().a();
826  aTrack.setFitted(1);
827  //if(m_mDouble["charge"]<0)
828  if (m_mDouble["charge2D"] < 0) {
829  a[1] = fmod(m_mDouble["phi0"] + m_mConstD["Trg_PI"], 2 * m_mConstD["Trg_PI"]);
830  } else {
831  a[1] = m_mDouble["phi0"];
832  }
833  //a[2] = 1/m_mDouble["pt"]*m_mDouble["charge"];
834  a[2] = 1 / m_mDouble["pt"] * m_mDouble["charge2D"];
835  a[3] = m_mDouble["z0"];
836  a[4] = m_mDouble["cot"];
837  helix.a(a);
838  aTrack.setHelix(helix);
839  aTrack.set2DFitChi2(m_mDouble["fit2DChi2"]);
840  aTrack.set3DFitChi2(m_mDouble["zChi2"]);
841 
843  // Save values
844  if (m_mBool["fRootFile"]) {
845  if (m_fileFitter3D == 0) {
846  m_fileFitter3D = new TFile(m_rootFitter3DFileName.c_str(), "RECREATE");
847  HandleRoot::initializeRoot("fitter3D", &m_treeConstantsFitter3D, &m_treeTrackFitter3D,
848  m_mTVectorD, m_mTClonesArray,
849  m_mConstD, m_mConstV,
850  m_mDouble, m_mVector
851  );
852  }
853  HandleRoot::saveTrackValues("fitter3D",
854  m_mTClonesArray, m_mDouble, m_mVector
855  );
856  }
857 
858  } // End track loop
859 
860  // Save values to file
861  // To prevent crash when there are no tracks in first event.
862  // If there is no track in first evnet then the HandleRoot saving functions will fail.
863  // This is due to bad HandleRoot software design.
864  // The first event that have tracks will be event 0 in ROOT file.
865  if (m_mBool["fRootFile"]) {
866  if (m_fileFitter3D != 0) m_treeTrackFitter3D->Fill();
867  }
868 
869  if (m_mBool["debugFitted"]) {
870  for (unsigned iTrack = 0; iTrack < trackList.size(); iTrack++) {
871  TCTrack& aTrack = * trackList[iTrack];
872  if (aTrack.fitted() == 0) aTrack.setDebugValue(TRGCDCTrack::EDebugValueType::fit3D, 1);
873  }
874  }
875  if (m_mBool["debugLargeZ0"]) {
876  // If 3D fit z0 is larger or smaller than expected.
877  for (unsigned iTrack = 0; iTrack < trackList.size(); iTrack++) {
878  TCTrack& aTrack = *trackList[iTrack];
879  if (aTrack.fitted()) {
880  double fitZ0 = aTrack.helix().dz();
881  if (fitZ0 > 20 || fitZ0 < -20) aTrack.setDebugValue(TRGCDCTrack::EDebugValueType::fit3D, 1);
882  }
883  }
884  }
885 
886  TRGDebug::leaveStage("Fitter 3D");
887 
888  return 1;
889 
890  }
891 
892  double TRGCDCFitter3D::calPhi(TRGCDCSegmentHit const* segmentHit, double eventTime)
893  {
894  const CDC::CDCGeometryPar& cdcp = CDC::CDCGeometryPar::Instance();
895  unsigned localId = segmentHit->segment().center().localId();
896  unsigned layerId = segmentHit->segment().center().layerId();
897  int nWires = cdcp.nWiresInLayer(layerId) * 2;
898  double rr = cdcp.senseWireR(layerId);
899  double tdc = segmentHit->segment().priorityTime();
900  int lr = segmentHit->segment().LUT()->getValue(segmentHit->segment().lutPattern());
901  return Fitter3DUtility::calPhi(localId, nWires, tdc, eventTime, rr, lr);
902  }
903 
904 
908  // Functions for saving.
909  void TRGCDCFitter3D::saveVhdlAndCoe()
910  {
911  // Print Vhdl
912  if (m_commonData->getPrintedToFile() == 0) {
913  if (m_commonData->getPrintVhdl() == 0) {
914  m_commonData->setVhdlOutputFile("Fitter3D.vhd");
915  m_commonData->setPrintVhdl(1);
916  } else {
917  m_commonData->setPrintVhdl(0);
918  m_commonData->entryVhdlCode();
919  m_commonData->signalsVhdlCode();
920  m_commonData->buffersVhdlCode();
921  m_commonData->printToFile();
922  // Print LUTs.
923  for (map<string, TRGCDCJLUT*>::iterator it = m_mLutStorage.begin(); it != m_mLutStorage.end(); ++it) {
924  //it->second->makeCOE("./VHDL/LutData/"+it->first+".coe");
925  it->second->makeCOE(it->first + ".coe");
926  }
927  }
928  }
929  }
930 
931  void TRGCDCFitter3D::saveAllSignals()
932  {
933  // Save all signals to m_mSavedSignals. Limit to 10 times.
934  if ((*m_mSavedSignals.begin()).second.size() < 10 || m_mSavedSignals.empty()) {
935  for (auto it = m_mSignalStorage.begin(); it != m_mSignalStorage.end(); ++it) {
936  m_mSavedSignals[(*it).first].push_back((*it).second.getInt());
937  }
938  }
939  }
940 
941  void TRGCDCFitter3D::saveIoSignals()
942  {
943  // Save input output signals. Limit to 1024.
944  if ((*m_mSavedIoSignals.begin()).second.size() < 1025 || m_mSavedIoSignals.empty()) {
945  m_mSavedIoSignals["phi0"].push_back(m_mSignalStorage["phi0"].getInt());
946  for (unsigned iSt = 0; iSt < 4;
947  iSt++) m_mSavedIoSignals["phi_" + to_string(iSt)].push_back(m_mSignalStorage["phi_" + to_string(iSt)].getInt());
948  m_mSavedIoSignals["rho"].push_back(m_mSignalStorage["rho"].getInt());
949  for (unsigned iSt = 0; iSt < 4;
950  iSt++) m_mSavedIoSignals["iZError2_" + to_string(iSt)].push_back(m_mSignalStorage["iZError2_" + to_string(iSt)].getInt());
951  //m_mSavedIoSignals["charge"].push_back(m_mSignalStorage["charge"].getInt());
952  m_mSavedIoSignals["charge"].push_back(m_mSignalStorage["charge"].getInt());
953  m_mSavedIoSignals["z0"].push_back(m_mSignalStorage["z0"].getInt());
954  m_mSavedIoSignals["cot"].push_back(m_mSignalStorage["cot"].getInt());
955  m_mSavedIoSignals["chi2Sum"].push_back(m_mSignalStorage["chi2Sum"].getInt());
956  }
957  }
958 
959 
960  void TRGCDCFitter3D::getMCValues(const TRGCDC& m_cdc_in, TRGCDCTrack* aTrack, const std::map<std::string, double>& m_mConstD_in,
961  std::map<std::string, double>& m_mDouble_in, std::map<std::string, std::vector<double> >& m_mVector_in)
962  {
963  // Access to track's MC particle.
964  const TCRelation& trackRelation = aTrack->relation();
965  // Biggest contibutor is 0. Next is 1 and so on.
966  const MCParticle& trackMCParticle = trackRelation.mcParticle(0);
967 
968  // Calculated impact position
969  TVector3 vertex = trackMCParticle.getVertex();
970  TLorentzVector vector4 = trackMCParticle.get4Vector();
971  TVector2 helixCenter;
972  TVector3 impactPosition;
973  Fitter3DUtility::findImpactPosition(&vertex, &vector4, int(m_mDouble_in["mcCharge"]), helixCenter, impactPosition);
974  m_mVector_in["mcVertex"] = vector<double> ({vertex.X(), vertex.Y(), vertex.Z()});
975  m_mVector_in["mcMomentum"] = vector<double> ({vector4.Px(), vector4.Py(), vector4.Pz()});
976  m_mVector_in["helixCenter"] = vector<double> ({helixCenter.X(), helixCenter.Y()});
977  m_mVector_in["impactPosition"] = vector<double> ({impactPosition.X(), impactPosition.Y(), impactPosition.Z()});
978 
979  // Access track's particle parameters
980  m_mDouble_in["mcPt"] = trackMCParticle.getMomentum().Pt();
981  m_mDouble_in["mcPhi0"] = 0;
982  if (trackMCParticle.getCharge() > 0) m_mDouble_in["mcPhi0"] = trackMCParticle.getMomentum().Phi() - m_mConstD_in.at("Trg_PI") / 2;
983  if (trackMCParticle.getCharge() < 0) m_mDouble_in["mcPhi0"] = trackMCParticle.getMomentum().Phi() + m_mConstD_in.at("Trg_PI") / 2;
984  // Change range to [0,2pi]
985  if (m_mDouble_in["mcPhi0"] < 0) m_mDouble_in["mcPhi0"] += 2 * m_mConstD_in.at("Trg_PI");
986  //m_mDouble["mcZ0"] = trackMCParticle.getVertex().Z();
987  m_mDouble_in["mcZ0"] = impactPosition.Z();
988  m_mDouble_in["mcCot"] = trackMCParticle.getMomentum().Pz() / trackMCParticle.getMomentum().Pt();
989  m_mDouble_in["mcCharge"] = trackMCParticle.getCharge();
990 
991  // mcStatus[0]: statusbit, mcStatus[1]: pdg, mcStatus[2]: charge
992  TVectorD mcStatus(3);
993  m_mDouble_in["mcStatus"] = trackMCParticle.getStatus();
994  m_mDouble_in["pdgId"] = trackMCParticle.getPDG();
995 
996  // Find position of track for each super layer
997  //...G4 trackID...
998  unsigned id = trackRelation.contributor(0);
999  vector<const TCSHit*> mcAllTSList[9];
1000  vector<const TCSHit*> mcTSList(9);
1001  //...Segment loop...
1002  const vector<const TCSHit*> hits = m_cdc_in.segmentHits();
1003  for (unsigned i = 0; i < hits.size(); i++) {
1004  const TCSHit& ts = * hits[i];
1005  if (! ts.signal().active()) continue;
1006  const TCWHit* wh = ts.segment().center().hit();
1007  if (! wh) continue;
1008  const unsigned trackId = wh->iMCParticle();
1009  if (id == trackId)
1010  mcAllTSList[wh->wire().superLayerId()].push_back(& ts);
1011  }
1012  //...Select best one in each super layer...
1013  for (unsigned i = 0; i < 9; i++) {
1014  const TCSHit* best = 0;
1015  if (mcAllTSList[i].size() == 0) {
1016  } else if (mcAllTSList[i].size() == 1) {
1017  best = mcAllTSList[i][0];
1018  } else {
1019  int timeMin = 99999;
1020  for (unsigned k = 0; k < mcAllTSList[i].size(); k++) {
1021  const TRGSignal& timing = mcAllTSList[i][k]->signal();
1022  const TRGTime& t = * timing[0];
1023  if (t.time() < timeMin) {
1024  timeMin = t.time();
1025  best = mcAllTSList[i][k];
1026  }
1027  }
1028  }
1029  mcTSList[i] = best;
1030  }
1031 
1032 
1033  // Get mc track positions. Unit is cm.
1034  m_mVector_in["mcPosX"] = vector<double> ({9999, 9999, 9999, 9999, 9999, 9999, 9999, 9999, 9999});
1035  m_mVector_in["mcPosY"] = vector<double> ({9999, 9999, 9999, 9999, 9999, 9999, 9999, 9999, 9999});
1036  m_mVector_in["mcPosZ"] = vector<double> ({9999, 9999, 9999, 9999, 9999, 9999, 9999, 9999, 9999});
1037  for (unsigned iSL = 0; iSL < 9; iSL++) {
1038  if (mcTSList[iSL] != 0) {
1039  TVector3 posTrack = mcTSList[iSL]->simHit()->getPosTrack();
1040  m_mVector_in["mcPosX"][iSL] = posTrack.X();
1041  m_mVector_in["mcPosY"][iSL] = posTrack.Y();
1042  m_mVector_in["mcPosZ"][iSL] = posTrack.Z();
1043  }
1044  }
1045  // Get mc LR
1046  m_mVector_in["simMcLR"] = vector<double> (9);
1047  for (unsigned iSL = 0; iSL < 9; iSL++) {
1048  if (mcTSList[iSL] != 0) {
1049  m_mVector_in["simMcLR"][iSL] = mcTSList[iSL]->simHit()->getPosFlag();
1050  }
1051  }
1052 
1054  //for(unsigned iSL=0; iSL<9; iSL++){
1055  // if(mcTSList[iSL]!=0) {
1056  // TVectorD t_helixParameters;
1057  // TVector3 t_positionAtR;
1058  // TVector3 t_momentumAtR;
1059  // Fitter3DUtility::calHelixParameters(mcTSList[iSL]->simHit()->getPosIn(), mcTSList[iSL]->simHit()->getMomentum(),trackMCParticle.getCharge(),t_helixParameters);
1060  // //cout<<"dr: "<<helixParameters[0]<<" phi0: "<<helixParameters[1]<<" R: "<<1/helixParameters[2]/0.3/1.5*100<<" dz: "<<helixParameters[3]<<" tanLambda: "<<helixParameters[4]<<endl;
1061  // //calVectorsAtR(t_helixParameters, trackMCParticle.getCharge(), m_rr[iSL]*100, t_positionAtR, t_momentumAtR);
1062  // //cout<<" x: "<<t_positionAtR.X()<<" y: "<<t_positionAtR.Y()<<" z: "<<t_positionAtR.Z()<<endl;
1063  // //cout<<"Px: "<<t_momentumAtR.X()<<" Py: "<<t_momentumAtR.Y()<<" Pz: "<<t_momentumAtR.Z()<<endl;
1064  // }
1065  //}
1066  }
1067 
1068  bool TRGCDCFitter3D::isAxialTrackFull(const TRGCDCTrack& aTrack)
1069  {
1070  bool trackFull = 1;
1071  for (unsigned iAx = 0; iAx < 5; iAx++) {
1072  // Check if all superlayers have one TS
1073  const vector<TCLink*>& links = aTrack.links(iAx * 2);
1074  const unsigned nSegments = links.size();
1075  // Find if there is a TS with a priority hit.
1076  // Loop over all TS in same superlayer.
1077  bool priorityHitTS = 0;
1078  for (unsigned iTS = 0; iTS < nSegments; iTS++) {
1079  const TCSegment* t_segment = dynamic_cast<const TCSegment*>(& links[iTS]->hit()->cell());
1080  if (t_segment->center().hit() != 0) priorityHitTS = 1;
1081  }
1082  if (nSegments != 1) {
1083  if (nSegments == 0) {
1084  trackFull = 0;
1085  if (m_mBool["fIsPrintError"]) cout << "Fitter3D::isAxialTrackFull() => There are no TS." << endl;
1086  } else {
1087  if (m_mBool["fIsPrintError"]) cout << "Fitter3D::isAxialTrackFull() => multiple TS are assigned." << endl;
1088  }
1089  } else {
1090  if (priorityHitTS == 0) {
1091  trackFull = 0;
1092  if (m_mBool["fIsPrintError"]) cout << "Fitter3D::isAxialTrackFull() => There are no priority hit TS" << endl;
1093  }
1094  }
1095  } // End superlayer loop
1096  return trackFull;
1097  }
1098 
1099  bool TRGCDCFitter3D::isStereoTrackFull(const TRGCDCTrack& aTrack)
1100  {
1101  bool trackFull = 1;
1102  for (unsigned iSt = 0; iSt < 4; iSt++) {
1103  // Check if all superlayers have one TS
1104  const vector<TCLink*>& links = aTrack.links(iSt * 2 + 1);
1105  const unsigned nSegments = links.size();
1106  // Find if there is a TS with a priority hit.
1107  // Loop over all TS in same superlayer.
1108  bool priorityHitTS = 0;
1109  for (unsigned iTS = 0; iTS < nSegments; iTS++) {
1110  const TCSegment* t_segment = dynamic_cast<const TCSegment*>(& links[iTS]->hit()->cell());
1111  if (t_segment->center().hit() != 0) priorityHitTS = 1;
1112  }
1113  if (nSegments != 1) {
1114  if (nSegments == 0) {
1115  trackFull = 0;
1116  if (m_mBool["fIsPrintError"]) cout << "Fitter3D::isStereoTrackFull() => There are no TS." << endl;
1117  } else {
1118  if (m_mBool["fIsPrintError"]) cout << "Fitter3D::isStereoTrackFull() => multiple TS are assigned." << endl;
1119  }
1120  } else {
1121  if (priorityHitTS == 0) {
1122  trackFull = 0;
1123  if (m_mBool["fIsPrintError"]) cout << "Fitter3D::isStereoTrackFull() => There are no priority hit TS" << endl;
1124  }
1125  }
1126  } // End superlayer loop
1127  return trackFull;
1128  }
1129 
1130  void TRGCDCFitter3D::findHitAxialSuperlayers(const TRGCDCTrack& aTrack, vector<double>& useAxSl, bool printError)
1131  {
1132  useAxSl.assign(5, 1);
1133  for (unsigned iAx = 0; iAx < 5; iAx++) {
1134  // Check if all superlayers have one TS
1135  const vector<TCLink*>& links = aTrack.links(iAx * 2);
1136  const unsigned nSegments = links.size();
1137  // Find if there is a TS with a priority hit.
1138  // Loop over all TS in same superlayer.
1139  bool priorityHitTS = 0;
1140  //cout<<"iAx:"<<iAx<<" nSegments:"<<nSegments<<endl;
1141  for (unsigned iTS = 0; iTS < nSegments; iTS++) {
1142  const TCSegment* t_segment = dynamic_cast<const TCSegment*>(& links[iTS]->hit()->cell());
1143  if (t_segment->center().hit() != 0) priorityHitTS = 1;
1144  //cout<<" iTS:"<<iTS<<" priority:"<<t_segment->center().hit()<<" tsId:"<<t_segment->localId()<<endl;
1145  }
1146  if (nSegments != 1) {
1147  if (nSegments == 0) {
1148  useAxSl[iAx] = 0;
1149  } else {
1150  if (printError) cout << "Fitter3D::findHitAxialSuperlayers() => multiple TS are assigned." << endl;
1151  }
1152  } else {
1153  if (priorityHitTS == 0) {
1154  useAxSl[iAx] = 0;
1155  if (printError) cout << "Fitter3D::findHitAxialSuperlayers() => There are no priority hit TS" << endl;
1156  }
1157  }
1158  } // End superlayer loop
1159  }
1160 
1161  void TRGCDCFitter3D::findHitStereoSuperlayers(const TRGCDCTrack& aTrack, vector<double>& useStSl, bool printError)
1162  {
1163  useStSl.assign(4, 1);
1164  for (unsigned iSt = 0; iSt < 4; iSt++) {
1165  // Check if all superlayers have one TS
1166  const vector<TCLink*>& links = aTrack.links(iSt * 2 + 1);
1167  const unsigned nSegments = links.size();
1168  // Find if there is a TS with a priority hit.
1169  // Loop over all TS in same superlayer.
1170  bool priorityHitTS = 0;
1171  for (unsigned iTS = 0; iTS < nSegments; iTS++) {
1172  const TCSegment* t_segment = dynamic_cast<const TCSegment*>(& links[iTS]->hit()->cell());
1173  if (t_segment->center().hit() != 0) priorityHitTS = 1;
1174  }
1175  if (nSegments != 1) {
1176  if (nSegments == 0) {
1177  useStSl[iSt] = 0;
1178  } else {
1179  if (printError) cout << "Fitter3D::findHitStereoSuperlayers() => multiple TS are assigned." << endl;
1180  }
1181  } else {
1182  if (priorityHitTS == 0) {
1183  useStSl[iSt] = 0;
1184  if (printError) cout << "Fitter3D::findHitStereoSuperlayers() => There are no priority hit TS" << endl;
1185  }
1186  }
1187  } // End superlayer loop
1188  }
1189 
1190  void TRGCDCFitter3D::removeImpossibleStereoSuperlayers(vector<double>& useStSl)
1191  {
1192  for (unsigned iSt = 0; iSt < 4; iSt++) {
1193  // Check if rho is large enough for stereo super layer.
1194  if (2 * m_mDouble["rho"] < m_mConstV["rr3D"][iSt]) {
1195  useStSl[iSt] = 0;
1196  if (m_mBool["fIsPrintError"]) cout << "Fitter3D::removeImpossibleStereoSuperlayers() => pt is too low for SL." << endl;
1197  }
1198  } // End superlayer loop
1199  }
1200 
1201  void TRGCDCFitter3D::selectAxialTSs(const TRGCDCTrack& aTrack, vector<int>& bestTSIndex)
1202  {
1203  bestTSIndex.resize(5);
1204  std::fill_n(bestTSIndex.begin(), 5, -1);
1205  for (unsigned iAx = 0; iAx < 5; iAx++) {
1206  // Check if all superlayers have one TS
1207  const vector<TCLink*>& links = aTrack.links(iAx * 2);
1208  const unsigned nSegments = links.size();
1209  //cout<<"iAx:"<<iAx<<" nSegments:"<<nSegments<<endl;
1210  // tsCandiateInfo[iTS] = 1st priority, tdc
1211  vector<tuple<int, double> > tsCandiateInfo(nSegments);
1212  // Get information from candidates.
1213  for (unsigned iTS = 0; iTS < nSegments; iTS++) {
1214  const TCSegment* t_segment = dynamic_cast<const TCSegment*>(& links[iTS]->hit()->cell());
1215  int firstPriority = 0;
1216  if (t_segment->center().hit() != 0) firstPriority = 1;
1217  double tdc = t_segment->priorityTime();
1218  std::get<0>(tsCandiateInfo[iTS]) = firstPriority;
1219  std::get<1>(tsCandiateInfo[iTS]) = tdc;
1220  //cout<<"Ax:"<<iAx<<" iTS:"<<iTS<<" firstPriority:"<<firstPriority<<" tdc:"<<tdc<<endl;
1221  }
1222  // Select best candidate.
1223  // bestTS => tsIndex, { 1st priority, tdc }
1224  pair<int, tuple<int, double> > bestTS = make_pair(-1, make_tuple(-1, 9999));
1225  // If there is a candidate.
1226  if (tsCandiateInfo.size() != 0) {
1228  //if (std::get<0>(tsCandiateInfo[0]) == 1) {
1229  // bestTS.first = 0;
1230  // bestTS.second = tsCandiateInfo[0];
1231  //}
1233  //int randomIndex = gRandom->Integer(nSegments);
1234  //if (std::get<0>(tsCandiateInfo[randomIndex]) == 1) {
1235  // bestTS.first = randomIndex;
1236  // bestTS.second = tsCandiateInfo[randomIndex];
1237  //}
1238  // Selection type 1st priority.
1239  for (unsigned iTS = 0; iTS < nSegments; iTS++) {
1240  if (std::get<0>(tsCandiateInfo[iTS]) == 1) {
1241  bool select = 0;
1242  if (bestTS.first == -1) select = 1;
1243  else if (std::get<1>(bestTS.second) > std::get<1>(tsCandiateInfo[iTS])) select = 1;
1244  if (select) {
1245  bestTS.first = iTS;
1246  bestTS.second = tsCandiateInfo[iTS];
1247  }
1248  }
1249  }
1250  }
1251  //cout<<"Ax:"<<iAx<<" best: iTS:"<<bestTS.first<<" firstPriority:"<<std::get<0>(bestTS.second)<<" tdc:"<<std::get<1>(bestTS.second)<<endl;
1252  bestTSIndex[iAx] = bestTS.first;
1253  } // End superlayer loop
1254  }
1255 
1256  int TRGCDCFitter3D::do2DFit(TRGCDCTrack& aTrack, const std::map<std::string, bool>& m_mBool_in,
1257  const std::map<std::string, double>& m_mConstD_in,
1258  std::map<std::string, std::vector<double> >& m_mConstV_in, std::map<std::string, double>& m_mDouble_in,
1259  std::map<std::string, std::vector<double> >& m_mVector_in)
1260  {
1261  m_mVector_in["useAxSl"] = vector<double> (5);
1262  //findHitAxialSuperlayers(aTrack, useAxSl, m_mBool_in["fIsPrintError"]);
1263 
1264  // Find best TS between links for each SL.
1265  vector<int> bestTSIndex(5);
1266  selectAxialTSs(aTrack, bestTSIndex);
1267  for (unsigned iAx = 0; iAx < 5; iAx++) {
1268  if (bestTSIndex[iAx] != -1) m_mVector_in["useAxSl"][iAx] = 1;
1269  //cout<<"useAxSl["<<iAx<<"]:"<<useAxSl[iAx]<<endl;
1270  }
1271 
1272  // Check if number of axial super layer hits is smaller or equal to 1.
1273  m_mDouble_in["nHitAx"] = m_mVector_in["useAxSl"][0] + m_mVector_in["useAxSl"][1] + m_mVector_in["useAxSl"][2] +
1274  m_mVector_in["useAxSl"][3] +
1275  m_mVector_in["useAxSl"][4];
1276  if (m_mDouble_in["nHitAx"] <= 1) {
1277  if (m_mBool_in.at("fVerbose") == 1) cout << "[2DFit] Exiting because nHitAx is " << m_mDouble_in["nHitAx"] << endl;
1278  aTrack.setFitted(0);
1279  return 1;
1280  }
1281 
1282  // Fill information for axial layers
1283  m_mVector_in["tsId"] = vector<double> (9);
1284  m_mVector_in["tsId2D"] = vector<double> (5);
1285  m_mVector_in["wirePhi"] = vector<double> (9);
1286  m_mVector_in["lutLR"] = vector<double> (9);
1287  m_mVector_in["LR"] = vector<double> (9);
1288  m_mVector_in["driftLength"] = vector<double> (9);
1289  m_mVector_in["tdc"] = vector<double> (9);
1290  if (!m_mVector_in.count("mcLR")) m_mVector_in["mcLR"] = vector<double> (9);
1291  for (unsigned iAx = 0; iAx < 5; iAx++) {
1292  if (m_mVector_in["useAxSl"][iAx] == 1) {
1293  const vector<TCLink*>& links = aTrack.links(iAx * 2);
1294  //const TCSegment * t_segment = dynamic_cast<const TCSegment *>(& links[0]->hit()->cell());
1295  const TCSegment* t_segment = dynamic_cast<const TCSegment*>(& links[bestTSIndex[iAx]]->hit()->cell());
1296  m_mVector_in["tsId"][iAx * 2] = t_segment->localId();
1297  m_mVector_in["tsId2D"][iAx] = m_mVector_in["tsId"][iAx * 2];
1298  m_mVector_in["wirePhi"][iAx * 2] = (double) t_segment->localId() / m_mConstV_in["nWires"][iAx * 2] * 4 * m_mConstD_in.at("Trg_PI");
1299  m_mVector_in["lutLR"][iAx * 2] = t_segment->LUT()->getValue(t_segment->lutPattern());
1300  // mcLR should be removed.
1301  if (m_mBool_in.at("fMc")) m_mVector_in["mcLR"][iAx * 2] = t_segment->hit()->mcLR() + 1;
1302  m_mVector_in["driftLength"][iAx * 2] = t_segment->hit()->drift();
1303  m_mVector_in["tdc"][iAx * 2] = t_segment->priorityTime();
1304  if (m_mBool_in.at("fmcLR") == 1) m_mVector_in["LR"][iAx * 2] = m_mVector_in["mcLR"][iAx * 2];
1305  else if (m_mBool_in.at("fLRLUT") == 1) m_mVector_in["LR"][iAx * 2] = m_mVector_in["lutLR"][iAx * 2];
1306  else m_mVector_in["LR"][iAx * 2] = 3;
1307  } else {
1308  m_mVector_in["tsId"][iAx * 2] = 9999;
1309  m_mVector_in["wirePhi"][iAx * 2] = 9999;
1310  m_mVector_in["lutLR"][iAx * 2] = 0;
1311  // mcLR should be removed.
1312  if (m_mBool_in.at("fMc")) m_mVector_in["mcLR"][iAx * 2] = 9999;
1313  m_mVector_in["driftLength"][iAx * 2] = 9999;
1314  m_mVector_in["tdc"][iAx * 2] = 9999;
1315  if (m_mBool_in.at("fmcLR") == 1) m_mVector_in["LR"][iAx * 2] = 9999;
1316  else if (m_mBool_in.at("fLRLUT") == 1) m_mVector_in["LR"][iAx * 2] = 9999;
1317  else m_mVector_in["LR"][iAx * 2] = 9999;
1318  }
1319  } // End superlayer loop
1321  //int minTSTdc = 9999;
1322  //for(unsigned iAx=0; iAx<9; iAx++){
1323  // if (minTSTdc > m_mVector_in["tdc"][2*iAx]) minTSTdc = m_mVector_in["tdc"][2*iAx];
1324  //}
1325  //m_mDouble_in["eventTime"] = minTSTdc;
1326 
1328  // Get 2D fit values
1329  // Get 2D fit values from IW 2D fitter
1330  m_mDouble_in["phi02D"] = aTrack.helix().phi0();
1331  m_mDouble_in["pt2D"] = aTrack.pt();
1332  if (aTrack.charge() < 0) {
1333  m_mDouble_in["phi02D"] -= m_mConstD_in.at("Trg_PI");
1334  if (m_mDouble_in["phi02D"] < 0) m_mDouble_in["phi02D"] += 2 * m_mConstD_in.at("Trg_PI");
1335  }
1336  m_mDouble_in["dr2D"] = aTrack.helix().dr() * 0.01;
1337  // Get 2D fit values from JB 2D fitter
1338  // Currently using JB fitter for 3D fitting
1339  m_mDouble_in["charge"] = double(aTrack.charge());
1340  // Set phi2DError for 2D fit
1341  m_mVector_in["phi2DError"] = vector<double> (5);
1342  for (unsigned iAx = 0; iAx < 5; iAx++) {
1343  if (m_mVector_in["useAxSl"][iAx] == 1) {
1344  // Check LR.
1345  if (m_mVector_in["LR"][2 * iAx] != 3) m_mVector_in["phi2DError"][iAx] = m_mConstV_in["driftPhi2DError"][iAx];
1346  else m_mVector_in["phi2DError"][iAx] = m_mConstV_in["wirePhi2DError"][iAx];
1347  // Check event time.
1348  if (m_mDouble_in["eventTime"] == 9999) m_mVector_in["phi2DError"][iAx] = m_mConstV_in["wirePhi2DError"][iAx];
1349  } else {
1350  m_mVector_in["phi2DError"][iAx] = 9999;
1351  }
1352  }
1353  // Set invPhi2DError for 2D fit
1354  m_mVector_in["phi2DInvError"] = vector<double> (5);
1355  for (unsigned iAx = 0; iAx < 5; iAx++) {
1356  if (m_mVector_in["useAxSl"][iAx] == 1) {
1357  m_mVector_in["phi2DInvError"][iAx] = 1 / m_mVector_in["phi2DError"][iAx];
1358  } else {
1359  m_mVector_in["phi2DInvError"][iAx] = 0;
1360  }
1361  }
1362  // Calculate phi2D.
1363  m_mVector_in["phi2D"] = vector<double> (5);
1364  if (m_mBool_in.at("f2DFitDrift") == 0 || m_mDouble_in["eventTime"] == 9999) {
1365  for (unsigned iAx = 0; iAx < 5; iAx++) {
1366  m_mVector_in["phi2D"][iAx] = m_mVector_in["wirePhi"][iAx * 2];
1367  }
1368  } else {
1369  for (unsigned iAx = 0; iAx < 5; iAx++) {
1370  if (m_mVector_in["useAxSl"][iAx] == 1) {
1371  // Get drift length from table.
1372  string tableName = "driftLengthTableSL" + to_string(iAx * 2);
1373  double t_driftTime = m_mVector_in["tdc"][iAx * 2] - m_mDouble_in["eventTime"];
1374  if (t_driftTime < 0) t_driftTime = 0;
1375  if (t_driftTime > 511) t_driftTime = 511;
1376  double t_driftLength = m_mConstV_in[tableName][(unsigned)t_driftTime];
1377  m_mVector_in["phi2D"][iAx] = Fitter3DUtility::calPhi(m_mVector_in["wirePhi"][iAx * 2], t_driftLength, m_mConstV_in["rr"][iAx * 2],
1378  m_mVector_in["LR"][iAx * 2]);
1379  } else {
1380  m_mVector_in["phi2D"][iAx] = 9999;
1381  }
1382  }
1383  }
1384  // Fit2D
1385  if (m_mBool_in.at("f2DFit") == 0) {
1386  m_mDouble_in["rho"] = m_mDouble_in["pt2D"] / 0.01 / 1.5 / 0.299792458;
1387  m_mDouble_in["pt"] = 0.299792458 * 1.5 * m_mDouble_in["rho"] / 100;
1388  m_mDouble_in["phi0"] = m_mDouble_in["phi02D"];
1389  m_mDouble_in["fit2DChi2"] = 9999;
1390  } else {
1391  m_mDouble_in["rho"] = 0;
1392  m_mDouble_in["phi0"] = 0;
1393  m_mDouble_in["fit2DChi2"] = 0;
1394  Fitter3DUtility::rPhiFitter(&m_mConstV_in["rr2D"][0], &m_mVector_in["phi2D"][0], &m_mVector_in["phi2DInvError"][0],
1395  m_mDouble_in["rho"],
1396  m_mDouble_in["phi0"], m_mDouble_in["fit2DChi2"]);
1397  m_mDouble_in["pt"] = 0.3 * 1.5 * m_mDouble_in["rho"] / 100;
1398  }
1399 
1400  // Find charge of particle.
1401  Fitter3DUtility::chargeFinder(&m_mConstV_in["nTSs2D"][0], &m_mVector_in["tsId2D"][0], &m_mVector_in["useAxSl"][0],
1402  m_mDouble_in["phi0"],
1403  m_mDouble_in["charge"], m_mDouble_in["charge2D"]);
1404 
1405  if (m_mBool_in.at("fVerbose")) {
1406  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]f2DFit: " <<
1407  m_mBool_in.at("f2DFit") <<
1408  endl;
1409  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]evtTime: " <<
1410  m_mDouble_in["eventTime"]
1411  << endl;
1412  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]wirePhi: " <<
1413  m_mVector_in["wirePhi"][0]
1414  << " " << m_mVector_in["wirePhi"][1] << " " << m_mVector_in["wirePhi"][2] << " " << m_mVector_in["wirePhi"][3] << " " <<
1415  m_mVector_in["wirePhi"][4] << " " << m_mVector_in["wirePhi"][5] << " " << m_mVector_in["wirePhi"][6] << " " <<
1416  m_mVector_in["wirePhi"][7] << " "
1417  << m_mVector_in["wirePhi"][8] << endl;
1418  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]LR: " << int(
1419  m_mVector_in["LR"][0]) << " " << int(m_mVector_in["LR"][1]) << " " << int(m_mVector_in["LR"][2]) << " " << int(
1420  m_mVector_in["LR"][3]) << " " << int(m_mVector_in["LR"][4]) << " " << int(m_mVector_in["LR"][5]) << " " << int(
1421  m_mVector_in["LR"][6]) << " " << int(m_mVector_in["LR"][7]) << " " << int(m_mVector_in["LR"][8]) << endl;
1422  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]drift: " <<
1423  m_mVector_in["driftLength"][0] << " " << m_mVector_in["driftLength"][1] << " " << m_mVector_in["driftLength"][2] << " " <<
1424  m_mVector_in["driftLength"][3] << " " << m_mVector_in["driftLength"][4] << " " << m_mVector_in["driftLength"][5] << " " <<
1425  m_mVector_in["driftLength"][6] << " " << m_mVector_in["driftLength"][7] << " " << m_mVector_in["driftLength"][8] << endl;
1426  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]tdc: " <<
1427  m_mVector_in["tdc"][0] <<
1428  " " << m_mVector_in["tdc"][1] << " " << m_mVector_in["tdc"][2] << " " << m_mVector_in["tdc"][3] << " " << m_mVector_in["tdc"][4] <<
1429  " " <<
1430  m_mVector_in["tdc"][5] << " " << m_mVector_in["tdc"][6] << " " << m_mVector_in["tdc"][7] << " " << m_mVector_in["tdc"][8] << endl;
1431  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]rr2D: " <<
1432  m_mConstV_in["rr2D"][0] <<
1433  " " << m_mConstV_in["rr2D"][1] << " " << m_mConstV_in["rr2D"][2] << " " << m_mConstV_in["rr2D"][3] << " " << m_mConstV_in["rr2D"][4]
1434  << endl;
1435  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]Phi2D: " <<
1436  m_mVector_in["phi2D"][0]
1437  << " " << m_mVector_in["phi2D"][1] << " " << m_mVector_in["phi2D"][2] << " " << m_mVector_in["phi2D"][3] << " " <<
1438  m_mVector_in["phi2D"][4] <<
1439  endl;
1440  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]Phi2DInvError: " <<
1441  m_mVector_in["phi2DInvError"][0] << " " << m_mVector_in["phi2DInvError"][1] << " " << m_mVector_in["phi2DInvError"][2] << " " <<
1442  m_mVector_in["phi2DInvError"][3] << " " << m_mVector_in["phi2DInvError"][4] << endl;
1443  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]charge: " << int(
1444  m_mDouble_in["charge"]) << endl;
1445  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]charge2D: " << int(
1446  m_mDouble_in["charge2D"]) << endl;
1447  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]pt: " <<
1448  m_mDouble_in["pt"] <<
1449  endl;
1450  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]rho: " <<
1451  m_mDouble_in["rho"] <<
1452  endl;
1453  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]phi0: " <<
1454  m_mDouble_in["phi0"] <<
1455  " " << m_mDouble_in["phi0"] / m_mConstD_in.at("Trg_PI") * 180 << endl;
1456  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]fit2DChi2: " <<
1457  m_mDouble_in["fit2DChi2"]
1458  << endl;
1459  cout << "[E" << int(m_mDouble_in["eventNumber"]) << "][T" << int(m_mDouble_in["trackId"]) << "]useAxSl: " << int(
1460  m_mVector_in["useAxSl"][0]) << " " << int(m_mVector_in["useAxSl"][1]) << " " << int(m_mVector_in["useAxSl"][2]) << " " << int(
1461  m_mVector_in["useAxSl"][3]) << endl;
1462  }
1463 
1464  if (std::isnan(m_mDouble_in["rho"]) || std::isnan(m_mDouble_in["phi0"])) {
1465  if (m_mBool_in.at("fVerbose") == 1) cout << "[2Dfit] Exiting because rho or phi0 is nan." << endl;
1466  return 2;
1467  }
1468  return 0;
1469  }
1470 
1471  void TRGCDCFitter3D::print3DInformation(int iTrack)
1472  {
1473  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]evtTime: " << m_mDouble["eventTime"] << endl;
1474  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]wirePhi: " << m_mVector["wirePhi"][0] << " " <<
1475  m_mVector["wirePhi"][1] << " " << m_mVector["wirePhi"][2] << " " << m_mVector["wirePhi"][3] << " " << m_mVector["wirePhi"][4] << " "
1476  << m_mVector["wirePhi"][5] << " " << m_mVector["wirePhi"][6] << " " << m_mVector["wirePhi"][7] << " " << m_mVector["wirePhi"][8] <<
1477  endl;
1478  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]LR: " << int(m_mVector["LR"][0]) << " " << int(
1479  m_mVector["LR"][1]) << " " << int(m_mVector["LR"][2]) << " " << int(m_mVector["LR"][3]) << " " << int(
1480  m_mVector["LR"][4]) << " " << int(m_mVector["LR"][5]) << " " << int(m_mVector["LR"][6]) << " " << int(
1481  m_mVector["LR"][7]) << " " << int(m_mVector["LR"][8]) << endl;
1482  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]lutLR: " << int(m_mVector["lutLR"][0]) << " " << int(
1483  m_mVector["lutLR"][1]) << " " << int(m_mVector["lutLR"][2]) << " " << int(m_mVector["lutLR"][3]) << " " << int(
1484  m_mVector["lutLR"][4]) << " " << int(m_mVector["lutLR"][5]) << " " << int(m_mVector["lutLR"][6]) << " " << int(
1485  m_mVector["lutLR"][7]) << " " << int(m_mVector["lutLR"][8]) << endl;
1486  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]useStSl: " << int(m_mVector["useStSl"][0]) << " " << int(
1487  m_mVector["useStSl"][1]) << " " << int(m_mVector["useStSl"][2]) << " " << int(m_mVector["useStSl"][3]) << endl;
1488  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]drift: " << m_mVector["driftLength"][0] << " " <<
1489  m_mVector["driftLength"][1] << " " << m_mVector["driftLength"][2] << " " << m_mVector["driftLength"][3] << " " <<
1490  m_mVector["driftLength"][4] << " " << m_mVector["driftLength"][5] << " " << m_mVector["driftLength"][6] << " " <<
1491  m_mVector["driftLength"][7] << " " << m_mVector["driftLength"][8] << endl;
1492  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]tdc: " << m_mVector["tdc"][0] << " " <<
1493  m_mVector["tdc"][1] << " " << m_mVector["tdc"][2] << " " << m_mVector["tdc"][3] << " " << m_mVector["tdc"][4] << " " <<
1494  m_mVector["tdc"][5] << " " << m_mVector["tdc"][6] << " " << m_mVector["tdc"][7] << " " << m_mVector["tdc"][8] << endl;
1495  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]Phi2D: " << m_mVector["phi2D"][0] << " " <<
1496  m_mVector["phi2D"][1] << " " << m_mVector["phi2D"][2] << " " << m_mVector["phi2D"][3] << " " << m_mVector["phi2D"][4] << endl;
1497  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]Phi3D: " << m_mVector["phi3D"][0] << " " <<
1498  m_mVector["phi3D"][1] << " " << m_mVector["phi3D"][2] << " " << m_mVector["phi3D"][3] << endl;
1499  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]zz: " << m_mVector["zz"][0] << " " << m_mVector["zz"][1]
1500  << " " << m_mVector["zz"][2] << " " << m_mVector["zz"][3] << endl;
1501  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]arcS: " << m_mVector["arcS"][0] << " " <<
1502  m_mVector["arcS"][1] << " " << m_mVector["arcS"][2] << " " << m_mVector["arcS"][3] << endl;
1503  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]zerror: " << m_mVector["zError"][0] << " " <<
1504  m_mVector["zError"][1] << " " << m_mVector["zError"][2] << " " << m_mVector["zError"][3] << endl;
1505  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]izerror: " << m_mVector["iZError2"][0] << " " <<
1506  m_mVector["iZError2"][1] << " " << m_mVector["iZError2"][2] << " " << m_mVector["iZError2"][3] << endl;
1507  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]charge: " << int(m_mDouble["charge"]) << endl;
1508  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]pt: " << m_mDouble["pt"] << endl;
1509  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]phi0: " << m_mDouble["phi0"] << endl;
1510  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]z0: " << m_mDouble["z0"] << endl;
1511  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]cot: " << m_mDouble["cot"] << endl;
1512  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]chi2: " << m_mDouble["zChi2"] << endl;
1513  if (m_mBool["fMc"]) {
1514  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]mcCharge: " << int(m_mDouble["mcCharge"]) << endl;
1515  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]mcPosZ: " << m_mVector["mcPosZ"][1] << " " <<
1516  m_mVector["mcPosZ"][3] << " " << m_mVector["mcPosZ"][5] << " " << m_mVector["mcPosZ"][7] << endl;
1517  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]mcPosZ: " << m_mVector["mcPosZ"][1] << " " <<
1518  m_mVector["mcPosZ"][3] << " " << m_mVector["mcPosZ"][5] << " " << m_mVector["mcPosZ"][7] << endl;
1519  cout << "[E" << int(m_mDouble["eventNumber"]) << "][T" << iTrack << "]mcLR: " << int(m_mVector["mcLR"][0]) << " " << int(
1520  m_mVector["mcLR"][1]) << " " << int(m_mVector["mcLR"][2]) << " " << int(m_mVector["mcLR"][3]) << " " << int(
1521  m_mVector["mcLR"][4]) << " " << int(m_mVector["mcLR"][5]) << " " << int(m_mVector["mcLR"][6]) << " " << int(
1522  m_mVector["mcLR"][7]) << " " << int(m_mVector["mcLR"][8]) << endl;
1523  }
1524  }
1525 
1526 
1527  void TRGCDCFitter3D::terminate()
1528  {
1529  if (m_mBool["fRootFile"]) {
1530  HandleRoot::writeRoot(m_fileFitter3D);
1531  HandleRoot::terminateRoot(m_mTVectorD, m_mTClonesArray);
1532  delete m_fileFitter3D;
1533  }
1534 
1535  if (m_mBool["fVHDLFile"]) {
1536  //if(m_mSavedIoSignals.size()!=0) FpgaUtility::multipleWriteCoe(10, m_mSavedIoSignals, "./VHDL/LutData/Io/");
1537  //if(m_mSavedSignals.size()!=0) FpgaUtility::writeSignals("./VHDL/signals",m_mSavedSignals);
1538  if (m_mSavedIoSignals.size() != 0) FpgaUtility::multipleWriteCoe(10, m_mSavedIoSignals, "./");
1539  if (m_mSavedSignals.size() != 0) FpgaUtility::writeSignals("signalValues", m_mSavedSignals);
1540  }
1541 
1542  if (m_commonData) delete m_commonData;
1543  }
1544 
1545  string TRGCDCFitter3D::version(void) const
1546  {
1547  return string("TRGCDCFitter3D 6.0");
1548  }
1549 
1550  std::string TRGCDCFitter3D::name(void) const
1551  {
1552  return m_name;
1553  }
1554 
1555  void TRGCDCFitter3D::getStereoGeometry(map<string, vector<double> >& stGeometry)
1556  {
1557  stGeometry["priorityLayer"] = {10, 22, 34, 46};
1558  stGeometry["nWires"] = vector<double> (4);
1559  stGeometry["cdcRadius"] = vector<double> (4);
1560  stGeometry["zToStraw"] = vector<double> (4);
1561  stGeometry["nShift"] = vector<double> (4);
1562  stGeometry["angleSt"] = vector<double> (4);
1563  const CDC::CDCGeometryPar& cdc = CDC::CDCGeometryPar::Instance();
1564  for (int iSt = 0; iSt < 4; ++iSt) {
1565  stGeometry["nWires"][iSt] = cdc.nWiresInLayer(stGeometry["priorityLayer"][iSt]) * 2;
1566  stGeometry["cdcRadius"][iSt] = cdc.senseWireR(stGeometry["priorityLayer"][iSt]);
1567  stGeometry["zToStraw"][iSt] = cdc.senseWireBZ(stGeometry["priorityLayer"][iSt]);
1568  stGeometry["nShift"][iSt] = cdc.nShifts(stGeometry["priorityLayer"][iSt]);
1569  stGeometry["angleSt"][iSt] = 2 * stGeometry["cdcRadius"][iSt] * sin(M_PI * stGeometry["nShift"][iSt] /
1570  (stGeometry["nWires"][iSt])) /
1571  (cdc.senseWireFZ(stGeometry["priorityLayer"][iSt]) - stGeometry["zToStraw"][iSt]);
1572  }
1573  }
1574 
1575  void TRGCDCFitter3D::getStereoXt(vector<double> const& stPriorityLayer, vector<vector<double> >& stXts, bool isSimple)
1576  {
1577  stXts.resize(stPriorityLayer.size(), vector<double> (512));
1578  const CDC::CDCGeometryPar& cdc = CDC::CDCGeometryPar::Instance();
1579  for (unsigned iSt = 0; iSt < stPriorityLayer.size(); ++iSt) {
1580  for (unsigned iTick = 0; iTick < stXts[iSt].size(); ++iTick) {
1581  double t = iTick * 2 * cdc.getTdcBinWidth();
1582  if (isSimple) {
1583  stXts[iSt][iTick] = cdc.getNominalDriftV() * t;
1584  } else {
1585  double driftLength_0 = cdc.getDriftLength(t, stPriorityLayer[iSt], 0);
1586  double driftLength_1 = cdc.getDriftLength(t, stPriorityLayer[iSt], 1);
1587  stXts[iSt][iTick] = (driftLength_0 + driftLength_1) / 2;
1588  }
1589  }
1590  }
1591  }
1592 
1593  void TRGCDCFitter3D::getConstants(std::map<std::string, double>& mConstD, std::map<std::string, std::vector<double> >& mConstV,
1594  bool isXtSimple)
1595  {
1596  const CDC::CDCGeometryPar& cdc = CDC::CDCGeometryPar::Instance();
1597  mConstD["Trg_PI"] = 3.141592653589793;
1598  mConstV["priorityLayer"] = {3, 10, 16, 22, 28, 34, 40, 46, 52};
1599  mConstV["rr"] = vector<double> (9);
1600  mConstV["nWires"] = vector<double> (9);
1601  mConstV["nTSs"] = vector<double> (9);
1602  for (unsigned iSL = 0; iSL < 9; iSL++) {
1603  unsigned t_layerId = mConstV["priorityLayer"][iSL];
1604  mConstV["rr"][iSL] = cdc.senseWireR(t_layerId);
1605  mConstV["nWires"][iSL] = cdc.nWiresInLayer(t_layerId) * 2;
1606  mConstV["nTSs"][iSL] = cdc.nWiresInLayer(t_layerId);
1607  }
1608  mConstV["nTSs2D"] = vector<double> (5);
1609  for (unsigned iAx = 0; iAx < 5; iAx++) {
1610  mConstV["nTSs2D"][iAx] = mConstV["nTSs"][2 * iAx];
1611  }
1612  mConstV["zToStraw"] = vector<double> (4);
1613  mConstV["zToOppositeStraw"] = vector<double> (4);
1614  mConstV["angleSt"] = vector<double> (4);
1615  mConstV["nShift"] = vector<double> (4);
1616  for (int iSt = 0; iSt < 4; iSt++) {
1617  unsigned t_layerId = mConstV["priorityLayer"][iSt * 2 + 1];
1618  mConstV["zToStraw"][iSt] = cdc.senseWireBZ(t_layerId);
1619  mConstV["zToOppositeStraw"][iSt] = cdc.senseWireFZ(t_layerId);
1620  mConstV["angleSt"][iSt] = 2 * mConstV["rr"][2 * iSt + 1] * sin(mConstD["Trg_PI"] * cdc.nShifts(t_layerId) /
1621  (2 * cdc.nWiresInLayer(t_layerId))) / (cdc.senseWireFZ(t_layerId) - cdc.senseWireBZ(t_layerId));
1622  mConstV["nShift"][iSt] = cdc.nShifts(t_layerId);
1623  }
1624 
1625  mConstV["rr2D"] = vector<double> (5);
1626  mConstV["rr3D"] = vector<double> (4);
1627  for (int iAx = 0; iAx < 5; iAx++) mConstV["rr2D"][iAx] = mConstV["rr"][iAx * 2];
1628  for (int iSt = 0; iSt < 4; iSt++) mConstV["rr3D"][iSt] = mConstV["rr"][iSt * 2 + 1];
1629 
1630  mConstV["wirePhi2DError"] = vector<double> (5);
1631  mConstV["driftPhi2DError"] = vector<double> (5);
1632  mConstV["wirePhi2DError"][0] = 0.00085106;
1633  mConstV["wirePhi2DError"][1] = 0.00039841;
1634  mConstV["wirePhi2DError"][2] = 0.00025806;
1635  mConstV["wirePhi2DError"][3] = 0.00019084;
1636  mConstV["wirePhi2DError"][4] = 0.0001514;
1637  mConstV["driftPhi2DError"][0] = 0.00085106;
1638  mConstV["driftPhi2DError"][1] = 0.00039841;
1639  mConstV["driftPhi2DError"][2] = 0.00025806;
1640  mConstV["driftPhi2DError"][3] = 0.00019084;
1641  mConstV["driftPhi2DError"][4] = 0.0001514;
1642  mConstV["driftZError"] = vector<double> ({0.7676, 0.9753, 1.029, 1.372});
1643  mConstV["wireZError"] = vector<double> ({0.7676, 0.9753, 1.029, 1.372});
1644 
1645  // Make driftLength table for each superlayer. Up to 511 clock ticks.
1646  // driftLengthTableSLX[ tdcCount (~2ns unit) ] = drift length (cm)
1647  for (unsigned iSl = 0; iSl < 9; iSl++) {
1648  string tableName = "driftLengthTableSL" + to_string(iSl);
1649  unsigned tableSize = 512;
1650  mConstV[tableName] = vector<double> (tableSize);
1651  unsigned t_layer = mConstV["priorityLayer"][iSl];
1652  for (unsigned iTick = 0; iTick < tableSize; iTick++) {
1653  double t_driftTime = iTick * 2 * cdc.getTdcBinWidth();
1654  double avgDriftLength = 0;
1655  if (isXtSimple == 1) {
1656  avgDriftLength = cdc.getNominalDriftV() * t_driftTime;
1657  } else {
1658  double driftLength_0 = cdc.getDriftLength(t_driftTime, t_layer, 0);
1659  double driftLength_1 = cdc.getDriftLength(t_driftTime, t_layer, 1);
1660  avgDriftLength = (driftLength_0 + driftLength_1) / 2;
1661  }
1662  mConstV[tableName][iTick] = avgDriftLength;
1663  }
1664  }
1665 
1666  mConstD["tdcBitSize"] = 9;
1667  mConstD["rhoBitSize"] = 11;
1668  mConstD["iError2BitSize"] = 8;
1669  mConstD["iError2Max"] = 1 / pow(mConstV["wireZError"][0], 2);
1670  // LUT values
1671  mConstD["JB"] = 0;
1672  mConstD["phiMax"] = mConstD["Trg_PI"];
1673  mConstD["phiMin"] = -mConstD["Trg_PI"];
1674  mConstD["rhoMin"] = 20;
1675  mConstD["rhoMax"] = 2500;
1676  mConstD["phiBitSize"] = 13;
1677  mConstD["driftPhiLUTOutBitSize"] = mConstD["phiBitSize"] - 1;
1678  mConstD["driftPhiLUTInBitSize"] = mConstD["tdcBitSize"];
1679  mConstD["acosLUTOutBitSize"] = mConstD["phiBitSize"] - 1;
1680  mConstD["acosLUTInBitSize"] = mConstD["rhoBitSize"];
1681  mConstD["zLUTInBitSize"] = mConstD["phiBitSize"];
1682  mConstD["zLUTOutBitSize"] = 9;
1683  mConstD["iDenLUTInBitSize"] = 13;
1684  mConstD["iDenLUTOutBitSize"] = 11; // To increase wireZError = 2.5*driftZError
1685  }
1686 
1688 } // namespace Belle2
Belle2::CDC::CDCGeometryPar
The Class for CDC Geometry Parameters.
Definition: CDCGeometryPar.h:73
Belle2::CDC::CDCGeometryPar::nShifts
int nShifts(int layerId) const
Returns number shift.
Definition: CDCGeometryPar.h:1212
Belle2::CDC::CDCGeometryPar::getTdcBinWidth
double getTdcBinWidth() const
Return TDC bin width (nsec).
Definition: CDCGeometryPar.h:744
Belle2::CDC::CDCGeometryPar::senseWireBZ
double senseWireBZ(int layerId) const
Returns backward z position of sense wire in each layer.
Definition: CDCGeometryPar.h:1252
Belle2::CDC::CDCGeometryPar::nWiresInLayer
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
Definition: CDCGeometryPar.h:1222
Belle2::CDC::CDCGeometryPar::senseWireFZ
double senseWireFZ(int layerId) const
Returns forward z position of sense wire in each layer.
Definition: CDCGeometryPar.h:1247
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:754
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:2161
Belle2::CDC::CDCGeometryPar::Instance
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
Definition: CDCGeometryPar.cc:39
Belle2::CDC::CDCGeometryPar::senseWireR
double senseWireR(int layerId) const
Returns radius of sense wire in each layer.
Definition: CDCGeometryPar.h:1242
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
Belle2::MCParticle::getVertex
TVector3 getVertex() const
Return production vertex position, shorthand for getProductionVertex().
Definition: MCParticle.h:183
Belle2::MCParticle::getMomentum
TVector3 getMomentum() const
Return momentum.
Definition: MCParticle.h:198
Belle2::MCParticle::getStatus
unsigned int getStatus(unsigned short int bitmask=USHRT_MAX) const
Return status code of particle.
Definition: MCParticle.h:122
Belle2::MCParticle::getCharge
float getCharge() const
Return the particle charge defined in TDatabasePDG.
Definition: MCParticle.cc:34
Belle2::MCParticle::getPDG
int getPDG() const
Return PDG code of particle.
Definition: MCParticle.h:112
Belle2::MCParticle::get4Vector
TLorentzVector get4Vector() const
Return 4Vector of particle.
Definition: MCParticle.h:207
Belle2::StoreAccessorBase::isRequired
bool isRequired(const std::string &name="")
Ensure this array/object has been registered previously.
Definition: StoreAccessorBase.h:78
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:95
Belle2::TRGCDCFitter3D::m_mTVectorD
std::map< std::string, TVectorD * > m_mTVectorD
TVectorD map for saving values to root file.
Definition: Fitter3D.h:159
Belle2::TRGCDCFitter3D::m_mVector
std::map< std::string, std::vector< double > > m_mVector
Map to hold vector values for Fitter3D.
Definition: Fitter3D.h:131
Belle2::TRGCDCFitter3D::m_mBool
std::map< std::string, bool > m_mBool
Map to hold input options.
Definition: Fitter3D.h:137
Belle2::TRGCDCFitter3D::m_mConstD
std::map< std::string, double > m_mConstD
Map to hold constant values for Fitter3D.
Definition: Fitter3D.h:133
Belle2::TRGCDCFitter3D::m_mSavedSignals
std::map< std::string, std::vector< signed long long > > m_mSavedSignals
Array of saved signals.
Definition: Fitter3D.h:163
Belle2::TRGCDCFitter3D::m_mLutStorage
std::map< std::string, TRGCDCJLUT * > m_mLutStorage
Map to hold JLuts.
Definition: Fitter3D.h:142
Belle2::TRGCDCFitter3D::m_mConstV
std::map< std::string, std::vector< double > > m_mConstV
Map to hold constant vectcors for Fitter3D.
Definition: Fitter3D.h:135
Belle2::TRGCDCFitter3D::m_mDouble
std::map< std::string, double > m_mDouble
Map to hold double values for Fitter3D.
Definition: Fitter3D.h:129
Belle2::TRGCDCFitter3D::m_name
const std::string m_name
Name.
Definition: Fitter3D.h:123
Belle2::TRGCDCFitter3D::m_mSignalStorage
std::map< std::string, TRGCDCJSignal > m_mSignalStorage
Map to hold JSignals.
Definition: Fitter3D.h:140
Belle2::TRGCDCFitter3D::m_mTClonesArray
std::map< std::string, TClonesArray * > m_mTClonesArray
TClonesArray map for saving values to root file.
Definition: Fitter3D.h:161
Belle2::TRGCDCFitter3D::m_fileFitter3D
TFile * m_fileFitter3D
Tfile for Fitter3D root file.
Definition: Fitter3D.h:150
Belle2::TRGCDCFitter3D::m_commonData
TRGCDCJSignalData * m_commonData
For VHDL code.
Definition: Fitter3D.h:144
Belle2::TRGCDCFitter3D::m_treeTrackFitter3D
TTree * m_treeTrackFitter3D
TTree for tracks of fitter3D.
Definition: Fitter3D.h:153
Belle2::TRGCDCFitter3D::m_mSavedIoSignals
std::map< std::string, std::vector< signed long long > > m_mSavedIoSignals
Array of I/O signals.
Definition: Fitter3D.h:165
Belle2::TRGCDCFitter3D::m_treeConstantsFitter3D
TTree * m_treeConstantsFitter3D
TTree for constants of fitter3D.
Definition: Fitter3D.h:156
Belle2::TRGCDCFitter3D::m_rootFitter3DFileName
std::string m_rootFitter3DFileName
Members for saving.
Definition: Fitter3D.h:147
Belle2::TRGCDCFitter3D::m_cdc
const TRGCDC & m_cdc
CDCTRG.
Definition: Fitter3D.h:126
Belle2::TRGCDCHelix
TRGCDCHelix parameter class.
Definition: Helix.h:34
Belle2::TRGCDCJSignalData
A class to hold common data for JSignals.
Definition: JSignalData.h:33
Belle2::TRGCDCSegmentHit
A class to represent a track segment hit in CDC.
Definition: SegmentHit.h:31
Belle2::TRGCDCSegment::priorityTime
float priorityTime(void) const
return priority time in TSHit.
Belle2::TRGCDCTrack
A class to represent a reconstructed charged track in TRGCDC.
Definition: TRGCDCTrack.h:38
Belle2::TRGCDC
The instance of TRGCDC is a singleton.
Definition: TRGCDC.h:69
Belle2::TRGSignal
A class to represent a digitized signal. Unit is nano second.
Definition: Signal.h:23
Belle2::TRGTime
A class to represent a signal timing in the trigger system.
Definition: Time.h:25
Fitter3DUtility::rotatePhi
static double rotatePhi(double value, double refPhi)
Rotates to range [-pi, pi].
Definition: Fitter3DUtility.cc:616
Fitter3DUtility::findQuadrant
static int findQuadrant(double value)
Finds quadrant of angle. Angle is in rad.
Definition: Fitter3DUtility.cc:649
Fitter3DUtility::rPhiFitter
static void rPhiFitter(double *rr, double *phi2, double *invphierror, double &rho, double &myphi0)
A circle fitter with invPhiError without fit chi2 output.
Definition: Fitter3DUtility.cc:76
Fitter3DUtility::calZ
static double calZ(int charge, double anglest, double ztostraw, double rr, double phi, double rho, double phi0)
Calculates z.
Definition: Fitter3DUtility.cc:903
Fitter3DUtility::setError
static void setError(std::map< std::string, double > const &mConstD, std::map< std::string, std::vector< double > > const &mConstV, std::map< std::string, Belle2::TRGCDCJSignal > &mSignalStorage)
Sets error using JSignal class.
Definition: Fitter3DUtility.cc:690
Fitter3DUtility::calS
static double calS(double rho, double rr)
Calculates arc length.
Definition: Fitter3DUtility.cc:1216
Fitter3DUtility::calPhi
static double calPhi(double wirePhi, double driftLength, double rr, int lr)
Pre 3D fitter functions. rr is in cm scale. driftLength is in cm scale.
Definition: Fitter3DUtility.cc:246
Fitter3DUtility::chargeFinder
static void chargeFinder(double *nTSs, double *tsIds, double *tsHitmap, double phi0, double inCharge, double &outCharge)
Charge finder using circle fitter output and axial TSs.
Definition: Fitter3DUtility.cc:147
Fitter3DUtility::rotateTsId
static int rotateTsId(int value, int refId, int nTSs)
Rotates to range [0, nTSs-1].
Definition: Fitter3DUtility.cc:637
Fitter3DUtility::toUnsignedTdc
static unsigned toUnsignedTdc(int tdc, int nBits)
Changes tdc and event time to unsigned value that has # bits.
Definition: Fitter3DUtility.cc:233
Fitter3DUtility::rSFit
static void rSFit(double *iezz2, double *arcS, double *zz, double &z0, double &cot, double &zchi2)
3D fitter functions Fits z and arc S.
Definition: Fitter3DUtility.cc:1315
Fitter3DUtility::findImpactPosition
static void findImpactPosition(TVector3 *mcPosition, TLorentzVector *mcMomentum, int charge, TVector2 &helixCenter, TVector3 &impactPosition)
MC calculation functions Calculates the impact position of track.
Definition: Fitter3DUtility.cc:1934
Fitter3DUtility::constrainRPerStSl
static void constrainRPerStSl(std::map< std::string, std::vector< double > > const &mConstV, std::map< std::string, Belle2::TRGCDCJSignal > &mSignalStorage)
Constrains R for each SL differently using JSignal and multiple constants.
Definition: Fitter3DUtility.cc:836
FpgaUtility::writeSignals
static void writeSignals(std::string outFilePath, std::map< std::string, std::vector< signed long long > > &data)
COE file functions.
Definition: FpgaUtility.cc:202
FpgaUtility::multipleWriteCoe
static void multipleWriteCoe(int lutInBitsize, std::map< std::string, std::vector< signed long long > > &data, const std::string &fileDirectory)
Writes multiple signal values to a file in coe format.
Definition: FpgaUtility.cc:220
Belle2::TRGCDCTrack::helix
const TRGCDCHelix & helix(void) const
returns helix parameter.
Definition: TRGCDCTrack.h:143
Belle2::TRGCDCTrackBase::relation
const TRGCDCRelation relation(void) const
returns MC information.
Definition: TrackBase.cc:142
Belle2::TRGCDC::segmentHits
std::vector< const TRGCDCSegmentHit * > segmentHits(void) const
returns a list of TRGCDCSegmentHit.
Definition: TRGCDC.h:996
Belle2::TRGCDCSegment::center
const TRGCDCWire & center(void) const
returns a center wire.
Definition: Segment.h:250
Belle2::TRGCDCJSignalData::printToFile
void printToFile()
Utilities Function to print VHDL code.
Definition: JSignalData.cc:106
Belle2::TRGCDCFitter3D::findHitAxialSuperlayers
static void findHitAxialSuperlayers(const TRGCDCTrack &aTrack, std::vector< double > &useAxSL, bool printError)
Finds which axial superlayers has TSs. useAxSL array indicating hit superlayers.
Definition: Fitter3D.cc:1130
Belle2::TRGCDCFitter3D::getStereoGeometry
static void getStereoGeometry(std::map< std::string, std::vector< double > > &stGeometry)
Get stereo geometry.
Definition: Fitter3D.cc:1555
Belle2::TRGCDCFitter3D::initialize
void initialize()
Initialization.
Definition: Fitter3D.cc:75
Belle2::TRGCDCJSignal::valuesToMapSignals
static void valuesToMapSignals(std::vector< std::tuple< std::string, double, int, double, double, int > > const &inValues, Belle2::TRGCDCJSignalData *inCommonData, std::map< std::string, Belle2::TRGCDCJSignal > &outMap)
Values => [name, value, bitwidth, min, max, clock] Changes values to signals.
Definition: JSignal.cc:2206
Belle2::TRGCDCTrackBase::charge
double charge(void) const
returns charge.
Definition: TrackBase.h:248
Belle2::TRGCDCFitter3D::calPhi
static double calPhi(TRGCDCSegmentHit const *segmentHit, double eventTime)
Utility functions.
Definition: Fitter3D.cc:892
Belle2::TRGCDCFitter3D::isStereoTrackFull
bool isStereoTrackFull(const TRGCDCTrack &aTrack)
Checks if stereo track has 4 TSs. One per each superlayer.
Definition: Fitter3D.cc:1099
Belle2::TRGCDCFitter3D::saveVhdlAndCoe
void saveVhdlAndCoe()
Functions for saving.
Definition: Fitter3D.cc:909
Belle2::TRGCDCJSignalData::getPrintVhdl
bool getPrintVhdl() const
Gets the status of m_printVhdl.
Definition: JSignalData.cc:75
Belle2::TRGCDCFitter3D::doit
int doit(std::vector< TRGCDCTrack * > &trackList)
Does track fitting.
Definition: Fitter3D.cc:210
Belle2::TRGCDCFitter3D::getStereoXt
static void getStereoXt(std::vector< double > const &stPriorityLayer, std::vector< std::vector< double > > &stXts, bool isSimple=0)
Get stereo Xt.
Definition: Fitter3D.cc:1575
Belle2::TRGCDCFitter3D::removeImpossibleStereoSuperlayers
void removeImpossibleStereoSuperlayers(std::vector< double > &useStSL)
Removes TSs that are not possible with track Pt.
Definition: Fitter3D.cc:1190
Belle2::TRGCDCCell::layerId
unsigned layerId(void) const
returns layer id.
Definition: Cell.h:211
Belle2::TRGCDC::stereoSegment
const TRGCDCSegment & stereoSegment(unsigned lyrId, unsigned id) const
returns a track segment in stereo layers. (lyrId is stereo #)
Definition: TRGCDC.h:947
Belle2::TRGCDCLUT::getValue
int getValue(unsigned) const
get LUT Values
Definition: LUT.cc:47
Belle2::TRGCDCJSignalData::getPrintedToFile
bool getPrintedToFile() const
Gets the status of m_printedToFile.
Definition: JSignalData.cc:80
Belle2::TRGCDCFitter3D::selectAxialTSs
static void selectAxialTSs(const TRGCDCTrack &aTrack, std::vector< int > &bestTSIndex)
Selects priority TSs when there are multiple candidate TSs for a superlayer.
Definition: Fitter3D.cc:1201
Belle2::TRGDebug::enterStage
static void enterStage(const std::string &stageName)
Declare that you enter new stage.
Definition: Debug.cc:24
Belle2::TRGCDCFitter3D::terminate
void terminate()
Termination.
Definition: Fitter3D.cc:1527
Belle2::TRGCDCFitter3D::saveIoSignals
void saveIoSignals()
Saves all I/O signals for debugging.
Definition: Fitter3D.cc:941
Belle2::TRGCDCHelix::phi0
double phi0(void) const
returns phi0.
Definition: Helix.h:278
Belle2::TRGCDCHelix::a
const CLHEP::HepVector & a(void) const
returns helix parameters.
Definition: Helix.h:313
Belle2::TRGCDCHelix::dr
double dr(void) const
returns dr.
Definition: Helix.h:271
Belle2::TRGCDCFitter3D::print3DInformation
void print3DInformation(int iTrack)
Print's information for debugging 3D.
Definition: Fitter3D.cc:1471
Belle2::TRGCDCSegment::lutPattern
unsigned lutPattern(void) const
hit pattern containing bit for priority position
Definition: Segment.cc:559
Belle2::TRGCDCFitter3D::~TRGCDCFitter3D
virtual ~TRGCDCFitter3D()
Destructor.
Definition: Fitter3D.cc:71
Belle2::TRGCDCFitter3D::findHitStereoSuperlayers
static void findHitStereoSuperlayers(const TRGCDCTrack &aTrack, std::vector< double > &useStSL, bool printError)
Finds which stereo superlayers has TSs. useStSL array indicating hit superlayers.
Definition: Fitter3D.cc:1161
Belle2::TRGCDC::segment
const TRGCDCSegment & segment(unsigned id) const
returns a track segment.
Definition: TRGCDC.h:933
Belle2::TRGCDCSegment::LUT
const TRGCDCLUT * LUT(void) const
returns LUT
Definition: Segment.h:243
Belle2::TRGCDC::getEventTime
int getEventTime(void) const
returns bad hits(finding invalid hits).
Definition: TRGCDC.cc:2741
Belle2::TRGCDCFitter3D::do2DFit
static int do2DFit(TRGCDCTrack &aTrack, const std::map< std::string, bool > &m_mBool_in, const 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:1256
Belle2::ORIGIN
const HepGeom::Point3D< double > ORIGIN
Origin 3D point.
Belle2::TRGCDCFitter3D::name
std::string name(void) const
Gets name of class.
Definition: Fitter3D.cc:1550
Belle2::TRGDebug::leaveStage
static void leaveStage(const std::string &stageName)
Declare that you leave a stage.
Definition: Debug.cc:34
Belle2::TRGCDCJSignalData::entryVhdlCode
void entryVhdlCode()
Function to print entry VHDL code.
Definition: JSignalData.cc:195
Belle2::TRGCDCFitter3D::getConstants
static void getConstants(std::map< std::string, double > &mConstD, std::map< std::string, std::vector< double > > &mConstV, bool isXtSimple=0)
Get constants for firmwareFit.
Definition: Fitter3D.cc:1593
Belle2::TRGCDCFitter3D::version
std::string version(void) const
Gets version of class.
Definition: Fitter3D.cc:1545
Belle2::TRGCDCFitter3D::getMCValues
static void getMCValues(const TRGCDC &m_cdc_in, TRGCDCTrack *aTrack, const 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:960
Belle2::TRGCDCTrackBase::links
const std::vector< TRGCDCLink * > & links(void) const
returns a vector to track segments.
Definition: TrackBase.cc:123
Belle2::TRGCDCSegmentHit::segment
const TRGCDCSegment & segment(void) const
returns a pointer to a track segment.
Definition: SegmentHit.cc:78
Belle2::TRGCDCTrack::pt
virtual double pt(void) const override
returns Pt.
Definition: TRGCDCTrack.h:165
Belle2::TRGCDCFitter3D::saveAllSignals
void saveAllSignals()
Saves all signals for debugging.
Definition: Fitter3D.cc:931
Belle2::TRGCDCJSignalData::signalsVhdlCode
void signalsVhdlCode()
Function to print definition of signal VHDL code.
Definition: JSignalData.cc:178
Belle2::TRGCDCJSignalData::setVhdlOutputFile
void setVhdlOutputFile(const std::string &)
Sets the filename for VHDL output.
Definition: JSignalData.cc:45
Belle2::TRGCDCFitter3D::doitComplex
int doitComplex(std::vector< TRGCDCTrack * > &trackList)
Does track fitting using JSignals.
Definition: Fitter3D.cc:457
Belle2::TRGCDCJSignalData::buffersVhdlCode
void buffersVhdlCode()
Function to print buffer VHDL code.
Definition: JSignalData.cc:150
Belle2::TRGCDCTrackBase::setFitted
void setFitted(bool fitted)
set fit status
Definition: TrackBase.h:234
Belle2::TRGCDCFitter3D::isAxialTrackFull
bool isAxialTrackFull(const TRGCDCTrack &aTrack)
Checks if axial track has 5 TSs. One per each superlayer.
Definition: Fitter3D.cc:1068
Belle2::TRGCDCJSignalData::setPrintVhdl
void setPrintVhdl(bool)
Sets if to print VHDL output.
Definition: JSignalData.cc:50
Belle2::TRGCDCCell::localId
unsigned localId(void) const
returns local id in a layer.
Definition: Cell.h:204
Belle2::TRGCDCJSignal::mapSignalsToValues
static void mapSignalsToValues(std::map< std::string, Belle2::TRGCDCJSignal >const &inMap, std::vector< std::pair< std::string, int > > const &inChoose, std::vector< std::tuple< std::string, double, int, double, double, int > > &outValues)
Choose => [signalName, FpgaEffects(=1)/NoFpgaEffects(=0)] Values => [name, value, bitwidth,...
Definition: JSignal.cc:2224
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17