Belle II Software development
Fitter3DUtility.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#ifndef __EXTERNAL__
9#include "trg/cdc/Fitter3DUtility.h"
10#include "trg/cdc/FpgaUtility.h"
11#include "trg/cdc/JLUT.h"
12#include "trg/cdc/JSignal.h"
13#include "trg/cdc/JSignalData.h"
14#include <cmath>
15#else
16#include "Fitter3DUtility.h"
17#include "JLUT.h"
18#include "JSignal.h"
19#include "JSignalData.h"
20#endif
21#include <utility>
22#include <iostream>
23#include <fstream>
24
25using std::cout;
26using std::endl;
27using std::pair;
28using std::make_pair;
29using std::tuple;
30using std::vector;
31using std::map;
32using std::string;
33using std::to_string;
34using std::get;
35using std::function;
36using std::ofstream;
37using std::ifstream;
38using std::make_tuple;
39
40int Fitter3DUtility::findSign(double* phi2)
41{
42 int mysign;
43 double sign_phi[2];
44
45 if ((phi2[0] - phi2[4]) > M_PI || (phi2[0] - phi2[4]) < -M_PI) {
46 if (phi2[0] > M_PI) {sign_phi[0] = phi2[0] - 2 * M_PI;}
47 else {sign_phi[0] = phi2[0];}
48 if (phi2[4] > M_PI) {sign_phi[1] = phi2[4] - 2 * M_PI;}
49 else {sign_phi[1] = phi2[4];}
50 } else {
51 sign_phi[0] = phi2[0];
52 sign_phi[1] = phi2[4];
53 }
54 if ((sign_phi[1] - sign_phi[0]) > 0) {mysign = 0;}
55 else {mysign = 1;}
56
57 return mysign;
58}
59
60void Fitter3DUtility::rPhiFit(double* rr, double* phi2, double* phierror, double& rho, double& myphi0)
61{
62
63 cout << "Fitter3DUtility::rPhiFit() will be deprecated. Please use Fitter3DUtility::rPhiFitter(). phierror was changed to inv phi error."
64 << endl;
65 double invphierror[5];
66 for (unsigned i = 0; i < 5; i++) {
67 invphierror[i] = 1 / phierror[i];
68 }
69 rPhiFitter(rr, phi2, invphierror, rho, myphi0);
70}
71
72void Fitter3DUtility::rPhiFitter(double* rr, double* phi2, double* invphierror, double& rho, double& myphi0)
73{
74
75 double chi2;
76 rPhiFitter(rr, phi2, invphierror, rho, myphi0, chi2);
77
78}
79
80void Fitter3DUtility::rPhiFitter(double* rr, double* phi2, double* invphierror, double& rho, double& myphi0, double& chi2)
81{
82
83 // Print input values
84 //for(unsigned iSL=0; iSL<5; iSL++){
85 // cout<<"SL["<<iSL<<"] rr: "<<rr[iSL]<<" phi: "<<phi2[iSL]<<" phiError: "<<phierror[iSL]<<endl;
86 //}
87
88 double A, B, C, D, E, hcx, hcy;
89 double invFiterror[5];
90 //double G;
91 //Calculate fit error
92 for (unsigned i = 0; i < 5; i++) {
93 // Sometimes SL8 and SL4 will not be hit. So cannot use below calculation.
94 //invFiterror[i]=1/sqrt((rr[4]*rr[4]-2*rr[4]*rr[2]*cos(phi2[4]-phi2[2])+rr[2]*rr[2])/(sin(phi2[4]-phi2[2])*sin(phi2[4]-phi2[2]))-rr[i]*rr[i])*invphierror[i];
95 //invFiterror[i]=invphierror[i];
96 //invFiterror[i]=invphierror[i]*rr[i];
97 invFiterror[i] = invphierror[i] / rr[i];
98 }
99
100 //r-phi fitter(2D Fitter) ->calculate pt and radius of track-> input for 3D fitter.
101 A = 0, B = 0, C = 0, D = 0, E = 0;
102 //G=0;
103 for (unsigned i = 0; i < 5; i++) {
104 A += cos(phi2[i]) * cos(phi2[i]) * (invFiterror[i] * invFiterror[i]);
105 B += sin(phi2[i]) * sin(phi2[i]) * (invFiterror[i] * invFiterror[i]);
106 C += cos(phi2[i]) * sin(phi2[i]) * (invFiterror[i] * invFiterror[i]);
107 D += rr[i] * cos(phi2[i]) * (invFiterror[i] * invFiterror[i]);
108 E += rr[i] * sin(phi2[i]) * (invFiterror[i] * invFiterror[i]);
109 //G+=rr[i]*rr[i]/(fiterror[i]*fiterror[i]);
110 }
111 hcx = D * B - E * C; //helix center x
112 hcx /= 2 * (A * B - C * C);
113 hcy = E * A - D * C; //helix center y
114 hcy /= 2 * (A * B - C * C);
115 rho = sqrt(hcx * hcx + hcy * hcy); //radius of helix
116 myphi0 = atan2(hcy, hcx);
117 if (myphi0 < 0) myphi0 += 2 * M_PI;
118 //myphi0=atan(hcy/hcx);
119 //if(hcx<0 && hcy>0) myphi0 += M_PI;
120 //if(hcx<0 && hcy<0) myphi0 += M_PI;
121 //if(hcx>0 && hcy<0) myphi0 += M_PI*2.0;
122
124 // double pchi2 = -2*hcx*D-2*hcy*E+G;
125 // pchi2/=3;
126
127 // Count number of total TS hits.
128 int nTSHits = 0;
129 for (unsigned iAx = 0; iAx < 5; iAx++) {
130 if (invphierror[iAx] != 0) nTSHits++;
131 }
132 // Another way to calculate chi2
133 chi2 = 0.;
134 for (unsigned i = 0; i < 5; i++) {
135 chi2 += (2 * (hcx * cos(phi2[i]) + hcy * sin(phi2[i])) - rr[i]) * (2 * (hcx * cos(phi2[i]) + hcy * sin(phi2[i])) - rr[i]) /
136 (invFiterror[i] * invFiterror[i]);
137 }
138 chi2 /= nTSHits - 2;
139
140}
141
142void Fitter3DUtility::chargeFinder(double* nTSs, double* tsIds, double* tsHitmap, double phi0, double inCharge,
143 double& foundCharge)
144{
145 vector<double> tsPhi(5);
146 vector<double> dPhi(5);
147 vector<double> dPhi_c(5);
148 vector<double> charge(5);
149 for (unsigned iAx = 0; iAx < 5; iAx++) {
150 //cout<<"iAx:"<<iAx<<" nTSs:"<<nTSs[iAx]<<" tsId:"<<tsIds[iAx]<<" hitmap:"<<tsHitmap[iAx]<<" phi0:"<<phi0<<" inCharge:"<<inCharge<<endl;
151 // Change tsId to rad unit.
152 tsPhi[iAx] = tsIds[iAx] * 2 * M_PI / nTSs[iAx];
153 dPhi[iAx] = tsPhi[iAx] - phi0;
154 //cout<<"iAx:"<<iAx<<" dPhi:"<<dPhi[iAx]<<endl;
155 // Constrain dPhi to [-pi, pi]
156 if (dPhi[iAx] < -M_PI) dPhi_c[iAx] = dPhi[iAx] + 2 * M_PI;
157 else if (dPhi[iAx] > M_PI) dPhi_c[iAx] = dPhi[iAx] - 2 * M_PI;
158 else dPhi_c[iAx] = dPhi[iAx];
159 //cout<<"iAx:"<<iAx<<" dPhi_c:"<<dPhi_c[iAx]<<endl;
160 // Calculate charge
161 if (tsHitmap[iAx] == 0) charge[iAx] = 0;
162 else if (dPhi_c[iAx] == 0) charge[iAx] = 0;
163 else if (dPhi_c[iAx] > 0) charge[iAx] = 1;
164 else charge[iAx] = -1;
165 //cout<<"iAx:"<<iAx<<" charge:"<<charge[iAx]<<endl;
166 }
167 // Sum up charge
168 double sumCharge = charge[0] + charge[1] + charge[2] + charge[3] + charge[4];
169 //cout<<"sumCharge:"<<sumCharge<<endl;
170 // Calculate foundCharge
171 if (sumCharge == 0) foundCharge = inCharge;
172 else if (sumCharge > 0) foundCharge = 1;
173 else foundCharge = -1;
174 //cout<<"foundCharge:"<<foundCharge<<endl;
175 //cout<<"inCharge:"<<inCharge<<endl;
176}
177
178
179void Fitter3DUtility::rPhiFit2(double* rr, double* phi2, double* phierror, double& rho, double& myphi0, int nTS)
180{
181
182 double A, B, C, D, E, hcx, hcy;
183 //double G;
184 double fiterror[5];
185 //Calculate fit error
186 for (int i = 0; i < nTS; i++) {
187 //fiterror[i]=sqrt((rr[4]*rr[4]-2*rr[4]*rr[2]*cos(phi2[4]-phi2[2])+rr[2]*rr[2])/(sin(phi2[4]-phi2[2])*sin(phi2[4]-phi2[2]))-rr[i]*rr[i])*phierror[i];
188 fiterror[i] = 1 + 0 * phierror[i];
189 }
190
191 //r-phi fitter(2D Fitter) ->calculate pt and radius of track-> input for 3D fitter.
192 A = 0, B = 0, C = 0, D = 0, E = 0;
193 //G=0;
194 for (int i = 0; i < nTS; i++) {
195 A += cos(phi2[i]) * cos(phi2[i]) / (fiterror[i] * fiterror[i]);
196 B += sin(phi2[i]) * sin(phi2[i]) / (fiterror[i] * fiterror[i]);
197 C += cos(phi2[i]) * sin(phi2[i]) / (fiterror[i] * fiterror[i]);
198 D += rr[i] * cos(phi2[i]) / (fiterror[i] * fiterror[i]);
199 E += rr[i] * sin(phi2[i]) / (fiterror[i] * fiterror[i]);
200 //G+=rr[i]*rr[i]/(fiterror[i]*fiterror[i]);
201 }
202 hcx = D * B - E * C; //helix center x
203 hcx /= 2 * (A * B - C * C);
204 hcy = E * A - D * C; //helix center y
205 hcy /= 2 * (A * B - C * C);
206 rho = sqrt(hcx * hcx + hcy * hcy); //radius of helix
207 myphi0 = atan2(hcy, hcx);
208 if (myphi0 < 0) myphi0 += 2 * M_PI;
209 //myphi0=atan(hcy/hcx);
210 //if(hcx<0 && hcy>0) myphi0 += M_PI;
211 //if(hcx<0 && hcy<0) myphi0 += M_PI;
212 //if(hcx>0 && hcy<0) myphi0 += M_PI*2.0;
213
215 // double pchi2 = -2*hcx*D-2*hcy*E+G;
216 // pchi2/=nTS-2;
217
219 //double pchi3 = 0.; //iw
220 //for(int i=0;i<nTS;i++){
221 // pchi3+=(2*(hcx*cos(phi2[i])+hcy*sin(phi2[i]))-rr[i])*(2*(hcx*cos(phi2[i])+hcy*sin(phi2[i]))-rr[i])/(fiterror[i]*fiterror[i]);
222 //}
223 //pchi3/=3;
224}
225
226unsigned Fitter3DUtility::toUnsignedTdc(int tdc, int nBits)
227{
228 bool notDone = 1;
229 int tdcMax = pow(2, nBits) - 1;
230 int resultTdc = tdc;
231 while (notDone) {
232 if (resultTdc > tdcMax) resultTdc -= (tdcMax + 1);
233 else if (resultTdc < 0) resultTdc += (tdcMax + 1);
234 else notDone = 0;
235 }
236 return (unsigned) resultTdc;
237}
238
239double Fitter3DUtility::calPhi(double wirePhi, double driftLength, double rr, int lr)
240{
241 //cout<<"rr:"<<rr<<" lr:"<<lr;
242 double result = wirePhi;
243 //cout<<" wirePhi:"<<result;
244 // time is in 2ns rms clock.
245 double t_dPhi = driftLength;
246 //cout<<" driftLength:"<<t_dPhi;
247 // Change to radian
248 // rr is cm scale.
249 t_dPhi = atan(t_dPhi / rr);
250 //cout<<" driftPhi:"<<t_dPhi;
251 // Use LR to add dPhi
252 if (lr == 1) result -= t_dPhi;
253 else if (lr == 2) result += t_dPhi;
254 //cout<<" phi:"<<result<<endl;
255 return result;
256}
257
258double Fitter3DUtility::calPhi(double wirePhi, double driftTime, double eventTime, double rr, int lr)
259{
260 // time is in 2ns rms clock.
261 double t_driftLength = (driftTime - eventTime) * 1000 / 1017.774 * 2 * 40 / 1000 / 10;
262 return calPhi(wirePhi, t_driftLength, rr, lr);
263}
264
265double Fitter3DUtility::calPhi(int localId, int nWires, double driftTime, double eventTime, double rr, int lr)
266{
267 double wirePhi = (double)localId / nWires * 4 * M_PI;
268 return Fitter3DUtility::calPhi(wirePhi, driftTime, eventTime, rr, lr);
269}
270
271void Fitter3DUtility::calPhi(std::map<std::string, double> const& mConstD,
272 std::map<std::string, std::vector<double> > const& mConstV, std::map<std::string, Belle2::TRGCDCJSignal>& mSignalStorage,
273 std::map<std::string, Belle2::TRGCDCJLUT* >& mLutStorage)
274{
275 Belle2::TRGCDCJSignalData* commonData = mSignalStorage.begin()->second.getCommonData();
276 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<tsId_"<<iSt<<">>>"<<endl; mSignalStorage["tsId_"+to_string(iSt)].dump();}
277 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<tdc_"<<iSt<<">>>"<<endl; mSignalStorage["tdc_"+to_string(iSt)].dump();}
278 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<lr_"<<iSt<<">>>"<<endl; mSignalStorage["lr_"+to_string(iSt)].dump();}
279 // Make phiFactor constants
280 if (mSignalStorage.find("phiFactor_0") == mSignalStorage.end()) {
281 for (unsigned iSt = 0; iSt < 4; iSt++) {
282 int nShiftBits = int(log(pow(2, 24) / 2 / mConstD.at("M_PI") * mConstV.at("nTSs")[2 * iSt + 1] *
283 mSignalStorage["phi0"].getToReal()) / log(2));
284 string t_name;
285 t_name = "phiFactor_" + to_string(iSt);
286 mSignalStorage[t_name] = Belle2::TRGCDCJSignal(2 * mConstD.at("M_PI") / mConstV.at("nTSs")[2 * iSt + 1],
287 mSignalStorage["phi0"].getToReal() / pow(2, nShiftBits), commonData);
288 }
289 }
290 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<phiFactor_"<<iSt<<">>>"<<endl; mSignalStorage["phiFactor_"+to_string(iSt)].dump();}
291 // wirePhi <= tsId * phiFactor
292 for (unsigned iSt = 0; iSt < 4;
293 iSt++) mSignalStorage["wirePhi_" + to_string(iSt)] <= mSignalStorage["tsId_" + to_string(iSt)] * mSignalStorage["phiFactor_" +
294 to_string(iSt)];
295 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<wirePhi_"<<iSt<<">>>"<<endl; mSignalStorage["wirePhi_"+to_string(iSt)].dump();}
296
297 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<tdc_"<<iSt<<">>>"<<endl; mSignalStorage["tdc_"+to_string(iSt)].dump();}
298 //cout<<"<<<eventTime>>>"<<endl; mSignalStorage["eventTime"].dump();
299
300 // Calculate driftTimeRaw (= tdc - eventTime)
301 for (unsigned iSt = 0; iSt < 4;
302 iSt++) mSignalStorage["driftTimeRaw_" + to_string(iSt)] <= mSignalStorage["tdc_" + to_string(iSt)] - mSignalStorage["eventTime"];
303
304 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<driftTimeRaw_"<<iSt<<">>>"<<endl; mSignalStorage["driftTimeRaw_"+to_string(iSt)].dump();}
305
306 // Confine drift time
307 // Constants for driftTime
308 if (mSignalStorage.find("maxDriftTimeLow") == mSignalStorage.end()) {
309 int maxDriftTime = 287;
310 mSignalStorage["maxDriftTimeLow"] = Belle2::TRGCDCJSignal(-512 + maxDriftTime + 1, mSignalStorage["eventTime"].getToReal(),
311 commonData);
312 mSignalStorage["maxDriftTimeHigh"] = Belle2::TRGCDCJSignal(maxDriftTime, mSignalStorage["eventTime"].getToReal(), commonData);
313 mSignalStorage["maxDriftTimeMax"] = Belle2::TRGCDCJSignal(512, mSignalStorage["eventTime"].getToReal(), commonData);
314 }
315 for (unsigned iSt = 0; iSt < 4; iSt++) {
316 string t_in1Name = "driftTimeRaw_" + to_string(iSt);
317 string t_valueName = "driftTime_c_" + to_string(iSt);
318 // Create data for ifElse
319 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
320 // Compare (eventTimeValid == 0)
321 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["eventTimeValid"], "=", Belle2::TRGCDCJSignal(0,
322 mSignalStorage["eventTimeValid"].getToReal(), commonData));
323 // Assignments (driftTime = 0)
324 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
325 make_pair(&mSignalStorage[t_valueName], Belle2::TRGCDCJSignal(0, mSignalStorage["eventTime"].getToReal(), commonData))
326 };
327 // Push to data.
328 t_data.push_back(make_pair(t_compare, t_assigns));
329 // Compare (driftTime < maxDriftTimeLow)
330 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], "<", mSignalStorage["maxDriftTimeLow"]);
331 // Assignments
332 t_assigns = {
333 make_pair(&mSignalStorage[t_valueName], (mSignalStorage[t_in1Name] + mSignalStorage["maxDriftTimeMax"]).limit(mSignalStorage["maxDriftTimeLow"], mSignalStorage["maxDriftTimeHigh"]))
334 };
335 // Push to data.
336 t_data.push_back(make_pair(t_compare, t_assigns));
337 // Compare (driftTime < maxDriftTimeHigh)
338 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], "<=", mSignalStorage["maxDriftTimeHigh"]);
339 // Assignments
340 t_assigns = {
341 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_in1Name].limit(mSignalStorage["maxDriftTimeLow"], mSignalStorage["maxDriftTimeHigh"]))
342 };
343 // Push to data.
344 t_data.push_back(make_pair(t_compare, t_assigns));
345 // Compare (else)
346 t_compare = Belle2::TRGCDCJSignal();
347 // Assignments
348 t_assigns = {
349 make_pair(&mSignalStorage[t_valueName], (mSignalStorage[t_in1Name] - mSignalStorage["maxDriftTimeMax"]).limit(mSignalStorage["maxDriftTimeLow"], mSignalStorage["maxDriftTimeHigh"]))
350 };
351 // Push to data.
352 t_data.push_back(make_pair(t_compare, t_assigns));
353 // Process ifElse data.
354 Belle2::TRGCDCJSignal::ifElse(t_data);
355 }
356 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<driftTime_c_"<<iSt<<">>>"<<endl; mSignalStorage["driftTime_c_"+to_string(iSt)].dump();}
357
359 //for(unsigned iSt=0; iSt<4; iSt++) {
360 // if(mSignalStorage["driftTimeRaw_"+to_string(iSt)].getInt() > maxDriftTime) cout<<"driftTimeRaw_"<<iSt<<":"<<mSignalStorage["driftTimeRaw_"+to_string(iSt)].getInt()<<" is larger than maxDriftTime:"<<maxDriftTime<<" changing to "<<mSignalStorage["driftTimeRaw_"+to_string(iSt)].getInt()-512<<endl;
361 // if(mSignalStorage["driftTimeRaw_"+to_string(iSt)].getInt() < -512+maxDriftTime+1) cout<<"driftTimeRaw_"<<iSt<<":"<<mSignalStorage["driftTimeRaw_"+to_string(iSt)].getInt()<<" is smaller than -512+maxDriftTime+1:"<<-512+maxDriftTime+1<<" changed to "<<mSignalStorage["driftTimeRaw_"+to_string(iSt)].getInt()+512<<endl;
362 //}
363
364 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<invErrorRho_"<<iSt<<">>>"<<endl; mSignalStorage["invErrorRho_"+to_string(iSt)].dump();}
365 // Remove negative driftTime
366 // Constants for driftTime
367 if (mSignalStorage.find("minDriftTime") == mSignalStorage.end()) {
368 int minDriftTime = -8;
369 mSignalStorage["minDriftTime"] = Belle2::TRGCDCJSignal(minDriftTime, mSignalStorage["eventTime"].getToReal(), commonData);
370 mSignalStorage["minDriftTimeLow"] = Belle2::TRGCDCJSignal(0, mSignalStorage["eventTime"].getToReal(), commonData);
371 }
372 for (unsigned iSt = 0; iSt < 4; iSt++) {
373 string t_in1Name = "driftTime_c_" + to_string(iSt);
374 string t_valueName = "driftTime_" + to_string(iSt);
375 string t_in2Name = "invErrorRho_" + to_string(iSt);
376 string t_value2Name = "iZError2_" + to_string(iSt);
377 string t_noneErrorName = "iZNoneError2_" + to_string(iSt);
378 // Create data for ifElse
379 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
380 // Compare (driftTime_c[i]>0)
381 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">=", mSignalStorage["minDriftTimeLow"]);
382 // Assignments
383 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
384 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_in1Name].limit(mSignalStorage["minDriftTimeLow"], mSignalStorage["maxDriftTimeHigh"])),
385 make_pair(&mSignalStorage[t_value2Name], mSignalStorage[t_in2Name])
386 };
387 // Push to data.
388 t_data.push_back(make_pair(t_compare, t_assigns));
389 // Compare (driftTime_c[i]>=minDriftTime)
390 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">=", mSignalStorage["minDriftTime"]);
391 // Assignments
392 t_assigns = {
393 make_pair(&mSignalStorage[t_valueName], mSignalStorage["minDriftTimeLow"]),
394 make_pair(&mSignalStorage[t_value2Name], mSignalStorage[t_in2Name])
395 };
396 // Push to data.
397 t_data.push_back(make_pair(t_compare, t_assigns));
398 // Compare (else)
399 t_compare = Belle2::TRGCDCJSignal();
400 // Assignments
401 t_assigns = {
402 make_pair(&mSignalStorage[t_valueName], mSignalStorage["minDriftTimeLow"]),
403 make_pair(&mSignalStorage[t_value2Name], mSignalStorage[t_noneErrorName])
404 };
405 // Push to data.
406 t_data.push_back(make_pair(t_compare, t_assigns));
407 // Process ifElse data.
408 Belle2::TRGCDCJSignal::ifElse(t_data);
409 }
410 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<driftTime_"<<iSt<<">>>"<<endl; mSignalStorage["driftTime_"+to_string(iSt)].dump();}
411 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<iZError2_"<<iSt<<">>>"<<endl; mSignalStorage["iZError2_"+to_string(iSt)].dump();}
412
414 //for(unsigned iSt=0; iSt<4; iSt++) {
415 // if(mSignalStorage["driftTime_c_"+to_string(iSt)].getInt() < minDriftTime) cout<<"driftTime_c_"<<iSt<<":"<<mSignalStorage["driftTime_c_"+to_string(iSt)].getInt()<<" is smaller than minDriftTime:"<<minDriftTime<<endl;
416 // else if(mSignalStorage["driftTime_c_"+to_string(iSt)].getInt() < 0) cout<<"driftTime_c_"<<iSt<<":"<<mSignalStorage["driftTime_c_"+to_string(iSt)].getInt()<<" is smaller than 0."<<endl;
417 //}
418
419 // Calculate minInvValue, maxInvValue for LUTs
420 if (!mSignalStorage.count("invDriftPhiMin")) { //nothing found
421 mSignalStorage["invDriftPhiMin"] = Belle2::TRGCDCJSignal(0, mSignalStorage["driftTime_0"].getToReal(), commonData);
422 mSignalStorage["invDriftPhiMax"] = Belle2::TRGCDCJSignal(pow(2, mConstD.at("tdcBitSize")) - 1,
423 mSignalStorage["driftTime_0"].getToReal(), commonData);
424 }
425 //cout<<"<<<invDriftPhiMin>>>"<<endl; mSignalStorage["invDriftPhiMin"].dump();
426 //cout<<"<<<invDriftPhiMax>>>"<<endl; mSignalStorage["invDriftPhiMax"].dump();
427 // Generate LUT(driftPhi[i]=arctan(driftLengthFunction(driftTime))/r[i])
428 for (unsigned iSt = 0; iSt < 4; iSt++) {
429 string t_valueName = "driftPhi_" + to_string(iSt);
430 string t_inName = "driftTime_" + to_string(iSt);
431 if (!mLutStorage.count(t_valueName)) { //if nothing found
432 mLutStorage[t_valueName] = new Belle2::TRGCDCJLUT(t_valueName);
433 // Lambda can not copy maps.
434 double t_parameter = mConstV.at("rr3D")[iSt];
435 mLutStorage[t_valueName]->setFloatFunction(
436 //[=](double aValue) -> double{return mConstV.at("driftLengthTableSL"+to_string(2*iSt+1))[aValue]/t_parameter;},
437 [ = ](double aValue) -> double{return atan(mConstV.at("driftLengthTableSL" + to_string(2 * iSt + 1))[aValue] / t_parameter);},
438 mSignalStorage[t_inName],
439 mSignalStorage["invDriftPhiMin"], mSignalStorage["invDriftPhiMax"], mSignalStorage["phi0"].getToReal(),
440 (int)mConstD.at("driftPhiLUTInBitSize"), (int)mConstD.at("driftPhiLUTOutBitSize"));
441 //mLutStorage[t_valueName]->makeCOE("./LutData/"+t_valueName+".coe");
442 }
443 }
446 for (unsigned iSt = 0; iSt < 4; iSt++) {
447 //cout<<"driftTime:"<<mSignalStorage["driftTime_"+to_string(iSt)].getActual()<<endl;
448 //cout<<"actual:"<<mLutStorage["driftPhi_"+to_string(iSt)]->getFloatOutput(mSignalStorage["driftTime_"+to_string(iSt)].getActual())<<endl;
449 mLutStorage["driftPhi_" + to_string(iSt)]->operate(mSignalStorage["driftTime_" + to_string(iSt)],
450 mSignalStorage["driftPhi_" + to_string(iSt)]);
451 }
452 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<driftPhi_"<<iSt<<">>>"<<endl; mSignalStorage["driftPhi_"+to_string(iSt)].dump();}
453
454 // Calculate finePhi depending on lr(0:no hit, 1: left, 2: right, 3: not determined)
455 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<lr_"<<iSt<<">>>"<<endl; mSignalStorage["lr_"+to_string(iSt)].dump();}
456 for (unsigned iSt = 0; iSt < 4; iSt++) {
457 string t_in1Name = "lr_" + to_string(iSt);
458 string t_valueName = "finePhi_" + to_string(iSt);
459 // Create data for ifElse
460 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
461 // Compare (lr == 1)
462 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], "=", Belle2::TRGCDCJSignal(1,
463 mSignalStorage[t_in1Name].getToReal(), commonData));
464 // Assignments
465 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
466 make_pair(&mSignalStorage[t_valueName], mSignalStorage["wirePhi_" + to_string(iSt)] - mSignalStorage["driftPhi_" + to_string(iSt)])
467 };
468 // Push to data.
469 t_data.push_back(make_pair(t_compare, t_assigns));
470 // Compare (lr == 2)
471 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], "=", Belle2::TRGCDCJSignal(2,
472 mSignalStorage[t_in1Name].getToReal(), commonData));
473 // Assignments
474 t_assigns = {
475 make_pair(&mSignalStorage[t_valueName], mSignalStorage["wirePhi_" + to_string(iSt)] + mSignalStorage["driftPhi_" + to_string(iSt)])
476 };
477 // Push to data.
478 t_data.push_back(make_pair(t_compare, t_assigns));
479 // Compare (else)
480 t_compare = Belle2::TRGCDCJSignal();
481 // Assignments
482 t_assigns = {
483 make_pair(&mSignalStorage[t_valueName], mSignalStorage["wirePhi_" + to_string(iSt)])
484 };
485 // Push to data.
486 t_data.push_back(make_pair(t_compare, t_assigns));
487 // Process ifElse data.
488 Belle2::TRGCDCJSignal::ifElse(t_data);
489 }
490 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<finePhi_"<<iSt<<">>>"<<endl; mSignalStorage["finePhi_"+to_string(iSt)].dump();}
491
492 // Constrain phi[i] to [-pi, pi]
493 if (mSignalStorage.find("phiMax_0") == mSignalStorage.end()) {
494 for (unsigned iSt = 0; iSt < 4; iSt++) {
495 string t_valueName = "finePhi_" + to_string(iSt);
496 string t_maxName = "phiMax_" + to_string(iSt);
497 string t_minName = "phiMin_" + to_string(iSt);
498 string t_2PiName = "phi2Pi_" + to_string(iSt);
499 mSignalStorage[t_maxName] = Belle2::TRGCDCJSignal(mConstD.at("M_PI"), mSignalStorage[t_valueName].getToReal(), commonData);
500 mSignalStorage[t_minName] = Belle2::TRGCDCJSignal(-mConstD.at("M_PI"), mSignalStorage[t_valueName].getToReal(), commonData);
501 mSignalStorage[t_2PiName] = Belle2::TRGCDCJSignal(2 * mConstD.at("M_PI"), mSignalStorage[t_valueName].getToReal(), commonData);
502 }
503 }
504 for (unsigned iSt = 0; iSt < 4; iSt++) {
505 string t_in1Name = "finePhi_" + to_string(iSt);
506 string t_valueName = "phi_" + to_string(iSt);
507 string t_maxName = "phiMax_" + to_string(iSt);
508 string t_minName = "phiMin_" + to_string(iSt);
509 string t_2PiName = "phi2Pi_" + to_string(iSt);
510 // Create data for ifElse
511 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
512 // Compare (dPhi[i] > dPhiPiMax)
513 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_maxName]);
514 // Assignments (dPhiPi_c[i] = dPhi[i] - dPhiPi2Pi[i])
515 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
516 make_pair(&mSignalStorage[t_valueName], (mSignalStorage[t_in1Name] - mSignalStorage[t_2PiName]).limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
517 };
518 // Push to data.
519 t_data.push_back(make_pair(t_compare, t_assigns));
520 // Compare (dPhi[i] > dPhiPiMin)
521 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_minName]);
522 // Assignments (dPhiPi_c[i] = dPhi[i])
523 t_assigns = {
524 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_in1Name].limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
525 };
526 // Push to data.
527 t_data.push_back(make_pair(t_compare, t_assigns));
528 // Compare (else)
529 t_compare = Belle2::TRGCDCJSignal();
530 // Assignments (dPhiPi_c[i] = dPhi[i] + dPhiPi2Pi[i])
531 t_assigns = {
532 make_pair(&mSignalStorage[t_valueName], (mSignalStorage[t_in1Name] + mSignalStorage[t_2PiName]).limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
533 };
534 // Push to data.
535 t_data.push_back(make_pair(t_compare, t_assigns));
536 // Process ifElse data.
537 Belle2::TRGCDCJSignal::ifElse(t_data);
538 }
539 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<phi_"<<iSt<<">>>"<<endl; mSignalStorage["phi_"+to_string(iSt)].dump();}
540
541}
542
543void Fitter3DUtility::saveStereoXt(std::vector<std::vector<double> >& stXts, string const& filePrefix)
544{
545 for (unsigned iSt = 0; iSt < stXts.size(); ++iSt) {
546 string filename = filePrefix + ".st" + to_string(iSt);
547 ofstream xtFile;
548 cout << "Fitter3DUtility::saveStereoXt() Saving xt file to " << filename << endl;
549 xtFile.open(filename);
550 for (unsigned iTick = 0; iTick < stXts[iSt].size(); ++iTick) {
551 xtFile << stXts[iSt][iTick] << endl;
552 }
553 xtFile.close();
554 }
555}
556
557void Fitter3DUtility::loadStereoXt(string const& filePrefix, int nFiles, std::vector<std::vector<double> >& stXts)
558{
559 stXts = vector<vector<double> > (nFiles);
560 for (unsigned iSt = 0; iSt < stXts.size(); ++iSt) {
561 string filename = filePrefix + ".st" + to_string(iSt);
562 ifstream xtFile;
563 xtFile.open(filename.c_str());
564 if (xtFile.fail()) {
565 cout << "Fitter3DUtility::loadStereoXt() Cannot load " << filename << endl;
566 return;
567 }
568 string line;
569 while (getline(xtFile, line)) {
570 stXts[iSt].push_back(std::stof(line));
571 }
572 xtFile.close();
573 }
574}
575
576double Fitter3DUtility::stereoIdToPhi(vector<double>& stNWires, int iSt, int id)
577{
578 return (double)id / stNWires[iSt] * 4 * M_PI;
579}
580
581// rawStTSs[iSt] = [TS ID, TS LR, TS driftTime]
582// stTSs is 9999 when not valid
583void Fitter3DUtility::calPhiFast(map<string, vector<double> >& stGeometry, vector<vector<double> > const& stXts, int eventTimeValid,
584 int eventTime, vector<vector<int> > const& rawStTSs, vector<double>& stTSs)
585{
586 stTSs = vector<double> (4, 9999);
587 if (eventTimeValid == 0) {
588 for (unsigned iSt = 0; iSt < 4; iSt++) {
589 stTSs[iSt] = stereoIdToPhi(stGeometry["nWires"], iSt, rawStTSs[iSt][0]);
590 }
591 } else {
592 //cout<<"eventTime: "<<eventTime;
593 for (unsigned iSt = 0; iSt < 4; iSt++) {
594 if (rawStTSs[iSt][1] != 0) {
595 // Get drift length from table.
596 int t = rawStTSs[iSt][2] - eventTime;
597 //cout<<" iSt:"<<iSt<<" rt: "<<rawStTSs[iSt][2];
598 if (t < 0) t = 0;
599 if (t > 511) t = 511;
600 double driftLength = stXts[iSt][t];
601 stTSs[iSt] = Fitter3DUtility::calPhi(stereoIdToPhi(stGeometry["nWires"], iSt, rawStTSs[iSt][0]), driftLength,
602 stGeometry["cdcRadius"][iSt], rawStTSs[iSt][1]);
603 }
604 }
605 //cout<<endl;
606 }
607}
608
609double Fitter3DUtility::rotatePhi(double value, double refPhi)
610{
611 double phiMin = -M_PI;
612 double phiMax = M_PI;
613 double result = value - refPhi;
614 bool rangeOK = 0;
615 while (rangeOK == 0) {
616 if (result > phiMax) result -= 2 * M_PI;
617 else if (result < phiMin) result += 2 * M_PI;
618 else rangeOK = 1;
619 }
620 return result;
621}
622
623double Fitter3DUtility::rotatePhi(double value, int refId, int nTSs)
624{
625 double refPhi = (double)refId / nTSs * 2 * M_PI;
626 return rotatePhi(value, refPhi);
627}
628
629int Fitter3DUtility::rotateTsId(int value, int refId, int nTSs)
630{
631 int result = value - refId;
632 bool rangeOk = 0;
633 while (rangeOk == 0) {
634 if (result >= nTSs) result -= nTSs;
635 else if (result < 0) result += nTSs;
636 else rangeOk = 1;
637 }
638 return result;
639}
640
641int Fitter3DUtility::findQuadrant(double value)
642{
643 // Rotate to [-pi,pi] range.
644 double phiMin = -M_PI;
645 double phiMax = M_PI;
646 bool rangeOK = 0;
647 while (rangeOK == 0) {
648 if (value > phiMax) value -= 2 * M_PI;
649 else if (value < phiMin) value += 2 * M_PI;
650 else rangeOK = 1;
651 }
652 // Find quadrant.
653 int result = -1;
654 if (value > M_PI_2) result = 2;
655 else if (value > 0) result = 1;
656 else if (value > -M_PI_2) result = 4;
657 else if (value > -M_PI) result = 3;
658 return result;
659}
660
661
662// rawStTSs[iSt] = [TS ID, TS LR, TS driftTime]
663void Fitter3DUtility::setErrorFast(std::vector<std::vector<int> > const& rawStTSs, int eventTimeValid,
664 std::vector<double>& invZError2)
665{
666 vector<double> driftZError({0.7676, 0.9753, 1.029, 1.372});
667 vector<double> wireZError({0.7676, 0.9753, 1.029, 1.372});
668 vector<double> zError(4, 9999);
669 invZError2 = vector<double> (4, 0);
670 for (unsigned iSt = 0; iSt < 4; ++iSt) {
671 if (rawStTSs[iSt][1] != 0) {
672 // Check LR and eventTime
673 if (rawStTSs[iSt][1] != 3 && eventTimeValid != 0) zError[iSt] = driftZError[iSt];
674 else zError[iSt] = wireZError[iSt];
675 // Get inverse zerror ^ 2
676 invZError2[iSt] = 1 / pow(zError[iSt], 2);
677 }
678 }
679}
680
681void Fitter3DUtility::setError(std::map<std::string, double> const& mConstD,
682 std::map<std::string, std::vector<double> > const& mConstV, std::map<std::string, Belle2::TRGCDCJSignal>& mSignalStorage)
683{
684 Belle2::TRGCDCJSignalData* commonData = mSignalStorage.begin()->second.getCommonData();
685 // Make constants for wireError,driftError,noneError
686 if (mSignalStorage.find("iZWireError2_0") == mSignalStorage.end()) {
687 for (unsigned iSt = 0; iSt < 4; iSt++) {
688 string t_name;
689 t_name = "iZWireError2_" + to_string(iSt);
690 mSignalStorage[t_name] = Belle2::TRGCDCJSignal(mConstD.at("iError2BitSize"), 1 / pow(mConstV.at("wireZError")[iSt], 2), 0,
691 mConstD.at("iError2Max"), -1, commonData);
692 t_name = "iZDriftError2_" + to_string(iSt);
693 mSignalStorage[t_name] = Belle2::TRGCDCJSignal(mConstD.at("iError2BitSize"), 1 / pow(mConstV.at("driftZError")[iSt], 2), 0,
694 mConstD.at("iError2Max"), -1, commonData);
695 t_name = "iZNoneError2_" + to_string(iSt);
696 mSignalStorage[t_name] = Belle2::TRGCDCJSignal(mConstD.at("iError2BitSize"), 0, 0, mConstD.at("iError2Max"), -1, commonData);
697 }
698 }
699 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<iZWireError2_"<<iSt<<">>>"<<endl; mSignalStorage["iZWireError2_"+to_string(iSt)].dump();}
700 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<iZDriftError2_"<<iSt<<">>>"<<endl; mSignalStorage["iZDriftError2_"+to_string(iSt)].dump();}
701 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<iZNoneError2_"<<iSt<<">>>"<<endl; mSignalStorage["iZNoneError2_"+to_string(iSt)].dump();}
702
703 // Set error depending on lr(0:no hit, 1: left, 2: right, 3: not determined)
704 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<lr_"<<iSt<<">>>"<<endl; mSignalStorage["lr_"+to_string(iSt)].dump();}
705 for (unsigned iSt = 0; iSt < 4; iSt++) {
706 string t_in1Name = "lr_" + to_string(iSt);
707 string t_valueName = "invErrorLR_" + to_string(iSt);
708 string t_noneErrorName = "iZNoneError2_" + to_string(iSt);
709 string t_driftErrorName = "iZDriftError2_" + to_string(iSt);
710 string t_wireErrorName = "iZWireError2_" + to_string(iSt);
711 // Create data for ifElse
712 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
713 // Compare (lr == 0)
714 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], "=", Belle2::TRGCDCJSignal(0,
715 mSignalStorage[t_in1Name].getToReal(), commonData));
716 // Assignments (invErrorLR <= iZNoneError2)
717 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
718 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_noneErrorName])
719 };
720 // Push to data.
721 t_data.push_back(make_pair(t_compare, t_assigns));
722 // Compare (lr == 3)
723 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], "=", Belle2::TRGCDCJSignal(3,
724 mSignalStorage[t_in1Name].getToReal(), commonData));
725 // Assignments (invErrorLR <= iZWireError)
726 t_assigns = {
727 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_wireErrorName])
728 };
729 // Push to data.
730 t_data.push_back(make_pair(t_compare, t_assigns));
731 // Compare (else)
732 t_compare = Belle2::TRGCDCJSignal();
733 // Assignments (invErrorLR <= iZDriftError)
734 t_assigns = {
735 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_driftErrorName])
736 };
737 // Push to data.
738 t_data.push_back(make_pair(t_compare, t_assigns));
739 // Process ifElse data.
740 Belle2::TRGCDCJSignal::ifElse(t_data);
741 }
742 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<invErrorLR_"<<iSt<<">>>"<<endl; mSignalStorage["invErrorLR_"+to_string(iSt)].dump();}
743
744 // Make constants for half radius of SL
745 if (mSignalStorage.find("halfRadius_0") == mSignalStorage.end()) {
746 for (unsigned iSt = 0; iSt < 4; iSt++) {
747 string t_name;
748 t_name = "halfRadius_" + to_string(iSt);
749 mSignalStorage[t_name] = Belle2::TRGCDCJSignal(mConstV.at("rr")[iSt * 2 + 1] / 2, mSignalStorage["rho"].getToReal(), commonData);
750 }
751 }
752 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<halfRadius_"<<iSt<<">>>"<<endl; mSignalStorage["halfRadius_"+to_string(iSt)].dump();}
753 // Set error depending on R(helix radius)
754 //cout<<"<<<rho>>>"<<endl; mSignalStorage["rho"].dump();
755 for (unsigned iSt = 0; iSt < 4; iSt++) {
756 string t_compareName = "halfRadius_" + to_string(iSt);
757 //string t_valueName = "iZError2_" + to_string(iSt);
758 string t_valueName = "invErrorRho_" + to_string(iSt);
759 string t_noneErrorName = "iZNoneError2_" + to_string(iSt);
760 string t_in1Name = "invErrorLR_" + to_string(iSt);
761 // Create data for ifElse
762 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
763 // Compare (R < r/2)
764 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["rho"], "<", mSignalStorage[t_compareName]);
765 // Assignments (invError <= 0)
766 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
767 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_noneErrorName])
768 };
769 // Push to data.
770 t_data.push_back(make_pair(t_compare, t_assigns));
771 // Compare (else)
772 t_compare = Belle2::TRGCDCJSignal();
773 // Assignments (invError <= invErrorLR)
774 t_assigns = {
775 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_in1Name])
776 };
777 // Push to data.
778 t_data.push_back(make_pair(t_compare, t_assigns));
779 // Process ifElse data.
780 Belle2::TRGCDCJSignal::ifElse(t_data);
781 }
782 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<invError_"<<iSt<<">>>"<<endl; mSignalStorage["invError_"+to_string(iSt)].dump();}
783}
784
785double Fitter3DUtility::calStAxPhi(int mysign, double anglest, double ztostraw, double rr, double rho, double myphi0)
786{
787 if (false) cout << anglest << ztostraw << endl; // Removes warnings when compileing
788
789 double myphiz, acos_real;
790 //Find phifit-phist
791 double t_rho = rho;
792 if (rr > (2 * rho)) t_rho = rr / 2;
793 acos_real = acos(rr / (2 * t_rho));
794 if (mysign == 1) {
795 myphiz = +acos_real + myphi0;
796 } else {
797 myphiz = -acos_real + myphi0;
798 }
799 if (myphiz > M_PI) myphiz -= 2 * M_PI;
800 if (myphiz < -M_PI) myphiz += 2 * M_PI;
801
802 return myphiz;
803}
804
805double Fitter3DUtility::calDeltaPhi(int mysign, double anglest, double ztostraw, double rr, double phi2, double rho, double myphi0)
806{
807 if (false) cout << anglest << ztostraw << endl; // Removes warnings when compileing
808
809 double myphiz, acos_real;
810 //Find phifit-phist
811 double t_rho = rho;
812 if (rr > (2 * rho)) t_rho = rr / 2;
813 acos_real = acos(rr / (2 * t_rho));
814 if (mysign == 1) {
815 myphiz = -acos_real - myphi0 + phi2;
816 } else {
817 myphiz = +acos_real - myphi0 + phi2;
818 }
819 if (myphiz > M_PI) myphiz -= 2 * M_PI;
820 if (myphiz < -M_PI) myphiz += 2 * M_PI;
821
822 return myphiz;
823}
824
825void Fitter3DUtility::constrainRPerStSl(std::map<std::string, std::vector<double> > const& mConstV,
826 std::map<std::string, Belle2::TRGCDCJSignal>& mSignalStorage)
827{
828 Belle2::TRGCDCJSignalData* commonData = mSignalStorage.begin()->second.getCommonData();
829 // Make constants for min and max values for each super layer.
830 if (mSignalStorage.find("rhoMin_0") == mSignalStorage.end()) {
831 for (unsigned iSt = 0; iSt < 4; iSt++) {
832 // Min value
833 string t_minName = "rhoMin_" + to_string(iSt);
834 double t_actual = (mConstV.find("rr")->second)[2 * iSt + 1] / 2;
835 double t_toReal = mSignalStorage["rho"].getToReal();
836 mSignalStorage[t_minName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
837 // Add one to prevent arccos, arcsin becoming infinity.
838 if (mSignalStorage[t_minName].getRealInt() < t_actual) {
839 t_actual += 1 * t_toReal;
840 mSignalStorage[t_minName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
841 //cout<<"[Warning] "<<t_minName<<" is too small. Adding one unit to prevent arccos output being nan."<<endl;
842 //cout<<t_minName<<".getRealInt():"<<mSignalStorage[t_minName].getRealInt()<<" is new min."<<endl;
843 }
844 // Max value
845 string t_maxName = "rhoMax_" + to_string(iSt);
846 t_actual = mSignalStorage["rho"].getMaxActual();
847 mSignalStorage[t_maxName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
848 }
849 }
850 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<rhoMin_"<<iSt<<">>>"<<endl; mSignalStorage["rhoMin_"+to_string(iSt)].dump();}
851 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<rhoMax_"<<iSt<<">>>"<<endl; mSignalStorage["rhoMax_"+to_string(iSt)].dump();}
852
853 //cout<<"<<<rho>>>"<<endl; mSignalStorage["rho"].dump();
854 // Constrain rho to min and max per layer.
855 for (unsigned iSt = 0; iSt < 4; iSt++) {
856 string t_in1Name = "rho";
857 string t_valueName = "rho_c_" + to_string(iSt);
858 string t_maxName = "rhoMax_" + to_string(iSt);
859 string t_minName = "rhoMin_" + to_string(iSt);
860 // Create data for ifElse
861 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
862 // Compare
863 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_maxName]);
864 // Assignments
865 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
866 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_maxName]),
867 };
868 // Push to data.
869 t_data.push_back(make_pair(t_compare, t_assigns));
870 // Compare
871 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_minName]);
872 // Assignments
873 t_assigns = {
874 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_in1Name].limit(mSignalStorage[t_minName], mSignalStorage[t_maxName]))
875 };
876 // Push to data.
877 t_data.push_back(make_pair(t_compare, t_assigns));
878 // Compare
879 t_compare = Belle2::TRGCDCJSignal();
880 // Assignments
881 t_assigns = {
882 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_minName])
883 };
884 // Push to data.
885 t_data.push_back(make_pair(t_compare, t_assigns));
886 Belle2::TRGCDCJSignal::ifElse(t_data);
887 }
888 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<rho_c_"<<iSt<<">>>"<<endl; mSignalStorage["rho_c_"+to_string(iSt)].dump();}
889
890}
891
892double Fitter3DUtility::calZ(int mysign, double anglest, double ztostraw, double rr, double phi2, double rho, double myphi0)
893{
894 double myphiz = 0, acos_real = 0, dPhiAx = 0;
895 //Find phifit-phist
896 double t_rho = rho;
897 if (rr > (2 * rho)) t_rho = rr / 2;
898 acos_real = acos(rr / (2 * t_rho));
899 if (mysign == 1) {
900 dPhiAx = -acos_real - myphi0;
901 } else {
902 dPhiAx = acos_real - myphi0;
903 }
904 myphiz = dPhiAx + phi2;
905 if (myphiz > M_PI) myphiz -= 2 * M_PI;
906 if (myphiz < -M_PI) myphiz += 2 * M_PI;
907
908 return (ztostraw - rr * 2 * sin(myphiz / 2) / anglest);
909}
910
911void Fitter3DUtility::calZ(std::map<std::string, double> const& mConstD, std::map<std::string, std::vector<double> > const& mConstV,
912 std::map<std::string, Belle2::TRGCDCJSignal>& mSignalStorage, std::map<std::string, Belle2::TRGCDCJLUT* >& mLutStorage)
913{
914 Belle2::TRGCDCJSignalData* commonData = mSignalStorage.begin()->second.getCommonData();
915 // Calculate zz[i]
916 if (mSignalStorage.find("invPhiAxMin_0") == mSignalStorage.end()) {
917 //for(unsigned iSt=0; iSt<4; iSt++) {
918 // string t_invMinName = "invPhiAxMin_" + to_string(iSt);
919 // double t_actual = (mConstV.find("rr")->second)[2*iSt+1]/2;
920 // double t_toReal = mSignalStorage["rho"].getToReal();
921 // mSignalStorage[t_invMinName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
922 // // Add one to prevent arccos becoming infinity.
923 // if(mSignalStorage[t_invMinName].getRealInt() < t_actual) {
924 // t_actual += 1*t_toReal;
925 // mSignalStorage[t_invMinName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
926 // //cout<<"[Warning] "<<t_invMinName<<" is too small. Adding one unit to prevent arccos output being nan."<<endl;
927 // //cout<<t_invMinName<<".getRealInt():"<<mSignalStorage[t_invMinName].getRealInt()<<" is new min."<<endl;
928 // }
929 // string t_invMaxName = "invPhiAxMax_" + to_string(iSt);
930 // t_actual = mSignalStorage["rho"].getMaxActual();
931 // mSignalStorage[t_invMaxName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
932 //}
933 for (unsigned iSt = 0; iSt < 4; iSt++) {
934 string t_invMinName = "invPhiAxMin_" + to_string(iSt);
935 //double t_actual = (mConstV.find("rr")->second)[2*iSt+1]/2;
936 double t_actual = mSignalStorage["rho_c_" + to_string(iSt)].getMinActual();
937 double t_toReal = mSignalStorage["rho_c_" + to_string(iSt)].getToReal();
938 mSignalStorage[t_invMinName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
939 string t_invMaxName = "invPhiAxMax_" + to_string(iSt);
940 t_actual = mSignalStorage["rho_c_" + to_string(iSt)].getMaxActual();
941 mSignalStorage[t_invMaxName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
942 }
943 }
944 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<rho_c_"<<iSt<<">>>"<<endl; mSignalStorage["rho_c_"+to_string(iSt)].dump();}
945 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<invPhiAxMin_"<<iSt<<">>>"<<endl; mSignalStorage["invPhiAxMin_"+to_string(iSt)].dump();}
946 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<invPhiAxMax_"<<iSt<<">>>"<<endl; mSignalStorage["invPhiAxMax_"+to_string(iSt)].dump();}
947 // Generate LUT(phiAx[i]=acos(rr[i]/2/rho)).
948 for (unsigned iSt = 0; iSt < 4; iSt++) {
949 string t_valueName = "phiAx_" + to_string(iSt);
950 string t_minName = "phiAxMin_" + to_string(iSt);
951 string t_maxName = "phiAxMax_" + to_string(iSt);
952 string t_invMinName = "invPhiAxMin_" + to_string(iSt);
953 string t_invMaxName = "invPhiAxMax_" + to_string(iSt);
954 if (!mLutStorage.count(t_valueName)) { //if nothing found
955 mLutStorage[t_valueName] = new Belle2::TRGCDCJLUT(t_valueName);
956 // Lambda can not copy maps.
957 double t_parameter = mConstV.at("rr3D")[iSt];
958 mLutStorage[t_valueName]->setFloatFunction(
959 [ = ](double aValue) -> double{return acos(t_parameter / 2 / aValue);},
960 //mSignalStorage["rho"],
961 mSignalStorage["rho_c_" + to_string(iSt)],
962 mSignalStorage[t_invMinName], mSignalStorage[t_invMaxName], mSignalStorage["phi0"].getToReal(),
963 (int)mConstD.at("acosLUTInBitSize"), (int)mConstD.at("acosLUTOutBitSize"));
964 //mLutStorage[t_valueName]->makeCOE("./LutData/"+t_valueName+".coe");
965 }
966 }
967 // phiAx[i] = acos(rr[i]/2/rho).
968 // Operate using LUT(phiAx[i]).
969 for (unsigned iSt = 0; iSt < 4; iSt++) {
970 // Set output name
971 string t_valueName = "phiAx_" + to_string(iSt);
972 //mLutStorage[t_valueName]->operate(mSignalStorage["rho"], mSignalStorage[t_valueName]);
973 mLutStorage[t_valueName]->operate(mSignalStorage["rho_c_" + to_string(iSt)], mSignalStorage[t_valueName]);
974 }
975 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<phiAx_"<<iSt<<">>>"<<endl; mSignalStorage["phiAx_"+to_string(iSt)].dump();}
976 //cout<<"<<<phi0>>>"<<endl; mSignalStorage["phi0"].dump();
977 //cout<<"<<<charge>>>"<<endl; mSignalStorage["charge"].dump();
978 // dPhiAx[i] = -+phiAx[i] - phi0
979 {
980 // Create data for ifElse
981 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
982 // Compare
983 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(Belle2::TRGCDCJSignal(1, mSignalStorage["charge"].getToReal(),
984 commonData), "=", mSignalStorage["charge"]);
985 // Assignments
986 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
987 make_pair(&mSignalStorage["dPhiAx_0"], -mSignalStorage["phiAx_0"] - mSignalStorage["phi0"]),
988 make_pair(&mSignalStorage["dPhiAx_1"], -mSignalStorage["phiAx_1"] - mSignalStorage["phi0"]),
989 make_pair(&mSignalStorage["dPhiAx_2"], -mSignalStorage["phiAx_2"] - mSignalStorage["phi0"]),
990 make_pair(&mSignalStorage["dPhiAx_3"], -mSignalStorage["phiAx_3"] - mSignalStorage["phi0"])
991 };
992 // Push to data.
993 t_data.push_back(make_pair(t_compare, t_assigns));
994 // Compare
995 t_compare = Belle2::TRGCDCJSignal();
996 // Assignments
997 t_assigns = {
998 make_pair(&mSignalStorage["dPhiAx_0"], mSignalStorage["phiAx_0"] - mSignalStorage["phi0"]),
999 make_pair(&mSignalStorage["dPhiAx_1"], mSignalStorage["phiAx_1"] - mSignalStorage["phi0"]),
1000 make_pair(&mSignalStorage["dPhiAx_2"], mSignalStorage["phiAx_2"] - mSignalStorage["phi0"]),
1001 make_pair(&mSignalStorage["dPhiAx_3"], mSignalStorage["phiAx_3"] - mSignalStorage["phi0"])
1002 };
1003 // Push to data.
1004 t_data.push_back(make_pair(t_compare, t_assigns));
1005 // Process ifElse data.
1006 Belle2::TRGCDCJSignal::ifElse(t_data);
1007 }
1008 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<dPhiAx_"<<iSt<<">>>"<<endl; mSignalStorage["dPhiAx_"+to_string(iSt)].dump();}
1009 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<phi_"<<iSt<<">>>"<<endl; mSignalStorage["phi_"+to_string(iSt)].dump();}
1010 // dPhi[i] = dPhiAx[i] + phi[i]
1011 for (unsigned iSt = 0; iSt < 4; iSt++) {
1012 mSignalStorage["dPhi_" + to_string(iSt)] <= mSignalStorage["dPhiAx_" + to_string(iSt)] + mSignalStorage["phi_" + to_string(iSt)];
1013 }
1014 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<dPhi_"<<iSt<<">>>"<<endl; mSignalStorage["dPhi_"+to_string(iSt)].dump();}
1015
1016 // Rotate dPhi to [-pi, pi]
1017 if (mSignalStorage.find("dPhiPiMax_0") == mSignalStorage.end()) {
1018 for (unsigned iSt = 0; iSt < 4; iSt++) {
1019 string t_valueName = "dPhi_" + to_string(iSt);
1020 string t_maxName = "dPhiPiMax_" + to_string(iSt);
1021 string t_minName = "dPhiPiMin_" + to_string(iSt);
1022 string t_2PiName = "dPhiPi2Pi_" + to_string(iSt);
1023 mSignalStorage[t_maxName] = Belle2::TRGCDCJSignal(mConstD.at("M_PI"), mSignalStorage[t_valueName].getToReal(), commonData);
1024 mSignalStorage[t_minName] = Belle2::TRGCDCJSignal(-mConstD.at("M_PI"), mSignalStorage[t_valueName].getToReal(), commonData);
1025 mSignalStorage[t_2PiName] = Belle2::TRGCDCJSignal(2 * mConstD.at("M_PI"), mSignalStorage[t_valueName].getToReal(), commonData);
1026 }
1027 }
1028 for (unsigned iSt = 0; iSt < 4; iSt++) {
1029 string t_in1Name = "dPhi_" + to_string(iSt);
1030 string t_valueName = "dPhiPi_c_" + to_string(iSt);
1031 string t_maxName = "dPhiPiMax_" + to_string(iSt);
1032 string t_minName = "dPhiPiMin_" + to_string(iSt);
1033 string t_2PiName = "dPhiPi2Pi_" + to_string(iSt);
1034 // Create data for ifElse
1035 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1036 // Compare (dPhi[i] > dPhiPiMax)
1037 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_maxName]);
1038 // Assignments (dPhiPi_c[i] = dPhi[i] - dPhiPi2Pi[i])
1039 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1040 make_pair(&mSignalStorage[t_valueName], (mSignalStorage[t_in1Name] - mSignalStorage[t_2PiName]).limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
1041 };
1042 // Push to data.
1043 t_data.push_back(make_pair(t_compare, t_assigns));
1044 // Compare (dPhi[i] > dPhiPiMin)
1045 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_minName]);
1046 // Assignments (dPhiPi_c[i] = dPhi[i])
1047 t_assigns = {
1048 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_in1Name].limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
1049 };
1050 // Push to data.
1051 t_data.push_back(make_pair(t_compare, t_assigns));
1052 // Compare (else)
1053 t_compare = Belle2::TRGCDCJSignal();
1054 // Assignments (dPhiPi_c[i] = dPhi[i] + dPhiPi2Pi[i])
1055 t_assigns = {
1056 make_pair(&mSignalStorage[t_valueName], (mSignalStorage[t_in1Name] + mSignalStorage[t_2PiName]).limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
1057 };
1058 // Push to data.
1059 t_data.push_back(make_pair(t_compare, t_assigns));
1060 // Process ifElse data.
1061 Belle2::TRGCDCJSignal::ifElse(t_data);
1062 }
1063 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<dPhiPi_c_"<<iSt<<">>>"<<endl; mSignalStorage["dPhiPi_c_"+to_string(iSt)].dump();}
1064
1065 // Calculate dPhiMax[i], dPhiMin[i]=0
1066 if (mSignalStorage.find("dPhiMax_0") == mSignalStorage.end()) {
1067 for (unsigned iSt = 0; iSt < 4; iSt++) {
1068 string t_maxName = "dPhiMax_" + to_string(iSt);
1069 string t_minName = "dPhiMin_" + to_string(iSt);
1070 string t_valueName = "dPhi_" + to_string(iSt);
1071 double t_value = 2 * mConstD.at("M_PI") * fabs(mConstV.at("nShift")[iSt]) / (mConstV.at("nWires")[2 * iSt + 1]);
1072 mSignalStorage[t_maxName] = Belle2::TRGCDCJSignal(t_value, mSignalStorage[t_valueName].getToReal(), commonData);
1073 mSignalStorage[t_minName] = Belle2::TRGCDCJSignal(0, mSignalStorage[t_valueName].getToReal(), commonData);
1074 }
1075 }
1076 // Constrain dPhiPi to [0, 2*pi*#shift/#holes]
1077 for (unsigned iSt = 0; iSt < 4; iSt++) {
1078 string t_maxName = "dPhiMax_" + to_string(iSt);
1079 string t_minName = "dPhiMin_" + to_string(iSt);
1080 string t_valueName = "dPhi_c_" + to_string(iSt);
1081 string t_in1Name = "dPhiPi_c_" + to_string(iSt);
1082 // For SL1, SL5
1083 if (iSt % 2 == 0) {
1084 // Create data for ifElse
1085 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1086 // Compare
1087 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", Belle2::TRGCDCJSignal(0,
1088 mSignalStorage[t_in1Name].getToReal(), commonData));
1089 // Assignments
1090 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1091 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_minName]),
1092 };
1093 // Push to data.
1094 t_data.push_back(make_pair(t_compare, t_assigns));
1095 // Compare
1096 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", -mSignalStorage[t_maxName]);
1097 // Assignments
1098 t_assigns = {
1099 make_pair(&mSignalStorage[t_valueName], Belle2::TRGCDCJSignal::absolute(mSignalStorage[t_in1Name]).limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
1100 };
1101 // Push to data.
1102 t_data.push_back(make_pair(t_compare, t_assigns));
1103 // Compare
1104 t_compare = Belle2::TRGCDCJSignal();
1105 // Assignments
1106 t_assigns = {
1107 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_maxName]),
1108 };
1109 // Push to data.
1110 t_data.push_back(make_pair(t_compare, t_assigns));
1111 // Process ifElse data.
1112 Belle2::TRGCDCJSignal::ifElse(t_data);
1113 // For SL3, SL7
1114 } else {
1115 // Create data for ifElse
1116 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1117 // Compare (dPhi > dPhiMax)
1118 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_maxName]);
1119 // Assignments (dPhi_c = dPhiMax)
1120 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1121 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_maxName]),
1122 };
1123 // Push to data.
1124 t_data.push_back(make_pair(t_compare, t_assigns));
1125 // Compare (dPhi > 0)
1126 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", Belle2::TRGCDCJSignal(0,
1127 mSignalStorage[t_in1Name].getToReal(), commonData));
1128 // Assignments (dPhi_c = dPhi)
1129 t_assigns = {
1130 make_pair(&mSignalStorage[t_valueName], Belle2::TRGCDCJSignal::absolute(mSignalStorage[t_in1Name]).limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
1131 };
1132 // Push to data.
1133 t_data.push_back(make_pair(t_compare, t_assigns));
1134 // Compare (else) => (dPhi > -Pi)
1135 t_compare = Belle2::TRGCDCJSignal();
1136 // Assignments (dPhi_c = dPhiMin)
1137 t_assigns = {
1138 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_minName]),
1139 };
1140 // Push to data.
1141 t_data.push_back(make_pair(t_compare, t_assigns));
1142 // Process ifElse.
1143 Belle2::TRGCDCJSignal::ifElse(t_data);
1144 }
1145 }
1146
1147 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<dPhiMax_"<<iSt<<">>>"<<endl; mSignalStorage["dPhiMax_"+to_string(iSt)].dump();}
1148 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<dPhiMin_"<<iSt<<">>>"<<endl; mSignalStorage["dPhiMin_"+to_string(iSt)].dump();}
1149 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<dPhi_c_"<<iSt<<">>>"<<endl; mSignalStorage["dPhi_c_"+to_string(iSt)].dump();}
1150 // Calculate minInvValue, maxInvValue for LUTs
1151 if (mSignalStorage.find("invZMin_0") == mSignalStorage.end()) {
1152 for (unsigned iSt = 0; iSt < 4; iSt++) {
1153 string t_minName = "invZMin_" + to_string(iSt);
1154 string t_maxName = "invZMax_" + to_string(iSt);
1155 string t_valueName = "dPhi_c_" + to_string(iSt);
1156 unsigned long long t_int = mSignalStorage[t_valueName].getMinInt();
1157 double t_toReal = mSignalStorage[t_valueName].getToReal();
1158 double t_actual = mSignalStorage[t_valueName].getMinActual();
1159 mSignalStorage[t_minName] = Belle2::TRGCDCJSignal(t_int, t_toReal, t_int, t_int, t_actual, t_actual, t_actual, -1, commonData);
1160 t_int = mSignalStorage[t_valueName].getMaxInt();
1161 t_actual = mSignalStorage[t_valueName].getMaxActual();
1162 mSignalStorage[t_maxName] = Belle2::TRGCDCJSignal(t_int, t_toReal, t_int, t_int, t_actual, t_actual, t_actual, -1, commonData);
1163 }
1164 }
1165 // Generate LUT(zz[i] = ztostraw[0]-+rr[0]*2*sin(dPhi_c[i]/2)/angleSt[i], -+ depends on SL)
1166 for (unsigned iSt = 0; iSt < 4; iSt++) {
1167 string t_inputName = "dPhi_c_" + to_string(iSt);
1168 string t_outputName = "zz_" + to_string(iSt);
1169 string t_invMinName = "invZMin_" + to_string(iSt);
1170 string t_invMaxName = "invZMax_" + to_string(iSt);
1171
1172 if (!mLutStorage.count(t_outputName)) { //if nothing found
1173 mLutStorage[t_outputName] = new Belle2::TRGCDCJLUT(t_outputName);
1174 // Lambda can not copy maps.
1175 double t_zToStraw = mConstV.at("zToStraw")[iSt];
1176 double t_rr3D = mConstV.at("rr3D")[iSt];
1177 double t_angleSt = mConstV.at("angleSt")[iSt];
1178
1179 if (iSt % 2 == 0) {
1180 mLutStorage[t_outputName]->setFloatFunction(
1181 [ = ](double aValue) -> double{return t_zToStraw + t_rr3D * 2 * sin(aValue / 2) / t_angleSt;},
1182 mSignalStorage[t_inputName],
1183 mSignalStorage[t_invMinName], mSignalStorage[t_invMaxName], mSignalStorage["rho"].getToReal(),
1184 (int) mConstD.at("zLUTInBitSize"), (int) mConstD.at("zLUTOutBitSize"));
1185 } else {
1186 mLutStorage[t_outputName]->setFloatFunction(
1187 [ = ](double aValue) -> double{return t_zToStraw - t_rr3D * 2 * sin(aValue / 2) / t_angleSt;},
1188 mSignalStorage[t_inputName],
1189 mSignalStorage[t_invMinName], mSignalStorage[t_invMaxName], mSignalStorage["rho"].getToReal(),
1190 (int) mConstD.at("zLUTInBitSize"), (int) mConstD.at("zLUTOutBitSize"));
1191 }
1192 //mLutStorage[t_outputName]->makeCOE("./LutData/"+t_outputName+".coe");
1193 }
1194 }
1195 // zz[i] = ztostraw[0]-+rr[0]*2*sin(dPhi_c[i]/2)/angleSt[i], -+ depends on SL)
1196 // Operate using LUT(zz[i]).
1197 for (unsigned iSt = 0; iSt < 4; iSt++) {
1198 string t_outputName = "zz_" + to_string(iSt);
1199 string t_inputName = "dPhi_c_" + to_string(iSt);
1200 mLutStorage[t_outputName]->operate(mSignalStorage[t_inputName], mSignalStorage[t_outputName]);
1201 }
1202 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<zz_"<<iSt<<">>>"<<endl; mSignalStorage["zz_"+to_string(iSt)].dump();}
1203}
1204double Fitter3DUtility::calS(double rho, double rr)
1205{
1206 double result;
1207 result = rho * 2 * asin(rr / 2 / rho);
1208 return result;
1209}
1210
1211void Fitter3DUtility::calS(std::map<std::string, double> const& mConstD, std::map<std::string, std::vector<double> > const& mConstV,
1212 std::map<std::string, Belle2::TRGCDCJSignal>& mSignalStorage, std::map<std::string, Belle2::TRGCDCJLUT* >& mLutStorage)
1213{
1214 Belle2::TRGCDCJSignalData* commonData = mSignalStorage.begin()->second.getCommonData();
1215 // Calculate minInvValue, maxInvValue for LUTs
1216 if (mSignalStorage.find("invArcSMin_0") == mSignalStorage.end()) {
1217 //for(unsigned iSt=0; iSt<4; iSt++) {
1218 // string t_invMinName = "invArcSMin_" + to_string(iSt);
1219 // double t_actual = mSignalStorage["rho"].getMaxActual();
1220 // double t_toReal = mSignalStorage["rho"].getToReal();
1221 // mSignalStorage[t_invMinName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
1222 // string t_invMaxName = "invArcSMax_" + to_string(iSt);
1223 // t_actual = (mConstV.find("rr")->second)[2*iSt+1]/2;
1224 // mSignalStorage[t_invMaxName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
1225 // // Subtract one to prevent arcsin becoming infinity.
1226 // if(mSignalStorage[t_invMaxName].getRealInt() < t_actual) {
1227 // t_actual += 1*t_toReal;
1228 // mSignalStorage[t_invMaxName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
1229 // //cout<<"[Warning] "<<t_invMaxName<<" is too small. Adding one unit to prevent arcsin output being nan."<<endl;
1230 // //cout<<t_invMaxName<<".getRealInt():"<<mSignalStorage[t_invMaxName].getRealInt()<<" is new max."<<endl;
1231 // }
1232 //}
1233 for (unsigned iSt = 0; iSt < 4; iSt++) {
1234 string t_invMinName = "invArcSMin_" + to_string(iSt);
1235 //double t_actual = (mConstV.find("rr")->second)[2*iSt+1]/2;
1236 double t_actual = mSignalStorage["rho_c_" + to_string(iSt)].getMaxActual();
1237 double t_toReal = mSignalStorage["rho_c_" + to_string(iSt)].getToReal();
1238 mSignalStorage[t_invMinName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
1239 string t_invMaxName = "invArcSMax_" + to_string(iSt);
1240 t_actual = mSignalStorage["rho_c_" + to_string(iSt)].getMinActual();
1241 mSignalStorage[t_invMaxName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
1242 }
1243 }
1244 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<invArcSMin_"<<iSt<<">>>"<<endl; mSignalStorage["invArcSMin_"+to_string(iSt)].dump();}
1245 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<invArcSMax_"<<iSt<<">>>"<<endl; mSignalStorage["invArcSMax_"+to_string(iSt)].dump();}
1246
1247 // Generate LUT(arcS[i]=2*rho*asin(rr[i]/2/rho).
1248 for (unsigned iSt = 0; iSt < 4; iSt++) {
1249 string t_valueName = "arcS_" + to_string(iSt);
1250 string t_invMinName = "invArcSMin_" + to_string(iSt);
1251 string t_invMaxName = "invArcSMax_" + to_string(iSt);
1252 if (!mLutStorage.count(t_valueName)) { //if nothing found
1253 mLutStorage[t_valueName] = new Belle2::TRGCDCJLUT(t_valueName);
1254 // Lambda can not copy maps.
1255 double t_parameter = mConstV.at("rr3D")[iSt];
1256 mLutStorage[t_valueName]->setFloatFunction(
1257 [ = ](double aValue) -> double{return 2 * aValue * asin(t_parameter / 2 / aValue);},
1258 //mSignalStorage["rho"],
1259 mSignalStorage["rho_c_" + to_string(iSt)],
1260 mSignalStorage[t_invMinName], mSignalStorage[t_invMaxName], mSignalStorage["rho"].getToReal(),
1261 (int)mConstD.at("acosLUTInBitSize"), (int)mConstD.at("acosLUTOutBitSize"));
1262 //mLutStorage[t_valueName]->makeCOE("./LutData/"+t_valueName+".coe");
1263 }
1264 }
1265 // arcS[i]=2*rho*asin(rr[i]/2/rho).
1266 // Operate using LUT(arcS[i]).
1267 for (unsigned iSt = 0; iSt < 4; iSt++) {
1268 // Set output name
1269 string t_valueName = "arcS_" + to_string(iSt);
1270 //mLutStorage[t_valueName]->operate(mSignalStorage["rho"], mSignalStorage[t_valueName]);
1271 mLutStorage[t_valueName]->operate(mSignalStorage["rho_c_" + to_string(iSt)], mSignalStorage[t_valueName]);
1272 }
1273 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<arcS_"<<iSt<<">>>"<<endl; mSignalStorage["arcS_"+to_string(iSt)].dump();}
1274}
1275
1276double Fitter3DUtility::calDen(double* arcS, double* zError)
1277{
1278 double t_sum1 = 0;
1279 double t_sumSS = 0;
1280 double t_sumS = 0;
1281 for (unsigned iSt = 0; iSt < 4; iSt++) {
1282 t_sum1 += pow(1 / zError[iSt], 2);
1283 t_sumSS += pow(arcS[iSt] / zError[iSt], 2);
1284 t_sumS += arcS[iSt] / pow(zError[iSt], 2);
1285 }
1286 return t_sum1 * t_sumSS - pow(t_sumS, 2);
1287}
1288
1289double Fitter3DUtility::calDenWithIZError(double* arcS, double* iZError)
1290{
1291 double t_sum1 = 0;
1292 double t_sumSS = 0;
1293 double t_sumS = 0;
1294 for (unsigned iSt = 0; iSt < 4; iSt++) {
1295 t_sum1 += pow(iZError[iSt], 2);
1296 t_sumSS += pow(arcS[iSt] * iZError[iSt], 2);
1297 t_sumS += arcS[iSt] * pow(iZError[iSt], 2);
1298 }
1299 return t_sum1 * t_sumSS - pow(t_sumS, 2);
1300}
1301
1302
1303void Fitter3DUtility::rSFit(double* iezz2, double* arcS, double* zz, double& z0, double& cot, double& zchi2)
1304{
1305
1306 double ss = 0, sx = 0, sxx = 0, cotnum = 0, z0num = 0;
1307 double z0nump1[4], z0nump2[4];
1308 double z0den, iz0den;
1309
1310 for (unsigned i = 0; i < 4; i++) {
1311 ss += iezz2[i];
1312 sx += arcS[i] * iezz2[i];
1313 sxx += arcS[i] * arcS[i] * iezz2[i];
1314 }
1315
1316 for (unsigned i = 0; i < 4; i++) {
1317 cotnum += (ss * arcS[i] - sx) * iezz2[i] * zz[i];
1318 z0nump1[i] = sxx - sx * arcS[i];
1319 z0nump2[i] = z0nump1[i] * iezz2[i] * zz[i];
1320 z0num += z0nump2[i];
1321 }
1322 z0den = (ss * sxx) - (sx * sx);
1323 iz0den = 1. / z0den;
1324 z0num *= iz0den;
1325 cotnum *= iz0den;
1326 z0 = z0num;
1327 cot = cotnum;
1328
1329 // Count number of total TS hits.
1330 int nTSHits = 0;
1331 for (unsigned iSt = 0; iSt < 4; iSt++) {
1332 if (iezz2[iSt] != 0) nTSHits++;
1333 }
1334 // Calculate chi2 of z0
1335 zchi2 = 0.;
1336 for (unsigned i = 0; i < 4; i++) {
1337 zchi2 += (zz[i] - z0 - cot * arcS[i]) * (zz[i] - z0 - cot * arcS[i]) * iezz2[i];
1338 }
1339 zchi2 /= (nTSHits - 2);
1340
1341}
1342
1343void Fitter3DUtility::rSFit2(double* iezz21, double* iezz22, double* arcS, double* zz, int* lutv, double& z0, double& cot,
1344 double& zchi2)
1345{
1346// cout << "rs2" << endl;
1347
1348 double ss = 0, sx = 0, sxx = 0, cotnum = 0, z0num = 0;
1349 double z0nump1[4], z0nump2[4];
1350 double z0den, iz0den;
1351
1352 for (unsigned i = 0; i < 4; i++) {
1353 if (lutv[i] == 2) {
1354 ss += iezz21[i];
1355 sx += arcS[i] * iezz21[i];
1356 sxx += arcS[i] * arcS[i] * iezz21[i];
1357 } else {
1358 ss += iezz22[i];
1359 sx += arcS[i] * iezz22[i];
1360 sxx += arcS[i] * arcS[i] * iezz22[i];
1361 }
1362 }
1363
1364 for (unsigned i = 0; i < 4; i++) {
1365 if (lutv[i] == 2) {
1366 cotnum += (ss * arcS[i] - sx) * iezz21[i] * zz[i];
1367 z0nump1[i] = sxx - sx * arcS[i];
1368 z0nump2[i] = z0nump1[i] * iezz21[i] * zz[i];
1369 z0num += z0nump2[i];
1370 } else {
1371 cotnum += (ss * arcS[i] - sx) * iezz22[i] * zz[i];
1372 z0nump1[i] = sxx - sx * arcS[i];
1373 z0nump2[i] = z0nump1[i] * iezz22[i] * zz[i];
1374 z0num += z0nump2[i];
1375 }
1376 }
1377 z0den = (ss * sxx) - (sx * sx);
1378 iz0den = 1. / z0den;
1379 z0num *= iz0den;
1380 cotnum *= iz0den;
1381 z0 = z0num;
1382 cot = cotnum;
1383
1384 // Calculate chi2 of z0
1385// double zchi2 = 0.;
1386 for (unsigned i = 0; i < 4; i++) {
1387 if (lutv[i] == 2) {
1388 zchi2 += (zz[i] - z0 - cot * arcS[i]) * (zz[i] - z0 - cot * arcS[i]) * iezz21[i];
1389 } else {
1390 zchi2 += (zz[i] - z0 - cot * arcS[i]) * (zz[i] - z0 - cot * arcS[i]) * iezz22[i];
1391 }
1392 }
1393 zchi2 /= (4 - 2);
1394}
1395
1396void Fitter3DUtility::rSFit(std::map<std::string, double> const& mConstD,
1397 std::map<std::string, std::vector<double> > const& mConstV, std::map<std::string, Belle2::TRGCDCJSignal>& mSignalStorage,
1398 std::map<std::string, Belle2::TRGCDCJLUT* >& mLutStorage)
1399{
1400 Belle2::TRGCDCJSignalData* commonData = mSignalStorage.begin()->second.getCommonData();
1401 // sum1_p1_0 = iZError2[0] + iZError2[1]
1402 mSignalStorage["sum1_p1_0"] <= mSignalStorage["iZError2_0"] + mSignalStorage["iZError2_1"];
1403 // sum1_p1_1 = iZError2[2] + iZError2[3]
1404 mSignalStorage["sum1_p1_1"] <= mSignalStorage["iZError2_2"] + mSignalStorage["iZError2_3"];
1405 //for(int iSt=0; iSt<2; iSt++) {cout<<"<<<sum1_p1_"<<iSt<<">>>"<<endl; mSignalStorage["sum1_p1_"+to_string(iSt)].dump();}
1406 // sum1 = sum1_p1[0] + sum1_p1[1]
1407 mSignalStorage["sum1"] <= mSignalStorage["sum1_p1_0"] + mSignalStorage["sum1_p1_1"];
1408 //cout<<"<<<sum1>>>"<<endl; mSignalStorage["sum1"].dump();
1409 // sumS_p1[i] = arcS[i] * iZError2[i]
1410 for (unsigned iSt = 0; iSt < 4; iSt++) {
1411 mSignalStorage["sumS_p1_" + to_string(iSt)] <= mSignalStorage["arcS_" + to_string(iSt)] * mSignalStorage["iZError2_" + to_string(
1412 iSt)];
1413 }
1414 // sumSS_p1[i] = arcS[i] * arcS[i]
1415 for (unsigned iSt = 0; iSt < 4; iSt++) {
1416 mSignalStorage["sumSS_p1_" + to_string(iSt)] <= mSignalStorage["arcS_" + to_string(iSt)] * mSignalStorage["arcS_" + to_string(iSt)];
1417 }
1418 // cotNumS_p1[i] = sum1 * arcS[i]
1419 for (unsigned iSt = 0; iSt < 4; iSt++) {
1420 mSignalStorage["cotNumS_p1_" + to_string(iSt)] <= mSignalStorage["sum1"] * mSignalStorage["arcS_" + to_string(iSt)];
1421 }
1422 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<sumS_p1_"<<iSt<<">>>"<<endl; mSignalStorage["sumS_p1_"+to_string(iSt)].dump();}
1423 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<sumSS_p1_"<<iSt<<">>>"<<endl; mSignalStorage["sumSS_p1_"+to_string(iSt)].dump();}
1424 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<cotNumS_p1_"<<iSt<<">>>"<<endl; mSignalStorage["cotNumS_p1_"+to_string(iSt)].dump();}
1425 // sumS_p2[0] = sumS_p1[0] + sumS_p1[1]
1426 mSignalStorage["sumS_p2_0"] <= mSignalStorage["sumS_p1_0"] + mSignalStorage["sumS_p1_1"];
1427 // sumS_p2[1] = sumS_p1[2] + sumS_p1[3]
1428 mSignalStorage["sumS_p2_1"] <= mSignalStorage["sumS_p1_2"] + mSignalStorage["sumS_p1_3"];
1429 // sumSS_p2[i] = sumSS_p1[i] * iZError2[i]
1430 for (unsigned iSt = 0; iSt < 4; iSt++) {
1431 mSignalStorage["sumSS_p2_" + to_string(iSt)] <= mSignalStorage["sumSS_p1_" + to_string(iSt)] * mSignalStorage["iZError2_" +
1432 to_string(iSt)];
1433 }
1434 //for(int iSt=0; iSt<2; iSt++) {cout<<"<<<sumS_p2_"<<iSt<<">>>"<<endl; mSignalStorage["sumS_p2_"+to_string(iSt)].dump();}
1435 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<sumSS_p2_"<<iSt<<">>>"<<endl; mSignalStorage["sumSS_p2_"+to_string(iSt)].dump();}
1436 // sumS = sumS_p2[0] + sumS_p2[1]
1437 mSignalStorage["sumS"] <= mSignalStorage["sumS_p2_0"] + mSignalStorage["sumS_p2_1"];
1438 // sumSS_p3[0] = sumSS_p2[0] + sumSS_p2[1]
1439 mSignalStorage["sumSS_p3_0"] <= mSignalStorage["sumSS_p2_0"] + mSignalStorage["sumSS_p2_1"];
1440 // sumSS_p3[1] = sumSS_p2[2] + sumSS_p2[3]
1441 mSignalStorage["sumSS_p3_1"] <= mSignalStorage["sumSS_p2_2"] + mSignalStorage["sumSS_p2_3"];
1442 //cout<<"<<<sumS>>>"<<endl; mSignalStorage["sumS"].dump();
1443 //for(int iSt=0; iSt<2; iSt++) {cout<<"<<<sumSS_p3_"<<iSt<<">>>"<<endl; mSignalStorage["sumSS_p3_"+to_string(iSt)].dump();}
1444 //cout<<"===================PIPELINE7============================"<<endl;
1445 // z0NumS_p1[i] = arcS[i] * sumS
1446 for (unsigned iSt = 0; iSt < 4; iSt++) {
1447 mSignalStorage["z0NumS_p1_" + to_string(iSt)] <= mSignalStorage["arcS_" + to_string(iSt)] * mSignalStorage["sumS"];
1448 }
1449 // den_p1 = sumS * sumS
1450 mSignalStorage["den_p1"] <= mSignalStorage["sumS"] * mSignalStorage["sumS"];
1451 // cotNumS[i] = cotNumS_p1[i] - sumS
1452 for (unsigned iSt = 0; iSt < 4; iSt++) {
1453 mSignalStorage["cotNumS_" + to_string(iSt)] <= mSignalStorage["cotNumS_p1_" + to_string(iSt)] - mSignalStorage["sumS"];
1454 }
1455 // sumSS = sumSS_p3[0] + sumSS_p3[1]
1456 mSignalStorage["sumSS"] <= mSignalStorage["sumSS_p3_0"] + mSignalStorage["sumSS_p3_1"];
1457 // zxiZError[i] = zz[i] * iZError2[i]
1458 for (unsigned iSt = 0; iSt < 4; iSt++) {
1459 mSignalStorage["zxiZError_" + to_string(iSt)] <= mSignalStorage["zz_" + to_string(iSt)] * mSignalStorage["iZError2_" + to_string(
1460 iSt)];
1461 }
1462 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<z0NumS_p1_"<<iSt<<">>>"<<endl; mSignalStorage["z0NumS_p1_"+to_string(iSt)].dump();}
1463 //cout<<"<<<den_p1>>>"<<endl; mSignalStorage["den_p1"].dump();
1464 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<cotNumS_"<<iSt<<">>>"<<endl; mSignalStorage["cotNumS_"+to_string(iSt)].dump();}
1465 //cout<<"<<<sumSS>>>"<<endl; mSignalStorage["sumSS"].dump();
1466 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<zxiZError_"<<iSt<<">>>"<<endl; mSignalStorage["zxiZError_"+to_string(iSt)].dump();}
1467 // den_p2 = sum1 * sumSS
1468 mSignalStorage["den_p2"] <= mSignalStorage["sum1"] * mSignalStorage["sumSS"];
1469 // z0NumS[i] =sumSS - z0NumS_p1[i]
1470 for (unsigned iSt = 0; iSt < 4; iSt++) {
1471 mSignalStorage["z0NumS_" + to_string(iSt)] <= mSignalStorage["sumSS"] - mSignalStorage["z0NumS_p1_" + to_string(iSt)];
1472 }
1473 //cout<<"<<<den_p2>>>"<<endl; mSignalStorage["den_p2"].dump();
1474 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<z0NumS_"<<iSt<<">>>"<<endl; mSignalStorage["z0NumS_"+to_string(iSt)].dump();}
1475 // cotNumSZ[i] = cotNumS[i] * zxiZError[i]
1476 for (unsigned iSt = 0; iSt < 4; iSt++) {
1477 mSignalStorage["cotNumSZ_" + to_string(iSt)] <= mSignalStorage["cotNumS_" + to_string(iSt)] * mSignalStorage["zxiZError_" +
1478 to_string(iSt)];
1479 }
1480 // den = den_p2 - den_p1
1481 mSignalStorage["den"] <= mSignalStorage["den_p2"] - mSignalStorage["den_p1"];
1482 // Calculate denMax and denMin.
1483 if (mSignalStorage.find("denMin") == mSignalStorage.end()) {
1484 // Create input for calDen.
1485 double t_rr3D[4], t_wireZError[4];
1486 for (int iSt = 0; iSt < 2; iSt++) {
1487 t_rr3D[iSt] = mConstV.at("rr3D")[iSt];
1488 t_wireZError[iSt] = 1 / mConstV.at("wireZError")[iSt];
1489 }
1490 for (int iSt = 2; iSt < 4; iSt++) {
1491 t_rr3D[iSt] = mConstV.at("rr3D")[iSt];
1492 t_wireZError[iSt] = 0;
1493 }
1494 double t_denMin = Fitter3DUtility::calDenWithIZError(t_rr3D, t_wireZError);
1495 mSignalStorage["denMin"] = Belle2::TRGCDCJSignal(t_denMin, mSignalStorage["den"].getToReal(), commonData);
1496 double t_arcSMax[4], t_driftZError[4];
1497 for (unsigned iSt = 0; iSt < 4; iSt++) {
1498 t_arcSMax[iSt] = mConstD.at("M_PI") * mConstV.at("rr3D")[iSt] / 2;
1499 t_driftZError[iSt] = 1 / mConstV.at("driftZError")[iSt];
1500 }
1501 double t_denMax = Fitter3DUtility::calDenWithIZError(t_arcSMax, t_driftZError);
1502 mSignalStorage["denMax"] = Belle2::TRGCDCJSignal(t_denMax, mSignalStorage["den"].getToReal(), commonData);
1503 }
1504 // Constrain den.
1505 {
1506 // Make ifElse data.
1507 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1508 // Compare
1509 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["den"], ">", mSignalStorage["denMax"]);
1510 // Assignments
1511 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1512 make_pair(&mSignalStorage["den_c"], mSignalStorage["denMax"])
1513 };
1514 // Push to data.
1515 t_data.push_back(make_pair(t_compare, t_assigns));
1516 // Compare
1517 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["den"], ">", mSignalStorage["denMin"]);
1518 // Assignments
1519 t_assigns = {
1520 make_pair(&mSignalStorage["den_c"], mSignalStorage["den"].limit(mSignalStorage["denMin"], mSignalStorage["denMax"]))
1521 };
1522 // Push to data.
1523 t_data.push_back(make_pair(t_compare, t_assigns));
1524 // Compare
1525 t_compare = Belle2::TRGCDCJSignal();
1526 // Assignments
1527 t_assigns = {
1528 make_pair(&mSignalStorage["den_c"], mSignalStorage["denMin"])
1529 };
1530 // Push to data.
1531 t_data.push_back(make_pair(t_compare, t_assigns));
1532 Belle2::TRGCDCJSignal::ifElse(t_data);
1533 }
1534 // z0NumSZ[i] = z0NumS[i]*zxiZError[i]
1535 for (unsigned iSt = 0; iSt < 4; iSt++) {
1536 mSignalStorage["z0NumSZ_" + to_string(iSt)] <= mSignalStorage["z0NumS_" + to_string(iSt)] * mSignalStorage["zxiZError_" + to_string(
1537 iSt)];
1538 }
1539 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<cotNumSZ_"<<iSt<<">>>"<<endl; mSignalStorage["cotNumSZ_"+to_string(iSt)].dump();}
1540 //cout<<"<<<den>>>"<<endl; mSignalStorage["den"].dump();
1541 //cout<<"<<<den_c>>>"<<endl; mSignalStorage["den_c"].dump();
1542 //cout<<"<<<denMin>>>"<<endl; mSignalStorage["denMin"].dump();
1543 //cout<<"<<<denMax>>>"<<endl; mSignalStorage["denMax"].dump();
1544 //cout<<"<<<offsetDen>>>"<<endl; mSignalStorage["offsetDen"].dump();
1545 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<z0NumSZ_"<<iSt<<">>>"<<endl; mSignalStorage["z0NumSZ_"+to_string(iSt)].dump();}
1546 // cot_p1[0] = cotNumSZ[0] + cotNumSZ[1]
1547 mSignalStorage["cot_p1_0"] <= mSignalStorage["cotNumSZ_0"] + mSignalStorage["cotNumSZ_1"];
1548 // cot_p1[1] = cotNumSZ[2] + cotNumSZ[3]
1549 mSignalStorage["cot_p1_1"] <= mSignalStorage["cotNumSZ_2"] + mSignalStorage["cotNumSZ_3"];
1550 // z0_p1_0 = z0NumSZ_0 + z0NumSZ_1
1551 mSignalStorage["z0_p1_0"] <= mSignalStorage["z0NumSZ_0"] + mSignalStorage["z0NumSZ_1"];
1552 // z0_p1_1 = z0NumSZ_2 + z0NumSZ_3
1553 mSignalStorage["z0_p1_1"] <= mSignalStorage["z0NumSZ_2"] + mSignalStorage["z0NumSZ_3"];
1555 //mSignalStorage["iDenMin"] <= Belle2::TRGCDCJSignal(1/t_denMax, pow(1/mSignalStorage["rho"].getToReal()/mSignalStorage["iZError2_0"].getToReal(),2));
1556 //mSignalStorage["iDenMax"] <= Belle2::TRGCDCJSignal(1/t_denMin, pow(1/mSignalStorage["rho"].getToReal()/mSignalStorage["iZError2_0"].getToReal(),2));
1557 // Calculate minInvValue, maxInvValue for LUTs
1558 if (mSignalStorage.find("invIDenMin") == mSignalStorage.end()) {
1559 unsigned long long t_int = mSignalStorage["den_c"].getMaxInt();
1560 double t_toReal = mSignalStorage["den_c"].getToReal();
1561 double t_actual = mSignalStorage["den_c"].getMaxActual();
1562 mSignalStorage["invIDenMin"] = Belle2::TRGCDCJSignal(t_int, t_toReal, t_int, t_int, t_actual, t_actual, t_actual, -1, commonData);
1563 t_int = mSignalStorage["den_c"].getMinInt();
1564 t_actual = mSignalStorage["den_c"].getMinActual();
1565 mSignalStorage["invIDenMax"] = Belle2::TRGCDCJSignal(t_int, t_toReal, t_int, t_int, t_actual, t_actual, t_actual, -1, commonData);
1566 }
1567 //cout<<"<<<invIDenMin>>>"<<endl; mSignalStorage["invIDenMin"].dump();
1568 //cout<<"<<<invIDenMax>>>"<<endl; mSignalStorage["invIDenMax"].dump();
1569 // Generate LUT(iDen = 1/den)
1570 if (!mLutStorage.count("iDen")) {
1571 mLutStorage["iDen"] = new Belle2::TRGCDCJLUT("iDen");
1572 mLutStorage["iDen"]->setFloatFunction(
1573 [ = ](double aValue) -> double{return 1 / aValue;},
1574 mSignalStorage["den_c"],
1575 mSignalStorage["invIDenMin"], mSignalStorage["invIDenMax"],
1576 pow(1 / mSignalStorage["rho"].getToReal() / mSignalStorage["iZError2_0"].getToReal(), 2),
1577 (int)mConstD.at("iDenLUTInBitSize"), (int)mConstD.at("iDenLUTOutBitSize"));
1578 //mLutStorage["iDen"]->makeCOE("./LutData/iDen.coe");
1579 }
1580 // Operate using LUT(iDen = 1/den)
1581 mLutStorage["iDen"]->operate(mSignalStorage["den_c"], mSignalStorage["iDen"]);
1582 //for(int iSt=0; iSt<2; iSt++) {cout<<"<<<cot_p1_"<<iSt<<">>>"<<endl; mSignalStorage["cot_p1_"+to_string(iSt)].dump();}
1583 //for(int iSt=0; iSt<2; iSt++) {cout<<"<<<z0_p1_"<<iSt<<">>>"<<endl; mSignalStorage["z0_p1_"+to_string(iSt)].dump();}
1584 //cout<<"<<<offsetDen_c>>>"<<endl; mSignalStorage["offsetDen_c"].dump();
1585 //cout<<"<<<denOffset_c>>>"<<endl; mSignalStorage["denOffset_c"].dump();
1586 //cout<<"<<<iDen>>>"<<endl; mSignalStorage["iDen"].dump();
1587 // iDen = offsetIDen + iDenOffset (Already done in LUT->operate)
1588 // cot_p2 = cot_p1_0 + cot_p1_1
1589 mSignalStorage["cot_p2"] <= mSignalStorage["cot_p1_0"] + mSignalStorage["cot_p1_1"];
1590 // z0_p2 = z0_p1_0 + z0_p1_1
1591 mSignalStorage["z0_p2"] <= mSignalStorage["z0_p1_0"] + mSignalStorage["z0_p1_1"];
1592 //cout<<"<<<cot_p2>>>"<<endl; mSignalStorage["cot_p2"].dump();
1593 //cout<<"<<<z0_p2>>>"<<endl; mSignalStorage["z0_p2"].dump();
1594 // cot = cot_p2 * iDen
1595 mSignalStorage["cot"] <= mSignalStorage["cot_p2"] * mSignalStorage["iDen"];
1596 // z0 = z0_p2 * iDen
1597 mSignalStorage["z0"] <= mSignalStorage["z0_p2"] * mSignalStorage["iDen"];
1598 if (mConstD.at("JB") == 1) {
1599 cout << "<<<cot>>>" << endl; mSignalStorage["cot"].dump();
1600 cout << "<<<z0>>>" << endl; mSignalStorage["z0"].dump();
1601 }
1602
1603 // Constrain z0 to [-30,30]
1604 /* cppcheck-suppress variableScope */
1605 double z0Min = -30;
1606 /* cppcheck-suppress variableScope */
1607 double z0Max = 30;
1608 // Calculate z0Max and z0Min
1609 if (!mSignalStorage.count("z0Min")) {
1610 mSignalStorage["z0Min"] = Belle2::TRGCDCJSignal(z0Min, mSignalStorage["z0"].getToReal(), commonData);
1611 mSignalStorage["z0Max"] = Belle2::TRGCDCJSignal(z0Max, mSignalStorage["z0"].getToReal(), commonData);
1612 }
1613 {
1614 // Make ifElse data.
1615 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1616 // Compare
1617 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["z0"], ">", mSignalStorage["z0Max"]);
1618 // Assignments
1619 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1620 make_pair(&mSignalStorage["z0_c"], mSignalStorage["z0Max"])
1621 };
1622 // Push to data.
1623 t_data.push_back(make_pair(t_compare, t_assigns));
1624 // Compare
1625 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["z0"], ">", mSignalStorage["z0Min"]);
1626 // Assignments
1627 t_assigns = {
1628 make_pair(&mSignalStorage["z0_c"], mSignalStorage["z0"].limit(mSignalStorage["z0Min"], mSignalStorage["z0Max"]))
1629 };
1630 // Push to data.
1631 t_data.push_back(make_pair(t_compare, t_assigns));
1632 // Compare
1633 t_compare = Belle2::TRGCDCJSignal();
1634 // Assignments
1635 t_assigns = {
1636 make_pair(&mSignalStorage["z0_c"], mSignalStorage["z0Min"])
1637 };
1638 // Push to data.
1639 t_data.push_back(make_pair(t_compare, t_assigns));
1640 Belle2::TRGCDCJSignal::ifElse(t_data);
1641 }
1642 //cout<<"<<<z0Min>>>"<<endl; mSignalStorage["z0Min"].dump();
1643 //cout<<"<<<z0Max>>>"<<endl; mSignalStorage["z0Max"].dump();
1644 //cout<<"<<<z0_c>>>"<<endl; mSignalStorage["z0_c"].dump();
1645 // Constraint cot to [-0.753, 1.428]
1646 double cotMin = cos(127 * mConstD.at("M_PI") / 180) / sin(127 * mConstD.at("M_PI") / 180);
1647 double cotMax = cos(35 * mConstD.at("M_PI") / 180) / sin(35 * mConstD.at("M_PI") / 180);
1648 if (!mSignalStorage.count("cotMin")) {
1649 mSignalStorage["cotMin"] = Belle2::TRGCDCJSignal(cotMin, mSignalStorage["cot"].getToReal(), commonData);
1650 mSignalStorage["cotMax"] = Belle2::TRGCDCJSignal(cotMax, mSignalStorage["cot"].getToReal(), commonData);
1651 }
1652 {
1653 // Make ifElse data.
1654 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1655 // Compare
1656 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["cot"], ">", mSignalStorage["cotMax"]);
1657 // Assignments
1658 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1659 make_pair(&mSignalStorage["cot_c"], mSignalStorage["cotMax"])
1660 };
1661 // Push to data.
1662 t_data.push_back(make_pair(t_compare, t_assigns));
1663 // Compare
1664 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["cot"], ">", mSignalStorage["cotMin"]);
1665 // Assignments
1666 t_assigns = {
1667 make_pair(&mSignalStorage["cot_c"], mSignalStorage["cot"].limit(mSignalStorage["cotMin"], mSignalStorage["cotMax"]))
1668 };
1669 // Push to data.
1670 t_data.push_back(make_pair(t_compare, t_assigns));
1671 // Compare
1672 t_compare = Belle2::TRGCDCJSignal();
1673 // Assignments
1674 t_assigns = {
1675 make_pair(&mSignalStorage["cot_c"], mSignalStorage["cotMin"])
1676 };
1677 // Push to data.
1678 t_data.push_back(make_pair(t_compare, t_assigns));
1679 Belle2::TRGCDCJSignal::ifElse(t_data);
1680 }
1681 //cout<<"<<<cotMin>>>"<<endl; mSignalStorage["cotMin"].dump();
1682 //cout<<"<<<cotMax>>>"<<endl; mSignalStorage["cotMax"].dump();
1683 //cout<<"<<<cot_c>>>"<<endl; mSignalStorage["cot_c"].dump();
1684 // Reduce number of bits for z0 to 11 bits and cot to 11 bits
1685 /* cppcheck-suppress variableScope */
1686 int z0Bits = 11;
1687 /* cppcheck-suppress variableScope */
1688 int cotBits = 11;
1689 {
1690 int z0Shift = mSignalStorage["z0_c"].getBitsize() - z0Bits;
1691 mSignalStorage["z0_r"] <= mSignalStorage["z0_c"].shift(z0Shift, 0);
1692 int cotShift = mSignalStorage["cot_c"].getBitsize() - cotBits;
1693 mSignalStorage["cot_r"] <= mSignalStorage["cot_c"].shift(cotShift, 0);
1694 }
1695 //cout<<"<<<z0_c>>>"<<endl; mSignalStorage["z0_c"].dump();
1696 //cout<<"<<<cot_c>>>"<<endl; mSignalStorage["cot_c"].dump();
1697 //cout<<"<<<z0_r>>>"<<endl; mSignalStorage["z0_r"].dump();
1698 //cout<<"<<<cot_r>>>"<<endl; mSignalStorage["cot_r"].dump();
1699
1700
1701 // fitDistFromZ0[i] = cot*arcS[i]
1702 for (unsigned iSt = 0; iSt < 4; iSt++) {
1703 mSignalStorage["fitDistFromZ0_" + to_string(iSt)] <= mSignalStorage["cot"] * mSignalStorage["arcS_" + to_string(iSt)];
1704 }
1705 // distFromZ0[i] = z[i] - z0
1706 for (unsigned iSt = 0; iSt < 4; iSt++) {
1707 mSignalStorage["distFromZ0_" + to_string(iSt)] <= mSignalStorage["zz_" + to_string(iSt)] - mSignalStorage["z0"];
1708 }
1709 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<fitDistFromZ0_"<<iSt<<">>>"<<endl; mSignalStorage["fitDistFromZ0_"+to_string(iSt)].dump();}
1710 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<distFromZ0_"<<iSt<<">>>"<<endl; mSignalStorage["distFromZ0_"+to_string(iSt)].dump();}
1711 // chiNum[i] = distFromZ0[i] - fitDistFromZ0[i]
1712 for (unsigned iSt = 0; iSt < 4; iSt++) {
1713 mSignalStorage["chiNum_" + to_string(iSt)] <= mSignalStorage["distFromZ0_" + to_string(iSt)] - mSignalStorage["fitDistFromZ0_" +
1714 to_string(iSt)];
1715 }
1716 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<chiNum_"<<iSt<<">>>"<<endl; mSignalStorage["chiNum_"+to_string(iSt)].dump();}
1717 // Calculate chiNumMax[i], chiNumMin[i].
1718 if (mSignalStorage.find("chiNumMax_0") == mSignalStorage.end()) {
1719 for (unsigned iSt = 0; iSt < 4; iSt++) {
1720 string t_maxName = "chiNumMax_" + to_string(iSt);
1721 string t_minName = "chiNumMin_" + to_string(iSt);
1722 string t_valueName = "chiNum_" + to_string(iSt);
1723 double t_maxValue = 2 * mConstV.at("zToOppositeStraw")[iSt];
1724 double t_minValue = 2 * mConstV.at("zToStraw")[iSt];
1725 mSignalStorage[t_maxName] = Belle2::TRGCDCJSignal(t_maxValue, mSignalStorage[t_valueName].getToReal(), commonData);
1726 mSignalStorage[t_minName] = Belle2::TRGCDCJSignal(t_minValue, mSignalStorage[t_valueName].getToReal(), commonData);
1727 }
1728 }
1729 // Constrain chiNum[i] to [-2*z_min[i], 2*z_max[i]]
1730 for (unsigned iSt = 0; iSt < 4; iSt++) {
1731 string t_maxName = "chiNumMax_" + to_string(iSt);
1732 string t_minName = "chiNumMin_" + to_string(iSt);
1733 string t_valueName = "chiNum_c_" + to_string(iSt);
1734 string t_in1Name = "chiNum_" + to_string(iSt);
1735 // Create data for ifElse
1736 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1737 // Compare
1738 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_maxName]);
1739 // Assignments
1740 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1741 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_maxName]),
1742 };
1743 // Push to data.
1744 t_data.push_back(make_pair(t_compare, t_assigns));
1745 // Compare
1746 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_minName]);
1747 // Assignments
1748 t_assigns = {
1749 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_in1Name].limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
1750 };
1751 // Push to data.
1752 t_data.push_back(make_pair(t_compare, t_assigns));
1753 // Compare
1754 t_compare = Belle2::TRGCDCJSignal();
1755 // Assignments
1756 t_assigns = {
1757 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_minName]),
1758 };
1759 // Push to data.
1760 t_data.push_back(make_pair(t_compare, t_assigns));
1761 // Process ifElse data.
1762 Belle2::TRGCDCJSignal::ifElse(t_data);
1763 }
1764 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<chiNumMax_"<<iSt<<">>>"<<endl; mSignalStorage["chiNumMax_"+to_string(iSt)].dump();}
1765 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<chiNumMin_"<<iSt<<">>>"<<endl; mSignalStorage["chiNumMin_"+to_string(iSt)].dump();}
1766 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<chiNum_c_"<<iSt<<">>>"<<endl; mSignalStorage["chiNum_c_"+to_string(iSt)].dump();}
1767 // chi2_p1[i] = chiNum_c[i] * chiNum_c[i]
1768 for (unsigned iSt = 0; iSt < 4; iSt++) {
1769 mSignalStorage["chi2_p1_" + to_string(iSt)] <= mSignalStorage["chiNum_c_" + to_string(iSt)] * mSignalStorage["chiNum_c_" +
1770 to_string(iSt)];
1771 }
1772 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<chi2_p1_"<<iSt<<">>>"<<endl; mSignalStorage["chi2_p1_"+to_string(iSt)].dump();}
1773 // chi2[i] = chi2_p1[i] * iError2[i]
1774 for (unsigned iSt = 0; iSt < 4; iSt++) {
1775 mSignalStorage["chi2_" + to_string(iSt)] <= mSignalStorage["chi2_p1_" + to_string(iSt)] * mSignalStorage["iZError2_" + to_string(
1776 iSt)];
1777 }
1778 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<chi2_"<<iSt<<">>>"<<endl; mSignalStorage["chi2_"+to_string(iSt)].dump();}
1779 // chi2Sum_p1_0 = chi2_0 + chi2_1
1780 mSignalStorage["chi2Sum_p1_0"] <= mSignalStorage["chi2_0"] + mSignalStorage["chi2_1"];
1781 // chi2Sum_p1_1 = chi2_2 + chi2_3
1782 mSignalStorage["chi2Sum_p1_1"] <= mSignalStorage["chi2_2"] + mSignalStorage["chi2_3"];
1783 //cout<<"<<<chi2Sum_p1_0>>>"<<endl; mSignalStorage["chi2Sum_p1_0"].dump();
1784 //cout<<"<<<chi2Sum_p1_1>>>"<<endl; mSignalStorage["chi2Sum_p1_1"].dump();
1785 // chi2Sum = chi2Sum_p1_0 + chi2Sum_p1_1
1786 mSignalStorage["zChi2"] <= (mSignalStorage["chi2Sum_p1_0"] + mSignalStorage["chi2Sum_p1_1"]).shift(1);
1787 if (mConstD.at("JB") == 1) {cout << "<<<zChi2>>>" << endl; mSignalStorage["zChi2"].dump();}
1788
1789 // Constrain zChi2 to [0, 100]
1790 /* cppcheck-suppress variableScope */
1791 double zChi2Min = 0;
1792 /* cppcheck-suppress variableScope */
1793 double zChi2Max = 100;
1794 /* cppcheck-suppress variableScope */
1795 int zChi2Bits = 4;
1796 // Calculate zChi2Max and zChi2Min
1797 if (!mSignalStorage.count("zChi2Min")) {
1798 mSignalStorage["zChi2Min"] = Belle2::TRGCDCJSignal(zChi2Min, mSignalStorage["zChi2"].getToReal(), commonData);
1799 mSignalStorage["zChi2Max"] = Belle2::TRGCDCJSignal(zChi2Max, mSignalStorage["zChi2"].getToReal(), commonData);
1800 }
1801 {
1802 // Make ifElse data.
1803 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1804 // Compare
1805 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["zChi2"], ">", mSignalStorage["zChi2Max"]);
1806 // Assignments
1807 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1808 make_pair(&mSignalStorage["zChi2_c"], mSignalStorage["zChi2Max"])
1809 };
1810 // Push to data.
1811 t_data.push_back(make_pair(t_compare, t_assigns));
1812 // Compare
1813 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["zChi2"], ">", mSignalStorage["zChi2Min"]);
1814 // Assignments
1815 t_assigns = {
1816 make_pair(&mSignalStorage["zChi2_c"], mSignalStorage["zChi2"].limit(mSignalStorage["zChi2Min"], mSignalStorage["zChi2Max"]))
1817 };
1818 // Push to data.
1819 t_data.push_back(make_pair(t_compare, t_assigns));
1820 // Compare
1821 t_compare = Belle2::TRGCDCJSignal();
1822 // Assignments
1823 t_assigns = {
1824 make_pair(&mSignalStorage["zChi2_c"], mSignalStorage["zChi2Min"])
1825 };
1826 // Push to data.
1827 t_data.push_back(make_pair(t_compare, t_assigns));
1828 Belle2::TRGCDCJSignal::ifElse(t_data);
1829 }
1830 //cout<<"<<<zChi2Min>>>"<<endl; mSignalStorage["zChi2Min"].dump();
1831 //cout<<"<<<zChi2Max>>>"<<endl; mSignalStorage["zChi2Max"].dump();
1832 //cout<<"<<<zChi2_c>>>"<<endl; mSignalStorage["zChi2_c"].dump();
1833 {
1834 int zChi2Shift = mSignalStorage["zChi2_c"].getBitsize() - zChi2Bits;
1835 mSignalStorage["zChi2_r"] <= mSignalStorage["zChi2_c"].shift(zChi2Shift, 0);
1836 }
1837 //cout<<"<<<zChi2_c>>>"<<endl; mSignalStorage["zChi2_c"].dump();
1838 //cout<<"<<<zChi2_r>>>"<<endl; mSignalStorage["zChi2_r"].dump();
1839
1840}
1841
1842void Fitter3DUtility::fitter3D(std::map<std::string, std::vector<double> >& stGeometry,
1843 std::vector<std::vector<double> > const& stXts,
1844 int eventTimeValid, int eventTime,
1845 std::vector<std::vector<int> > const& rawStTSs,
1846 int charge, double radius, double phi_c, double& z0, double& cot, double& chi2)
1847{
1848 vector<double> zz;
1849 vector<double> arcS;
1850 vector<double> invZError2;
1851 fitter3D(stGeometry, stXts, eventTimeValid, eventTime, rawStTSs, charge, radius, phi_c, z0, cot, chi2, arcS, zz, invZError2);
1852}
1853
1854void Fitter3DUtility::fitter3D(std::map<std::string, std::vector<double> >& stGeometry,
1855 std::vector<std::vector<double> > const& stXts,
1856 int eventTimeValid, int eventTime,
1857 std::vector<std::vector<int> > const& rawStTSs,
1858 int charge, double radius, double phi_c, double& z0, double& cot, double& chi2, vector<double>& arcS, vector<double>& zz,
1859 vector<double>& invZError2)
1860{
1861 vector<double> stTSs(4);
1862 Fitter3DUtility::calPhiFast(stGeometry, stXts, eventTimeValid, eventTime, rawStTSs, stTSs);
1863 Fitter3DUtility::setErrorFast(rawStTSs, eventTimeValid, invZError2);
1864 zz = vector<double> (4, 0);
1865 arcS = vector<double> (4, 0);
1866 for (unsigned iSt = 0; iSt < 4; iSt++) {
1867 if (rawStTSs[iSt][1] != 0) {
1868 zz[iSt] = Fitter3DUtility::calZ(charge, stGeometry["angleSt"][iSt], stGeometry["zToStraw"][iSt],
1869 stGeometry["cdcRadius"][iSt], stTSs[iSt], radius, phi_c);
1870 arcS[iSt] = Fitter3DUtility::calS(radius, stGeometry["cdcRadius"][iSt]);
1871 }
1872 }
1873 Fitter3DUtility::rSFit(&invZError2[0], &arcS[0], &zz[0], z0, cot, chi2);
1874}
1875
1876void Fitter3DUtility::fitter3DFirm(std::map<std::string, double>& mConstD,
1877 const std::map<std::string, std::vector<double> >& mConstV,
1878 int eventTimeValid, int eventTime,
1879 std::vector<std::vector<int> > const& rawStTSs,
1880 int charge, double radius, double phi_c,
1881 Belle2::TRGCDCJSignalData* commonData, std::map<std::string, Belle2::TRGCDCJSignal>& mSignalStorage,
1882 std::map<std::string, Belle2::TRGCDCJLUT*>& mLutStorage)
1883{
1884 // Change to Signals.
1885 // rawStTSs[iSt] = [TS ID, TS LR, TS driftTime]
1886 vector<tuple<string, double, int, double, double, int> > t_values = {
1887 make_tuple("phi0", phi_c, mConstD["phiBitSize"], mConstD["phiMin"], mConstD["phiMax"], 0),
1888 make_tuple("rho", radius, mConstD["rhoBitSize"], mConstD["rhoMin"], mConstD["rhoMax"], 0),
1889 make_tuple("charge", (int)(charge == 1 ? 1 : 0), 1, 0, 1.5, 0),
1890 make_tuple("lr_0", rawStTSs[0][1], 2, 0, 3.5, 0),
1891 make_tuple("lr_1", rawStTSs[1][1], 2, 0, 3.5, 0),
1892 make_tuple("lr_2", rawStTSs[2][1], 2, 0, 3.5, 0),
1893 make_tuple("lr_3", rawStTSs[3][1], 2, 0, 3.5, 0),
1894 make_tuple("tsId_0", rawStTSs[0][0], ceil(log(mConstV.at("nTSs")[1]) / log(2)), 0, pow(2, ceil(log(mConstV.at("nTSs")[1]) / log(2))) - 0.5, 0),
1895 make_tuple("tsId_1", rawStTSs[1][0], ceil(log(mConstV.at("nTSs")[3]) / log(2)), 0, pow(2, ceil(log(mConstV.at("nTSs")[3]) / log(2))) - 0.5, 0),
1896 make_tuple("tsId_2", rawStTSs[2][0], ceil(log(mConstV.at("nTSs")[5]) / log(2)), 0, pow(2, ceil(log(mConstV.at("nTSs")[5]) / log(2))) - 0.5, 0),
1897 make_tuple("tsId_3", rawStTSs[3][0], ceil(log(mConstV.at("nTSs")[7]) / log(2)), 0, pow(2, ceil(log(mConstV.at("nTSs")[7]) / log(2))) - 0.5, 0),
1898 make_tuple("tdc_0", rawStTSs[0][2], mConstD["tdcBitSize"], 0, pow(2, mConstD["tdcBitSize"]) - 0.5, 0),
1899 make_tuple("tdc_1", rawStTSs[1][2], mConstD["tdcBitSize"], 0, pow(2, mConstD["tdcBitSize"]) - 0.5, 0),
1900 make_tuple("tdc_2", rawStTSs[2][2], mConstD["tdcBitSize"], 0, pow(2, mConstD["tdcBitSize"]) - 0.5, 0),
1901 make_tuple("tdc_3", rawStTSs[3][2], mConstD["tdcBitSize"], 0, pow(2, mConstD["tdcBitSize"]) - 0.5, 0),
1902 make_tuple("eventTime", eventTime, mConstD["tdcBitSize"], 0, pow(2, mConstD["tdcBitSize"]) - 0.5, 0),
1903 make_tuple("eventTimeValid", eventTimeValid, 1, 0, 1.5, 0),
1904 };
1905 Belle2::TRGCDCJSignal::valuesToMapSignals(t_values, commonData, mSignalStorage);
1906 // Calculate complex fit3D
1907 Fitter3DUtility::setError(mConstD, mConstV, mSignalStorage);
1908 Fitter3DUtility::calPhi(mConstD, mConstV, mSignalStorage, mLutStorage);
1909 Fitter3DUtility::constrainRPerStSl(mConstV, mSignalStorage);
1910 Fitter3DUtility::calZ(mConstD, mConstV, mSignalStorage, mLutStorage);
1911 Fitter3DUtility::calS(mConstD, mConstV, mSignalStorage, mLutStorage);
1912 Fitter3DUtility::rSFit(mConstD, mConstV, mSignalStorage, mLutStorage);
1913 // Post handling of signals.
1914 // Name all signals.
1915 if ((*mSignalStorage.begin()).second.getName() == "") {
1916 for (auto it = mSignalStorage.begin(); it != mSignalStorage.end(); ++it) {
1917 (*it).second.setName((*it).first);
1918 }
1919 }
1920}
1921
1922void Fitter3DUtility::findImpactPosition(ROOT::Math::XYZVector* mcPosition, ROOT::Math::PxPyPzEVector* mcMomentum, int charge,
1923 ROOT::Math::XYVector& helixCenter,
1924 ROOT::Math::XYZVector& impactPosition)
1925{
1926
1927 // Finds the impact position. Everything is in cm, and GeV.
1928 // Input: production vertex (mcx, mcy, mcz),
1929 // momentum at production vertex (px, py, pz)
1930 // charge of particle.
1931 // Output: helix center's coordinate (hcx, hcy)
1932 // impact position (impactX, impactY, impactZ)
1933
1934 double rho = sqrt(pow(mcMomentum->X(), 2) + pow(mcMomentum->Y(), 2)) / 0.3 / 1.5 * 100;
1935 double hcx = mcPosition->X() + rho * cos(atan2(mcMomentum->Y(), mcMomentum->X()) - charge * M_PI_2);
1936 double hcy = mcPosition->Y() + rho * sin(atan2(mcMomentum->Y(), mcMomentum->X()) - charge * M_PI_2);
1937 helixCenter.SetXY(hcx, hcy);
1938 double impactX = (helixCenter.R() - rho) / helixCenter.R() * helixCenter.X();
1939 double impactY = (helixCenter.R() - rho) / helixCenter.R() * helixCenter.Y();
1940 int signdS;
1941 if (atan2(impactY, impactX) < atan2(mcPosition->Y(), mcPosition->X())) signdS = -1;
1942 else signdS = 1;
1943 double dS = 2 * rho * asin(sqrt(pow(impactX - mcPosition->X(), 2) + pow(impactY - mcPosition->Y(), 2)) / 2 / rho);
1944 double impactZ = mcMomentum->Pz() / mcMomentum->Pt() * dS * signdS + mcPosition->Z();
1945 impactPosition.SetXYZ(impactX, impactY, impactZ);
1946
1947}
1948
1949void Fitter3DUtility::calHelixParameters(ROOT::Math::XYZVector position, ROOT::Math::XYZVector momentum, int charge,
1950 TVectorD& helixParameters)
1951{
1952 // HelixParameters: dR, phi0, keppa, dz, tanLambda
1953 double t_alpha = 1 / 0.3 / 1.5 * 100;
1954 double t_pT = momentum.Rho();
1955 double t_R = t_pT * t_alpha;
1956 helixParameters.Clear();
1957 helixParameters.ResizeTo(5);
1958 helixParameters[2] = t_alpha / t_R;
1959 helixParameters[1] = atan2(position.Y() - t_R * momentum.X() / t_pT * charge, position.X() + t_R * momentum.Y() / t_pT * charge);
1960 helixParameters[0] = (position.X() + t_R * momentum.Y() / t_pT * charge) / cos(helixParameters[1]) - t_R;
1961 double t_phi = atan2(-momentum.X() * charge, momentum.Y() * charge) - helixParameters[1];
1962 if (t_phi > M_PI) t_phi -= 2 * M_PI;
1963 if (t_phi < -M_PI) t_phi += 2 * M_PI;
1964 helixParameters[4] = momentum.Z() / t_pT * charge;
1965 helixParameters[3] = position.Z() + helixParameters[4] * t_R * t_phi;
1966}
1967void Fitter3DUtility::calVectorsAtR(const TVectorD& helixParameters, int charge, double cdcRadius, ROOT::Math::XYZVector& position,
1968 ROOT::Math::XYZVector& momentum)
1969{
1970 // HelixParameters: dR, phi0, keppa, dz, tanLambda
1971 double t_alpha = 1 / 0.3 / 1.5 * 100;
1972 double t_R = t_alpha / helixParameters[2];
1973 double t_pT = t_R / t_alpha;
1974 double t_phi = -1 * charge * acos((-pow(cdcRadius, 2) + pow(t_R + helixParameters[0], 2) + pow(t_R,
1975 2)) / (2 * t_R * (t_R + helixParameters[0])));
1976 double t_X = (helixParameters[0] + t_R) * cos(helixParameters[1]) - t_R * cos(helixParameters[1] + t_phi);
1977 double t_Y = (helixParameters[0] + t_R) * sin(helixParameters[1]) - t_R * sin(helixParameters[1] + t_phi);
1978 double t_Z = helixParameters[3] - helixParameters[4] * t_R * t_phi;
1979 double t_Px = -t_pT * sin(helixParameters[1] + t_phi) * charge;
1980 double t_Py = t_pT * cos(helixParameters[1] + t_phi) * charge;
1981 double t_Pz = t_pT * helixParameters[4] * charge;
1982 position.SetXYZ(t_X, t_Y, t_Z);
1983 momentum.SetXYZ(t_Px, t_Py, t_Pz);
1984}
1985
1986int Fitter3DUtility::bitSize(int numberBits, int mode)
1987{
1988 int bitsize = 1;
1989 if (mode == 1) {
1990 for (int i = 0; i < numberBits - 1; i++) {
1991 bitsize *= 2;
1992 }
1993 bitsize -= 1;
1994 } else if (mode == 0) {
1995 for (int i = 0; i < numberBits; i++) {
1996 bitsize *= 2;
1997 }
1998 }
1999 return bitsize;
2000}
2001
2002void Fitter3DUtility::changeInteger(int& integer, double real, double minValue, double maxValue, int bitSize)
2003{
2004 // Version 1.
2005 //double range = maxValue - minValue;
2006 //double convert = bitSize/range;
2008 //integer = int( (real - minValue) * convert +0.5);
2009 // Version 2.
2010 double range = maxValue - minValue;
2011 double convert = (bitSize - 0.5) / range;
2012 integer = FpgaUtility::roundInt((real - minValue) * convert);
2013}
2014
2015void Fitter3DUtility::changeReal(double& real, int integer, double minValue, double maxValue, int bitSize)
2016{
2018 //double range = maxValue - minValue;
2019 //double convert = bitSize/range;
2020 //real = (integer/convert) + minValue;
2021 // Version 2.
2022 double range = maxValue - minValue;
2023 double convert = (bitSize - 0.5) / range;
2024 real = (integer / convert) + minValue;
2025}
2026
2027
2028void Fitter3DUtility::findExtreme(double& m_max, double& m_min, double value)
2029{
2030 if (value > m_max) m_max = value;
2031 if (value < m_min) m_min = value;
2032}
E
R E
internal precision of FFTW codelets
Definition DiscreteCosineTransform_31points.cc:36
Belle2::TRGCDCJLUT
A class to use LUTs for TRGCDC.
Definition JLUT.h:35
Belle2::TRGCDCJSignalData
A class to hold common data for JSignals.
Definition JSignalData.h:33
Belle2::TRGCDCJSignal
A class to use Signals for TRGCDC 3D tracker.
Definition JSignal.h:35
Fitter3DUtility::calVectorsAtR
static void calVectorsAtR(const TVectorD &helixParameters, int charge, double radius, ROOT::Math::XYZVector &position, ROOT::Math::XYZVector &momentum)
Calculates position and momentum at a certain radius.
Definition Fitter3DUtility.cc:1967
Fitter3DUtility::loadStereoXt
static void loadStereoXt(std::string const &filePrefix, int nFiles, std::vector< std::vector< double > > &stXts)
Load stereo Xt file.
Definition Fitter3DUtility.cc:557
Fitter3DUtility::rotatePhi
static double rotatePhi(double value, double refPhi)
Rotates to range [-pi, pi].
Definition Fitter3DUtility.cc:609
Fitter3DUtility::calStAxPhi
static double calStAxPhi(int charge, double anglest, double ztostraw, double rr, double rho, double phi0)
Calculates the fitted axial phi for the stereo super layer.
Definition Fitter3DUtility.cc:785
Fitter3DUtility::changeReal
static void changeReal(double &real, int integer, double minValue, double maxValue, int bitSize)
Changes integer to float value.
Definition Fitter3DUtility.cc:2015
Fitter3DUtility::calHelixParameters
static void calHelixParameters(ROOT::Math::XYZVector position, ROOT::Math::XYZVector momentum, int charge, TVectorD &helixParameters)
HelixParameters: dR, phi0, keppa, dz, tanLambda Calculates the helix parameters of track.
Definition Fitter3DUtility.cc:1949
Fitter3DUtility::findQuadrant
static int findQuadrant(double value)
Finds quadrant of angle. Angle is in rad.
Definition Fitter3DUtility.cc:641
Fitter3DUtility::bitSize
static int bitSize(int numberBits, int mode)
Firmware convert functions.
Definition Fitter3DUtility.cc:1986
Fitter3DUtility::fitter3DFirm
static void fitter3DFirm(std::map< std::string, double > &mConstD, const std::map< std::string, std::vector< double > > &mConstV, int eventTimeValid, int eventTime, std::vector< std::vector< int > > const &rawStTSs, int charge, double radius, double phi_c, Belle2::TRGCDCJSignalData *commonData, std::map< std::string, Belle2::TRGCDCJSignal > &mSignalStorage, std::map< std::string, Belle2::TRGCDCJLUT * > &mLutStorage)
Combines several functions for fitter3D firmware.
Definition Fitter3DUtility.cc:1876
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::changeInteger
static void changeInteger(int &integer, double real, double minValue, double maxValue, int bitSize)
Changes float to integer value.
Definition Fitter3DUtility.cc:2002
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::calDenWithIZError
static double calDenWithIZError(double *arcS, double *iZError)
Calculates the denominator for fitting z and arc s.
Definition Fitter3DUtility.cc:1289
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::setErrorFast
static void setErrorFast(std::vector< std::vector< int > > const &rawStTSs, int eventTimeValid, std::vector< double > &invZError2)
Sets error for fitting.
Definition Fitter3DUtility.cc:663
Fitter3DUtility::calPhiFast
static void calPhiFast(std::map< std::string, std::vector< double > > &stGeometry, std::vector< std::vector< double > > const &stXts, int eventTimeValid, int eventTime, std::vector< std::vector< int > > const &rawStTSs, std::vector< double > &stTSs)
Calculates phi with fast simulation.
Definition Fitter3DUtility.cc:583
Fitter3DUtility::calDeltaPhi
static double calDeltaPhi(int charge, double anglest, double ztostraw, double rr, double phi, double rho, double phi0)
Calculates the phi difference between fitted axial phi and stereo phi.
Definition Fitter3DUtility.cc:805
Fitter3DUtility::calS
static double calS(double rho, double rr)
Calculates arc length.
Definition Fitter3DUtility.cc:1204
Fitter3DUtility::stereoIdToPhi
static double stereoIdToPhi(std::vector< double > &stNWires, int iSt, int id)
Converts stereo ts id to phi.
Definition Fitter3DUtility.cc:576
Fitter3DUtility::fitter3D
static void fitter3D(std::map< std::string, std::vector< double > > &stGeometry, std::vector< std::vector< double > > const &stXts, int eventTimeValid, int eventTime, std::vector< std::vector< int > > const &rawStTSs, int charge, double radius, double phi_c, double &z0, double &cot, double &chi2)
Combines several functions for fitter3D.
Definition Fitter3DUtility.cc:1842
Fitter3DUtility::findSign
static int findSign(double *phi2)
2D fitter functions
Definition Fitter3DUtility.cc:40
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::saveStereoXt
static void saveStereoXt(std::vector< std::vector< double > > &stXts, std::string const &filePrefix)
Saves stereo Xt to file.
Definition Fitter3DUtility.cc:543
Fitter3DUtility::rPhiFit
static void rPhiFit(double *rr, double *phi2, double *phierror, double &rho, double &myphi0)
A circle fitter.
Definition Fitter3DUtility.cc:60
Fitter3DUtility::rSFit2
static void rSFit2(double *iezz21, double *iezz22, double *arcS, double *zz, int *lutv, double &z0, double &cot, double &zchi2)
Fits z and arc S. (Will be deprecated.)
Definition Fitter3DUtility.cc:1343
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::calDen
static double calDen(double *arcS, double *zError)
Calculates the denominator for fitting z and arc s.
Definition Fitter3DUtility.cc:1276
Fitter3DUtility::rPhiFit2
static void rPhiFit2(double *rr, double *phi2, double *phierror, double &rho, double &myphi0, int nTS)
A circle fitter.
Definition Fitter3DUtility.cc:179
Fitter3DUtility::findExtreme
static void findExtreme(double &m_max, double &m_min, double value)
Finds maximum and minimum values.
Definition Fitter3DUtility.cc:2028
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::roundInt
static double roundInt(double value)
Round double value.
Definition FpgaUtility.cc:38
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::TRGCDCJSignal::ifElse
static void ifElse(std::vector< std::pair< TRGCDCJSignal, std::vector< std::pair< TRGCDCJSignal *, TRGCDCJSignal > > > > &data, int targetClock)
If else implementation with target clock.
Definition JSignal.cc:800
Belle2::TRGCDCJSignal::absolute
static TRGCDCJSignal const absolute(TRGCDCJSignal const &first)
Absolute TRGCDCJSignal. Removes 1 bit if signed or minus unsigned.
Definition JSignal.cc:1562
Belle2::TRGCDCJSignal::comp
static TRGCDCJSignal comp(TRGCDCJSignal const &lhs, const std::string &operate, TRGCDCJSignal const &rhs)
Compare two signals.
Definition JSignal.cc:1169