event_display/macro/track.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_zx;
0010   TF1 *track_zy;
0011
0012   TF1 *track_rz_inv;
0013
0014   double phi_in;
0015   double phi_out;
0016
0017   double theta_in;
0018   double theta_out;
0019
0020   double d_phi;
0021   double d_theta;
0022
0023   double dca[3];
0024   double dca_mean[3];
0025
0026   vector < int > cluster_ids;
0027   TVector3 p1;
0028   TVector3 p2;
0029
0030   double dca_2d;
0031   double phi;
0032   double theta;
0033   double phi_tracklet;
0034   double theta_tracklet;
0035
0036   double charge;
0037   vector < vector < double > > zline;
0038
0039   vector<double> x;
0040   vector<double> y;
0041   vector<double> z;
0042   vector<double> r;
0043   vector<double> x_err;
0044   vector<double> y_err;
0045   vector<double> z_err;
0046
0047   double pT;
0048   double rad;
0049   double cx, cy;
0050   double p_reco[3];
0051
0052   double track_theta; // in z-r plane, degree
0053
0054   double A[3][3];
0055   double B[1][3];
0056   double beta[1][3];
0057
0058   bool dca2d_range_out = false;
0059   bool dcaz_range_out = false;
0060   bool dca_range_out = false;
0061   bool is_deleted = false;
0062
0063   track() : track_xy(nullptr),
0064         track_rz(nullptr),
0065         track_zx(nullptr),
0066         track_zy(nullptr),
0067         track_rz_inv(nullptr),
0068         zline(2, vector<double>(0))
0069
0070   {
0071     p1 = TVector3();
0072     p2 = TVector3();
0073   };
0074
0075   virtual ~track()
0076   {
0077     if (track_xy != nullptr)
0078       delete track_xy;
0079     if (track_rz != nullptr)
0080       delete track_rz;
0081     if (track_zx != nullptr)
0082       delete track_zx;
0083     if (track_zy != nullptr)
0084       delete track_zy;
0085     if (track_rz_inv != nullptr)
0086       delete track_rz_inv;
0087   };
0088
0089   void calc_pT()
0090   {
0091     TVector3 p3 = {dca_mean[0], dca_mean[1], dca_mean[2]}; // vertex
0092
0093     TVector3 mid1 = TVector3((p1.x() + p2.x()) / 2, (p1.y() + p2.y()) / 2, (p1.z() + p2.z()) / 2);
0094     TVector3 mid2 = TVector3((p3.x() + p1.x()) / 2, (p3.y() + p1.y()) / 2, (p3.z() + p1.z()) / 2);
0095
0096     double alpha = atan2((p2.y() - p1.y()), (p2.x() - p1.x())) + M_PI / 2;
0097     double beta = atan2((p1.y() - p3.y()), (p1.x() - p3.x())) + M_PI / 2;
0098
0099     cx = (mid1.y() - mid2.y() - mid1.x() * tan(alpha) + mid2.x() * tan(beta)) / (tan(beta) - tan(alpha));
0100     cy = mid1.y() - (mid1.x() - cx) * tan(alpha);
0101
0102     rad = sqrt((p1.x() - cx) * (p1.x() - cx) + (p1.y() - cy) * (p1.y() - cy)); //[cm]
0103     pT = /*charge **/ 0.3 * 1.4 * rad / 100;                                   //[GeV/c]
0104
0105     double phi_p = atan2(p3.y() - cy, p3.x() - cx) + (M_PI / 2);
0106     double phi_p1 = atan2(p1.y() - cy, p1.x() - cx) + (M_PI / 2);
0107     // cout<<"p1 : "<<p1.x()<<"  "<<p1.y()<<endl;
0108     // cout<<"phi : "<<phi_p/M_PI * 180<<"  "<<phi_p1/M_PI * 180<<endl;
0109     double d_phi = phi_p1 - phi_p;
0110     // cout<<"d_phi : "<<d_phi/M_PI * 180<<endl;
0111     // d_phi -= -2 * (M_PI * (int)(d_phi / (M_PI)));
0112
0113     // if(fabs(d_phi) > M_PI)
0114     // {
0115     //     cout<<"fabs(d_phi) > M_PI  "<<d_phi<<endl;
0116     //     cout<<(M_PI * (int)(d_phi / (M_PI)))<<endl;
0117     // // d_phi -=  (M_PI * (int)(d_phi / (M_PI)));
0118     // d_phi -=  2 * (M_PI * (int)(d_phi / (M_PI)));
0119
0120     //     // d_phi -=d_phi * fabs(d_phi) *  M_PI * (int)(d_phi / M_PI);
0121     // }
0122     // cout<<" d_phi : "<<d_phi/M_PI * 180<<endl;
0123
0124     double phi_offset = 0;
0125     double cross = p3.x() * p1.y() - p3.y() * p1.x();
0126
0127     if ((d_phi) < 0)
0128       // if((phi_p1 - phi_p) < 0)
0129       phi_offset = M_PI;
0130     // double sign = d_phi / fabs(d_phi);
0131
0132     // double sign = ((phi_p1 -(M_PI/2))  - (phi_p - (M_PI/2)))/fabs((phi_p1 - (M_PI/2)) - (phi_p - (M_PI/2)));
0133     // cout<<sign<<endl;
0134     // double sign_rz = (p1.z() - dca_mean[2]) / fabs(p1.z() - dca_mean[2]);
0135     // double sign_x = (p1.x() - p3.x()) / fabs(p1.x() - p3.x());
0136     // double sign_y = (p1.y() - p3.y()) / fabs(p1.y() - p3.y());
0137     // cout<<"  sign_x,y : "<<sign_x<<"  "<<sign_y<<endl;
0138     double t_ = sqrt(p1.x() * p1.x() + p1.y() * p1.y());
0139     // cout<<"no_sign : "<<fabs(pT) * cos(phi_p)<<"  "<<fabs(pT) * sin(phi_p)<<endl;
0140     // cout<<"phi : "<<phi_p + phi_offset<<endl;
0141
0142     p_reco[0] = /*sign_x **/ fabs(pT) * cos(phi_p + phi_offset);
0143     p_reco[1] = /*sign_y **/ fabs(pT) * sin(phi_p + phi_offset);
0144     // cout<<"p1 : "<<p1.x() - dca_mean[0]<<"  "<<p1.y() - dca_mean[1]<<endl;
0145     // cout<<"p  : "<<p_reco[0]<<"  "<<p_reco[1]<<endl;
0146     p_reco[2] = fabs(pT) / t_ * p1.z();
0147     // cout<<endl;
0148   }
0149
0150   void calc_line()
0151   {
0152     x.push_back(p1.x());
0153     x.push_back(p2.x());
0154     x.push_back(dca_mean[0]);
0155
0156     x_err.push_back(1.);
0157     x_err.push_back(1.);
0158     x_err.push_back(1.);
0159
0160     y.push_back(p1.y());
0161     y.push_back(p2.y());
0162     y.push_back(dca_mean[1]);
0163
0164     y_err.push_back(1.);
0165     y_err.push_back(1.);
0166     y_err.push_back(1.);
0167
0168     LeastSquare least_square;
0169     least_square.Setdatax(x);
0170     least_square.Setdatay(y);
0171     least_square.Setdatayerr(y_err);
0172
0173     least_square.SetDebugMode(false);
0174     least_square.Calc(0);
0175
0176     if (track_xy != nullptr)
0177       {
0178     delete track_xy;
0179     track_xy = nullptr;
0180       }
0181
0182     track_xy = new TF1("track_xy", "pol1", -15, 15);
0183
0184     track_xy->SetParameter(1, least_square.GetSlope());
0185     track_xy->SetParameter(0, least_square.GetIntercept());
0186
0187     z.push_back(p1.z());
0188     z.push_back(p2.z());
0189     z.push_back(dca_mean[2]);
0190
0191     r.push_back(getpointr(1));
0192     r.push_back(getpointr(2));
0193     r.push_back(getpointr(3));
0194
0195     z_err.push_back(20. / sqrt(12.)); // [mm]
0196     z_err.push_back(20. / sqrt(12.)); // [mm]
0197     z_err.push_back(0.8);             // [mm]
0198
0199     ///////////////////////////////////////////
0200     // r-z plane                             //
0201     ///////////////////////////////////////////
0202     LeastSquare least_square_zr;
0203     least_square_zr.Setdatax(r);
0204     least_square_zr.Setdatay(z);
0205     least_square_zr.Setdatayerr(z_err);
0206
0207     // least_square_zr.SetDebugMode(false);
0208     least_square_zr.Calc(1);
0209
0210     if (track_rz != nullptr)
0211       {
0212     delete track_rz;
0213     track_rz = nullptr;
0214       }
0215     double slope_inv = least_square_zr.GetSlope();
0216     double intercept_inv = least_square_zr.GetIntercept();
0217
0218     track_rz = new TF1("track_rz", "pol1", -25, 25);
0219
0220     track_rz->SetParameter(1, getslope_inv(slope_inv, intercept_inv));
0221     track_rz->SetParameter(0, getintercept_inv(slope_inv, intercept_inv));
0222
0223     track_theta = 360.0 * TMath::ATan( track_rz->GetParameter(1) ) / TMath::Pi() / 2 ; // degree
0224
0225     // cout << "Track z-r: " << track_rz->GetParameter( 1 ) << "\t"
0226     //   << track_theta << endl; // values look OK
0227     //        track_rz->SetParameter(1, least_square_zr.GetSlope());
0228     // track_rz->SetParameter(0, least_square_zr.GetIntercept());
0229
0230     /*
0231 ///////////////////////////////////////////
0232 // z-x plane                             //
0233 ///////////////////////////////////////////
0234 LeastSquare least_square_zx;
0235 least_square_zx.Setdatax(z);
0236 least_square_zx.Setdatay(x);
0237 least_square_zx.Setdatayerr(x_err);
0238
0239 // least_square_zx.SetDebugMode(false);
0240 least_square_zx.Calc(1);
0241
0242 if (track_zx != nullptr)
0243 {
0244 delete track_zx;
0245 track_zx = nullptr;
0246 }
0247 double slope_inv_zx = least_square_zx.GetSlope();
0248 double intercept_inv_zx = least_square_zx.GetIntercept();
0249
0250 track_zx = new TF1("track_zx", "pol1", -25, 25);
0251
0252 track_zx->SetParameter(1, getslope_inv(slope_inv_zx, intercept_inv_zx));
0253 track_zx->SetParameter(0, getintercept_inv(slope_inv_zx, intercept_inv_zx));
0254
0255 ///////////////////////////////////////////
0256 // z-y plane                             //
0257 ///////////////////////////////////////////
0258 LeastSquare least_square_zy;
0259 least_square_zy.Setdatax(z);
0260 least_square_zy.Setdatay(y);
0261 least_square_zy.Setdatayerr(y_err);
0262
0263 // least_square_zy.SetDebugMode(false);
0264 least_square_zy.Calc(1);
0265
0266 if (track_zy != nullptr)
0267 {
0268 delete track_zy;
0269 track_zy = nullptr;
0270 }
0271 double slope_inv_zy = least_square_zy.GetSlope();
0272 double intercept_inv_zy = least_square_zy.GetIntercept();
0273
0274 track_zy = new TF1("track_zy", "pol1", -25, 25);
0275
0276 track_zy->SetParameter(1, getslope_inv(slope_inv_zy, intercept_inv_zy));
0277 track_zy->SetParameter(0, getintercept_inv(slope_inv_zy, intercept_inv_zy));
0278
0279 cout << "parameter calc" << endl;
0280 cout << "slope in x-y: "
0281 << track_xy->GetParameter(1) << "\t"
0282 << track_zy->GetParameter(1) / track_zx->GetParameter(1) << "\t"
0283 << 100 * (track_xy->GetParameter(1) -  track_zy->GetParameter(1) / track_zx->GetParameter(1) ) / track_xy->GetParameter(1) <<  "%"
0284 << endl;
0285     */
0286   }
0287
0288   void calc_inv()
0289   {
0290     if (track_rz_inv != nullptr)
0291       {
0292     delete track_rz_inv;
0293     track_rz_inv = nullptr;
0294       }
0295
0296     double slope_ = 1 / track_rz->GetParameter(1);
0297     double intercept_ = -track_rz->GetParameter(0) / track_rz->GetParameter(1);
0298
0299     track_rz_inv = new TF1("track_rz_inv", "pol1", -25, 25);
0300
0301     track_rz_inv->SetParameter(1, slope_);
0302     track_rz_inv->SetParameter(0, intercept_);
0303   }
0304
0305   double getslope_inv(double Slope, double Intercept)
0306   {
0307     return 1 / Slope;
0308   }
0309
0310   double getintercept_inv(double Slope, double Intercept)
0311   {
0312     return -Intercept / Slope;
0313   }
0314
0315   double getpointr(int mode)
0316   {
0317     double r = 0;
0318
0319     if (mode == 1)
0320       {
0321     r = (p1.y() / fabs(p1.y())) * sqrt(p1.x() * p1.x() + p1.y() * p1.y()); // returns +r or -r, the sign is the sign of y
0322       }
0323     else if (mode == 2)
0324       {
0325     r = (p2.y() / fabs(p2.y())) * sqrt(p2.x() * p2.x() + p2.y() * p2.y());
0326       }
0327     else if (mode == 3)
0328       {
0329     r = (dca_mean[1] / fabs(dca_mean[1])) * sqrt(dca_mean[0] * dca_mean[0] + dca_mean[1] * dca_mean[1]);
0330       }
0331
0332     return r;
0333   }
0334
0335   double getphi()
0336   {
0337     phi = atan(track_xy->GetParameter(1));
0338
0339     if (p1.x() < dca_mean[0])
0340       {
0341     if (p1.y() < dca_mean[1])
0342       phi -= M_PI;
0343     else
0344       phi += M_PI;
0345       }
0346
0347     return phi;
0348   }
0349
0350   double gettheta()
0351   {
0352     theta = atan(track_rz->GetParameter(1));
0353
0354     if (p1.z() < dca_mean[2])
0355       {
0356     if (getpointr(1) < getpointr(3))
0357       theta -= M_PI;
0358     else
0359       theta += M_PI;
0360       }
0361
0362     if ( getpointr(1) < getpointr(3) )
0363       theta = fabs(theta);
0364
0365     return theta;
0366   }
0367
0368   double getphi_tracklet()
0369   {
0370     phi_tracklet = atan2((p2.y() - p1.y()), (p2.x() - p1.x()));
0371     return phi_tracklet;
0372   }
0373
0374   double gettheta_tracklet()
0375   {
0376     theta_tracklet = atan2((fabs(getpointr(2)) - fabs(getpointr(1))), (p2.z() - p1.z()));
0377
0378     return theta_tracklet;
0379   }
0380
0381   vector<double> get_crossline(double x1, double y1, double r1)
0382   {
0383     vector<double> cross_co;
0384     // ax + by + c = 0
0385     cross_co.push_back(2 * (cx - x1));                                                           // a
0386     cross_co.push_back(2 * (cy - y1));                                                           // b
0387     cross_co.push_back((x1 + cx) * (x1 - cx) + (y1 + cy) * (y1 - cy) + (rad + r1) * (rad - r1)); // c
0388
0389     return cross_co;
0390   }
0391
0392   vector<double> get_crosspoint(const vector<double> &V) // ax + by + c = 0
0393   {
0394     vector<double> cross;
0395     double a = V[0];
0396     double b = V[1];
0397     double c = V[2];
0398
0399     double d = fabs(a * cx + b * cy + c) / sqrt(a * a + b * b);
0400     // cout << "d ; " << d << " " << rad << endl;
0401
0402     double theta = atan2(b, a);
0403
0404     if (d > rad)
0405       {
0406     cout << "d > rad" << endl;
0407       }
0408     else if (d == rad)
0409       {
0410     cout << "d == rad" << endl;
0411
0412     if (a * cx + b * cy + c > 0)
0413       theta += M_PI;
0414     cross.push_back(rad * cos(theta) + cx);
0415     cross.push_back(rad * sin(theta) + cy);
0416       }
0417     else
0418       {
0419     cout << "d < rad" << endl;
0420     double alpha, beta, phi;
0421     phi = acos(d / rad);
0422     alpha = theta - phi;
0423     beta = theta + phi;
0424     if (a * cx + b * cy + c > 0)
0425       {
0426         alpha += M_PI;
0427         beta += M_PI;
0428       }
0429     // double temp_cross[2][2];
0430     vector<vector<double>> temp_cross = {{rad * cos(alpha) + cx, rad * cos(beta) + cx}, {rad * sin(alpha) + cy, rad * sin(beta) + cy}};
0431     // for (unsigned int j = 0; j < temp_cross.at(0).size(); j++)
0432     // {
0433     //     cout << "temp_cross : ";
0434     //     for (unsigned int i = 0; i < temp_cross.size(); i++)
0435     //         cout << temp_cross.at(j).at(i) << "  ";
0436     //     cout << endl;
0437     // }
0438     cross = get_closestpoint(temp_cross);
0439       }
0440     // cout << "cross size : " << cross.size() << endl;
0441     // for (int i = 0; i < cross.size(); i++)
0442     //     cout << cross[i] << endl;
0443     return cross;
0444   }
0445
0446   vector<double> get_closestpoint(const vector<vector<double>> VV)
0447   {
0448     vector<double> closest;
0449     double pre_phi = 999;
0450     double phi_p1 = atan2(p1.y() - dca_mean[1], p1.x() - dca_mean[0]);
0451     for (unsigned int i = 0; i < VV.at(0).size(); i++)
0452       {
0453     cout << "dca_mean : " << dca_mean[0] << "  " << dca_mean[1] << endl;
0454     cout << "VV : " << VV[0][i] << "  " << VV[1][i] << endl;
0455     cout << "pVV : " << VV[0][i] - dca_mean[0] << "  " << VV[1][i] - dca_mean[1] << endl;
0456     cout << "p1  : " << p1.x() - dca_mean[0] << "  " << p1.y() - dca_mean[1] << endl;
0457     // cout << VV.at(i).at(0) << endl;
0458     double dot = (VV[0][i] - dca_mean[0]) * (p1.x() - dca_mean[0]) + (VV[1][i] - dca_mean[1]) * (p1.y() - dca_mean[1]);
0459     double temp_phi = atan2(VV[1][i] - dca_mean[1], VV[0][i] - dca_mean[0]);
0460     cout << "dot : " << dot << endl;
0461     if (dot > 0)
0462       {
0463         if (fabs(temp_phi - phi_p1) < fabs(pre_phi - phi_p1))
0464           {
0465         closest.erase(closest.begin(), closest.end());
0466         cout << "VV size : " << VV.size() << " " << VV.at(0).size() << endl;
0467         for (int j = 0; j < 2; j++)
0468           closest.push_back(VV[j][i]);
0469         pre_phi = temp_phi;
0470           }
0471       }
0472       }
0473     // closest = temp_closest;
0474     // for (unsigned int i = 0; i < VV.at(0).size(); i++)
0475     // {
0476     //     for (int j = 0; j < 2; j++)
0477     //         closest.push_back(VV[j][i]);
0478     // }
0479     // for (unsigned int i = 0; i < closest.size(); i++)
0480     //     cout << closest[i] << endl;
0481
0482     return closest;
0483   }
0484
0485   vector<double> get_crossframe(int Detect)
0486   {
0487     cout << endl;
0488     cout << "calc frame" << endl;
0489     double size;
0490     if (Detect == 0)
0491       size = 15;
0492     else if (Detect == 1)
0493       size = 120;
0494     vector<vector<double>> crossframe(2, vector<double>(0));
0495     for (int j = 0; j < 2; j++)
0496       {
0497     for (int i = 0; i < 2; i++)
0498       {
0499         double a_ = (j == 0) ? 1 : 0;
0500         double b_ = (j == 0) ? 0 : 1;
0501         double c_ = (2 * i - 1) * size;
0502         cout << "kesu : " << a_ << " " << b_ << " " << c_ << endl;
0503         vector<double> co_ = {a_, b_, c_};
0504         vector<double> temp_crossframe = get_crosspoint(co_);
0505         if (temp_crossframe.size() == 2)
0506           {
0507         cout << "temp_crossframe : " << temp_crossframe[0] << "  " << temp_crossframe[1] << endl;
0508         crossframe.at(0).push_back(temp_crossframe[0]);
0509         crossframe.at(1).push_back(temp_crossframe[1]);
0510           }
0511       }
0512       }
0513     cout << "crossframe" << endl;
0514     cout << "size : " << crossframe.size() << "  " << crossframe.at(0).size() << endl;
0515     for (int j = 0; j < crossframe.at(0).size(); j++)
0516       {
0517     for (int i = 0; i < crossframe.size(); i++)
0518       {
0519         cout << crossframe[i][j] << ", ";
0520       }
0521     cout << endl;
0522       }
0523     vector<double> result;
0524     if (crossframe.at(0).size() != 0)
0525       result = get_closestpoint(crossframe);
0526     cout << "closest frame" << endl;
0527     for (int i = 0; i < result.size(); i++)
0528       {
0529     cout << result[i] << ", ";
0530       }
0531     cout << endl;
0532     cout << "calc end" << endl;
0533
0534     return result;
0535   }
0536
0537   void print()
0538     {
0539       cout << "p1("
0540        << setprecision(4) << setw(8) << p1.x() << ", "
0541        << setprecision(4) << setw(8) << p1.y() << ", "
0542        << setprecision(4) << setw(8) << p1.z() << ""
0543        << ") --- "
0544        << "p2("
0545        << setprecision(4) << setw(8) << p2.x() << ", "
0546        << setprecision(4) << setw(8) << p2.y() << ", "
0547        << setprecision(4) << setw(8) << p2.z() << ""
0548        << ") --- "
0549        << "DCA("
0550        << setprecision(4) << setw(8) << dca_mean[0] << ", "
0551        << setprecision(4) << setw(8) << dca_mean[1] << ", "
0552        << setprecision(4) << setw(8) << dca_mean[2] << ""
0553        << ")"
0554        << dca_range_out << is_deleted << "  "
0555        << "slope = " << track_rz->GetParameter(1)
0556        << "  Intercept = " << track_rz->GetParameter(0) << " "
0557        << "theta = " << track_theta
0558        << endl;
0559     }
0560
0561
0562   double GetTrackTheta()
0563   {
0564     return track_theta;
0565   }
0566
0567   double GetTrackPhi( double r )
0568   {
0569     cerr << "double GetTrackPhi( double r ) is under construction" << endl;
0570     return 0;
0571   }
0572
0573   void AddAssociatedClusterId( int i )
0574   {
0575     cluster_ids.push_back( i );
0576   }
0577
0578 };