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