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