202409_performance/z_vtx/track_pT.hh

0001 #pragma once
0002 #include "least_square2.cc"
0003
0004 class track
0005 {
0006 public:
0007     TF1 *track_xy;
0008     TF1 *track_rz;
0009     TF1 *track_rz_inv;
0010     double phi_in;
0011     double phi_out;
0012     double theta_in;
0013     double theta_out;
0014     double d_phi;
0015     double d_theta;
0016     double dca[5];
0017     double dca_mean[3];
0018     TVector3 p1;
0019     TVector3 p2;
0020     double dca_2d;
0021     double phi;
0022     double theta;
0023     double phi_tracklet;
0024     double theta_tracklet;
0025     double charge;
0026
0027     int ladder1;
0028     int ladder2;
0029
0030     vector<double> x;
0031     vector<double> y;
0032     vector<double> z;
0033     vector<double> r;
0034     vector<double> y_err;
0035     vector<double> z_err;
0036
0037     double pT;
0038     double rad;
0039     double cx, cy;
0040     double p_reco[3];
0041
0042     double A[3][3];
0043     double B[1][3];
0044     double beta[1][3];
0045
0046     bool dca2d_range_out = false;
0047     bool dcaz_range_out = false;
0048     bool dca_range_out = false;
0049     bool is_deleted = false;
0050
0051     track() : track_xy(nullptr),
0052               track_rz(nullptr),
0053               track_rz_inv(nullptr)
0054
0055     {
0056         p1 = TVector3();
0057         p2 = TVector3();
0058     };
0059
0060     virtual ~track()
0061     {
0062         if (track_xy != nullptr)
0063             delete track_xy;
0064         if (track_rz != nullptr)
0065             delete track_rz;
0066         if (track_rz_inv != nullptr)
0067             delete track_rz_inv;
0068     };
0069
0070     void calc_pT()
0071     {
0072         TVector3 p3 = {dca_mean[0], dca_mean[1], dca_mean[2]}; // vertex
0073
0074         TVector3 mid1 = TVector3((p1.x() + p2.x()) / 2, (p1.y() + p2.y()) / 2, (p1.z() + p2.z()) / 2);
0075         TVector3 mid2 = TVector3((p3.x() + p1.x()) / 2, (p3.y() + p1.y()) / 2, (p3.z() + p1.z()) / 2);
0076
0077         double alpha = atan2((p2.y() - p1.y()), (p2.x() - p1.x())) + M_PI / 2;
0078         double beta = atan2((p1.y() - p3.y()), (p1.x() - p3.x())) + M_PI / 2;
0079
0080         cx = (mid1.y() - mid2.y() - mid1.x() * tan(alpha) + mid2.x() * tan(beta)) / (tan(beta) - tan(alpha));
0081         cy = mid1.y() - (mid1.x() - cx) * tan(alpha);
0082
0083         rad = sqrt((p1.x() - cx) * (p1.x() - cx) + (p1.y() - cy) * (p1.y() - cy)); //[cm]
0084         pT = /*charge **/ 0.3 * 1.4 * rad / 100;                                   //[GeV/c]
0085
0086         double phi_p = atan2(p3.y() - cy, p3.x() - cx) + (M_PI / 2);
0087         double phi_p1 = atan2(p1.y() - cy, p1.x() - cx) + (M_PI / 2);
0088         // cout<<"p1 : "<<p1.x()<<"  "<<p1.y()<<endl;
0089         // cout<<"phi : "<<phi_p/M_PI * 180<<"  "<<phi_p1/M_PI * 180<<endl;
0090         double d_phi = phi_p1 - phi_p;
0091         // cout << "d_phi : " << d_phi / M_PI * 180 << endl;
0092         // d_phi -= -2 * (M_PI * (int)(d_phi / (M_PI)));
0093
0094         if (fabs(d_phi) > M_PI)
0095         {
0096             // cout << "fabs(d_phi) > M_PI  " << d_phi << endl;
0097             // cout << (M_PI * (int)(d_phi / (M_PI))) << endl;
0098             // d_phi -=  (M_PI * (int)(d_phi / (M_PI)));
0099             d_phi -= 2 * (M_PI * (int)(d_phi / (M_PI)));
0100
0101             // d_phi -=d_phi * fabs(d_phi) *  M_PI * (int)(d_phi / M_PI);
0102         }
0103         // cout << " d_phi : " << d_phi / M_PI * 180 << endl;
0104
0105         double phi_offset = 0;
0106         double cross = p3.x() * p1.y() - p3.y() * p1.x();
0107
0108         if ((d_phi) < 0)
0109             // if((phi_p1 - phi_p) < 0)
0110             phi_offset = M_PI;
0111         // double sign = d_phi / fabs(d_phi);
0112
0113         // double sign = ((phi_p1 -(M_PI/2))  - (phi_p - (M_PI/2)))/fabs((phi_p1 - (M_PI/2)) - (phi_p - (M_PI/2)));
0114         // cout<<sign<<endl;
0115         // double sign_rz = (p1.z() - dca_mean[2]) / fabs(p1.z() - dca_mean[2]);
0116         // double sign_x = (p1.x() - p3.x()) / fabs(p1.x() - p3.x());
0117         // double sign_y = (p1.y() - p3.y()) / fabs(p1.y() - p3.y());
0118         // cout<<"  sign_x,y : "<<sign_x<<"  "<<sign_y<<endl;
0119         double t_ = sqrt(p1.x() * p1.x() + p1.y() * p1.y());
0120         // cout<<"no_sign : "<<fabs(pT) * cos(phi_p)<<"  "<<fabs(pT) * sin(phi_p)<<endl;
0121         // cout << "phi : " << phi_p + phi_offset << endl;
0122
0123         p_reco[0] = /*sign_x **/ fabs(pT) * cos(phi_p + phi_offset);
0124         p_reco[1] = /*sign_y **/ fabs(pT) * sin(phi_p + phi_offset);
0125         // cout << "p1 : " << p1.x() - dca_mean[0] << "  " << p1.y() - dca_mean[1] << endl;
0126         // cout << "p  : " << p_reco[0] << "  " << p_reco[1] << endl;
0127         p_reco[2] = fabs(pT) / t_ * p1.z();
0128         // cout << endl;
0129     }
0130
0131     void calc_line()
0132     {
0133         x.push_back(p1.x());
0134         x.push_back(p2.x());
0135         x.push_back(dca_mean[0]);
0136
0137         y.push_back(p1.y());
0138         y.push_back(p2.y());
0139         y.push_back(dca_mean[1]);
0140
0141         y_err.push_back(1.);
0142         y_err.push_back(1.);
0143         y_err.push_back(1.);
0144
0145         LeastSquare least_square;
0146         least_square.Setdatax(x);
0147         least_square.Setdatay(y);
0148         least_square.Setdatayerr(y_err);
0149
0150         least_square.SetDebugMode(false);
0151         least_square.Calc(0);
0152
0153         if (track_xy != nullptr)
0154         {
0155             delete track_xy;
0156             track_xy = nullptr;
0157         }
0158
0159         track_xy = new TF1("track_xy", "pol1", -15, 15);
0160
0161         track_xy->SetParameter(1, least_square.GetSlope());
0162         track_xy->SetParameter(0, least_square.GetIntercept());
0163
0164         z.push_back(p1.z());
0165         z.push_back(p2.z());
0166         z.push_back(dca_mean[2]);
0167
0168         r.push_back(getpointr(1));
0169         r.push_back(getpointr(2));
0170         r.push_back(getpointr(3));
0171
0172         z_err.push_back(20. / sqrt(12.)); // [mm]
0173         z_err.push_back(20. / sqrt(12.)); // [mm]
0174         z_err.push_back(0.8);             // [mm]
0175
0176         // r-z plane
0177         LeastSquare least_square_zr;
0178         least_square_zr.Setdatax(r);
0179         least_square_zr.Setdatay(z);
0180         least_square_zr.Setdatayerr(z_err);
0181
0182         // least_square_zr.SetDebugMode(false);
0183         least_square_zr.Calc(1);
0184
0185         if (track_rz != nullptr)
0186         {
0187             delete track_rz;
0188             track_rz = nullptr;
0189         }
0190         double slope_inv = least_square_zr.GetSlope();
0191         double intercept_inv = least_square_zr.GetIntercept();
0192
0193         track_rz = new TF1("track_rz", "pol1", -25, 25);
0194
0195         track_rz->SetParameter(1, getslope_inv(slope_inv, intercept_inv));
0196         track_rz->SetParameter(0, getintercept_inv(slope_inv, intercept_inv));
0197
0198         //        track_rz->SetParameter(1, least_square_zr.GetSlope());
0199         // track_rz->SetParameter(0, least_square_zr.GetIntercept());
0200     }
0201
0202     void calc_inv()
0203     {
0204         if (track_rz_inv != nullptr)
0205         {
0206             delete track_rz_inv;
0207             track_rz_inv = nullptr;
0208         }
0209
0210         double slope_ = 1 / track_rz->GetParameter(1);
0211         double intercept_ = -track_rz->GetParameter(0) / track_rz->GetParameter(1);
0212
0213         track_rz_inv = new TF1("track_rz_inv", "pol1", -25, 25);
0214
0215         track_rz_inv->SetParameter(1, slope_);
0216         track_rz_inv->SetParameter(0, intercept_);
0217     }
0218
0219     double getslope_inv(double Slope, double Intercept)
0220     {
0221         return 1 / Slope;
0222     }
0223
0224     double getintercept_inv(double Slope, double Intercept)
0225     {
0226         return -Intercept / Slope;
0227     }
0228
0229     double getpointr(int mode)
0230     {
0231         double r = 0;
0232
0233         if (mode == 1)
0234         {
0235             r = (p1.y() / fabs(p1.y())) * sqrt(p1.x() * p1.x() + p1.y() * p1.y());
0236         }
0237         else if (mode == 2)
0238         {
0239             r = (p2.y() / fabs(p2.y())) * sqrt(p2.x() * p2.x() + p2.y() * p2.y());
0240         }
0241         else if (mode == 3)
0242         {
0243             r = (dca_mean[1] / fabs(dca_mean[1])) * sqrt(dca_mean[0] * dca_mean[0] + dca_mean[1] * dca_mean[1]);
0244         }
0245
0246         return r;
0247     }
0248
0249     double getphi()
0250     {
0251         phi = atan(track_xy->GetParameter(1));
0252
0253         if (p1.x() < dca_mean[0])
0254         {
0255             if (p1.y() < dca_mean[1])
0256                 phi -= M_PI;
0257             else
0258                 phi += M_PI;
0259         }
0260
0261         return phi;
0262     }
0263
0264     double gettheta()
0265     {
0266         theta = atan(track_rz->GetParameter(1));
0267
0268         if (p1.z() < dca_mean[2])
0269         {
0270             if (getpointr(1) < getpointr(3))
0271                 theta -= M_PI;
0272             else
0273                 theta += M_PI;
0274         }
0275
0276         if (getpointr(1) < getpointr(3))
0277             theta = fabs(theta);
0278
0279         return theta;
0280     }
0281
0282     double getphi_tracklet()
0283     {
0284         phi_tracklet = atan2((p2.y() - p1.y()), (p2.x() - p1.x()));
0285         return phi_tracklet;
0286     }
0287
0288     double gettheta_tracklet()
0289     {
0290         theta_tracklet = atan2((fabs(getpointr(2)) - fabs(getpointr(1))), (p2.z() - p1.z()));
0291
0292         return theta_tracklet;
0293     }
0294
0295     /*void print()
0296     {
0297         cout << "p1("
0298              << setprecision(4) << setw(8) << p1.x() << ", "
0299              << setprecision(4) << setw(8) << p1.y() << ", "
0300              << setprecision(4) << setw(8) << p1.z() << ""
0301              << ") --- "
0302              << "p2("
0303              << setprecision(4) << setw(8) << p2.x() << ", "
0304              << setprecision(4) << setw(8) << p2.y() << ", "
0305              << setprecision(4) << setw(8) << p2.z() << ""
0306              << ") --- "
0307              << "DCA("
0308              << setprecision(4) << setw(8) << dca_mean[0] << ", "
0309              << setprecision(4) << setw(8) << dca_mean[1] << ", "
0310              << setprecision(4) << setw(8) << dca_mean[2] << ""
0311              << ")"
0312              << dca_range_out << is_deleted << "  "
0313              << "slope = " << track_rz->GetParameter(1) << "  Intercept = " << track_rz->GetParameter(0)
0314              << endl;
0315     }*/
0316 };
0317
0318 class truth
0319 {
0320 public:
0321     double m_truthpx;
0322     double m_truthpy;
0323     double m_truthpz;
0324     double m_truthpt;
0325     double m_trutheta;
0326     double m_truththeta;
0327     double m_truthphi;
0328     double m_truthz_out;
0329
0330     double m_xvtx;
0331     double m_yvtx;
0332     double m_zvtx;
0333
0334     int m_truthpid;
0335     int m_status;
0336
0337     TF1 *truth_xy = nullptr;
0338     TF1 *truth_rz = nullptr;
0339
0340     double dca_mean[3];
0341
0342     double center[2];
0343     double m_rad;
0344
0345     bool is_intt = false;
0346     bool is_charged = false;
0347     bool have_cluster = false;
0348     bool have_track = false;
0349
0350     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)
0351     {
0352         m_truthpx = M_truthpx;
0353         m_truthpy = M_truthpy;
0354         m_truthpz = M_truthpz;
0355         m_truthpt = M_truthpt;
0356         m_status = M_status;
0357         m_trutheta = M_trutheta;
0358         m_truthpid = M_truthpid;
0359         m_truthphi = M_truthphi;
0360         m_truththeta = M_truththeta;
0361         m_xvtx = M_xvtx;
0362         m_yvtx = M_yvtx;
0363         m_zvtx = M_zvtx;
0364     }
0365
0366     double getzout()
0367     {
0368         double y = m_truthpy / fabs(m_truthpy) * 10;
0369         double z_outer = (y - truth_rz->GetParameter(0)) / truth_rz->GetParameter(1);
0370
0371         return z_outer;
0372     }
0373
0374     double getpointr(int mode)
0375     {
0376         double r = 0;
0377         if (mode == 1) // truth
0378         {
0379             r = m_truthpy / fabs(m_truthpy) * m_truthpt;
0380         }
0381         else if (mode == 2) // dca_mean
0382         {
0383             r = dca_mean[1] / fabs(dca_mean[1]) * sqrt(dca_mean[0] * dca_mean[0] + dca_mean[1] * dca_mean[1]);
0384         }
0385         else if (mode == 3) // vertex
0386         {
0387             r = m_yvtx / fabs(m_yvtx) * sqrt(m_xvtx * m_xvtx + m_yvtx * m_yvtx);
0388         }
0389
0390         return r;
0391     }
0392
0393     void calc_line()
0394     {
0395         if (truth_xy != nullptr)
0396         {
0397             delete truth_xy;
0398             truth_xy = nullptr;
0399         }
0400
0401         truth_xy = new TF1("truth_xy", "pol1", -1000, 1000);
0402
0403         truth_xy->SetParameter(1, m_truthpy / m_truthpx);
0404         truth_xy->SetParameter(0, -m_yvtx / m_truthpx * m_xvtx + m_yvtx);
0405
0406         if (truth_rz != nullptr)
0407         {
0408             delete truth_rz;
0409             truth_rz = nullptr;
0410         }
0411
0412         truth_rz = new TF1("truth_rz", "pol1", -1000, 1000);
0413
0414         // from the truth vertex
0415         double slope = getpointr(1) / m_truthpz;
0416         double intercept = -slope * m_zvtx + getpointr(3);
0417
0418         truth_rz->SetParameter(1, slope);
0419         truth_rz->SetParameter(0, intercept);
0420
0421         m_truthz_out = getzout();
0422     }
0423
0424     void calc_center()
0425     {
0426         m_rad = m_truthpt / (0.3 * 1.4) * 100;
0427         double sign = m_truthpid / fabs(m_truthpid);
0428         center[0] = -sign * (m_xvtx + m_rad * (-sin(m_truthphi)));
0429         center[1] = -sign * (m_yvtx + m_rad * cos(m_truthphi));
0430     }
0431 };