Belle II Software development
Static Public Member Functions | List of all members
Fitter3DUtility Class Reference

A class that holds functions for 3D tracking. More...

#include <Fitter3DUtility.h>

Static Public Member Functions

static int findSign (double *phi2)
 2D fitter functions
 
static void rPhiFit (double *rr, double *phi2, double *phierror, double &rho, double &myphi0)
 A circle fitter.
 
static void rPhiFitter (double *rr, double *phi2, double *invphierror, double &rho, double &myphi0)
 A circle fitter with invPhiError without fit chi2 output.
 
static void rPhiFitter (double *rr, double *phi2, double *invphierror, double &rho, double &myphi0, double &chi2)
 A circle fitter with invPhiError with fit chi2 output.
 
static void chargeFinder (double *nTSs, double *tsIds, double *tsHitmap, double phi0, double inCharge, double &outCharge)
 Charge finder using circle fitter output and axial TSs.
 
static void rPhiFit2 (double *rr, double *phi2, double *phierror, double &rho, double &myphi0, int nTS)
 A circle fitter.
 
static unsigned toUnsignedTdc (int tdc, int nBits)
 Changes tdc and event time to unsigned value that has # bits.
 
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.
 
static double calPhi (double wirePhi, double driftTime, double eventTime, double rr, int lr)
 Pre 3D fitter functions.
 
static double calPhi (int localId, int nWires, double driftTime, double eventTime, double rr, int lr)
 Pre 3D fitter functions.
 
static void calPhi (std::map< std::string, double > const &mConstD, std::map< std::string, std::vector< double > > const &mConstV, std::map< std::string, Belle2::TRGCDCJSignal > &mSignalStorage, std::map< std::string, Belle2::TRGCDCJLUT * > &mLutStorage)
 Calculates phi.
 
static void saveStereoXt (std::vector< std::vector< double > > &stXts, std::string const &filePrefix)
 Saves stereo Xt to file.
 
static void loadStereoXt (std::string const &filePrefix, int nFiles, std::vector< std::vector< double > > &stXts)
 Load stereo Xt file.
 
static double stereoIdToPhi (std::vector< double > &stNWires, int iSt, int id)
 Converts stereo ts id to phi.
 
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.
 
static double rotatePhi (double value, double refPhi)
 Rotates to range [-pi, pi].
 
static double rotatePhi (double value, int refId, int nTSs)
 Rotates to range [-pi, pi]. Use tsId as reference.
 
static int rotateTsId (int value, int refId, int nTSs)
 Rotates to range [0, nTSs-1].
 
static int findQuadrant (double value)
 Finds quadrant of angle. Angle is in rad.
 
static void setErrorFast (std::vector< std::vector< int > > const &rawStTSs, int eventTimeValid, std::vector< double > &invZError2)
 Sets error for fitting.
 
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.
 
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.
 
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.
 
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.
 
static double calZ (int charge, double anglest, double ztostraw, double rr, double phi, double rho, double phi0)
 Calculates z.
 
static void calZ (std::map< std::string, double > const &mConstD, std::map< std::string, std::vector< double > > const &mConstV, std::map< std::string, Belle2::TRGCDCJSignal > &mSignalStorage, std::map< std::string, Belle2::TRGCDCJLUT * > &mLutStorage)
 Calculates z using JSignal and multiple constants.
 
static double calS (double rho, double rr)
 Calculates arc length.
 
static void calS (std::map< std::string, double > const &mConstD, std::map< std::string, std::vector< double > > const &mConstV, std::map< std::string, Belle2::TRGCDCJSignal > &mSignalStorage, std::map< std::string, Belle2::TRGCDCJLUT * > &mLutStorage)
 Calculates arc length using JSignal and multiple constants.
 
static double calIZError2 (int lr, double driftZError, double wireZError)
 Chooses and calculates inverse z error.
 
static double calDen (double *arcS, double *zError)
 Calculates the denominator for fitting z and arc s.
 
static double calDenWithIZError (double *arcS, double *iZError)
 Calculates the denominator for fitting z and arc s.
 
static void rSFit (double *iezz2, double *arcS, double *zz, double &z0, double &cot, double &zchi2)
 3D fitter functions Fits z and arc S.
 
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.)
 
static void rSFit (std::map< std::string, double > const &mConstD, std::map< std::string, std::vector< double > > const &mConstV, std::map< std::string, Belle2::TRGCDCJSignal > &mSignalStorage, std::map< std::string, Belle2::TRGCDCJLUT * > &mLutStorage)
 Fits z and arc S using JSingal with multiple constants.
 
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.
 
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, std::vector< double > &arcS, std::vector< double > &zz, std::vector< double > &invZError2)
 Combines several functions for fitter3D. Also outputs arcS, zz and invZError2.
 
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.
 
static void findImpactPosition (ROOT::Math::XYZVector *mcPosition, ROOT::Math::PxPyPzEVector *mcMomentum, int charge, ROOT::Math::XYVector &helixCenter, ROOT::Math::XYZVector &impactPosition)
 MC calculation functions Calculates the impact position of track.
 
static void calHelixParameters (ROOT::Math::XYZVector position, ROOT::Math::XYZVector momentum, int charge, TVectorD &helixParameters)
 HelixParameters: dR, phi0, keppa, dz, tanLambda Calculates the helix parameters of track.
 
static void calVectorsAtR (const TVectorD &helixParameters, int charge, double radius, ROOT::Math::XYZVector &position, ROOT::Math::XYZVector &momentum)
 Calculates position and momentum at a certain radius.
 
static int bitSize (int numberBits, int mode)
 Firmware convert functions.
 
static void changeInteger (int &integer, double real, double minValue, double maxValue, int bitSize)
 Changes float to integer value.
 
static void changeReal (double &real, int integer, double minValue, double maxValue, int bitSize)
 Changes integer to float value.
 
static void findExtreme (double &m_max, double &m_min, double value)
 Finds maximum and minium values.
 

Detailed Description

A class that holds functions for 3D tracking.

Definition at line 32 of file Fitter3DUtility.h.

Member Function Documentation

◆ bitSize()

int bitSize ( int  numberBits,
int  mode 
)
static

Firmware convert functions.

Old firmware functions. Currently used in Hough3DFinder. Calculates maximum for a given number of bits. (Will be deprecated.)

Definition at line 1986 of file Fitter3DUtility.cc.

1987{
1988 int bitsize = 1;
1989 if (mode == 1) {
1990 for (int i = 0; i < numberBits - 1; i++) {
1991 bitsize *= 2;
1992 }
1993 bitsize -= 1;
1994 } else if (mode == 0) {
1995 for (int i = 0; i < numberBits; i++) {
1996 bitsize *= 2;
1997 }
1998 }
1999 return bitsize;
2000}

◆ calDeltaPhi()

double calDeltaPhi ( int  charge,
double  anglest,
double  ztostraw,
double  rr,
double  phi,
double  rho,
double  phi0 
)
static

Calculates the phi difference between fitted axial phi and stereo phi.

Definition at line 805 of file Fitter3DUtility.cc.

806{
807 if (false) cout << anglest << ztostraw << endl; // Removes warnings when compileing
808
809 double myphiz, acos_real;
810 //Find phifit-phist
811 double t_rho = rho;
812 if (rr > (2 * rho)) t_rho = rr / 2;
813 acos_real = acos(rr / (2 * t_rho));
814 if (mysign == 1) {
815 myphiz = -acos_real - myphi0 + phi2;
816 } else {
817 myphiz = +acos_real - myphi0 + phi2;
818 }
819 if (myphiz > M_PI) myphiz -= 2 * M_PI;
820 if (myphiz < -M_PI) myphiz += 2 * M_PI;
821
822 return myphiz;
823}

◆ calDen()

double calDen ( double *  arcS,
double *  zError 
)
static

Calculates the denominator for fitting z and arc s.

Definition at line 1276 of file Fitter3DUtility.cc.

1277{
1278 double t_sum1 = 0;
1279 double t_sumSS = 0;
1280 double t_sumS = 0;
1281 for (unsigned iSt = 0; iSt < 4; iSt++) {
1282 t_sum1 += pow(1 / zError[iSt], 2);
1283 t_sumSS += pow(arcS[iSt] / zError[iSt], 2);
1284 t_sumS += arcS[iSt] / pow(zError[iSt], 2);
1285 }
1286 return t_sum1 * t_sumSS - pow(t_sumS, 2);
1287}

◆ calDenWithIZError()

double calDenWithIZError ( double *  arcS,
double *  iZError 
)
static

Calculates the denominator for fitting z and arc s.

Definition at line 1289 of file Fitter3DUtility.cc.

1290{
1291 double t_sum1 = 0;
1292 double t_sumSS = 0;
1293 double t_sumS = 0;
1294 for (unsigned iSt = 0; iSt < 4; iSt++) {
1295 t_sum1 += pow(iZError[iSt], 2);
1296 t_sumSS += pow(arcS[iSt] * iZError[iSt], 2);
1297 t_sumS += arcS[iSt] * pow(iZError[iSt], 2);
1298 }
1299 return t_sum1 * t_sumSS - pow(t_sumS, 2);
1300}

◆ calHelixParameters()

void calHelixParameters ( ROOT::Math::XYZVector  position,
ROOT::Math::XYZVector  momentum,
int  charge,
TVectorD &  helixParameters 
)
static

HelixParameters: dR, phi0, keppa, dz, tanLambda Calculates the helix parameters of track.

Definition at line 1949 of file Fitter3DUtility.cc.

1951{
1952 // HelixParameters: dR, phi0, keppa, dz, tanLambda
1953 double t_alpha = 1 / 0.3 / 1.5 * 100;
1954 double t_pT = momentum.Rho();
1955 double t_R = t_pT * t_alpha;
1956 helixParameters.Clear();
1957 helixParameters.ResizeTo(5);
1958 helixParameters[2] = t_alpha / t_R;
1959 helixParameters[1] = atan2(position.Y() - t_R * momentum.X() / t_pT * charge, position.X() + t_R * momentum.Y() / t_pT * charge);
1960 helixParameters[0] = (position.X() + t_R * momentum.Y() / t_pT * charge) / cos(helixParameters[1]) - t_R;
1961 double t_phi = atan2(-momentum.X() * charge, momentum.Y() * charge) - helixParameters[1];
1962 if (t_phi > M_PI) t_phi -= 2 * M_PI;
1963 if (t_phi < -M_PI) t_phi += 2 * M_PI;
1964 helixParameters[4] = momentum.Z() / t_pT * charge;
1965 helixParameters[3] = position.Z() + helixParameters[4] * t_R * t_phi;
1966}
Belle2::EvtPDLUtil::charge
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition: EvtPDLUtil.cc:44

◆ calPhi() [1/4]

double calPhi ( double  wirePhi,
double  driftLength,
double  rr,
int  lr 
)
static

Pre 3D fitter functions. rr is in cm scale. driftLength is in cm scale.

Definition at line 239 of file Fitter3DUtility.cc.

240{
241 //cout<<"rr:"<<rr<<" lr:"<<lr;
242 double result = wirePhi;
243 //cout<<" wirePhi:"<<result;
244 // time is in 2ns rms clock.
245 double t_dPhi = driftLength;
246 //cout<<" driftLength:"<<t_dPhi;
247 // Change to radian
248 // rr is cm scale.
249 t_dPhi = atan(t_dPhi / rr);
250 //cout<<" driftPhi:"<<t_dPhi;
251 // Use LR to add dPhi
252 if (lr == 1) result -= t_dPhi;
253 else if (lr == 2) result += t_dPhi;
254 //cout<<" phi:"<<result<<endl;
255 return result;
256}

◆ calPhi() [2/4]

double calPhi ( double  wirePhi,
double  driftTime,
double  eventTime,
double  rr,
int  lr 
)
static

Pre 3D fitter functions.

Definition at line 258 of file Fitter3DUtility.cc.

259{
260 // time is in 2ns rms clock.
261 double t_driftLength = (driftTime - eventTime) * 1000 / 1017.774 * 2 * 40 / 1000 / 10;
262 return calPhi(wirePhi, t_driftLength, rr, lr);
263}
Fitter3DUtility::calPhi
static double calPhi(double wirePhi, double driftLength, double rr, int lr)
Pre 3D fitter functions. rr is in cm scale. driftLength is in cm scale.
Definition: Fitter3DUtility.cc:239

◆ calPhi() [3/4]

double calPhi ( int  localId,
int  nWires,
double  driftTime,
double  eventTime,
double  rr,
int  lr 
)
static

Pre 3D fitter functions.

Definition at line 265 of file Fitter3DUtility.cc.

266{
267 double wirePhi = (double)localId / nWires * 4 * M_PI;
268 return Fitter3DUtility::calPhi(wirePhi, driftTime, eventTime, rr, lr);
269}

◆ calPhi() [4/4]

void calPhi ( std::map< std::string, double > const &  mConstD,
std::map< std::string, std::vector< double > > const &  mConstV,
std::map< std::string, Belle2::TRGCDCJSignal > &  mSignalStorage,
std::map< std::string, Belle2::TRGCDCJLUT * > &  mLutStorage 
)
static

Calculates phi.

Definition at line 271 of file Fitter3DUtility.cc.

