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