event_display/macro/track_pT.hh

0001 #pragma once
0002 #include "include_tracking/least_square2.cc"
0003
0004 class track
0005 {
0006 public:
0007     TF1 *track_xy;
0008     TF1 *track_rz;
0009     TF1 *track_rz_inv;
0010
0011     double phi_in;
0012     double phi_out;
0013
0014     double theta_in;
0015     double theta_out;
0016
0017     double d_phi;
0018     double d_theta;
0019
0020     double dca[3];
0021     double dca_mean[3];
0022
0023     TVector3 p1;
0024     TVector3 p2;
0025
0026     double dca_2d;
0027     double phi;
0028     double theta;
0029     double phi_tracklet;
0030     double theta_tracklet;
0031
0032     double charge;
0033     vector < vector < double > > zline;
0034
0035     vector < double > x;
0036     vector < double > y;
0037     vector < double > z;
0038     vector < double > r;
0039     vector < double > y_err;
0040     vector < double > z_err;
0041
0042     double pT;
0043     double rad;
0044     double cx, cy;
0045     double p_reco[3];
0046
0047     double A[3][3];
0048     double B[1][3];
0049     double beta[1][3];
0050
0051     bool dca2d_range_out = false;
0052     bool dcaz_range_out = false;
0053     bool dca_range_out = false;
0054     bool is_deleted = false;
0055
0056     track() : track_xy(nullptr),
0057               track_rz(nullptr),
0058               track_rz_inv(nullptr),
0059               zline(2, vector<double>(0))
0060
0061     {
0062         p1 = TVector3();
0063         p2 = TVector3();
0064     };
0065
0066     virtual ~track()
0067     {
0068         if (track_xy != nullptr)
0069             delete track_xy;
0070         if (track_rz != nullptr)
0071             delete track_rz;
0072         if (track_rz_inv != nullptr)
0073             delete track_rz_inv;
0074     };
0075
0076     void calc_pT()
0077     {
0078         TVector3 p3 = {dca_mean[0], dca_mean[1], dca_mean[2]}; // vertex
0079
0080         TVector3 mid1 = TVector3((p1.x() + p2.x()) / 2, (p1.y() + p2.y()) / 2, (p1.z() + p2.z()) / 2);
0081         TVector3 mid2 = TVector3((p3.x() + p1.x()) / 2, (p3.y() + p1.y()) / 2, (p3.z() + p1.z()) / 2);
0082
0083         double alpha = atan2((p2.y() - p1.y()), (p2.x() - p1.x())) + M_PI / 2;
0084         double beta = atan2((p1.y() - p3.y()), (p1.x() - p3.x())) + M_PI / 2;
0085
0086         cx = (mid1.y() - mid2.y() - mid1.x() * tan(alpha) + mid2.x() * tan(beta)) / (tan(beta) - tan(alpha));
0087         cy = mid1.y() - (mid1.x() - cx) * tan(alpha);
0088
0089         rad = sqrt((p1.x() - cx) * (p1.x() - cx) + (p1.y() - cy) * (p1.y() - cy)); //[cm]
0090         pT = /*charge **/ 0.3 * 1.4 * rad / 100;                                   //[GeV/c]
0091
0092         double phi_p = atan2(p3.y() - cy, p3.x() - cx) + (M_PI / 2);
0093         double phi_p1 = atan2(p1.y() - cy, p1.x() - cx) + (M_PI / 2);
0094         // cout<<"p1 : "<<p1.x()<<"  "<<p1.y()<<endl;
0095         // cout<<"phi : "<<phi_p/M_PI * 180<<"  "<<phi_p1/M_PI * 180<<endl;
0096         double d_phi = phi_p1 - phi_p;
0097         // cout<<"d_phi : "<<d_phi/M_PI * 180<<endl;
0098         // d_phi -= -2 * (M_PI * (int)(d_phi / (M_PI)));
0099
0100         // if(fabs(d_phi) > M_PI)
0101         // {
0102         //     cout<<"fabs(d_phi) > M_PI  "<<d_phi<<endl;
0103         //     cout<<(M_PI * (int)(d_phi / (M_PI)))<<endl;
0104         // // d_phi -=  (M_PI * (int)(d_phi / (M_PI)));
0105         // d_phi -=  2 * (M_PI * (int)(d_phi / (M_PI)));
0106
0107         //     // d_phi -=d_phi * fabs(d_phi) *  M_PI * (int)(d_phi / M_PI);
0108         // }
0109         // cout<<" d_phi : "<<d_phi/M_PI * 180<<endl;
0110
0111         double phi_offset = 0;
0112         double cross = p3.x() * p1.y() - p3.y() * p1.x();
0113
0114         if ((d_phi) < 0)
0115             // if((phi_p1 - phi_p) < 0)
0116             phi_offset = M_PI;
0117         // double sign = d_phi / fabs(d_phi);
0118
0119         // double sign = ((phi_p1 -(M_PI/2))  - (phi_p - (M_PI/2)))/fabs((phi_p1 - (M_PI/2)) - (phi_p - (M_PI/2)));
0120         // cout<<sign<<endl;
0121         // double sign_rz = (p1.z() - dca_mean[2]) / fabs(p1.z() - dca_mean[2]);
0122         // double sign_x = (p1.x() - p3.x()) / fabs(p1.x() - p3.x());
0123         // double sign_y = (p1.y() - p3.y()) / fabs(p1.y() - p3.y());
0124         // cout<<"  sign_x,y : "<<sign_x<<"  "<<sign_y<<endl;
0125         double t_ = sqrt(p1.x() * p1.x() + p1.y() * p1.y());
0126         // cout<<"no_sign : "<<fabs(pT) * cos(phi_p)<<"  "<<fabs(pT) * sin(phi_p)<<endl;
0127         // cout<<"phi : "<<phi_p + phi_offset<<endl;
0128
0129         p_reco[0] = /*sign_x **/ fabs(pT) * cos(phi_p + phi_offset);
0130         p_reco[1] = /*sign_y **/ fabs(pT) * sin(phi_p + phi_offset);
0131         // cout<<"p1 : "<<p1.x() - dca_mean[0]<<"  "<<p1.y() - dca_mean[1]<<endl;
0132         // cout<<"p  : "<<p_reco[0]<<"  "<<p_reco[1]<<endl;
0133         p_reco[2] = fabs(pT) / t_ * p1.z();
0134         // cout<<endl;
0135     }
0136
0137     void calc_line()
0138     {
0139         x.push_back(p1.x());
0140         x.push_back(p2.x());
0141         x.push_back(dca_mean[0]);
0142
0143         y.push_back(p1.y());
0144         y.push_back(p2.y());
0145         y.push_back(dca_mean[1]);
0146
0147         y_err.push_back(1.);
0148         y_err.push_back(1.);
0149         y_err.push_back(1.);
0150
0151         LeastSquare least_square;
0152         least_square.Setdatax(x);
0153         least_square.Setdatay(y);
0154         least_square.Setdatayerr(y_err);
0155
0156         least_square.SetDebugMode(false);
0157         least_square.Calc(0);
0158
0159         if (track_xy != nullptr)
0160         {
0161             delete track_xy;
0162             track_xy = nullptr;
0163         }
0164
0165         track_xy = new TF1("track_xy", "pol1", -15, 15);
0166
0167         track_xy->SetParameter(1, least_square.GetSlope());
0168         track_xy->SetParameter(0, least_square.GetIntercept());
0169
0170         z.push_back(p1.z());
0171         z.push_back(p2.z());
0172         z.push_back(dca_mean[2]);
0173
0174         r.push_back(getpointr(1));
0175         r.push_back(getpointr(2));
0176         r.push_back(getpointr(3));
0177
0178         z_err.push_back(20. / sqrt(12.)); // [mm]
0179         z_err.push_back(20. / sqrt(12.)); // [mm]
0180         z_err.push_back(0.8);             // [mm]
0181
0182         // r-z plane
0183         LeastSquare least_square_zr;
0184         least_square_zr.Setdatax(r);
0185         least_square_zr.Setdatay(z);
0186         least_square_zr.Setdatayerr(z_err);
0187
0188         // least_square_zr.SetDebugMode(false);
0189         least_square_zr.Calc(1);
0190
0191         if (track_rz != nullptr)
0192         {
0193             delete track_rz;
0194             track_rz = nullptr;
0195         }
0196         double slope_inv = least_square_zr.GetSlope();
0197         double intercept_inv = least_square_zr.GetIntercept();
0198
0199         track_rz = new TF1("track_rz", "pol1", -25, 25);
0200
0201         track_rz->SetParameter(1, getslope_inv(slope_inv, intercept_inv));
0202         track_rz->SetParameter(0, getintercept_inv(slope_inv, intercept_inv));
0203
0204         //        track_rz->SetParameter(1, least_square_zr.GetSlope());
0205         // track_rz->SetParameter(0, least_square_zr.GetIntercept());
0206     }
0207
0208     void calc_inv()
0209     {
0210         if (track_rz_inv != nullptr)
0211         {
0212             delete track_rz_inv;
0213             track_rz_inv = nullptr;
0214         }
0215
0216         double slope_ = 1 / track_rz->GetParameter(1);
0217         double intercept_ = -track_rz->GetParameter(0) / track_rz->GetParameter(1);
0218
0219         track_rz_inv = new TF1("track_rz_inv", "pol1", -25, 25);
0220
0221         track_rz_inv->SetParameter(1, slope_);
0222         track_rz_inv->SetParameter(0, intercept_);
0223     }
0224
0225     double getslope_inv(double Slope, double Intercept)
0226     {
0227         return 1 / Slope;
0228     }
0229
0230     double getintercept_inv(double Slope, double Intercept)
0231     {
0232         return -Intercept / Slope;
0233     }
0234
0235     double getpointr(int mode)
0236     {
0237         double r = 0;
0238
0239         if (mode == 1)
0240         {
0241             r = (p1.y() / fabs(p1.y())) * sqrt(p1.x() * p1.x() + p1.y() * p1.y());
0242         }
0243         else if (mode == 2)
0244         {
0245             r = (p2.y() / fabs(p2.y())) * sqrt(p2.x() * p2.x() + p2.y() * p2.y());
0246         }
0247         else if (mode == 3)
0248         {
0249             r = (dca_mean[1] / fabs(dca_mean[1])) * sqrt(dca_mean[0] * dca_mean[0] + dca_mean[1] * dca_mean[1]);
0250         }
0251
0252         return r;
0253     }
0254
0255     double getphi()
0256     {
0257         phi = atan(track_xy->GetParameter(1));
0258
0259         if (p1.x() < dca_mean[0])
0260         {
0261             if (p1.y() < dca_mean[1])
0262                 phi -= M_PI;
0263             else
0264                 phi += M_PI;
0265         }
0266
0267         return phi;
0268     }
0269
0270     double gettheta()
0271     {
0272         theta = atan(track_rz->GetParameter(1));
0273
0274         if (p1.z() < dca_mean[2])
0275         {
0276             if (getpointr(1) < getpointr(3))
0277                 theta -= M_PI;
0278             else
0279                 theta += M_PI;
0280         }
0281
0282         if (getpointr(1) < getpointr(3))
0283             theta = fabs(theta);
0284
0285         return theta;
0286     }
0287
0288     double getphi_tracklet()
0289     {
0290         phi_tracklet = atan2((p2.y() - p1.y()), (p2.x() - p1.x()));
0291         return phi_tracklet;
0292     }
0293
0294     double gettheta_tracklet()
0295     {
0296         theta_tracklet = atan2((fabs(getpointr(2)) - fabs(getpointr(1))), (p2.z() - p1.z()));
0297
0298         return theta_tracklet;
0299     }
0300
0301     vector<double> get_crossline(double x1, double y1, double r1)
0302     {
0303         vector<double> cross_co;
0304         // ax + by + c = 0
0305         cross_co.push_back(2 * (cx - x1));                                                           // a
0306         cross_co.push_back(2 * (cy - y1));                                                           // b
0307         cross_co.push_back((x1 + cx) * (x1 - cx) + (y1 + cy) * (y1 - cy) + (rad + r1) * (rad - r1)); // c
0308
0309         return cross_co;
0310     }
0311
0312     vector<double> get_crosspoint(const vector<double> &V) // ax + by + c = 0
0313     {
0314         vector<double> cross;
0315         double a = V[0];
0316         double b = V[1];
0317         double c = V[2];
0318
0319         double d = fabs(a * cx + b * cy + c) / sqrt(a * a + b * b);
0320         // cout << "d ; " << d << " " << rad << endl;
0321
0322         double theta = atan2(b, a);
0323
0324         if (d > rad)
0325         {
0326             cout << "d > rad" << endl;
0327         }
0328         else if (d == rad)
0329         {
0330             cout << "d == rad" << endl;
0331
0332             if (a * cx + b * cy + c > 0)
0333                 theta += M_PI;
0334             cross.push_back(rad * cos(theta) + cx);
0335             cross.push_back(rad * sin(theta) + cy);
0336         }
0337         else
0338         {
0339             cout << "d < rad" << endl;
0340             double alpha, beta, phi;
0341             phi = acos(d / rad);
0342             alpha = theta - phi;
0343             beta = theta + phi;
0344             if (a * cx + b * cy + c > 0)
0345             {
0346                 alpha += M_PI;
0347                 beta += M_PI;
0348             }
0349             // double temp_cross[2][2];
0350             vector<vector<double>> temp_cross = {{rad * cos(alpha) + cx, rad * cos(beta) + cx}, {rad * sin(alpha) + cy, rad * sin(beta) + cy}};
0351             // for (unsigned int j = 0; j < temp_cross.at(0).size(); j++)
0352             // {
0353             //     cout << "temp_cross : ";
0354             //     for (unsigned int i = 0; i < temp_cross.size(); i++)
0355             //         cout << temp_cross.at(j).at(i) << "  ";
0356             //     cout << endl;
0357             // }
0358             cross = get_closestpoint(temp_cross);
0359         }
0360         // cout << "cross size : " << cross.size() << endl;
0361         // for (int i = 0; i < cross.size(); i++)
0362         //     cout << cross[i] << endl;
0363         return cross;
0364     }
0365
0366     vector<double> get_closestpoint(const vector<vector<double>> VV)
0367     {
0368         vector<double> closest;
0369         double pre_phi = 999;
0370         double phi_p1 = atan2(p1.y() - dca_mean[1], p1.x() - dca_mean[0]);
0371         for (unsigned int i = 0; i < VV.at(0).size(); i++)
0372         {
0373             cout << "dca_mean : " << dca_mean[0] << "  " << dca_mean[1] << endl;
0374             cout << "VV : " << VV[0][i] << "  " << VV[1][i] << endl;
0375             cout << "pVV : " << VV[0][i] - dca_mean[0] << "  " << VV[1][i] - dca_mean[1] << endl;
0376             cout << "p1  : " << p1.x() - dca_mean[0] << "  " << p1.y() - dca_mean[1] << endl;
0377             // cout << VV.at(i).at(0) << endl;
0378             double dot = (VV[0][i] - dca_mean[0]) * (p1.x() - dca_mean[0]) + (VV[1][i] - dca_mean[1]) * (p1.y() - dca_mean[1]);
0379             double temp_phi = atan2(VV[1][i] - dca_mean[1], VV[0][i] - dca_mean[0]);
0380             cout << "dot : " << dot << endl;
0381             if (dot > 0)
0382             {
0383                 if (fabs(temp_phi - phi_p1) < fabs(pre_phi - phi_p1))
0384                 {
0385                     closest.erase(closest.begin(), closest.end());
0386                     cout << "VV size : " << VV.size() << " " << VV.at(0).size() << endl;
0387                     for (int j = 0; j < 2; j++)
0388                         closest.push_back(VV[j][i]);
0389                     pre_phi = temp_phi;
0390                 }
0391             }
0392         }
0393         // closest = temp_closest;
0394         // for (unsigned int i = 0; i < VV.at(0).size(); i++)
0395         // {
0396         //     for (int j = 0; j < 2; j++)
0397         //         closest.push_back(VV[j][i]);
0398         // }
0399         // for (unsigned int i = 0; i < closest.size(); i++)
0400         //     cout << closest[i] << endl;
0401
0402         return closest;
0403     }
0404
0405     vector<double> get_crossframe(int Detect)
0406     {
0407         cout << endl;
0408         cout << "calc frame" << endl;
0409         double size;
0410         if (Detect == 0)
0411             size = 15;
0412         else if (Detect == 1)
0413             size = 120;
0414         vector<vector<double>> crossframe(2, vector<double>(0));
0415         for (int j = 0; j < 2; j++)
0416         {
0417             for (int i = 0; i < 2; i++)
0418             {
0419                 double a_ = (j == 0) ? 1 : 0;
0420                 double b_ = (j == 0) ? 0 : 1;
0421                 double c_ = (2 * i - 1) * size;
0422                 cout << "kesu : " << a_ << " " << b_ << " " << c_ << endl;
0423                 vector<double> co_ = {a_, b_, c_};
0424                 vector<double> temp_crossframe = get_crosspoint(co_);
0425                 if (temp_crossframe.size() == 2)
0426                 {
0427                     cout << "temp_crossframe : " << temp_crossframe[0] << "  " << temp_crossframe[1] << endl;
0428                     crossframe.at(0).push_back(temp_crossframe[0]);
0429                     crossframe.at(1).push_back(temp_crossframe[1]);
0430                 }
0431             }
0432         }
0433         cout << "crossframe" << endl;
0434         cout << "size : " << crossframe.size() << "  " << crossframe.at(0).size() << endl;
0435         for (int j = 0; j < crossframe.at(0).size(); j++)
0436         {
0437             for (int i = 0; i < crossframe.size(); i++)
0438             {
0439                 cout << crossframe[i][j] << ", ";
0440             }
0441             cout << endl;
0442         }
0443         vector<double> result;
0444         if (crossframe.at(0).size() != 0)
0445             result = get_closestpoint(crossframe);
0446         cout << "closest frame" << endl;
0447         for (int i = 0; i < result.size(); i++)
0448         {
0449             cout << result[i] << ", ";
0450         }
0451         cout << endl;
0452         cout << "calc end" << endl;
0453
0454         return result;
0455     }
0456
0457     /*void print()
0458     {
0459         cout << "p1("
0460              << setprecision(4) << setw(8) << p1.x() << ", "
0461              << setprecision(4) << setw(8) << p1.y() << ", "
0462              << setprecision(4) << setw(8) << p1.z() << ""
0463              << ") --- "
0464              << "p2("
0465              << setprecision(4) << setw(8) << p2.x() << ", "
0466              << setprecision(4) << setw(8) << p2.y() << ", "
0467              << setprecision(4) << setw(8) << p2.z() << ""
0468              << ") --- "
0469              << "DCA("
0470              << setprecision(4) << setw(8) << dca_mean[0] << ", "
0471              << setprecision(4) << setw(8) << dca_mean[1] << ", "
0472              << setprecision(4) << setw(8) << dca_mean[2] << ""
0473              << ")"
0474              << dca_range_out << is_deleted << "  "
0475              << "slope = " << track_rz->GetParameter(1) << "  Intercept = " << track_rz->GetParameter(0)
0476              << endl;
0477     }*/
0478 };
0479
0480 class truth
0481 {
0482 public:
0483     double m_truthpx;
0484     double m_truthpy;
0485     double m_truthpz;
0486     double m_truthpt;
0487     double m_trutheta;
0488     double m_truththeta;
0489     double m_truthphi;
0490     double m_truthz_out;
0491
0492     double m_xvtx;
0493     double m_yvtx;
0494     double m_zvtx;
0495
0496     int m_truthpid;
0497     int m_status;
0498
0499     TF1 *truth_xy = nullptr;
0500     TF1 *truth_rz = nullptr;
0501
0502     double dca_mean[3];
0503
0504     double center[2];
0505     double m_rad;
0506
0507     bool is_intt = false;
0508     bool is_charged = false;
0509     bool have_cluster = false;
0510     bool have_track = false;
0511
0512     void set_truth(double M_truthpx, double M_truthpy, double M_truthpz, double M_truthpt, int M_status, double M_trutheta, int M_truthpid, double M_truththeta, double M_truthphi, double M_xvtx, double M_yvtx, double M_zvtx)
0513     {
0514         m_truthpx = M_truthpx;
0515         m_truthpy = M_truthpy;
0516         m_truthpz = M_truthpz;
0517         m_truthpt = M_truthpt;
0518         m_status = M_status;
0519         m_trutheta = M_trutheta;
0520         m_truthpid = M_truthpid;
0521         m_truthphi = M_truthphi;
0522         m_truththeta = M_truththeta;
0523         m_xvtx = M_xvtx;
0524         m_yvtx = M_yvtx;
0525         m_zvtx = M_zvtx;
0526     }
0527
0528     double getzout()
0529     {
0530         double y = m_truthpy / fabs(m_truthpy) * 10;
0531         // cout << "truth_rz->GetParameter(0) : " << truth_rz->GetParameter(0) << endl;
0532         // cout<<"y : "<<y<<endl;
0533         double z_outer = (y - truth_rz->GetParameter(0)) / truth_rz->GetParameter(1);
0534         // cout << "z_outer : " << z_outer << endl;
0535
0536         return z_outer;
0537     }
0538
0539     double getpointr(int mode)
0540     {
0541         double r = 0;
0542         if (mode == 1) // truth
0543         {
0544             r = m_truthpy / fabs(m_truthpy) * m_truthpt;
0545         }
0546         else if (mode == 2) // dca_mean
0547         {
0548             r = dca_mean[1] / fabs(dca_mean[1]) * sqrt(dca_mean[0] * dca_mean[0] + dca_mean[1] * dca_mean[1]);
0549         }
0550         else if (mode == 3) // vertex
0551         {
0552             r = /*m_yvtx / fabs(m_yvtx) **/ sqrt(m_xvtx * m_xvtx + m_yvtx * m_yvtx);
0553         }
0554
0555         return r;
0556     }
0557
0558     void calc_line()
0559     {
0560         if (truth_xy != nullptr)
0561         {
0562             delete truth_xy;
0563             truth_xy = nullptr;
0564         }
0565
0566         truth_xy = new TF1("truth_xy", "pol1", -1000, 1000);
0567
0568         truth_xy->SetParameter(1, m_truthpy / m_truthpx);
0569         truth_xy->SetParameter(0, -m_yvtx / m_truthpx * m_xvtx + m_yvtx);
0570
0571         if (truth_rz != nullptr)
0572         {
0573             delete truth_rz;
0574             truth_rz = nullptr;
0575         }
0576
0577         truth_rz = new TF1("truth_rz", "pol1", -1000, 1000);
0578
0579         // from the truth vertex
0580         double slope = getpointr(1) / m_truthpz;
0581         double intercept = -slope * m_zvtx + getpointr(3);
0582         // cout << slope << "  " << m_zvtx << "  " << getpointr(3) << endl;
0583
0584         truth_rz->SetParameter(1, slope);
0585         truth_rz->SetParameter(0, intercept);
0586
0587         m_truthz_out = getzout();
0588     }
0589
0590     void calc_center()
0591     {
0592         m_rad = m_truthpt / (0.3 * 1.4) * 100;
0593         double sign = m_truthpid / fabs(m_truthpid);
0594         center[0] = -sign * (m_xvtx + m_rad * (-sin(m_truthphi)));
0595         center[1] = -sign * (m_yvtx + m_rad * cos(m_truthphi));
0596     }
0597
0598     /////////////////////////////////////////
0599
0600     // vector<double> get_crossline(double x1, double y1, double r1)
0601     // {
0602     //     vector<double> cross_co;
0603     //     // ax + by + c = 0
0604     //     cross_co.push_back(2 * (cx - x1));                                                           // a
0605     //     cross_co.push_back(2 * (cy - y1));                                                           // b
0606     //     cross_co.push_back((x1 + cx) * (x1 - cx) + (y1 + cy) * (y1 - cy) + (rad + r1) * (rad - r1)); // c
0607
0608     //     return cross_co;
0609     // }
0610
0611     // vector<double> get_crosspoint(const vector<double> &V) // ax + by + c = 0
0612     // {
0613     //     vector<double> cross;
0614     //     double a = V[0];
0615     //     double b = V[1];
0616     //     double c = V[2];
0617
0618     //     double d = fabs(a * cx + b * cy + c) / sqrt(a * a + b * b);
0619     //     cout << "d ; " << d << " " << rad << endl;
0620
0621     //     double theta = atan2(b, a);
0622
0623     //     if (d > rad)
0624     //     {
0625     //         cout << "d > rad" << endl;
0626     //     }
0627     //     else if (d == rad)
0628     //     {
0629     //         cout << "d == rad" << endl;
0630
0631     //         if (a * cx + b * cy + c > 0)
0632     //             theta += M_PI;
0633     //         cross.push_back(rad * cos(theta) + cx);
0634     //         cross.push_back(rad * sin(theta) + cy);
0635     //     }
0636     //     else
0637     //     {
0638     //         cout << "d < rad" << endl;
0639     //         double alpha, beta, phi;
0640     //         phi = acos(d / rad);
0641     //         alpha = theta - phi;
0642     //         beta = theta + phi;
0643     //         if (a * cx + b * cy + c > 0)
0644     //         {
0645     //             alpha += M_PI;
0646     //             beta += M_PI;
0647     //         }
0648     //         // double temp_cross[2][2];
0649     //         vector<vector<double>> temp_cross = {{rad * cos(alpha) + cx, rad * cos(beta) + cx}, {rad * sin(alpha) + cy, rad * sin(beta) + cy}};
0650     //         // for (unsigned int j = 0; j < temp_cross.at(0).size(); j++)
0651     //         // {
0652     //         //     cout << "temp_cross : ";
0653     //         //     for (unsigned int i = 0; i < temp_cross.size(); i++)
0654     //         //         cout << temp_cross.at(j).at(i) << "  ";
0655     //         //     cout << endl;
0656     //         // }
0657     //         cross = get_closestpoint(temp_cross);
0658     //     }
0659     //     // cout << "cross size : " << cross.size() << endl;
0660     //     // for (int i = 0; i < cross.size(); i++)
0661     //     //     cout << cross[i] << endl;
0662     //     return cross;
0663     // }
0664
0665     // vector<double> get_closestpoint(const vector<vector<double>> VV)
0666     // {
0667     //     vector<double> closest;
0668     //     double pre_phi = 999;
0669     //     double phi_p1 = atan2(p1.y() - dca_mean[1], p1.x() - dca_mean[0]);
0670     //     for (unsigned int i = 0; i < VV.at(0).size(); i++)
0671     //     {
0672     //         cout << "dca_mean : " << dca_mean[0] << "  " << dca_mean[1] << endl;
0673     //         cout << "VV : " << VV[0][i] << "  " << VV[1][i] << endl;
0674     //         cout << "pVV : " << VV[0][i] - dca_mean[0] << "  " << VV[1][i] - dca_mean[1] << endl;
0675     //         cout << "p1  : " << p1.x() - dca_mean[0] << "  " << p1.y() - dca_mean[1] << endl;
0676     //         // cout << VV.at(i).at(0) << endl;
0677     //         double dot = (VV[0][i] - dca_mean[0]) * (p1.x() - dca_mean[0]) + (VV[1][i] - dca_mean[1]) * (p1.y() - dca_mean[1]);
0678     //         double temp_phi = atan2(VV[1][i] - dca_mean[1], VV[0][i] - dca_mean[0]);
0679     //         cout << "dot : " << dot << endl;
0680     //         if (dot > 0)
0681     //         {
0682     //             if (fabs(temp_phi - phi_p1) < fabs(pre_phi - phi_p1))
0683     //             {
0684     //                 closest.erase(closest.begin(), closest.end());
0685     //                 cout << "VV size : " << VV.size() << " " << VV.at(0).size() << endl;
0686     //                 for (int j = 0; j < 2; j++)
0687     //                     closest.push_back(VV[j][i]);
0688     //                 pre_phi = temp_phi;
0689     //             }
0690     //         }
0691     //     }
0692     //     // closest = temp_closest;
0693     //     // for (unsigned int i = 0; i < VV.at(0).size(); i++)
0694     //     // {
0695     //     //     for (int j = 0; j < 2; j++)
0696     //     //         closest.push_back(VV[j][i]);
0697     //     // }
0698     //     // for (unsigned int i = 0; i < closest.size(); i++)
0699     //     //     cout << closest[i] << endl;
0700
0701     //     return closest;
0702     // }
0703
0704     // vector<double> get_crossframe(int Detect)
0705     // {
0706     //     cout << endl;
0707     //     cout << "calc frame" << endl;
0708     //     double size;
0709     //     if (Detect == 0)
0710     //         size = 15;
0711     //     else if (Detect == 1)
0712     //         size = 120;
0713     //     vector<vector<double>> crossframe(2, vector<double>(0));
0714     //     for (int j = 0; j < 2; j++)
0715     //     {
0716     //         for (int i = 0; i < 2; i++)
0717     //         {
0718     //             double a_ = (j == 0) ? 1 : 0;
0719     //             double b_ = (j == 0) ? 0 : 1;
0720     //             double c_ = (2 * i - 1) * size;
0721     //             cout << "kesu : " << a_ << " " << b_ << " " << c_ << endl;
0722     //             vector<double> co_ = {a_, b_, c_};
0723     //             vector<double> temp_crossframe = get_crosspoint(co_);
0724     //             if (temp_crossframe.size() == 2)
0725     //             {
0726     //                 cout << "temp_crossframe : " << temp_crossframe[0] << "  " << temp_crossframe[1] << endl;
0727     //                 crossframe.at(0).push_back(temp_crossframe[0]);
0728     //                 crossframe.at(1).push_back(temp_crossframe[1]);
0729     //             }
0730     //         }
0731     //     }
0732     //     cout << "crossframe" << endl;
0733     //     cout << "size : " << crossframe.size() << "  " << crossframe.at(0).size() << endl;
0734     //     for (int j = 0; j < crossframe.at(0).size(); j++)
0735     //     {
0736     //         for (int i = 0; i < crossframe.size(); i++)
0737     //         {
0738     //             cout << crossframe[i][j] << ", ";
0739     //         }
0740     //         cout << endl;
0741     //     }
0742     //     vector<double> result;
0743     //     if (crossframe.at(0).size() != 0)
0744     //         result = get_closestpoint(crossframe);
0745     //     cout << "closest frame" << endl;
0746     //     for (int i = 0; i < result.size(); i++)
0747     //     {
0748     //         cout << result[i] << ", ";
0749     //     }
0750     //     cout << endl;
0751     //     cout << "calc end" << endl;
0752
0753     //     return result;
0754     // }
0755 };