274{
275 Belle2::TRGCDCJSignalData* commonData = mSignalStorage.begin()->second.getCommonData();
276 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<tsId_"<<iSt<<">>>"<<endl; mSignalStorage["tsId_"+to_string(iSt)].dump();}
277 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<tdc_"<<iSt<<">>>"<<endl; mSignalStorage["tdc_"+to_string(iSt)].dump();}
278 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<lr_"<<iSt<<">>>"<<endl; mSignalStorage["lr_"+to_string(iSt)].dump();}
279 // Make phiFactor constants
280 if (mSignalStorage.find("phiFactor_0") == mSignalStorage.end()) {
281 for (unsigned iSt = 0; iSt < 4; iSt++) {
282 int nShiftBits = int(log(pow(2, 24) / 2 / mConstD.at("M_PI") * mConstV.at("nTSs")[2 * iSt + 1] *
283 mSignalStorage["phi0"].getToReal()) / log(2));
284 string t_name;
285 t_name = "phiFactor_" + to_string(iSt);
286 mSignalStorage[t_name] = Belle2::TRGCDCJSignal(2 * mConstD.at("M_PI") / mConstV.at("nTSs")[2 * iSt + 1],
287 mSignalStorage["phi0"].getToReal() / pow(2, nShiftBits), commonData);
288 }
289 }
290 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<phiFactor_"<<iSt<<">>>"<<endl; mSignalStorage["phiFactor_"+to_string(iSt)].dump();}
291 // wirePhi <= tsId * phiFactor
292 for (unsigned iSt = 0; iSt < 4;
293 iSt++) mSignalStorage["wirePhi_" + to_string(iSt)] <= mSignalStorage["tsId_" + to_string(iSt)] * mSignalStorage["phiFactor_" +
294 to_string(iSt)];
295 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<wirePhi_"<<iSt<<">>>"<<endl; mSignalStorage["wirePhi_"+to_string(iSt)].dump();}
296
297 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<tdc_"<<iSt<<">>>"<<endl; mSignalStorage["tdc_"+to_string(iSt)].dump();}
298 //cout<<"<<<eventTime>>>"<<endl; mSignalStorage["eventTime"].dump();
299
300 // Calculate driftTimeRaw (= tdc - eventTime)
301 for (unsigned iSt = 0; iSt < 4;
302 iSt++) mSignalStorage["driftTimeRaw_" + to_string(iSt)] <= mSignalStorage["tdc_" + to_string(iSt)] - mSignalStorage["eventTime"];
303
304 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<driftTimeRaw_"<<iSt<<">>>"<<endl; mSignalStorage["driftTimeRaw_"+to_string(iSt)].dump();}
305
306 // Confine drift time
307 // Constants for driftTime
308 if (mSignalStorage.find("maxDriftTimeLow") == mSignalStorage.end()) {
309 int maxDriftTime = 287;
310 mSignalStorage["maxDriftTimeLow"] = Belle2::TRGCDCJSignal(-512 + maxDriftTime + 1, mSignalStorage["eventTime"].getToReal(),
311 commonData);
312 mSignalStorage["maxDriftTimeHigh"] = Belle2::TRGCDCJSignal(maxDriftTime, mSignalStorage["eventTime"].getToReal(), commonData);
313 mSignalStorage["maxDriftTimeMax"] = Belle2::TRGCDCJSignal(512, mSignalStorage["eventTime"].getToReal(), commonData);
314 }
315 for (unsigned iSt = 0; iSt < 4; iSt++) {
316 string t_in1Name = "driftTimeRaw_" + to_string(iSt);
317 string t_valueName = "driftTime_c_" + to_string(iSt);
318 // Create data for ifElse
319 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
320 // Compare (eventTimeValid == 0)
321 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["eventTimeValid"], "=", Belle2::TRGCDCJSignal(0,
322 mSignalStorage["eventTimeValid"].getToReal(), commonData));
323 // Assignments (driftTime = 0)
324 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
325 make_pair(&mSignalStorage[t_valueName], Belle2::TRGCDCJSignal(0, mSignalStorage["eventTime"].getToReal(), commonData))
326 };
327 // Push to data.
328 t_data.push_back(make_pair(t_compare, t_assigns));
329 // Compare (driftTime < maxDriftTimeLow)
330 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], "<", mSignalStorage["maxDriftTimeLow"]);
331 // Assignments
332 t_assigns = {
333 make_pair(&mSignalStorage[t_valueName], (mSignalStorage[t_in1Name] + mSignalStorage["maxDriftTimeMax"]).limit(mSignalStorage["maxDriftTimeLow"], mSignalStorage["maxDriftTimeHigh"]))
334 };
335 // Push to data.
336 t_data.push_back(make_pair(t_compare, t_assigns));
337 // Compare (driftTime < maxDriftTimeHigh)
338 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], "<=", mSignalStorage["maxDriftTimeHigh"]);
339 // Assignments
340 t_assigns = {
341 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_in1Name].limit(mSignalStorage["maxDriftTimeLow"], mSignalStorage["maxDriftTimeHigh"]))
342 };
343 // Push to data.
344 t_data.push_back(make_pair(t_compare, t_assigns));
345 // Compare (else)
346 t_compare = Belle2::TRGCDCJSignal();
347 // Assignments
348 t_assigns = {
349 make_pair(&mSignalStorage[t_valueName], (mSignalStorage[t_in1Name] - mSignalStorage["maxDriftTimeMax"]).limit(mSignalStorage["maxDriftTimeLow"], mSignalStorage["maxDriftTimeHigh"]))
350 };
351 // Push to data.
352 t_data.push_back(make_pair(t_compare, t_assigns));
353 // Process ifElse data.
354 Belle2::TRGCDCJSignal::ifElse(t_data);
355 }
356 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<driftTime_c_"<<iSt<<">>>"<<endl; mSignalStorage["driftTime_c_"+to_string(iSt)].dump();}
357
359 //for(unsigned iSt=0; iSt<4; iSt++) {
360 // if(mSignalStorage["driftTimeRaw_"+to_string(iSt)].getInt() > maxDriftTime) cout<<"driftTimeRaw_"<<iSt<<":"<<mSignalStorage["driftTimeRaw_"+to_string(iSt)].getInt()<<" is larger than maxDriftTime:"<<maxDriftTime<<" changing to "<<mSignalStorage["driftTimeRaw_"+to_string(iSt)].getInt()-512<<endl;
361 // if(mSignalStorage["driftTimeRaw_"+to_string(iSt)].getInt() < -512+maxDriftTime+1) cout<<"driftTimeRaw_"<<iSt<<":"<<mSignalStorage["driftTimeRaw_"+to_string(iSt)].getInt()<<" is smaller than -512+maxDriftTime+1:"<<-512+maxDriftTime+1<<" changed to "<<mSignalStorage["driftTimeRaw_"+to_string(iSt)].getInt()+512<<endl;
362 //}
363
364 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<invErrorRho_"<<iSt<<">>>"<<endl; mSignalStorage["invErrorRho_"+to_string(iSt)].dump();}
365 // Remove negative driftTime
366 // Constants for driftTime
367 if (mSignalStorage.find("minDriftTime") == mSignalStorage.end()) {
368 int minDriftTime = -8;
369 mSignalStorage["minDriftTime"] = Belle2::TRGCDCJSignal(minDriftTime, mSignalStorage["eventTime"].getToReal(), commonData);
370 mSignalStorage["minDriftTimeLow"] = Belle2::TRGCDCJSignal(0, mSignalStorage["eventTime"].getToReal(), commonData);
371 }
372 for (unsigned iSt = 0; iSt < 4; iSt++) {
373 string t_in1Name = "driftTime_c_" + to_string(iSt);
374 string t_valueName = "driftTime_" + to_string(iSt);
375 string t_in2Name = "invErrorRho_" + to_string(iSt);
376 string t_value2Name = "iZError2_" + to_string(iSt);
377 string t_noneErrorName = "iZNoneError2_" + to_string(iSt);
378 // Create data for ifElse
379 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
380 // Compare (driftTime_c[i]>0)
381 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">=", mSignalStorage["minDriftTimeLow"]);
382 // Assignments
383 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
384 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_in1Name].limit(mSignalStorage["minDriftTimeLow"], mSignalStorage["maxDriftTimeHigh"])),
385 make_pair(&mSignalStorage[t_value2Name], mSignalStorage[t_in2Name])
386 };
387 // Push to data.
388 t_data.push_back(make_pair(t_compare, t_assigns));
389 // Compare (driftTime_c[i]>=minDriftTime)
390 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">=", mSignalStorage["minDriftTime"]);
391 // Assignments
392 t_assigns = {
393 make_pair(&mSignalStorage[t_valueName], mSignalStorage["minDriftTimeLow"]),
394 make_pair(&mSignalStorage[t_value2Name], mSignalStorage[t_in2Name])
395 };
396 // Push to data.
397 t_data.push_back(make_pair(t_compare, t_assigns));
398 // Compare (else)
399 t_compare = Belle2::TRGCDCJSignal();
400 // Assignments
401 t_assigns = {
402 make_pair(&mSignalStorage[t_valueName], mSignalStorage["minDriftTimeLow"]),
403 make_pair(&mSignalStorage[t_value2Name], mSignalStorage[t_noneErrorName])
404 };
405 // Push to data.
406 t_data.push_back(make_pair(t_compare, t_assigns));
407 // Process ifElse data.
408 Belle2::TRGCDCJSignal::ifElse(t_data);
409 }
410 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<driftTime_"<<iSt<<">>>"<<endl; mSignalStorage["driftTime_"+to_string(iSt)].dump();}
411 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<iZError2_"<<iSt<<">>>"<<endl; mSignalStorage["iZError2_"+to_string(iSt)].dump();}
412
414 //for(unsigned iSt=0; iSt<4; iSt++) {
415 // if(mSignalStorage["driftTime_c_"+to_string(iSt)].getInt() < minDriftTime) cout<<"driftTime_c_"<<iSt<<":"<<mSignalStorage["driftTime_c_"+to_string(iSt)].getInt()<<" is smaller than minDriftTime:"<<minDriftTime<<endl;
416 // else if(mSignalStorage["driftTime_c_"+to_string(iSt)].getInt() < 0) cout<<"driftTime_c_"<<iSt<<":"<<mSignalStorage["driftTime_c_"+to_string(iSt)].getInt()<<" is smaller than 0."<<endl;
417 //}
418
419 // Calculate minInvValue, maxInvValue for LUTs
420 if (!mSignalStorage.count("invDriftPhiMin")) { //nothing found
421 mSignalStorage["invDriftPhiMin"] = Belle2::TRGCDCJSignal(0, mSignalStorage["driftTime_0"].getToReal(), commonData);
422 mSignalStorage["invDriftPhiMax"] = Belle2::TRGCDCJSignal(pow(2, mConstD.at("tdcBitSize")) - 1,
423 mSignalStorage["driftTime_0"].getToReal(), commonData);
424 }
425 //cout<<"<<<invDriftPhiMin>>>"<<endl; mSignalStorage["invDriftPhiMin"].dump();
426 //cout<<"<<<invDriftPhiMax>>>"<<endl; mSignalStorage["invDriftPhiMax"].dump();
427 // Generate LUT(driftPhi[i]=arctan(driftLengthFunction(driftTime))/r[i])
428 for (unsigned iSt = 0; iSt < 4; iSt++) {
429 string t_valueName = "driftPhi_" + to_string(iSt);
430 string t_inName = "driftTime_" + to_string(iSt);
431 if (!mLutStorage.count(t_valueName)) { //if nothing found
432 mLutStorage[t_valueName] = new Belle2::TRGCDCJLUT(t_valueName);
433 // Lambda can not copy maps.
434 double t_parameter = mConstV.at("rr3D")[iSt];
435 mLutStorage[t_valueName]->setFloatFunction(
436 //[=](double aValue) -> double{return mConstV.at("driftLengthTableSL"+to_string(2*iSt+1))[aValue]/t_parameter;},
437 [ = ](double aValue) -> double{return atan(mConstV.at("driftLengthTableSL" + to_string(2 * iSt + 1))[aValue] / t_parameter);},
438 mSignalStorage[t_inName],
439 mSignalStorage["invDriftPhiMin"], mSignalStorage["invDriftPhiMax"], mSignalStorage["phi0"].getToReal(),
440 (int)mConstD.at("driftPhiLUTInBitSize"), (int)mConstD.at("driftPhiLUTOutBitSize"));
441 //mLutStorage[t_valueName]->makeCOE("./LutData/"+t_valueName+".coe");
442 }
443 }
446 for (unsigned iSt = 0; iSt < 4; iSt++) {
447 //cout<<"driftTime:"<<mSignalStorage["driftTime_"+to_string(iSt)].getActual()<<endl;
448 //cout<<"actual:"<<mLutStorage["driftPhi_"+to_string(iSt)]->getFloatOutput(mSignalStorage["driftTime_"+to_string(iSt)].getActual())<<endl;
449 mLutStorage["driftPhi_" + to_string(iSt)]->operate(mSignalStorage["driftTime_" + to_string(iSt)],
450 mSignalStorage["driftPhi_" + to_string(iSt)]);
451 }
452 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<driftPhi_"<<iSt<<">>>"<<endl; mSignalStorage["driftPhi_"+to_string(iSt)].dump();}
453
454 // Calculate finePhi depending on lr(0:no hit, 1: left, 2: right, 3: not determined)
455 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<lr_"<<iSt<<">>>"<<endl; mSignalStorage["lr_"+to_string(iSt)].dump();}
456 for (unsigned iSt = 0; iSt < 4; iSt++) {
457 string t_in1Name = "lr_" + to_string(iSt);
458 string t_valueName = "finePhi_" + to_string(iSt);
459 // Create data for ifElse
460 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
461 // Compare (lr == 1)
462 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], "=", Belle2::TRGCDCJSignal(1,
463 mSignalStorage[t_in1Name].getToReal(), commonData));
464 // Assignments
465 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
466 make_pair(&mSignalStorage[t_valueName], mSignalStorage["wirePhi_" + to_string(iSt)] - mSignalStorage["driftPhi_" + to_string(iSt)])
467 };
468 // Push to data.
469 t_data.push_back(make_pair(t_compare, t_assigns));
470 // Compare (lr == 2)
471 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], "=", Belle2::TRGCDCJSignal(2,
472 mSignalStorage[t_in1Name].getToReal(), commonData));
473 // Assignments
474 t_assigns = {
475 make_pair(&mSignalStorage[t_valueName], mSignalStorage["wirePhi_" + to_string(iSt)] + mSignalStorage["driftPhi_" + to_string(iSt)])
476 };
477 // Push to data.
478 t_data.push_back(make_pair(t_compare, t_assigns));
479 // Compare (else)
480 t_compare = Belle2::TRGCDCJSignal();
481 // Assignments
482 t_assigns = {
483 make_pair(&mSignalStorage[t_valueName], mSignalStorage["wirePhi_" + to_string(iSt)])
484 };
485 // Push to data.
486 t_data.push_back(make_pair(t_compare, t_assigns));
487 // Process ifElse data.
488 Belle2::TRGCDCJSignal::ifElse(t_data);
489 }
490 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<finePhi_"<<iSt<<">>>"<<endl; mSignalStorage["finePhi_"+to_string(iSt)].dump();}
491
492 // Constrain phi[i] to [-pi, pi]
493 if (mSignalStorage.find("phiMax_0") == mSignalStorage.end()) {
494 for (unsigned iSt = 0; iSt < 4; iSt++) {
495 string t_valueName = "finePhi_" + to_string(iSt);
496 string t_maxName = "phiMax_" + to_string(iSt);
497 string t_minName = "phiMin_" + to_string(iSt);
498 string t_2PiName = "phi2Pi_" + to_string(iSt);
499 mSignalStorage[t_maxName] = Belle2::TRGCDCJSignal(mConstD.at("M_PI"), mSignalStorage[t_valueName].getToReal(), commonData);
500 mSignalStorage[t_minName] = Belle2::TRGCDCJSignal(-mConstD.at("M_PI"), mSignalStorage[t_valueName].getToReal(), commonData);
501 mSignalStorage[t_2PiName] = Belle2::TRGCDCJSignal(2 * mConstD.at("M_PI"), mSignalStorage[t_valueName].getToReal(), commonData);
502 }
503 }
504 for (unsigned iSt = 0; iSt < 4; iSt++) {
505 string t_in1Name = "finePhi_" + to_string(iSt);
506 string t_valueName = "phi_" + to_string(iSt);
507 string t_maxName = "phiMax_" + to_string(iSt);
508 string t_minName = "phiMin_" + to_string(iSt);
509 string t_2PiName = "phi2Pi_" + to_string(iSt);
510 // Create data for ifElse
511 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
512 // Compare (dPhi[i] > dPhiPiMax)
513 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_maxName]);
514 // Assignments (dPhiPi_c[i] = dPhi[i] - dPhiPi2Pi[i])
515 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
516 make_pair(&mSignalStorage[t_valueName], (mSignalStorage[t_in1Name] - mSignalStorage[t_2PiName]).limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
517 };
518 // Push to data.
519 t_data.push_back(make_pair(t_compare, t_assigns));
520 // Compare (dPhi[i] > dPhiPiMin)
521 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_minName]);
522 // Assignments (dPhiPi_c[i] = dPhi[i])
523 t_assigns = {
524 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_in1Name].limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
525 };
526 // Push to data.
527 t_data.push_back(make_pair(t_compare, t_assigns));
528 // Compare (else)
529 t_compare = Belle2::TRGCDCJSignal();
530 // Assignments (dPhiPi_c[i] = dPhi[i] + dPhiPi2Pi[i])
531 t_assigns = {
532 make_pair(&mSignalStorage[t_valueName], (mSignalStorage[t_in1Name] + mSignalStorage[t_2PiName]).limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
533 };
534 // Push to data.
535 t_data.push_back(make_pair(t_compare, t_assigns));
536 // Process ifElse data.
537 Belle2::TRGCDCJSignal::ifElse(t_data);
538 }
539 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<phi_"<<iSt<<">>>"<<endl; mSignalStorage["phi_"+to_string(iSt)].dump();}
540
541}
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
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::comp
static TRGCDCJSignal comp(TRGCDCJSignal const &lhs, const std::string &operate, TRGCDCJSignal const &rhs)
Compare two signals.
Definition: JSignal.cc:1169

◆ calPhiFast()

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 
)
static

Calculates phi with fast simulation.

Definition at line 583 of file Fitter3DUtility.cc.

585{
586 stTSs = vector<double> (4, 9999);
587 if (eventTimeValid == 0) {
588 for (unsigned iSt = 0; iSt < 4; iSt++) {
589 stTSs[iSt] = stereoIdToPhi(stGeometry["nWires"], iSt, rawStTSs[iSt][0]);
590 }
591 } else {
592 //cout<<"eventTime: "<<eventTime;
593 for (unsigned iSt = 0; iSt < 4; iSt++) {
594 if (rawStTSs[iSt][1] != 0) {
595 // Get drift length from table.
596 int t = rawStTSs[iSt][2] - eventTime;
597 //cout<<" iSt:"<<iSt<<" rt: "<<rawStTSs[iSt][2];
598 if (t < 0) t = 0;
599 if (t > 511) t = 511;
600 double driftLength = stXts[iSt][t];
601 stTSs[iSt] = Fitter3DUtility::calPhi(stereoIdToPhi(stGeometry["nWires"], iSt, rawStTSs[iSt][0]), driftLength,
602 stGeometry["cdcRadius"][iSt], rawStTSs[iSt][1]);
603 }
604 }
605 //cout<<endl;
606 }
607}
Fitter3DUtility::stereoIdToPhi
static double stereoIdToPhi(std::vector< double > &stNWires, int iSt, int id)
Converts stereo ts id to phi.
Definition: Fitter3DUtility.cc:576

◆ calS() [1/2]

double calS ( double  rho,
double  rr 
)
static

Calculates arc length.

Definition at line 1204 of file Fitter3DUtility.cc.

1205{
1206 double result;
1207 result = rho * 2 * asin(rr / 2 / rho);
1208 return result;
1209}

◆ calS() [2/2]

void calS ( std::map< std::string, double > const &  mConstD,
std::map< std::string, std::vector< double > > const &  mConstV,
std::map< std::string, Belle2::TRGCDCJSignal > &  mSignalStorage,
std::map< std::string, Belle2::TRGCDCJLUT * > &  mLutStorage 
)
static

Calculates arc length using JSignal and multiple constants.

Definition at line 1211 of file Fitter3DUtility.cc.

