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