1213{
1214 Belle2::TRGCDCJSignalData* commonData = mSignalStorage.begin()->second.getCommonData();
1215 // Calculate minInvValue, maxInvValue for LUTs
1216 if (mSignalStorage.find("invArcSMin_0") == mSignalStorage.end()) {
1217 //for(unsigned iSt=0; iSt<4; iSt++) {
1218 // string t_invMinName = "invArcSMin_" + to_string(iSt);
1219 // double t_actual = mSignalStorage["rho"].getMaxActual();
1220 // double t_toReal = mSignalStorage["rho"].getToReal();
1221 // mSignalStorage[t_invMinName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
1222 // string t_invMaxName = "invArcSMax_" + to_string(iSt);
1223 // t_actual = (mConstV.find("rr")->second)[2*iSt+1]/2;
1224 // mSignalStorage[t_invMaxName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
1225 // // Subtract one to prevent arcsin becoming infinity.
1226 // if(mSignalStorage[t_invMaxName].getRealInt() < t_actual) {
1227 // t_actual += 1*t_toReal;
1228 // mSignalStorage[t_invMaxName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
1229 // //cout<<"[Warning] "<<t_invMaxName<<" is too small. Adding one unit to prevent arcsin output being nan."<<endl;
1230 // //cout<<t_invMaxName<<".getRealInt():"<<mSignalStorage[t_invMaxName].getRealInt()<<" is new max."<<endl;
1231 // }
1232 //}
1233 for (unsigned iSt = 0; iSt < 4; iSt++) {
1234 string t_invMinName = "invArcSMin_" + to_string(iSt);
1235 //double t_actual = (mConstV.find("rr")->second)[2*iSt+1]/2;
1236 double t_actual = mSignalStorage["rho_c_" + to_string(iSt)].getMaxActual();
1237 double t_toReal = mSignalStorage["rho_c_" + to_string(iSt)].getToReal();
1238 mSignalStorage[t_invMinName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
1239 string t_invMaxName = "invArcSMax_" + to_string(iSt);
1240 t_actual = mSignalStorage["rho_c_" + to_string(iSt)].getMinActual();
1241 mSignalStorage[t_invMaxName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
1242 }
1243 }
1244 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<invArcSMin_"<<iSt<<">>>"<<endl; mSignalStorage["invArcSMin_"+to_string(iSt)].dump();}
1245 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<invArcSMax_"<<iSt<<">>>"<<endl; mSignalStorage["invArcSMax_"+to_string(iSt)].dump();}
1246
1247 // Generate LUT(arcS[i]=2*rho*asin(rr[i]/2/rho).
1248 for (unsigned iSt = 0; iSt < 4; iSt++) {
1249 string t_valueName = "arcS_" + to_string(iSt);
1250 string t_invMinName = "invArcSMin_" + to_string(iSt);
1251 string t_invMaxName = "invArcSMax_" + to_string(iSt);
1252 if (!mLutStorage.count(t_valueName)) { //if nothing found
1253 mLutStorage[t_valueName] = new Belle2::TRGCDCJLUT(t_valueName);
1254 // Lambda can not copy maps.
1255 double t_parameter = mConstV.at("rr3D")[iSt];
1256 mLutStorage[t_valueName]->setFloatFunction(
1257 [ = ](double aValue) -> double{return 2 * aValue * asin(t_parameter / 2 / aValue);},
1258 //mSignalStorage["rho"],
1259 mSignalStorage["rho_c_" + to_string(iSt)],
1260 mSignalStorage[t_invMinName], mSignalStorage[t_invMaxName], mSignalStorage["rho"].getToReal(),
1261 (int)mConstD.at("acosLUTInBitSize"), (int)mConstD.at("acosLUTOutBitSize"));
1262 //mLutStorage[t_valueName]->makeCOE("./LutData/"+t_valueName+".coe");
1263 }
1264 }
1265 // arcS[i]=2*rho*asin(rr[i]/2/rho).
1266 // Operate using LUT(arcS[i]).
1267 for (unsigned iSt = 0; iSt < 4; iSt++) {
1268 // Set output name
1269 string t_valueName = "arcS_" + to_string(iSt);
1270 //mLutStorage[t_valueName]->operate(mSignalStorage["rho"], mSignalStorage[t_valueName]);
1271 mLutStorage[t_valueName]->operate(mSignalStorage["rho_c_" + to_string(iSt)], mSignalStorage[t_valueName]);
1272 }
1273 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<arcS_"<<iSt<<">>>"<<endl; mSignalStorage["arcS_"+to_string(iSt)].dump();}
1274}

◆ calStAxPhi()

double calStAxPhi ( int  charge,
double  anglest,
double  ztostraw,
double  rr,
double  rho,
double  phi0 
)
static

Calculates the fitted axial phi for the stereo super layer.

Definition at line 785 of file Fitter3DUtility.cc.

786{
787 if (false) cout << anglest << ztostraw << endl; // Removes warnings when compileing
788
789 double myphiz, acos_real;
790 //Find phifit-phist
791 double t_rho = rho;
792 if (rr > (2 * rho)) t_rho = rr / 2;
793 acos_real = acos(rr / (2 * t_rho));
794 if (mysign == 1) {
795 myphiz = +acos_real + myphi0;
796 } else {
797 myphiz = -acos_real + myphi0;
798 }
799 if (myphiz > M_PI) myphiz -= 2 * M_PI;
800 if (myphiz < -M_PI) myphiz += 2 * M_PI;
801
802 return myphiz;
803}

◆ calVectorsAtR()

void calVectorsAtR ( const TVectorD &  helixParameters,
int  charge,
double  radius,
ROOT::Math::XYZVector &  position,
ROOT::Math::XYZVector &  momentum 
)
static

Calculates position and momentum at a certain radius.

Definition at line 1967 of file Fitter3DUtility.cc.

1969{
1970 // HelixParameters: dR, phi0, keppa, dz, tanLambda
1971 double t_alpha = 1 / 0.3 / 1.5 * 100;
1972 double t_R = t_alpha / helixParameters[2];
1973 double t_pT = t_R / t_alpha;
1974 double t_phi = -1 * charge * acos((-pow(cdcRadius, 2) + pow(t_R + helixParameters[0], 2) + pow(t_R,
1975 2)) / (2 * t_R * (t_R + helixParameters[0])));
1976 double t_X = (helixParameters[0] + t_R) * cos(helixParameters[1]) - t_R * cos(helixParameters[1] + t_phi);
1977 double t_Y = (helixParameters[0] + t_R) * sin(helixParameters[1]) - t_R * sin(helixParameters[1] + t_phi);
1978 double t_Z = helixParameters[3] - helixParameters[4] * t_R * t_phi;
1979 double t_Px = -t_pT * sin(helixParameters[1] + t_phi) * charge;
1980 double t_Py = t_pT * cos(helixParameters[1] + t_phi) * charge;
1981 double t_Pz = t_pT * helixParameters[4] * charge;
1982 position.SetXYZ(t_X, t_Y, t_Z);
1983 momentum.SetXYZ(t_Px, t_Py, t_Pz);
1984}

◆ calZ() [1/2]

double calZ ( int  charge,
double  anglest,
double  ztostraw,
double  rr,
double  phi,
double  rho,
double  phi0 
)
static

Calculates z.

Definition at line 892 of file Fitter3DUtility.cc.

893{
894 double myphiz = 0, acos_real = 0, dPhiAx = 0;
895 //Find phifit-phist
896 double t_rho = rho;
897 if (rr > (2 * rho)) t_rho = rr / 2;
898 acos_real = acos(rr / (2 * t_rho));
899 if (mysign == 1) {
900 dPhiAx = -acos_real - myphi0;
901 } else {
902 dPhiAx = acos_real - myphi0;
903 }
904 myphiz = dPhiAx + phi2;
905 if (myphiz > M_PI) myphiz -= 2 * M_PI;
906 if (myphiz < -M_PI) myphiz += 2 * M_PI;
907
908 return (ztostraw - rr * 2 * sin(myphiz / 2) / anglest);
909}

◆ calZ() [2/2]

void calZ ( std::map< std::string, double > const &  mConstD,
std::map< std::string, std::vector< double > > const &  mConstV,
std::map< std::string, Belle2::TRGCDCJSignal > &  mSignalStorage,
std::map< std::string, Belle2::TRGCDCJLUT * > &  mLutStorage 
)
static

Calculates z using JSignal and multiple constants.

Definition at line 911 of file Fitter3DUtility.cc.

913{
914 Belle2::TRGCDCJSignalData* commonData = mSignalStorage.begin()->second.getCommonData();
915 // Calculate zz[i]
916 if (mSignalStorage.find("invPhiAxMin_0") == mSignalStorage.end()) {
917 //for(unsigned iSt=0; iSt<4; iSt++) {
918 // string t_invMinName = "invPhiAxMin_" + to_string(iSt);
919 // double t_actual = (mConstV.find("rr")->second)[2*iSt+1]/2;
920 // double t_toReal = mSignalStorage["rho"].getToReal();
921 // mSignalStorage[t_invMinName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
922 // // Add one to prevent arccos becoming infinity.
923 // if(mSignalStorage[t_invMinName].getRealInt() < t_actual) {
924 // t_actual += 1*t_toReal;
925 // mSignalStorage[t_invMinName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
926 // //cout<<"[Warning] "<<t_invMinName<<" is too small. Adding one unit to prevent arccos output being nan."<<endl;
927 // //cout<<t_invMinName<<".getRealInt():"<<mSignalStorage[t_invMinName].getRealInt()<<" is new min."<<endl;
928 // }
929 // string t_invMaxName = "invPhiAxMax_" + to_string(iSt);
930 // t_actual = mSignalStorage["rho"].getMaxActual();
931 // mSignalStorage[t_invMaxName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
932 //}
933 for (unsigned iSt = 0; iSt < 4; iSt++) {
934 string t_invMinName = "invPhiAxMin_" + to_string(iSt);
935 //double t_actual = (mConstV.find("rr")->second)[2*iSt+1]/2;
936 double t_actual = mSignalStorage["rho_c_" + to_string(iSt)].getMinActual();
937 double t_toReal = mSignalStorage["rho_c_" + to_string(iSt)].getToReal();
938 mSignalStorage[t_invMinName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
939 string t_invMaxName = "invPhiAxMax_" + to_string(iSt);
940 t_actual = mSignalStorage["rho_c_" + to_string(iSt)].getMaxActual();
941 mSignalStorage[t_invMaxName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
942 }
943 }
944 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<rho_c_"<<iSt<<">>>"<<endl; mSignalStorage["rho_c_"+to_string(iSt)].dump();}
945 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<invPhiAxMin_"<<iSt<<">>>"<<endl; mSignalStorage["invPhiAxMin_"+to_string(iSt)].dump();}
946 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<invPhiAxMax_"<<iSt<<">>>"<<endl; mSignalStorage["invPhiAxMax_"+to_string(iSt)].dump();}
947 // Generate LUT(phiAx[i]=acos(rr[i]/2/rho)).
948 for (unsigned iSt = 0; iSt < 4; iSt++) {
949 string t_valueName = "phiAx_" + to_string(iSt);
950 string t_minName = "phiAxMin_" + to_string(iSt);
951 string t_maxName = "phiAxMax_" + to_string(iSt);
952 string t_invMinName = "invPhiAxMin_" + to_string(iSt);
953 string t_invMaxName = "invPhiAxMax_" + to_string(iSt);
954 if (!mLutStorage.count(t_valueName)) { //if nothing found
955 mLutStorage[t_valueName] = new Belle2::TRGCDCJLUT(t_valueName);
956 // Lambda can not copy maps.
957 double t_parameter = mConstV.at("rr3D")[iSt];
958 mLutStorage[t_valueName]->setFloatFunction(
959 [ = ](double aValue) -> double{return acos(t_parameter / 2 / aValue);},
960 //mSignalStorage["rho"],
961 mSignalStorage["rho_c_" + to_string(iSt)],
962 mSignalStorage[t_invMinName], mSignalStorage[t_invMaxName], mSignalStorage["phi0"].getToReal(),
963 (int)mConstD.at("acosLUTInBitSize"), (int)mConstD.at("acosLUTOutBitSize"));
964 //mLutStorage[t_valueName]->makeCOE("./LutData/"+t_valueName+".coe");
965 }
966 }
967 // phiAx[i] = acos(rr[i]/2/rho).
968 // Operate using LUT(phiAx[i]).
969 for (unsigned iSt = 0; iSt < 4; iSt++) {
970 // Set output name
971 string t_valueName = "phiAx_" + to_string(iSt);
972 //mLutStorage[t_valueName]->operate(mSignalStorage["rho"], mSignalStorage[t_valueName]);
973 mLutStorage[t_valueName]->operate(mSignalStorage["rho_c_" + to_string(iSt)], mSignalStorage[t_valueName]);
974 }
975 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<phiAx_"<<iSt<<">>>"<<endl; mSignalStorage["phiAx_"+to_string(iSt)].dump();}
976 //cout<<"<<<phi0>>>"<<endl; mSignalStorage["phi0"].dump();
977 //cout<<"<<<charge>>>"<<endl; mSignalStorage["charge"].dump();
978 // dPhiAx[i] = -+phiAx[i] - phi0
979 {
980 // Create data for ifElse
981 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
982 // Compare
983 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(Belle2::TRGCDCJSignal(1, mSignalStorage["charge"].getToReal(),
984 commonData), "=", mSignalStorage["charge"]);
985 // Assignments
986 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
987 make_pair(&mSignalStorage["dPhiAx_0"], -mSignalStorage["phiAx_0"] - mSignalStorage["phi0"]),
988 make_pair(&mSignalStorage["dPhiAx_1"], -mSignalStorage["phiAx_1"] - mSignalStorage["phi0"]),
989 make_pair(&mSignalStorage["dPhiAx_2"], -mSignalStorage["phiAx_2"] - mSignalStorage["phi0"]),
990 make_pair(&mSignalStorage["dPhiAx_3"], -mSignalStorage["phiAx_3"] - mSignalStorage["phi0"])
991 };
992 // Push to data.
993 t_data.push_back(make_pair(t_compare, t_assigns));
994 // Compare
995 t_compare = Belle2::TRGCDCJSignal();
996 // Assignments
997 t_assigns = {
998 make_pair(&mSignalStorage["dPhiAx_0"], mSignalStorage["phiAx_0"] - mSignalStorage["phi0"]),
999 make_pair(&mSignalStorage["dPhiAx_1"], mSignalStorage["phiAx_1"] - mSignalStorage["phi0"]),
1000 make_pair(&mSignalStorage["dPhiAx_2"], mSignalStorage["phiAx_2"] - mSignalStorage["phi0"]),
1001 make_pair(&mSignalStorage["dPhiAx_3"], mSignalStorage["phiAx_3"] - mSignalStorage["phi0"])
1002 };
1003 // Push to data.
1004 t_data.push_back(make_pair(t_compare, t_assigns));
1005 // Process ifElse data.
1006 Belle2::TRGCDCJSignal::ifElse(t_data);
1007 }
1008 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<dPhiAx_"<<iSt<<">>>"<<endl; mSignalStorage["dPhiAx_"+to_string(iSt)].dump();}
1009 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<phi_"<<iSt<<">>>"<<endl; mSignalStorage["phi_"+to_string(iSt)].dump();}
1010 // dPhi[i] = dPhiAx[i] + phi[i]
1011 for (unsigned iSt = 0; iSt < 4; iSt++) {
1012 mSignalStorage["dPhi_" + to_string(iSt)] <= mSignalStorage["dPhiAx_" + to_string(iSt)] + mSignalStorage["phi_" + to_string(iSt)];
1013 }
1014 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<dPhi_"<<iSt<<">>>"<<endl; mSignalStorage["dPhi_"+to_string(iSt)].dump();}
1015
1016 // Rotate dPhi to [-pi, pi]
1017 if (mSignalStorage.find("dPhiPiMax_0") == mSignalStorage.end()) {
1018 for (unsigned iSt = 0; iSt < 4; iSt++) {
1019 string t_valueName = "dPhi_" + to_string(iSt);
1020 string t_maxName = "dPhiPiMax_" + to_string(iSt);
1021 string t_minName = "dPhiPiMin_" + to_string(iSt);
1022 string t_2PiName = "dPhiPi2Pi_" + to_string(iSt);
1023 mSignalStorage[t_maxName] = Belle2::TRGCDCJSignal(mConstD.at("M_PI"), mSignalStorage[t_valueName].getToReal(), commonData);
1024 mSignalStorage[t_minName] = Belle2::TRGCDCJSignal(-mConstD.at("M_PI"), mSignalStorage[t_valueName].getToReal(), commonData);
1025 mSignalStorage[t_2PiName] = Belle2::TRGCDCJSignal(2 * mConstD.at("M_PI"), mSignalStorage[t_valueName].getToReal(), commonData);
1026 }
1027 }
1028 for (unsigned iSt = 0; iSt < 4; iSt++) {
1029 string t_in1Name = "dPhi_" + to_string(iSt);
1030 string t_valueName = "dPhiPi_c_" + to_string(iSt);
1031 string t_maxName = "dPhiPiMax_" + to_string(iSt);
1032 string t_minName = "dPhiPiMin_" + to_string(iSt);
1033 string t_2PiName = "dPhiPi2Pi_" + to_string(iSt);
1034 // Create data for ifElse
1035 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1036 // Compare (dPhi[i] > dPhiPiMax)
1037 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_maxName]);
1038 // Assignments (dPhiPi_c[i] = dPhi[i] - dPhiPi2Pi[i])
1039 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1040 make_pair(&mSignalStorage[t_valueName], (mSignalStorage[t_in1Name] - mSignalStorage[t_2PiName]).limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
1041 };
1042 // Push to data.
1043 t_data.push_back(make_pair(t_compare, t_assigns));
1044 // Compare (dPhi[i] > dPhiPiMin)
1045 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_minName]);
1046 // Assignments (dPhiPi_c[i] = dPhi[i])
1047 t_assigns = {
1048 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_in1Name].limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
1049 };
1050 // Push to data.
1051 t_data.push_back(make_pair(t_compare, t_assigns));
1052 // Compare (else)
1053 t_compare = Belle2::TRGCDCJSignal();
1054 // Assignments (dPhiPi_c[i] = dPhi[i] + dPhiPi2Pi[i])
1055 t_assigns = {
1056 make_pair(&mSignalStorage[t_valueName], (mSignalStorage[t_in1Name] + mSignalStorage[t_2PiName]).limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
1057 };
1058 // Push to data.
1059 t_data.push_back(make_pair(t_compare, t_assigns));
1060 // Process ifElse data.
1061 Belle2::TRGCDCJSignal::ifElse(t_data);
1062 }
1063 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<dPhiPi_c_"<<iSt<<">>>"<<endl; mSignalStorage["dPhiPi_c_"+to_string(iSt)].dump();}
1064
1065 // Calculate dPhiMax[i], dPhiMin[i]=0
1066 if (mSignalStorage.find("dPhiMax_0") == mSignalStorage.end()) {
1067 for (unsigned iSt = 0; iSt < 4; iSt++) {
1068 string t_maxName = "dPhiMax_" + to_string(iSt);
1069 string t_minName = "dPhiMin_" + to_string(iSt);
1070 string t_valueName = "dPhi_" + to_string(iSt);
1071 double t_value = 2 * mConstD.at("M_PI") * fabs(mConstV.at("nShift")[iSt]) / (mConstV.at("nWires")[2 * iSt + 1]);
1072 mSignalStorage[t_maxName] = Belle2::TRGCDCJSignal(t_value, mSignalStorage[t_valueName].getToReal(), commonData);
1073 mSignalStorage[t_minName] = Belle2::TRGCDCJSignal(0, mSignalStorage[t_valueName].getToReal(), commonData);
1074 }
1075 }
1076 // Constrain dPhiPi to [0, 2*pi*#shift/#holes]
1077 for (unsigned iSt = 0; iSt < 4; iSt++) {
1078 string t_maxName = "dPhiMax_" + to_string(iSt);
1079 string t_minName = "dPhiMin_" + to_string(iSt);
1080 string t_valueName = "dPhi_c_" + to_string(iSt);
1081 string t_in1Name = "dPhiPi_c_" + to_string(iSt);
1082 // For SL1, SL5
1083 if (iSt % 2 == 0) {
1084 // Create data for ifElse
1085 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1086 // Compare
1087 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", Belle2::TRGCDCJSignal(0,
1088 mSignalStorage[t_in1Name].getToReal(), commonData));
1089 // Assignments
1090 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1091 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_minName]),
1092 };
1093 // Push to data.
1094 t_data.push_back(make_pair(t_compare, t_assigns));
1095 // Compare
1096 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", -mSignalStorage[t_maxName]);
1097 // Assignments
1098 t_assigns = {
1099 make_pair(&mSignalStorage[t_valueName], Belle2::TRGCDCJSignal::absolute(mSignalStorage[t_in1Name]).limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
1100 };
1101 // Push to data.
1102 t_data.push_back(make_pair(t_compare, t_assigns));
1103 // Compare
1104 t_compare = Belle2::TRGCDCJSignal();
1105 // Assignments
1106 t_assigns = {
1107 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_maxName]),
1108 };
1109 // Push to data.
1110 t_data.push_back(make_pair(t_compare, t_assigns));
1111 // Process ifElse data.
1112 Belle2::TRGCDCJSignal::ifElse(t_data);
1113 // For SL3, SL7
1114 } else {
1115 // Create data for ifElse
1116 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1117 // Compare (dPhi > dPhiMax)
1118 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_maxName]);
1119 // Assignments (dPhi_c = dPhiMax)
1120 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1121 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_maxName]),
1122 };
1123 // Push to data.
1124 t_data.push_back(make_pair(t_compare, t_assigns));
1125 // Compare (dPhi > 0)
1126 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", Belle2::TRGCDCJSignal(0,
1127 mSignalStorage[t_in1Name].getToReal(), commonData));
1128 // Assignments (dPhi_c = dPhi)
1129 t_assigns = {
1130 make_pair(&mSignalStorage[t_valueName], Belle2::TRGCDCJSignal::absolute(mSignalStorage[t_in1Name]).limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
1131 };
1132 // Push to data.
1133 t_data.push_back(make_pair(t_compare, t_assigns));
1134 // Compare (else) => (dPhi > -Pi)
1135 t_compare = Belle2::TRGCDCJSignal();
1136 // Assignments (dPhi_c = dPhiMin)
1137 t_assigns = {
1138 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_minName]),
1139 };
1140 // Push to data.
1141 t_data.push_back(make_pair(t_compare, t_assigns));
1142 // Process ifElse.
1143 Belle2::TRGCDCJSignal::ifElse(t_data);
1144 }
1145 }
1146
1147 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<dPhiMax_"<<iSt<<">>>"<<endl; mSignalStorage["dPhiMax_"+to_string(iSt)].dump();}
1148 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<dPhiMin_"<<iSt<<">>>"<<endl; mSignalStorage["dPhiMin_"+to_string(iSt)].dump();}
1149 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<dPhi_c_"<<iSt<<">>>"<<endl; mSignalStorage["dPhi_c_"+to_string(iSt)].dump();}
1150 // Calculate minInvValue, maxInvValue for LUTs
1151 if (mSignalStorage.find("invZMin_0") == mSignalStorage.end()) {
1152 for (unsigned iSt = 0; iSt < 4; iSt++) {
1153 string t_minName = "invZMin_" + to_string(iSt);
1154 string t_maxName = "invZMax_" + to_string(iSt);
1155 string t_valueName = "dPhi_c_" + to_string(iSt);
1156 unsigned long long t_int = mSignalStorage[t_valueName].getMinInt();
1157 double t_toReal = mSignalStorage[t_valueName].getToReal();
1158 double t_actual = mSignalStorage[t_valueName].getMinActual();
1159 mSignalStorage[t_minName] = Belle2::TRGCDCJSignal(t_int, t_toReal, t_int, t_int, t_actual, t_actual, t_actual, -1, commonData);
1160 t_int = mSignalStorage[t_valueName].getMaxInt();
1161 t_actual = mSignalStorage[t_valueName].getMaxActual();
1162 mSignalStorage[t_maxName] = Belle2::TRGCDCJSignal(t_int, t_toReal, t_int, t_int, t_actual, t_actual, t_actual, -1, commonData);
1163 }
1164 }
1165 // Generate LUT(zz[i] = ztostraw[0]-+rr[0]*2*sin(dPhi_c[i]/2)/angleSt[i], -+ depends on SL)
1166 for (unsigned iSt = 0; iSt < 4; iSt++) {
1167 string t_inputName = "dPhi_c_" + to_string(iSt);
1168 string t_outputName = "zz_" + to_string(iSt);
1169 string t_invMinName = "invZMin_" + to_string(iSt);
1170 string t_invMaxName = "invZMax_" + to_string(iSt);
1171
1172 if (!mLutStorage.count(t_outputName)) { //if nothing found
1173 mLutStorage[t_outputName] = new Belle2::TRGCDCJLUT(t_outputName);
1174 // Lambda can not copy maps.
1175 double t_zToStraw = mConstV.at("zToStraw")[iSt];
1176 double t_rr3D = mConstV.at("rr3D")[iSt];
1177 double t_angleSt = mConstV.at("angleSt")[iSt];
1178
1179 if (iSt % 2 == 0) {
1180 mLutStorage[t_outputName]->setFloatFunction(
1181 [ = ](double aValue) -> double{return t_zToStraw + t_rr3D * 2 * sin(aValue / 2) / t_angleSt;},
1182 mSignalStorage[t_inputName],
1183 mSignalStorage[t_invMinName], mSignalStorage[t_invMaxName], mSignalStorage["rho"].getToReal(),
1184 (int) mConstD.at("zLUTInBitSize"), (int) mConstD.at("zLUTOutBitSize"));
1185 } else {
1186 mLutStorage[t_outputName]->setFloatFunction(
1187 [ = ](double aValue) -> double{return t_zToStraw - t_rr3D * 2 * sin(aValue / 2) / t_angleSt;},
1188 mSignalStorage[t_inputName],
1189 mSignalStorage[t_invMinName], mSignalStorage[t_invMaxName], mSignalStorage["rho"].getToReal(),
1190 (int) mConstD.at("zLUTInBitSize"), (int) mConstD.at("zLUTOutBitSize"));
1191 }
1192 //mLutStorage[t_outputName]->makeCOE("./LutData/"+t_outputName+".coe");
1193 }
1194 }
1195 // zz[i] = ztostraw[0]-+rr[0]*2*sin(dPhi_c[i]/2)/angleSt[i], -+ depends on SL)
1196 // Operate using LUT(zz[i]).
1197 for (unsigned iSt = 0; iSt < 4; iSt++) {
1198 string t_outputName = "zz_" + to_string(iSt);
1199 string t_inputName = "dPhi_c_" + to_string(iSt);
1200 mLutStorage[t_outputName]->operate(mSignalStorage[t_inputName], mSignalStorage[t_outputName]);
1201 }
1202 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<zz_"<<iSt<<">>>"<<endl; mSignalStorage["zz_"+to_string(iSt)].dump();}
1203}
Belle2::TRGCDCJSignal::absolute
static TRGCDCJSignal const absolute(TRGCDCJSignal const &first)
Absolute TRGCDCJSignal. Removes 1 bit if signed or minus unsigned.
Definition: JSignal.cc:1562

◆ changeInteger()

void changeInteger ( int &  integer,
double  real,
double  minValue,
double  maxValue,
int  bitSize 
)
static

Changes float to integer value.

Definition at line 2002 of file Fitter3DUtility.cc.

2003{
2004 // Version 1.
2005 //double range = maxValue - minValue;
2006 //double convert = bitSize/range;
2008 //integer = int( (real - minValue) * convert +0.5);
2009 // Version 2.
2010 double range = maxValue - minValue;
2011 double convert = (bitSize - 0.5) / range;
2012 integer = FpgaUtility::roundInt((real - minValue) * convert);
2013}
Fitter3DUtility::bitSize
static int bitSize(int numberBits, int mode)
Firmware convert functions.
Definition: Fitter3DUtility.cc:1986
FpgaUtility::roundInt
static double roundInt(double value)
Round double value.
Definition: FpgaUtility.cc:38

◆ changeReal()

void changeReal ( double &  real,
int  integer,
double  minValue,
double  maxValue,
int  bitSize 
)
static

Changes integer to float value.

Definition at line 2015 of file Fitter3DUtility.cc.

2016{
2018 //double range = maxValue - minValue;
2019 //double convert = bitSize/range;
2020 //real = (integer/convert) + minValue;
2021 // Version 2.
2022 double range = maxValue - minValue;
2023 double convert = (bitSize - 0.5) / range;
2024 real = (integer / convert) + minValue;
2025}

◆ chargeFinder()

void chargeFinder ( double *  nTSs,
double *  tsIds,
double *  tsHitmap,
double  phi0,
double  inCharge,
double &  outCharge 
)
static

Charge finder using circle fitter output and axial TSs.

Definition at line 142 of file Fitter3DUtility.cc.

144{
145 vector<double> tsPhi(5);
146 vector<double> dPhi(5);
147 vector<double> dPhi_c(5);
148 vector<double> charge(5);
149 for (unsigned iAx = 0; iAx < 5; iAx++) {
150 //cout<<"iAx:"<<iAx<<" nTSs:"<<nTSs[iAx]<<" tsId:"<<tsIds[iAx]<<" hitmap:"<<tsHitmap[iAx]<<" phi0:"<<phi0<<" inCharge:"<<inCharge<<endl;
151 // Change tsId to rad unit.
152 tsPhi[iAx] = tsIds[iAx] * 2 * M_PI / nTSs[iAx];
153 dPhi[iAx] = tsPhi[iAx] - phi0;
154 //cout<<"iAx:"<<iAx<<" dPhi:"<<dPhi[iAx]<<endl;
155 // Constrain dPhi to [-pi, pi]
156 if (dPhi[iAx] < -M_PI) dPhi_c[iAx] = dPhi[iAx] + 2 * M_PI;
157 else if (dPhi[iAx] > M_PI) dPhi_c[iAx] = dPhi[iAx] - 2 * M_PI;
158 else dPhi_c[iAx] = dPhi[iAx];
159 //cout<<"iAx:"<<iAx<<" dPhi_c:"<<dPhi_c[iAx]<<endl;
160 // Calculate charge
161 if (tsHitmap[iAx] == 0) charge[iAx] = 0;
162 else if (dPhi_c[iAx] == 0) charge[iAx] = 0;
163 else if (dPhi_c[iAx] > 0) charge[iAx] = 1;
164 else charge[iAx] = -1;
165 //cout<<"iAx:"<<iAx<<" charge:"<<charge[iAx]<<endl;
166 }
167 // Sum up charge
168 double sumCharge = charge[0] + charge[1] + charge[2] + charge[3] + charge[4];
169 //cout<<"sumCharge:"<<sumCharge<<endl;
170 // Calculate foundCharge
171 if (sumCharge == 0) foundCharge = inCharge;
172 else if (sumCharge > 0) foundCharge = 1;
173 else foundCharge = -1;
174 //cout<<"foundCharge:"<<foundCharge<<endl;
175 //cout<<"inCharge:"<<inCharge<<endl;
176}

◆ constrainRPerStSl()

void constrainRPerStSl ( std::map< std::string, std::vector< double > > const &  mConstV,
std::map< std::string, Belle2::TRGCDCJSignal > &  mSignalStorage 
)
static

Constrains R for each SL differently using JSignal and multiple constants.

Definition at line 825 of file Fitter3DUtility.cc.

827{
828 Belle2::TRGCDCJSignalData* commonData = mSignalStorage.begin()->second.getCommonData();
829 // Make constants for min and max values for each super layer.
830 if (mSignalStorage.find("rhoMin_0") == mSignalStorage.end()) {
831 for (unsigned iSt = 0; iSt < 4; iSt++) {
832 // Min value
833 string t_minName = "rhoMin_" + to_string(iSt);
834 double t_actual = (mConstV.find("rr")->second)[2 * iSt + 1] / 2;
835 double t_toReal = mSignalStorage["rho"].getToReal();
836 mSignalStorage[t_minName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
837 // Add one to prevent arccos, arcsin becoming infinity.
838 if (mSignalStorage[t_minName].getRealInt() < t_actual) {
839 t_actual += 1 * t_toReal;
840 mSignalStorage[t_minName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
841 //cout<<"[Warning] "<<t_minName<<" is too small. Adding one unit to prevent arccos output being nan."<<endl;
842 //cout<<t_minName<<".getRealInt():"<<mSignalStorage[t_minName].getRealInt()<<" is new min."<<endl;
843 }
844 // Max value
845 string t_maxName = "rhoMax_" + to_string(iSt);
846 t_actual = mSignalStorage["rho"].getMaxActual();
847 mSignalStorage[t_maxName] = Belle2::TRGCDCJSignal(t_actual, t_toReal, commonData);
848 }
849 }
850 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<rhoMin_"<<iSt<<">>>"<<endl; mSignalStorage["rhoMin_"+to_string(iSt)].dump();}
851 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<rhoMax_"<<iSt<<">>>"<<endl; mSignalStorage["rhoMax_"+to_string(iSt)].dump();}
852
853 //cout<<"<<<rho>>>"<<endl; mSignalStorage["rho"].dump();
854 // Constrain rho to min and max per layer.
855 for (unsigned iSt = 0; iSt < 4; iSt++) {
856 string t_in1Name = "rho";
857 string t_valueName = "rho_c_" + to_string(iSt);
858 string t_maxName = "rhoMax_" + to_string(iSt);
859 string t_minName = "rhoMin_" + to_string(iSt);
860 // Create data for ifElse
861 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
862 // Compare
863 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_maxName]);
864 // Assignments
865 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
866 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_maxName]),
867 };
868 // Push to data.
869 t_data.push_back(make_pair(t_compare, t_assigns));
870 // Compare
871 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_minName]);
872 // Assignments
873 t_assigns = {
874 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_in1Name].limit(mSignalStorage[t_minName], mSignalStorage[t_maxName]))
875 };
876 // Push to data.
877 t_data.push_back(make_pair(t_compare, t_assigns));
878 // Compare
879 t_compare = Belle2::TRGCDCJSignal();
880 // Assignments
881 t_assigns = {
882 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_minName])
883 };
884 // Push to data.
885 t_data.push_back(make_pair(t_compare, t_assigns));
886 Belle2::TRGCDCJSignal::ifElse(t_data);
887 }
888 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<rho_c_"<<iSt<<">>>"<<endl; mSignalStorage["rho_c_"+to_string(iSt)].dump();}
889
890}

◆ findExtreme()

void findExtreme ( double &  m_max,
double &  m_min,
double  value 
)
static

Finds maximum and minium values.

Definition at line 2028 of file Fitter3DUtility.cc.

2029{
2030 if (value > m_max) m_max = value;
2031 if (value < m_min) m_min = value;
2032}

◆ findImpactPosition()

void findImpactPosition ( ROOT::Math::XYZVector *  mcPosition,
ROOT::Math::PxPyPzEVector *  mcMomentum,
int  charge,
ROOT::Math::XYVector &  helixCenter,
ROOT::Math::XYZVector &  impactPosition 
)
static

MC calculation functions Calculates the impact position of track.

Definition at line 1922 of file Fitter3DUtility.cc.

1925{
1926
1927 // Finds the impact position. Everything is in cm, and GeV.
1928 // Input: production vertex (mcx, mcy, mcz),
1929 // momentum at production vertex (px, py, pz)
1930 // charge of particle.
1931 // Output: helix center's coordiante (hcx, hcy)
1932 // impact position (impactX, impactY, impactZ)
1933
1934 double rho = sqrt(pow(mcMomentum->X(), 2) + pow(mcMomentum->Y(), 2)) / 0.3 / 1.5 * 100;
1935 double hcx = mcPosition->X() + rho * cos(atan2(mcMomentum->Y(), mcMomentum->X()) - charge * M_PI_2);
1936 double hcy = mcPosition->Y() + rho * sin(atan2(mcMomentum->Y(), mcMomentum->X()) - charge * M_PI_2);
1937 helixCenter.SetXY(hcx, hcy);
1938 double impactX = (helixCenter.R() - rho) / helixCenter.R() * helixCenter.X();
1939 double impactY = (helixCenter.R() - rho) / helixCenter.R() * helixCenter.Y();
1940 int signdS;
1941 if (atan2(impactY, impactX) < atan2(mcPosition->Y(), mcPosition->X())) signdS = -1;
1942 else signdS = 1;
1943 double dS = 2 * rho * asin(sqrt(pow(impactX - mcPosition->X(), 2) + pow(impactY - mcPosition->Y(), 2)) / 2 / rho);
1944 double impactZ = mcMomentum->Pz() / mcMomentum->Pt() * dS * signdS + mcPosition->Z();
1945 impactPosition.SetXYZ(impactX, impactY, impactZ);
1946
1947}
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28

◆ findQuadrant()

int findQuadrant ( double  value)
static

Finds quadrant of angle. Angle is in rad.

Definition at line 641 of file Fitter3DUtility.cc.

642{
643 // Rotate to [-pi,pi] range.
644 double phiMin = -M_PI;
645 double phiMax = M_PI;
646 bool rangeOK = 0;
647 while (rangeOK == 0) {
648 if (value > phiMax) value -= 2 * M_PI;
649 else if (value < phiMin) value += 2 * M_PI;
650 else rangeOK = 1;
651 }
652 // Find quadrant.
653 int result = -1;
654 if (value > M_PI_2) result = 2;
655 else if (value > 0) result = 1;
656 else if (value > -M_PI_2) result = 4;
657 else if (value > -M_PI) result = 3;
658 return result;
659}

◆ findSign()

int findSign ( double *  phi2)
static

2D fitter functions

Definition at line 40 of file Fitter3DUtility.cc.

41{
42 int mysign;
43 double sign_phi[2];
44
45 if ((phi2[0] - phi2[4]) > M_PI || (phi2[0] - phi2[4]) < -M_PI) {
46 if (phi2[0] > M_PI) {sign_phi[0] = phi2[0] - 2 * M_PI;}
47 else {sign_phi[0] = phi2[0];}
48 if (phi2[4] > M_PI) {sign_phi[1] = phi2[4] - 2 * M_PI;}
49 else {sign_phi[1] = phi2[4];}
50 } else {
51 sign_phi[0] = phi2[0];
52 sign_phi[1] = phi2[4];
53 }
54 if ((sign_phi[1] - sign_phi[0]) > 0) {mysign = 0;}
55 else {mysign = 1;}
56
57 return mysign;
58}

◆ fitter3D() [1/2]

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 
)
static

Combines several functions for fitter3D.

Definition at line 1842 of file Fitter3DUtility.cc.

1847{
1848 vector<double> zz;
1849 vector<double> arcS;
1850 vector<double> invZError2;
1851 fitter3D(stGeometry, stXts, eventTimeValid, eventTime, rawStTSs, charge, radius, phi_c, z0, cot, chi2, arcS, zz, invZError2);
1852}
Fitter3DUtility::fitter3D
static void fitter3D(std::map< std::string, std::vector< double > > &stGeometry, std::vector< std::vector< double > > const &stXts, int eventTimeValid, int eventTime, std::vector< std::vector< int > > const &rawStTSs, int charge, double radius, double phi_c, double &z0, double &cot, double &chi2)
Combines several functions for fitter3D.
Definition: Fitter3DUtility.cc:1842

◆ fitter3D() [2/2]

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,
std::vector< double > &  arcS,
std::vector< double > &  zz,
std::vector< double > &  invZError2 
)
static

Combines several functions for fitter3D. Also outputs arcS, zz and invZError2.

Definition at line 1854 of file Fitter3DUtility.cc.

1860{
1861 vector<double> stTSs(4);
1862 Fitter3DUtility::calPhiFast(stGeometry, stXts, eventTimeValid, eventTime, rawStTSs, stTSs);
1863 Fitter3DUtility::setErrorFast(rawStTSs, eventTimeValid, invZError2);
1864 zz = vector<double> (4, 0);
1865 arcS = vector<double> (4, 0);
1866 for (unsigned iSt = 0; iSt < 4; iSt++) {
1867 if (rawStTSs[iSt][1] != 0) {
1868 zz[iSt] = Fitter3DUtility::calZ(charge, stGeometry["angleSt"][iSt], stGeometry["zToStraw"][iSt],
1869 stGeometry["cdcRadius"][iSt], stTSs[iSt], radius, phi_c);
1870 arcS[iSt] = Fitter3DUtility::calS(radius, stGeometry["cdcRadius"][iSt]);
1871 }
1872 }
1873 Fitter3DUtility::rSFit(&invZError2[0], &arcS[0], &zz[0], z0, cot, chi2);
1874}
Fitter3DUtility::calZ
static double calZ(int charge, double anglest, double ztostraw, double rr, double phi, double rho, double phi0)
Calculates z.
Definition: Fitter3DUtility.cc:892
Fitter3DUtility::setErrorFast
static void setErrorFast(std::vector< std::vector< int > > const &rawStTSs, int eventTimeValid, std::vector< double > &invZError2)
Sets error for fitting.
Definition: Fitter3DUtility.cc:663
Fitter3DUtility::calPhiFast
static void calPhiFast(std::map< std::string, std::vector< double > > &stGeometry, std::vector< std::vector< double > > const &stXts, int eventTimeValid, int eventTime, std::vector< std::vector< int > > const &rawStTSs, std::vector< double > &stTSs)
Calculates phi with fast simulation.
Definition: Fitter3DUtility.cc:583
Fitter3DUtility::calS
static double calS(double rho, double rr)
Calculates arc length.
Definition: Fitter3DUtility.cc:1204
Fitter3DUtility::rSFit
static void rSFit(double *iezz2, double *arcS, double *zz, double &z0, double &cot, double &zchi2)
3D fitter functions Fits z and arc S.
Definition: Fitter3DUtility.cc:1303

◆ fitter3DFirm()

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 
)
static

Combines several functions for fitter3D firmware.

Definition at line 1876 of file Fitter3DUtility.cc.

1883{
1884 // Change to Signals.
1885 // rawStTSs[iSt] = [TS ID, TS LR, TS driftTime]
1886 vector<tuple<string, double, int, double, double, int> > t_values = {
1887 make_tuple("phi0", phi_c, mConstD["phiBitSize"], mConstD["phiMin"], mConstD["phiMax"], 0),
1888 make_tuple("rho", radius, mConstD["rhoBitSize"], mConstD["rhoMin"], mConstD["rhoMax"], 0),
1889 make_tuple("charge", (int)(charge == 1 ? 1 : 0), 1, 0, 1.5, 0),
1890 make_tuple("lr_0", rawStTSs[0][1], 2, 0, 3.5, 0),
1891 make_tuple("lr_1", rawStTSs[1][1], 2, 0, 3.5, 0),
1892 make_tuple("lr_2", rawStTSs[2][1], 2, 0, 3.5, 0),
1893 make_tuple("lr_3", rawStTSs[3][1], 2, 0, 3.5, 0),
1894 make_tuple("tsId_0", rawStTSs[0][0], ceil(log(mConstV.at("nTSs")[1]) / log(2)), 0, pow(2, ceil(log(mConstV.at("nTSs")[1]) / log(2))) - 0.5, 0),
1895 make_tuple("tsId_1", rawStTSs[1][0], ceil(log(mConstV.at("nTSs")[3]) / log(2)), 0, pow(2, ceil(log(mConstV.at("nTSs")[3]) / log(2))) - 0.5, 0),
1896 make_tuple("tsId_2", rawStTSs[2][0], ceil(log(mConstV.at("nTSs")[5]) / log(2)), 0, pow(2, ceil(log(mConstV.at("nTSs")[5]) / log(2))) - 0.5, 0),
1897 make_tuple("tsId_3", rawStTSs[3][0], ceil(log(mConstV.at("nTSs")[7]) / log(2)), 0, pow(2, ceil(log(mConstV.at("nTSs")[7]) / log(2))) - 0.5, 0),
1898 make_tuple("tdc_0", rawStTSs[0][2], mConstD["tdcBitSize"], 0, pow(2, mConstD["tdcBitSize"]) - 0.5, 0),
1899 make_tuple("tdc_1", rawStTSs[1][2], mConstD["tdcBitSize"], 0, pow(2, mConstD["tdcBitSize"]) - 0.5, 0),
1900 make_tuple("tdc_2", rawStTSs[2][2], mConstD["tdcBitSize"], 0, pow(2, mConstD["tdcBitSize"]) - 0.5, 0),
1901 make_tuple("tdc_3", rawStTSs[3][2], mConstD["tdcBitSize"], 0, pow(2, mConstD["tdcBitSize"]) - 0.5, 0),
1902 make_tuple("eventTime", eventTime, mConstD["tdcBitSize"], 0, pow(2, mConstD["tdcBitSize"]) - 0.5, 0),
1903 make_tuple("eventTimeValid", eventTimeValid, 1, 0, 1.5, 0),
1904 };
1905 Belle2::TRGCDCJSignal::valuesToMapSignals(t_values, commonData, mSignalStorage);
1906 // Calculate complex fit3D
1907 Fitter3DUtility::setError(mConstD, mConstV, mSignalStorage);
1908 Fitter3DUtility::calPhi(mConstD, mConstV, mSignalStorage, mLutStorage);
1909 Fitter3DUtility::constrainRPerStSl(mConstV, mSignalStorage);
1910 Fitter3DUtility::calZ(mConstD, mConstV, mSignalStorage, mLutStorage);
1911 Fitter3DUtility::calS(mConstD, mConstV, mSignalStorage, mLutStorage);
1912 Fitter3DUtility::rSFit(mConstD, mConstV, mSignalStorage, mLutStorage);
1913 // Post handling of signals.
1914 // Name all signals.
1915 if ((*mSignalStorage.begin()).second.getName() == "") {
1916 for (auto it = mSignalStorage.begin(); it != mSignalStorage.end(); ++it) {
1917 (*it).second.setName((*it).first);
1918 }
1919 }
1920}
Fitter3DUtility::setError
static void setError(std::map< std::string, double > const &mConstD, std::map< std::string, std::vector< double > > const &mConstV, std::map< std::string, Belle2::TRGCDCJSignal > &mSignalStorage)
Sets error using JSignal class.
Definition: Fitter3DUtility.cc:681
Fitter3DUtility::constrainRPerStSl
static void constrainRPerStSl(std::map< std::string, std::vector< double > > const &mConstV, std::map< std::string, Belle2::TRGCDCJSignal > &mSignalStorage)
Constrains R for each SL differently using JSignal and multiple constants.
Definition: Fitter3DUtility.cc:825
Belle2::TRGCDCJSignal::valuesToMapSignals
static void valuesToMapSignals(std::vector< std::tuple< std::string, double, int, double, double, int > > const &inValues, Belle2::TRGCDCJSignalData *inCommonData, std::map< std::string, Belle2::TRGCDCJSignal > &outMap)
Values => [name, value, bitwidth, min, max, clock] Changes values to signals.
Definition: JSignal.cc:2208

◆ loadStereoXt()

void loadStereoXt ( std::string const &  filePrefix,
int  nFiles,
std::vector< std::vector< double > > &  stXts 
)
static

Load stereo Xt file.

Definition at line 557 of file Fitter3DUtility.cc.

558{
559 stXts = vector<vector<double> > (nFiles);
560 for (unsigned iSt = 0; iSt < stXts.size(); ++iSt) {
561 string filename = filePrefix + ".st" + to_string(iSt);
562 ifstream xtFile;
563 xtFile.open(filename.c_str());
564 if (xtFile.fail()) {
565 cout << "Fitter3DUtility::loadStereoXt() Cannot load " << filename << endl;
566 return;
567 }
568 string line;
569 while (getline(xtFile, line)) {
570 stXts[iSt].push_back(std::stof(line));
571 }
572 xtFile.close();
573 }
574}

◆ rotatePhi() [1/2]

double rotatePhi ( double  value,
double  refPhi 
)
static

Rotates to range [-pi, pi].

Definition at line 609 of file Fitter3DUtility.cc.

610{
611 double phiMin = -M_PI;
612 double phiMax = M_PI;
613 double result = value - refPhi;
614 bool rangeOK = 0;
615 while (rangeOK == 0) {
616 if (result > phiMax) result -= 2 * M_PI;
617 else if (result < phiMin) result += 2 * M_PI;
618 else rangeOK = 1;
619 }
620 return result;
621}

◆ rotatePhi() [2/2]

double rotatePhi ( double  value,
int  refId,
int  nTSs 
)
static

Rotates to range [-pi, pi]. Use tsId as reference.

Definition at line 623 of file Fitter3DUtility.cc.

624{
625 double refPhi = (double)refId / nTSs * 2 * M_PI;
626 return rotatePhi(value, refPhi);
627}
Fitter3DUtility::rotatePhi
static double rotatePhi(double value, double refPhi)
Rotates to range [-pi, pi].
Definition: Fitter3DUtility.cc:609

◆ rotateTsId()

int rotateTsId ( int  value,
int  refId,
int  nTSs 
)
static

Rotates to range [0, nTSs-1].

Definition at line 629 of file Fitter3DUtility.cc.

630{
631 int result = value - refId;
632 bool rangeOk = 0;
633 while (rangeOk == 0) {
634 if (result >= nTSs) result -= nTSs;
635 else if (result < 0) result += nTSs;
636 else rangeOk = 1;
637 }
638 return result;
639}

◆ rPhiFit()

void rPhiFit ( double *  rr,
double *  phi2,
double *  phierror,
double &  rho,
double &  myphi0 
)
static

A circle fitter.

Definition at line 60 of file Fitter3DUtility.cc.

61{
62
63 cout << "Fitter3DUtility::rPhiFit() will be deprecated. Please use Fitter3DUtility::rPhiFitter(). phierror was changed to inv phi error."
64 << endl;
65 double invphierror[5];
66 for (unsigned i = 0; i < 5; i++) {
67 invphierror[i] = 1 / phierror[i];
68 }
69 rPhiFitter(rr, phi2, invphierror, rho, myphi0);
70}
Fitter3DUtility::rPhiFitter
static void rPhiFitter(double *rr, double *phi2, double *invphierror, double &rho, double &myphi0)
A circle fitter with invPhiError without fit chi2 output.
Definition: Fitter3DUtility.cc:72

◆ rPhiFit2()

void rPhiFit2 ( double *  rr,
double *  phi2,
double *  phierror,
double &  rho,
double &  myphi0,
int  nTS 
)
static

A circle fitter.

Definition at line 179 of file Fitter3DUtility.cc.

180{
181
182 double A, B, C, D, E, hcx, hcy;
183 //double G;
184 double fiterror[5];
185 //Calculate fit error
186 for (int i = 0; i < nTS; i++) {
187 //fiterror[i]=sqrt((rr[4]*rr[4]-2*rr[4]*rr[2]*cos(phi2[4]-phi2[2])+rr[2]*rr[2])/(sin(phi2[4]-phi2[2])*sin(phi2[4]-phi2[2]))-rr[i]*rr[i])*phierror[i];
188 fiterror[i] = 1 + 0 * phierror[i];
189 }
190
191 //r-phi fitter(2D Fitter) ->calculate pt and radius of track-> input for 3D fitter.
192 A = 0, B = 0, C = 0, D = 0, E = 0;
193 //G=0;
194 for (int i = 0; i < nTS; i++) {
195 A += cos(phi2[i]) * cos(phi2[i]) / (fiterror[i] * fiterror[i]);
196 B += sin(phi2[i]) * sin(phi2[i]) / (fiterror[i] * fiterror[i]);
197 C += cos(phi2[i]) * sin(phi2[i]) / (fiterror[i] * fiterror[i]);
198 D += rr[i] * cos(phi2[i]) / (fiterror[i] * fiterror[i]);
199 E += rr[i] * sin(phi2[i]) / (fiterror[i] * fiterror[i]);
200 //G+=rr[i]*rr[i]/(fiterror[i]*fiterror[i]);
201 }
202 hcx = D * B - E * C; //helix center x
203 hcx /= 2 * (A * B - C * C);
204 hcy = E * A - D * C; //helix center y
205 hcy /= 2 * (A * B - C * C);
206 rho = sqrt(hcx * hcx + hcy * hcy); //radius of helix
207 myphi0 = atan2(hcy, hcx);
208 if (myphi0 < 0) myphi0 += 2 * M_PI;
209 //myphi0=atan(hcy/hcx);
210 //if(hcx<0 && hcy>0) myphi0 += M_PI;
211 //if(hcx<0 && hcy<0) myphi0 += M_PI;
212 //if(hcx>0 && hcy<0) myphi0 += M_PI*2.0;
213
215 // double pchi2 = -2*hcx*D-2*hcy*E+G;
216 // pchi2/=nTS-2;
217
219 //double pchi3 = 0.; //iw
220 //for(int i=0;i<nTS;i++){
221 // pchi3+=(2*(hcx*cos(phi2[i])+hcy*sin(phi2[i]))-rr[i])*(2*(hcx*cos(phi2[i])+hcy*sin(phi2[i]))-rr[i])/(fiterror[i]*fiterror[i]);
222 //}
223 //pchi3/=3;
224}
E
R E
internal precision of FFTW codelets
Definition: DiscreteCosineTransform_31points.cc:36

◆ rPhiFitter() [1/2]

void rPhiFitter ( double *  rr,
double *  phi2,
double *  invphierror,
double &  rho,
double &  myphi0 
)
static

A circle fitter with invPhiError without fit chi2 output.

Definition at line 72 of file Fitter3DUtility.cc.

73{
74
75 double chi2;
76 rPhiFitter(rr, phi2, invphierror, rho, myphi0, chi2);
77
78}

◆ rPhiFitter() [2/2]

void rPhiFitter ( double *  rr,
double *  phi2,
double *  invphierror,
double &  rho,
double &  myphi0,
double &  chi2 
)
static

A circle fitter with invPhiError with fit chi2 output.

Definition at line 80 of file Fitter3DUtility.cc.

81{
82
83 // Print input values
84 //for(unsigned iSL=0; iSL<5; iSL++){
85 // cout<<"SL["<<iSL<<"] rr: "<<rr[iSL]<<" phi: "<<phi2[iSL]<<" phiError: "<<phierror[iSL]<<endl;
86 //}
87
88 double A, B, C, D, E, hcx, hcy;
89 double invFiterror[5];
90 //double G;
91 //Calculate fit error
92 for (unsigned i = 0; i < 5; i++) {
93 // Sometimes SL8 and SL4 will not be hit. So cannot use below calculation.
94 //invFiterror[i]=1/sqrt((rr[4]*rr[4]-2*rr[4]*rr[2]*cos(phi2[4]-phi2[2])+rr[2]*rr[2])/(sin(phi2[4]-phi2[2])*sin(phi2[4]-phi2[2]))-rr[i]*rr[i])*invphierror[i];
95 //invFiterror[i]=invphierror[i];
96 //invFiterror[i]=invphierror[i]*rr[i];
97 invFiterror[i] = invphierror[i] / rr[i];
98 }
99
100 //r-phi fitter(2D Fitter) ->calculate pt and radius of track-> input for 3D fitter.
101 A = 0, B = 0, C = 0, D = 0, E = 0;
102 //G=0;
103 for (unsigned i = 0; i < 5; i++) {
104 A += cos(phi2[i]) * cos(phi2[i]) * (invFiterror[i] * invFiterror[i]);
105 B += sin(phi2[i]) * sin(phi2[i]) * (invFiterror[i] * invFiterror[i]);
106 C += cos(phi2[i]) * sin(phi2[i]) * (invFiterror[i] * invFiterror[i]);
107 D += rr[i] * cos(phi2[i]) * (invFiterror[i] * invFiterror[i]);
108 E += rr[i] * sin(phi2[i]) * (invFiterror[i] * invFiterror[i]);
109 //G+=rr[i]*rr[i]/(fiterror[i]*fiterror[i]);
110 }
111 hcx = D * B - E * C; //helix center x
112 hcx /= 2 * (A * B - C * C);
113 hcy = E * A - D * C; //helix center y
114 hcy /= 2 * (A * B - C * C);
115 rho = sqrt(hcx * hcx + hcy * hcy); //radius of helix
116 myphi0 = atan2(hcy, hcx);
117 if (myphi0 < 0) myphi0 += 2 * M_PI;
118 //myphi0=atan(hcy/hcx);
119 //if(hcx<0 && hcy>0) myphi0 += M_PI;
120 //if(hcx<0 && hcy<0) myphi0 += M_PI;
121 //if(hcx>0 && hcy<0) myphi0 += M_PI*2.0;
122
124 // double pchi2 = -2*hcx*D-2*hcy*E+G;
125 // pchi2/=3;
126
127 // Count number of total TS hits.
128 int nTSHits = 0;
129 for (unsigned iAx = 0; iAx < 5; iAx++) {
130 if (invphierror[iAx] != 0) nTSHits++;
131 }
132 // Another way to calculate chi2
133 chi2 = 0.;
134 for (unsigned i = 0; i < 5; i++) {
135 chi2 += (2 * (hcx * cos(phi2[i]) + hcy * sin(phi2[i])) - rr[i]) * (2 * (hcx * cos(phi2[i]) + hcy * sin(phi2[i])) - rr[i]) /
136 (invFiterror[i] * invFiterror[i]);
137 }
138 chi2 /= nTSHits - 2;
139
140}

◆ rSFit() [1/2]

void rSFit ( double *  iezz2,
double *  arcS,
double *  zz,
double &  z0,
double &  cot,
double &  zchi2 
)
static

3D fitter functions Fits z and arc S.

Definition at line 1303 of file Fitter3DUtility.cc.

1304{
1305
1306 double ss = 0, sx = 0, sxx = 0, cotnum = 0, z0num = 0;
1307 double z0nump1[4], z0nump2[4];
1308 double z0den, iz0den;
1309
1310 for (unsigned i = 0; i < 4; i++) {
1311 ss += iezz2[i];
1312 sx += arcS[i] * iezz2[i];
1313 sxx += arcS[i] * arcS[i] * iezz2[i];
1314 }
1315
1316 for (unsigned i = 0; i < 4; i++) {
1317 cotnum += (ss * arcS[i] - sx) * iezz2[i] * zz[i];
1318 z0nump1[i] = sxx - sx * arcS[i];
1319 z0nump2[i] = z0nump1[i] * iezz2[i] * zz[i];
1320 z0num += z0nump2[i];
1321 }
1322 z0den = (ss * sxx) - (sx * sx);
1323 iz0den = 1. / z0den;
1324 z0num *= iz0den;
1325 cotnum *= iz0den;
1326 z0 = z0num;
1327 cot = cotnum;
1328
1329 // Count number of total TS hits.
1330 int nTSHits = 0;
1331 for (unsigned iSt = 0; iSt < 4; iSt++) {
1332 if (iezz2[iSt] != 0) nTSHits++;
1333 }
1334 // Calculate chi2 of z0
1335 zchi2 = 0.;
1336 for (unsigned i = 0; i < 4; i++) {
1337 zchi2 += (zz[i] - z0 - cot * arcS[i]) * (zz[i] - z0 - cot * arcS[i]) * iezz2[i];
1338 }
1339 zchi2 /= (nTSHits - 2);
1340
1341}

◆ rSFit() [2/2]

void rSFit ( std::map< std::string, double > const &  mConstD,
std::map< std::string, std::vector< double > > const &  mConstV,
std::map< std::string, Belle2::TRGCDCJSignal > &  mSignalStorage,
std::map< std::string, Belle2::TRGCDCJLUT * > &  mLutStorage 
)
static

Fits z and arc S using JSingal with multiple constants.

Definition at line 1396 of file Fitter3DUtility.cc.

1399{
1400 Belle2::TRGCDCJSignalData* commonData = mSignalStorage.begin()->second.getCommonData();
1401 // sum1_p1_0 = iZError2[0] + iZError2[1]
1402 mSignalStorage["sum1_p1_0"] <= mSignalStorage["iZError2_0"] + mSignalStorage["iZError2_1"];
1403 // sum1_p1_1 = iZError2[2] + iZError2[3]
1404 mSignalStorage["sum1_p1_1"] <= mSignalStorage["iZError2_2"] + mSignalStorage["iZError2_3"];
1405 //for(int iSt=0; iSt<2; iSt++) {cout<<"<<<sum1_p1_"<<iSt<<">>>"<<endl; mSignalStorage["sum1_p1_"+to_string(iSt)].dump();}
1406 // sum1 = sum1_p1[0] + sum1_p1[1]
1407 mSignalStorage["sum1"] <= mSignalStorage["sum1_p1_0"] + mSignalStorage["sum1_p1_1"];
1408 //cout<<"<<<sum1>>>"<<endl; mSignalStorage["sum1"].dump();
1409 // sumS_p1[i] = arcS[i] * iZError2[i]
1410 for (unsigned iSt = 0; iSt < 4; iSt++) {
1411 mSignalStorage["sumS_p1_" + to_string(iSt)] <= mSignalStorage["arcS_" + to_string(iSt)] * mSignalStorage["iZError2_" + to_string(
1412 iSt)];
1413 }
1414 // sumSS_p1[i] = arcS[i] * arcS[i]
1415 for (unsigned iSt = 0; iSt < 4; iSt++) {
1416 mSignalStorage["sumSS_p1_" + to_string(iSt)] <= mSignalStorage["arcS_" + to_string(iSt)] * mSignalStorage["arcS_" + to_string(iSt)];
1417 }
1418 // cotNumS_p1[i] = sum1 * arcS[i]
1419 for (unsigned iSt = 0; iSt < 4; iSt++) {
1420 mSignalStorage["cotNumS_p1_" + to_string(iSt)] <= mSignalStorage["sum1"] * mSignalStorage["arcS_" + to_string(iSt)];
1421 }
1422 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<sumS_p1_"<<iSt<<">>>"<<endl; mSignalStorage["sumS_p1_"+to_string(iSt)].dump();}
1423 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<sumSS_p1_"<<iSt<<">>>"<<endl; mSignalStorage["sumSS_p1_"+to_string(iSt)].dump();}
1424 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<cotNumS_p1_"<<iSt<<">>>"<<endl; mSignalStorage["cotNumS_p1_"+to_string(iSt)].dump();}
1425 // sumS_p2[0] = sumS_p1[0] + sumS_p1[1]
1426 mSignalStorage["sumS_p2_0"] <= mSignalStorage["sumS_p1_0"] + mSignalStorage["sumS_p1_1"];
1427 // sumS_p2[1] = sumS_p1[2] + sumS_p1[3]
1428 mSignalStorage["sumS_p2_1"] <= mSignalStorage["sumS_p1_2"] + mSignalStorage["sumS_p1_3"];
1429 // sumSS_p2[i] = sumSS_p1[i] * iZError2[i]
1430 for (unsigned iSt = 0; iSt < 4; iSt++) {
1431 mSignalStorage["sumSS_p2_" + to_string(iSt)] <= mSignalStorage["sumSS_p1_" + to_string(iSt)] * mSignalStorage["iZError2_" +
1432 to_string(iSt)];
1433 }
1434 //for(int iSt=0; iSt<2; iSt++) {cout<<"<<<sumS_p2_"<<iSt<<">>>"<<endl; mSignalStorage["sumS_p2_"+to_string(iSt)].dump();}
1435 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<sumSS_p2_"<<iSt<<">>>"<<endl; mSignalStorage["sumSS_p2_"+to_string(iSt)].dump();}
1436 // sumS = sumS_p2[0] + sumS_p2[1]
1437 mSignalStorage["sumS"] <= mSignalStorage["sumS_p2_0"] + mSignalStorage["sumS_p2_1"];
1438 // sumSS_p3[0] = sumSS_p2[0] + sumSS_p2[1]
1439 mSignalStorage["sumSS_p3_0"] <= mSignalStorage["sumSS_p2_0"] + mSignalStorage["sumSS_p2_1"];
1440 // sumSS_p3[1] = sumSS_p2[2] + sumSS_p2[3]
1441 mSignalStorage["sumSS_p3_1"] <= mSignalStorage["sumSS_p2_2"] + mSignalStorage["sumSS_p2_3"];
1442 //cout<<"<<<sumS>>>"<<endl; mSignalStorage["sumS"].dump();
1443 //for(int iSt=0; iSt<2; iSt++) {cout<<"<<<sumSS_p3_"<<iSt<<">>>"<<endl; mSignalStorage["sumSS_p3_"+to_string(iSt)].dump();}
1444 //cout<<"===================PIPELINE7============================"<<endl;
1445 // z0NumS_p1[i] = arcS[i] * sumS
1446 for (unsigned iSt = 0; iSt < 4; iSt++) {
1447 mSignalStorage["z0NumS_p1_" + to_string(iSt)] <= mSignalStorage["arcS_" + to_string(iSt)] * mSignalStorage["sumS"];
1448 }
1449 // den_p1 = sumS * sumS
1450 mSignalStorage["den_p1"] <= mSignalStorage["sumS"] * mSignalStorage["sumS"];
1451 // cotNumS[i] = cotNumS_p1[i] - sumS
1452 for (unsigned iSt = 0; iSt < 4; iSt++) {
1453 mSignalStorage["cotNumS_" + to_string(iSt)] <= mSignalStorage["cotNumS_p1_" + to_string(iSt)] - mSignalStorage["sumS"];
1454 }
1455 // sumSS = sumSS_p3[0] + sumSS_p3[1]
1456 mSignalStorage["sumSS"] <= mSignalStorage["sumSS_p3_0"] + mSignalStorage["sumSS_p3_1"];
1457 // zxiZError[i] = zz[i] * iZError2[i]
1458 for (unsigned iSt = 0; iSt < 4; iSt++) {
1459 mSignalStorage["zxiZError_" + to_string(iSt)] <= mSignalStorage["zz_" + to_string(iSt)] * mSignalStorage["iZError2_" + to_string(
1460 iSt)];
1461 }
1462 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<z0NumS_p1_"<<iSt<<">>>"<<endl; mSignalStorage["z0NumS_p1_"+to_string(iSt)].dump();}
1463 //cout<<"<<<den_p1>>>"<<endl; mSignalStorage["den_p1"].dump();
1464 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<cotNumS_"<<iSt<<">>>"<<endl; mSignalStorage["cotNumS_"+to_string(iSt)].dump();}
1465 //cout<<"<<<sumSS>>>"<<endl; mSignalStorage["sumSS"].dump();
1466 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<zxiZError_"<<iSt<<">>>"<<endl; mSignalStorage["zxiZError_"+to_string(iSt)].dump();}
1467 // den_p2 = sum1 * sumSS
1468 mSignalStorage["den_p2"] <= mSignalStorage["sum1"] * mSignalStorage["sumSS"];
1469 // z0NumS[i] =sumSS - z0NumS_p1[i]
1470 for (unsigned iSt = 0; iSt < 4; iSt++) {
1471 mSignalStorage["z0NumS_" + to_string(iSt)] <= mSignalStorage["sumSS"] - mSignalStorage["z0NumS_p1_" + to_string(iSt)];
1472 }
1473 //cout<<"<<<den_p2>>>"<<endl; mSignalStorage["den_p2"].dump();
1474 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<z0NumS_"<<iSt<<">>>"<<endl; mSignalStorage["z0NumS_"+to_string(iSt)].dump();}
1475 // cotNumSZ[i] = cotNumS[i] * zxiZError[i]
1476 for (unsigned iSt = 0; iSt < 4; iSt++) {
1477 mSignalStorage["cotNumSZ_" + to_string(iSt)] <= mSignalStorage["cotNumS_" + to_string(iSt)] * mSignalStorage["zxiZError_" +
1478 to_string(iSt)];
1479 }
1480 // den = den_p2 - den_p1
1481 mSignalStorage["den"] <= mSignalStorage["den_p2"] - mSignalStorage["den_p1"];
1482 // Calculate denMax and denMin.
1483 if (mSignalStorage.find("denMin") == mSignalStorage.end()) {
1484 // Create input for calDen.
1485 double t_rr3D[4], t_wireZError[4];
1486 for (int iSt = 0; iSt < 2; iSt++) {
1487 t_rr3D[iSt] = mConstV.at("rr3D")[iSt];
1488 t_wireZError[iSt] = 1 / mConstV.at("wireZError")[iSt];
1489 }
1490 for (int iSt = 2; iSt < 4; iSt++) {
1491 t_rr3D[iSt] = mConstV.at("rr3D")[iSt];
1492 t_wireZError[iSt] = 0;
1493 }
1494 double t_denMin = Fitter3DUtility::calDenWithIZError(t_rr3D, t_wireZError);
1495 mSignalStorage["denMin"] = Belle2::TRGCDCJSignal(t_denMin, mSignalStorage["den"].getToReal(), commonData);
1496 double t_arcSMax[4], t_driftZError[4];
1497 for (unsigned iSt = 0; iSt < 4; iSt++) {
1498 t_arcSMax[iSt] = mConstD.at("M_PI") * mConstV.at("rr3D")[iSt] / 2;
1499 t_driftZError[iSt] = 1 / mConstV.at("driftZError")[iSt];
1500 }
1501 double t_denMax = Fitter3DUtility::calDenWithIZError(t_arcSMax, t_driftZError);
1502 mSignalStorage["denMax"] = Belle2::TRGCDCJSignal(t_denMax, mSignalStorage["den"].getToReal(), commonData);
1503 }
1504 // Constrain den.
1505 {
1506 // Make ifElse data.
1507 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1508 // Compare
1509 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["den"], ">", mSignalStorage["denMax"]);
1510 // Assignments
1511 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1512 make_pair(&mSignalStorage["den_c"], mSignalStorage["denMax"])
1513 };
1514 // Push to data.
1515 t_data.push_back(make_pair(t_compare, t_assigns));
1516 // Compare
1517 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["den"], ">", mSignalStorage["denMin"]);
1518 // Assignments
1519 t_assigns = {
1520 make_pair(&mSignalStorage["den_c"], mSignalStorage["den"].limit(mSignalStorage["denMin"], mSignalStorage["denMax"]))
1521 };
1522 // Push to data.
1523 t_data.push_back(make_pair(t_compare, t_assigns));
1524 // Compare
1525 t_compare = Belle2::TRGCDCJSignal();
1526 // Assignments
1527 t_assigns = {
1528 make_pair(&mSignalStorage["den_c"], mSignalStorage["denMin"])
1529 };
1530 // Push to data.
1531 t_data.push_back(make_pair(t_compare, t_assigns));
1532 Belle2::TRGCDCJSignal::ifElse(t_data);
1533 }
1534 // z0NumSZ[i] = z0NumS[i]*zxiZError[i]
1535 for (unsigned iSt = 0; iSt < 4; iSt++) {
1536 mSignalStorage["z0NumSZ_" + to_string(iSt)] <= mSignalStorage["z0NumS_" + to_string(iSt)] * mSignalStorage["zxiZError_" + to_string(
1537 iSt)];
1538 }
1539 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<cotNumSZ_"<<iSt<<">>>"<<endl; mSignalStorage["cotNumSZ_"+to_string(iSt)].dump();}
1540 //cout<<"<<<den>>>"<<endl; mSignalStorage["den"].dump();
1541 //cout<<"<<<den_c>>>"<<endl; mSignalStorage["den_c"].dump();
1542 //cout<<"<<<denMin>>>"<<endl; mSignalStorage["denMin"].dump();
1543 //cout<<"<<<denMax>>>"<<endl; mSignalStorage["denMax"].dump();
1544 //cout<<"<<<offsetDen>>>"<<endl; mSignalStorage["offsetDen"].dump();
1545 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<z0NumSZ_"<<iSt<<">>>"<<endl; mSignalStorage["z0NumSZ_"+to_string(iSt)].dump();}
1546 // cot_p1[0] = cotNumSZ[0] + cotNumSZ[1]
1547 mSignalStorage["cot_p1_0"] <= mSignalStorage["cotNumSZ_0"] + mSignalStorage["cotNumSZ_1"];
1548 // cot_p1[1] = cotNumSZ[2] + cotNumSZ[3]
1549 mSignalStorage["cot_p1_1"] <= mSignalStorage["cotNumSZ_2"] + mSignalStorage["cotNumSZ_3"];
1550 // z0_p1_0 = z0NumSZ_0 + z0NumSZ_1
1551 mSignalStorage["z0_p1_0"] <= mSignalStorage["z0NumSZ_0"] + mSignalStorage["z0NumSZ_1"];
1552 // z0_p1_1 = z0NumSZ_2 + z0NumSZ_3
1553 mSignalStorage["z0_p1_1"] <= mSignalStorage["z0NumSZ_2"] + mSignalStorage["z0NumSZ_3"];
1555 //mSignalStorage["iDenMin"] <= Belle2::TRGCDCJSignal(1/t_denMax, pow(1/mSignalStorage["rho"].getToReal()/mSignalStorage["iZError2_0"].getToReal(),2));
1556 //mSignalStorage["iDenMax"] <= Belle2::TRGCDCJSignal(1/t_denMin, pow(1/mSignalStorage["rho"].getToReal()/mSignalStorage["iZError2_0"].getToReal(),2));
1557 // Calculate minInvValue, maxInvValue for LUTs
1558 if (mSignalStorage.find("invIDenMin") == mSignalStorage.end()) {
1559 unsigned long long t_int = mSignalStorage["den_c"].getMaxInt();
1560 double t_toReal = mSignalStorage["den_c"].getToReal();
1561 double t_actual = mSignalStorage["den_c"].getMaxActual();
1562 mSignalStorage["invIDenMin"] = Belle2::TRGCDCJSignal(t_int, t_toReal, t_int, t_int, t_actual, t_actual, t_actual, -1, commonData);
1563 t_int = mSignalStorage["den_c"].getMinInt();
1564 t_actual = mSignalStorage["den_c"].getMinActual();
1565 mSignalStorage["invIDenMax"] = Belle2::TRGCDCJSignal(t_int, t_toReal, t_int, t_int, t_actual, t_actual, t_actual, -1, commonData);
1566 }
1567 //cout<<"<<<invIDenMin>>>"<<endl; mSignalStorage["invIDenMin"].dump();
1568 //cout<<"<<<invIDenMax>>>"<<endl; mSignalStorage["invIDenMax"].dump();
1569 // Generate LUT(iDen = 1/den)
1570 if (!mLutStorage.count("iDen")) {
1571 mLutStorage["iDen"] = new Belle2::TRGCDCJLUT("iDen");
1572 mLutStorage["iDen"]->setFloatFunction(
1573 [ = ](double aValue) -> double{return 1 / aValue;},
1574 mSignalStorage["den_c"],
1575 mSignalStorage["invIDenMin"], mSignalStorage["invIDenMax"],
1576 pow(1 / mSignalStorage["rho"].getToReal() / mSignalStorage["iZError2_0"].getToReal(), 2),
1577 (int)mConstD.at("iDenLUTInBitSize"), (int)mConstD.at("iDenLUTOutBitSize"));
1578 //mLutStorage["iDen"]->makeCOE("./LutData/iDen.coe");
1579 }
1580 // Operate using LUT(iDen = 1/den)
1581 mLutStorage["iDen"]->operate(mSignalStorage["den_c"], mSignalStorage["iDen"]);
1582 //for(int iSt=0; iSt<2; iSt++) {cout<<"<<<cot_p1_"<<iSt<<">>>"<<endl; mSignalStorage["cot_p1_"+to_string(iSt)].dump();}
1583 //for(int iSt=0; iSt<2; iSt++) {cout<<"<<<z0_p1_"<<iSt<<">>>"<<endl; mSignalStorage["z0_p1_"+to_string(iSt)].dump();}
1584 //cout<<"<<<offsetDen_c>>>"<<endl; mSignalStorage["offsetDen_c"].dump();
1585 //cout<<"<<<denOffset_c>>>"<<endl; mSignalStorage["denOffset_c"].dump();
1586 //cout<<"<<<iDen>>>"<<endl; mSignalStorage["iDen"].dump();
1587 // iDen = offsetIDen + iDenOffset (Already done in LUT->operate)
1588 // cot_p2 = cot_p1_0 + cot_p1_1
1589 mSignalStorage["cot_p2"] <= mSignalStorage["cot_p1_0"] + mSignalStorage["cot_p1_1"];
1590 // z0_p2 = z0_p1_0 + z0_p1_1
1591 mSignalStorage["z0_p2"] <= mSignalStorage["z0_p1_0"] + mSignalStorage["z0_p1_1"];
1592 //cout<<"<<<cot_p2>>>"<<endl; mSignalStorage["cot_p2"].dump();
1593 //cout<<"<<<z0_p2>>>"<<endl; mSignalStorage["z0_p2"].dump();
1594 // cot = cot_p2 * iDen
1595 mSignalStorage["cot"] <= mSignalStorage["cot_p2"] * mSignalStorage["iDen"];
1596 // z0 = z0_p2 * iDen
1597 mSignalStorage["z0"] <= mSignalStorage["z0_p2"] * mSignalStorage["iDen"];
1598 if (mConstD.at("JB") == 1) {
1599 cout << "<<<cot>>>" << endl; mSignalStorage["cot"].dump();
1600 cout << "<<<z0>>>" << endl; mSignalStorage["z0"].dump();
1601 }
1602
1603 // Constrain z0 to [-30,30]
1604 /* cppcheck-suppress variableScope */
1605 double z0Min = -30;
1606 /* cppcheck-suppress variableScope */
1607 double z0Max = 30;
1608 // Calculate z0Max and z0Min
1609 if (!mSignalStorage.count("z0Min")) {
1610 mSignalStorage["z0Min"] = Belle2::TRGCDCJSignal(z0Min, mSignalStorage["z0"].getToReal(), commonData);
1611 mSignalStorage["z0Max"] = Belle2::TRGCDCJSignal(z0Max, mSignalStorage["z0"].getToReal(), commonData);
1612 }
1613 {
1614 // Make ifElse data.
1615 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1616 // Compare
1617 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["z0"], ">", mSignalStorage["z0Max"]);
1618 // Assignments
1619 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1620 make_pair(&mSignalStorage["z0_c"], mSignalStorage["z0Max"])
1621 };
1622 // Push to data.
1623 t_data.push_back(make_pair(t_compare, t_assigns));
1624 // Compare
1625 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["z0"], ">", mSignalStorage["z0Min"]);
1626 // Assignments
1627 t_assigns = {
1628 make_pair(&mSignalStorage["z0_c"], mSignalStorage["z0"].limit(mSignalStorage["z0Min"], mSignalStorage["z0Max"]))
1629 };
1630 // Push to data.
1631 t_data.push_back(make_pair(t_compare, t_assigns));
1632 // Compare
1633 t_compare = Belle2::TRGCDCJSignal();
1634 // Assignments
1635 t_assigns = {
1636 make_pair(&mSignalStorage["z0_c"], mSignalStorage["z0Min"])
1637 };
1638 // Push to data.
1639 t_data.push_back(make_pair(t_compare, t_assigns));
1640 Belle2::TRGCDCJSignal::ifElse(t_data);
1641 }
1642 //cout<<"<<<z0Min>>>"<<endl; mSignalStorage["z0Min"].dump();
1643 //cout<<"<<<z0Max>>>"<<endl; mSignalStorage["z0Max"].dump();
1644 //cout<<"<<<z0_c>>>"<<endl; mSignalStorage["z0_c"].dump();
1645 // Constraint cot to [-0.753, 1.428]
1646 double cotMin = cos(127 * mConstD.at("M_PI") / 180) / sin(127 * mConstD.at("M_PI") / 180);
1647 double cotMax = cos(35 * mConstD.at("M_PI") / 180) / sin(35 * mConstD.at("M_PI") / 180);
1648 if (!mSignalStorage.count("cotMin")) {
1649 mSignalStorage["cotMin"] = Belle2::TRGCDCJSignal(cotMin, mSignalStorage["cot"].getToReal(), commonData);
1650 mSignalStorage["cotMax"] = Belle2::TRGCDCJSignal(cotMax, mSignalStorage["cot"].getToReal(), commonData);
1651 }
1652 {
1653 // Make ifElse data.
1654 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1655 // Compare
1656 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["cot"], ">", mSignalStorage["cotMax"]);
1657 // Assignments
1658 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1659 make_pair(&mSignalStorage["cot_c"], mSignalStorage["cotMax"])
1660 };
1661 // Push to data.
1662 t_data.push_back(make_pair(t_compare, t_assigns));
1663 // Compare
1664 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["cot"], ">", mSignalStorage["cotMin"]);
1665 // Assignments
1666 t_assigns = {
1667 make_pair(&mSignalStorage["cot_c"], mSignalStorage["cot"].limit(mSignalStorage["cotMin"], mSignalStorage["cotMax"]))
1668 };
1669 // Push to data.
1670 t_data.push_back(make_pair(t_compare, t_assigns));
1671 // Compare
1672 t_compare = Belle2::TRGCDCJSignal();
1673 // Assignments
1674 t_assigns = {
1675 make_pair(&mSignalStorage["cot_c"], mSignalStorage["cotMin"])
1676 };
1677 // Push to data.
1678 t_data.push_back(make_pair(t_compare, t_assigns));
1679 Belle2::TRGCDCJSignal::ifElse(t_data);
1680 }
1681 //cout<<"<<<cotMin>>>"<<endl; mSignalStorage["cotMin"].dump();
1682 //cout<<"<<<cotMax>>>"<<endl; mSignalStorage["cotMax"].dump();
1683 //cout<<"<<<cot_c>>>"<<endl; mSignalStorage["cot_c"].dump();
1684 // Reduce number of bits for z0 to 11 bits and cot to 11 bits
1685 /* cppcheck-suppress variableScope */
1686 int z0Bits = 11;
1687 /* cppcheck-suppress variableScope */
1688 int cotBits = 11;
1689 {
1690 int z0Shift = mSignalStorage["z0_c"].getBitsize() - z0Bits;
1691 mSignalStorage["z0_r"] <= mSignalStorage["z0_c"].shift(z0Shift, 0);
1692 int cotShift = mSignalStorage["cot_c"].getBitsize() - cotBits;
1693 mSignalStorage["cot_r"] <= mSignalStorage["cot_c"].shift(cotShift, 0);
1694 }
1695 //cout<<"<<<z0_c>>>"<<endl; mSignalStorage["z0_c"].dump();
1696 //cout<<"<<<cot_c>>>"<<endl; mSignalStorage["cot_c"].dump();
1697 //cout<<"<<<z0_r>>>"<<endl; mSignalStorage["z0_r"].dump();
1698 //cout<<"<<<cot_r>>>"<<endl; mSignalStorage["cot_r"].dump();
1699
1700
1701 // fitDistFromZ0[i] = cot*arcS[i]
1702 for (unsigned iSt = 0; iSt < 4; iSt++) {
1703 mSignalStorage["fitDistFromZ0_" + to_string(iSt)] <= mSignalStorage["cot"] * mSignalStorage["arcS_" + to_string(iSt)];
1704 }
1705 // distFromZ0[i] = z[i] - z0
1706 for (unsigned iSt = 0; iSt < 4; iSt++) {
1707 mSignalStorage["distFromZ0_" + to_string(iSt)] <= mSignalStorage["zz_" + to_string(iSt)] - mSignalStorage["z0"];
1708 }
1709 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<fitDistFromZ0_"<<iSt<<">>>"<<endl; mSignalStorage["fitDistFromZ0_"+to_string(iSt)].dump();}
1710 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<distFromZ0_"<<iSt<<">>>"<<endl; mSignalStorage["distFromZ0_"+to_string(iSt)].dump();}
1711 // chiNum[i] = distFromZ0[i] - fitDistFromZ0[i]
1712 for (unsigned iSt = 0; iSt < 4; iSt++) {
1713 mSignalStorage["chiNum_" + to_string(iSt)] <= mSignalStorage["distFromZ0_" + to_string(iSt)] - mSignalStorage["fitDistFromZ0_" +
1714 to_string(iSt)];
1715 }
1716 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<chiNum_"<<iSt<<">>>"<<endl; mSignalStorage["chiNum_"+to_string(iSt)].dump();}
1717 // Calculate chiNumMax[i], chiNumMin[i].
1718 if (mSignalStorage.find("chiNumMax_0") == mSignalStorage.end()) {
1719 for (unsigned iSt = 0; iSt < 4; iSt++) {
1720 string t_maxName = "chiNumMax_" + to_string(iSt);
1721 string t_minName = "chiNumMin_" + to_string(iSt);
1722 string t_valueName = "chiNum_" + to_string(iSt);
1723 double t_maxValue = 2 * mConstV.at("zToOppositeStraw")[iSt];
1724 double t_minValue = 2 * mConstV.at("zToStraw")[iSt];
1725 mSignalStorage[t_maxName] = Belle2::TRGCDCJSignal(t_maxValue, mSignalStorage[t_valueName].getToReal(), commonData);
1726 mSignalStorage[t_minName] = Belle2::TRGCDCJSignal(t_minValue, mSignalStorage[t_valueName].getToReal(), commonData);
1727 }
1728 }
1729 // Constrain chiNum[i] to [-2*z_min[i], 2*z_max[i]]
1730 for (unsigned iSt = 0; iSt < 4; iSt++) {
1731 string t_maxName = "chiNumMax_" + to_string(iSt);
1732 string t_minName = "chiNumMin_" + to_string(iSt);
1733 string t_valueName = "chiNum_c_" + to_string(iSt);
1734 string t_in1Name = "chiNum_" + to_string(iSt);
1735 // Create data for ifElse
1736 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1737 // Compare
1738 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_maxName]);
1739 // Assignments
1740 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1741 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_maxName]),
1742 };
1743 // Push to data.
1744 t_data.push_back(make_pair(t_compare, t_assigns));
1745 // Compare
1746 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], ">", mSignalStorage[t_minName]);
1747 // Assignments
1748 t_assigns = {
1749 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_in1Name].limit(mSignalStorage[t_minName], mSignalStorage[t_maxName])),
1750 };
1751 // Push to data.
1752 t_data.push_back(make_pair(t_compare, t_assigns));
1753 // Compare
1754 t_compare = Belle2::TRGCDCJSignal();
1755 // Assignments
1756 t_assigns = {
1757 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_minName]),
1758 };
1759 // Push to data.
1760 t_data.push_back(make_pair(t_compare, t_assigns));
1761 // Process ifElse data.
1762 Belle2::TRGCDCJSignal::ifElse(t_data);
1763 }
1764 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<chiNumMax_"<<iSt<<">>>"<<endl; mSignalStorage["chiNumMax_"+to_string(iSt)].dump();}
1765 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<chiNumMin_"<<iSt<<">>>"<<endl; mSignalStorage["chiNumMin_"+to_string(iSt)].dump();}
1766 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<chiNum_c_"<<iSt<<">>>"<<endl; mSignalStorage["chiNum_c_"+to_string(iSt)].dump();}
1767 // chi2_p1[i] = chiNum_c[i] * chiNum_c[i]
1768 for (unsigned iSt = 0; iSt < 4; iSt++) {
1769 mSignalStorage["chi2_p1_" + to_string(iSt)] <= mSignalStorage["chiNum_c_" + to_string(iSt)] * mSignalStorage["chiNum_c_" +
1770 to_string(iSt)];
1771 }
1772 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<chi2_p1_"<<iSt<<">>>"<<endl; mSignalStorage["chi2_p1_"+to_string(iSt)].dump();}
1773 // chi2[i] = chi2_p1[i] * iError2[i]
1774 for (unsigned iSt = 0; iSt < 4; iSt++) {
1775 mSignalStorage["chi2_" + to_string(iSt)] <= mSignalStorage["chi2_p1_" + to_string(iSt)] * mSignalStorage["iZError2_" + to_string(
1776 iSt)];
1777 }
1778 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<chi2_"<<iSt<<">>>"<<endl; mSignalStorage["chi2_"+to_string(iSt)].dump();}
1779 // chi2Sum_p1_0 = chi2_0 + chi2_1
1780 mSignalStorage["chi2Sum_p1_0"] <= mSignalStorage["chi2_0"] + mSignalStorage["chi2_1"];
1781 // chi2Sum_p1_1 = chi2_2 + chi2_3
1782 mSignalStorage["chi2Sum_p1_1"] <= mSignalStorage["chi2_2"] + mSignalStorage["chi2_3"];
1783 //cout<<"<<<chi2Sum_p1_0>>>"<<endl; mSignalStorage["chi2Sum_p1_0"].dump();
1784 //cout<<"<<<chi2Sum_p1_1>>>"<<endl; mSignalStorage["chi2Sum_p1_1"].dump();
1785 // chi2Sum = chi2Sum_p1_0 + chi2Sum_p1_1
1786 mSignalStorage["zChi2"] <= (mSignalStorage["chi2Sum_p1_0"] + mSignalStorage["chi2Sum_p1_1"]).shift(1);
1787 if (mConstD.at("JB") == 1) {cout << "<<<zChi2>>>" << endl; mSignalStorage["zChi2"].dump();}
1788
1789 // Constrain zChi2 to [0, 100]
1790 /* cppcheck-suppress variableScope */
1791 double zChi2Min = 0;
1792 /* cppcheck-suppress variableScope */
1793 double zChi2Max = 100;
1794 /* cppcheck-suppress variableScope */
1795 int zChi2Bits = 4;
1796 // Calculate zChi2Max and zChi2Min
1797 if (!mSignalStorage.count("zChi2Min")) {
1798 mSignalStorage["zChi2Min"] = Belle2::TRGCDCJSignal(zChi2Min, mSignalStorage["zChi2"].getToReal(), commonData);
1799 mSignalStorage["zChi2Max"] = Belle2::TRGCDCJSignal(zChi2Max, mSignalStorage["zChi2"].getToReal(), commonData);
1800 }
1801 {
1802 // Make ifElse data.
1803 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
1804 // Compare
1805 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["zChi2"], ">", mSignalStorage["zChi2Max"]);
1806 // Assignments
1807 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
1808 make_pair(&mSignalStorage["zChi2_c"], mSignalStorage["zChi2Max"])
1809 };
1810 // Push to data.
1811 t_data.push_back(make_pair(t_compare, t_assigns));
1812 // Compare
1813 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["zChi2"], ">", mSignalStorage["zChi2Min"]);
1814 // Assignments
1815 t_assigns = {
1816 make_pair(&mSignalStorage["zChi2_c"], mSignalStorage["zChi2"].limit(mSignalStorage["zChi2Min"], mSignalStorage["zChi2Max"]))
1817 };
1818 // Push to data.
1819 t_data.push_back(make_pair(t_compare, t_assigns));
1820 // Compare
1821 t_compare = Belle2::TRGCDCJSignal();
1822 // Assignments
1823 t_assigns = {
1824 make_pair(&mSignalStorage["zChi2_c"], mSignalStorage["zChi2Min"])
1825 };
1826 // Push to data.
1827 t_data.push_back(make_pair(t_compare, t_assigns));
1828 Belle2::TRGCDCJSignal::ifElse(t_data);
1829 }
1830 //cout<<"<<<zChi2Min>>>"<<endl; mSignalStorage["zChi2Min"].dump();
1831 //cout<<"<<<zChi2Max>>>"<<endl; mSignalStorage["zChi2Max"].dump();
1832 //cout<<"<<<zChi2_c>>>"<<endl; mSignalStorage["zChi2_c"].dump();
1833 {
1834 int zChi2Shift = mSignalStorage["zChi2_c"].getBitsize() - zChi2Bits;
1835 mSignalStorage["zChi2_r"] <= mSignalStorage["zChi2_c"].shift(zChi2Shift, 0);
1836 }
1837 //cout<<"<<<zChi2_c>>>"<<endl; mSignalStorage["zChi2_c"].dump();
1838 //cout<<"<<<zChi2_r>>>"<<endl; mSignalStorage["zChi2_r"].dump();
1839
1840}
Fitter3DUtility::calDenWithIZError
static double calDenWithIZError(double *arcS, double *iZError)
Calculates the denominator for fitting z and arc s.
Definition: Fitter3DUtility.cc:1289

◆ rSFit2()

void rSFit2 ( double *  iezz21,
double *  iezz22,
double *  arcS,
double *  zz,
int *  lutv,
double &  z0,
double &  cot,
double &  zchi2 
)
static

Fits z and arc S. (Will be deprecated.)

Definition at line 1343 of file Fitter3DUtility.cc.

1345{
1346// cout << "rs2" << endl;
1347
1348 double ss = 0, sx = 0, sxx = 0, cotnum = 0, z0num = 0;
1349 double z0nump1[4], z0nump2[4];
1350 double z0den, iz0den;
1351
1352 for (unsigned i = 0; i < 4; i++) {
1353 if (lutv[i] == 2) {
1354 ss += iezz21[i];
1355 sx += arcS[i] * iezz21[i];
1356 sxx += arcS[i] * arcS[i] * iezz21[i];
1357 } else {
1358 ss += iezz22[i];
1359 sx += arcS[i] * iezz22[i];
1360 sxx += arcS[i] * arcS[i] * iezz22[i];
1361 }
1362 }
1363
1364 for (unsigned i = 0; i < 4; i++) {
1365 if (lutv[i] == 2) {
1366 cotnum += (ss * arcS[i] - sx) * iezz21[i] * zz[i];
1367 z0nump1[i] = sxx - sx * arcS[i];
1368 z0nump2[i] = z0nump1[i] * iezz21[i] * zz[i];
1369 z0num += z0nump2[i];
1370 } else {
1371 cotnum += (ss * arcS[i] - sx) * iezz22[i] * zz[i];
1372 z0nump1[i] = sxx - sx * arcS[i];
1373 z0nump2[i] = z0nump1[i] * iezz22[i] * zz[i];
1374 z0num += z0nump2[i];
1375 }
1376 }
1377 z0den = (ss * sxx) - (sx * sx);
1378 iz0den = 1. / z0den;
1379 z0num *= iz0den;
1380 cotnum *= iz0den;
1381 z0 = z0num;
1382 cot = cotnum;
1383
1384 // Calculate chi2 of z0
1385// double zchi2 = 0.;
1386 for (unsigned i = 0; i < 4; i++) {
1387 if (lutv[i] == 2) {
1388 zchi2 += (zz[i] - z0 - cot * arcS[i]) * (zz[i] - z0 - cot * arcS[i]) * iezz21[i];
1389 } else {
1390 zchi2 += (zz[i] - z0 - cot * arcS[i]) * (zz[i] - z0 - cot * arcS[i]) * iezz22[i];
1391 }
1392 }
1393 zchi2 /= (4 - 2);
1394}

◆ saveStereoXt()

void saveStereoXt ( std::vector< std::vector< double > > &  stXts,
std::string const &  filePrefix 
)
static

Saves stereo Xt to file.

Definition at line 543 of file Fitter3DUtility.cc.

544{
545 for (unsigned iSt = 0; iSt < stXts.size(); ++iSt) {
546 string filename = filePrefix + ".st" + to_string(iSt);
547 ofstream xtFile;
548 cout << "Fitter3DUtility::saveStereoXt() Saving xt file to " << filename << endl;
549 xtFile.open(filename);
550 for (unsigned iTick = 0; iTick < stXts[iSt].size(); ++iTick) {
551 xtFile << stXts[iSt][iTick] << endl;
552 }
553 xtFile.close();
554 }
555}

◆ setError()

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 
)
static

Sets error using JSignal class.

Definition at line 681 of file Fitter3DUtility.cc.

683{
684 Belle2::TRGCDCJSignalData* commonData = mSignalStorage.begin()->second.getCommonData();
685 // Make constants for wireError,driftError,noneError
686 if (mSignalStorage.find("iZWireError2_0") == mSignalStorage.end()) {
687 for (unsigned iSt = 0; iSt < 4; iSt++) {
688 string t_name;
689 t_name = "iZWireError2_" + to_string(iSt);
690 mSignalStorage[t_name] = Belle2::TRGCDCJSignal(mConstD.at("iError2BitSize"), 1 / pow(mConstV.at("wireZError")[iSt], 2), 0,
691 mConstD.at("iError2Max"), -1, commonData);
692 t_name = "iZDriftError2_" + to_string(iSt);
693 mSignalStorage[t_name] = Belle2::TRGCDCJSignal(mConstD.at("iError2BitSize"), 1 / pow(mConstV.at("driftZError")[iSt], 2), 0,
694 mConstD.at("iError2Max"), -1, commonData);
695 t_name = "iZNoneError2_" + to_string(iSt);
696 mSignalStorage[t_name] = Belle2::TRGCDCJSignal(mConstD.at("iError2BitSize"), 0, 0, mConstD.at("iError2Max"), -1, commonData);
697 }
698 }
699 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<iZWireError2_"<<iSt<<">>>"<<endl; mSignalStorage["iZWireError2_"+to_string(iSt)].dump();}
700 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<iZDriftError2_"<<iSt<<">>>"<<endl; mSignalStorage["iZDriftError2_"+to_string(iSt)].dump();}
701 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<iZNoneError2_"<<iSt<<">>>"<<endl; mSignalStorage["iZNoneError2_"+to_string(iSt)].dump();}
702
703 // Set error depending on lr(0:no hit, 1: left, 2: right, 3: not determined)
704 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<lr_"<<iSt<<">>>"<<endl; mSignalStorage["lr_"+to_string(iSt)].dump();}
705 for (unsigned iSt = 0; iSt < 4; iSt++) {
706 string t_in1Name = "lr_" + to_string(iSt);
707 string t_valueName = "invErrorLR_" + to_string(iSt);
708 string t_noneErrorName = "iZNoneError2_" + to_string(iSt);
709 string t_driftErrorName = "iZDriftError2_" + to_string(iSt);
710 string t_wireErrorName = "iZWireError2_" + to_string(iSt);
711 // Create data for ifElse
712 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
713 // Compare (lr == 0)
714 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], "=", Belle2::TRGCDCJSignal(0,
715 mSignalStorage[t_in1Name].getToReal(), commonData));
716 // Assignments (invErrorLR <= iZNoneError2)
717 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
718 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_noneErrorName])
719 };
720 // Push to data.
721 t_data.push_back(make_pair(t_compare, t_assigns));
722 // Compare (lr == 3)
723 t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage[t_in1Name], "=", Belle2::TRGCDCJSignal(3,
724 mSignalStorage[t_in1Name].getToReal(), commonData));
725 // Assignments (invErrorLR <= iZWireError)
726 t_assigns = {
727 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_wireErrorName])
728 };
729 // Push to data.
730 t_data.push_back(make_pair(t_compare, t_assigns));
731 // Compare (else)
732 t_compare = Belle2::TRGCDCJSignal();
733 // Assignments (invErrorLR <= iZDriftError)
734 t_assigns = {
735 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_driftErrorName])
736 };
737 // Push to data.
738 t_data.push_back(make_pair(t_compare, t_assigns));
739 // Process ifElse data.
740 Belle2::TRGCDCJSignal::ifElse(t_data);
741 }
742 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<invErrorLR_"<<iSt<<">>>"<<endl; mSignalStorage["invErrorLR_"+to_string(iSt)].dump();}
743
744 // Make constants for half radius of SL
745 if (mSignalStorage.find("halfRadius_0") == mSignalStorage.end()) {
746 for (unsigned iSt = 0; iSt < 4; iSt++) {
747 string t_name;
748 t_name = "halfRadius_" + to_string(iSt);
749 mSignalStorage[t_name] = Belle2::TRGCDCJSignal(mConstV.at("rr")[iSt * 2 + 1] / 2, mSignalStorage["rho"].getToReal(), commonData);
750 }
751 }
752 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<halfRadius_"<<iSt<<">>>"<<endl; mSignalStorage["halfRadius_"+to_string(iSt)].dump();}
753 // Set error depending on R(helix radius)
754 //cout<<"<<<rho>>>"<<endl; mSignalStorage["rho"].dump();
755 for (unsigned iSt = 0; iSt < 4; iSt++) {
756 string t_compareName = "halfRadius_" + to_string(iSt);
757 //string t_valueName = "iZError2_" + to_string(iSt);
758 string t_valueName = "invErrorRho_" + to_string(iSt);
759 string t_noneErrorName = "iZNoneError2_" + to_string(iSt);
760 string t_in1Name = "invErrorLR_" + to_string(iSt);
761 // Create data for ifElse
762 vector<pair<Belle2::TRGCDCJSignal, vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > > > t_data;
763 // Compare (R < r/2)
764 Belle2::TRGCDCJSignal t_compare = Belle2::TRGCDCJSignal::comp(mSignalStorage["rho"], "<", mSignalStorage[t_compareName]);
765 // Assignments (invError <= 0)
766 vector<pair<Belle2::TRGCDCJSignal*, Belle2::TRGCDCJSignal> > t_assigns = {
767 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_noneErrorName])
768 };
769 // Push to data.
770 t_data.push_back(make_pair(t_compare, t_assigns));
771 // Compare (else)
772 t_compare = Belle2::TRGCDCJSignal();
773 // Assignments (invError <= invErrorLR)
774 t_assigns = {
775 make_pair(&mSignalStorage[t_valueName], mSignalStorage[t_in1Name])
776 };
777 // Push to data.
778 t_data.push_back(make_pair(t_compare, t_assigns));
779 // Process ifElse data.
780 Belle2::TRGCDCJSignal::ifElse(t_data);
781 }
782 //for(int iSt=0; iSt<4; iSt++) {cout<<"<<<invError_"<<iSt<<">>>"<<endl; mSignalStorage["invError_"+to_string(iSt)].dump();}
783}

◆ setErrorFast()

void setErrorFast ( std::vector< std::vector< int > > const &  rawStTSs,
int  eventTimeValid,
std::vector< double > &  invZError2 
)
static

Sets error for fitting.

Definition at line 663 of file Fitter3DUtility.cc.

665{
666 vector<double> driftZError({0.7676, 0.9753, 1.029, 1.372});
667 vector<double> wireZError({0.7676, 0.9753, 1.029, 1.372});
668 vector<double> zError(4, 9999);
669 invZError2 = vector<double> (4, 0);
670 for (unsigned iSt = 0; iSt < 4; ++iSt) {
671 if (rawStTSs[iSt][1] != 0) {
672 // Check LR and eventTime
673 if (rawStTSs[iSt][1] != 3 && eventTimeValid != 0) zError[iSt] = driftZError[iSt];
674 else zError[iSt] = wireZError[iSt];
675 // Get inverse zerror ^ 2
676 invZError2[iSt] = 1 / pow(zError[iSt], 2);
677 }
678 }
679}

◆ stereoIdToPhi()

double stereoIdToPhi ( std::vector< double > &  stNWires,
int  iSt,
int  id 
)
static

Converts stereo ts id to phi.

Definition at line 576 of file Fitter3DUtility.cc.

577{
578 return (double)id / stNWires[iSt] * 4 * M_PI;
579}

◆ toUnsignedTdc()

unsigned toUnsignedTdc ( int  tdc,
int  nBits 
)
static

Changes tdc and event time to unsigned value that has # bits.

Definition at line 226 of file Fitter3DUtility.cc.

227{
228 bool notDone = 1;
229 int tdcMax = pow(2, nBits) - 1;
230 int resultTdc = tdc;
231 while (notDone) {
232 if (resultTdc > tdcMax) resultTdc -= (tdcMax + 1);
233 else if (resultTdc < 0) resultTdc += (tdcMax + 1);
234 else notDone = 0;
235 }
236 return (unsigned) resultTdc;
237}
The documentation for this class was generated from the following files: