event_display/macro/DoTracking.hh

0001 #pragma once
0002
0003 #include "cluster.hh"
0004 #include "clustEvent.hh"
0005 #include "analysis_tracking.hh"
0006
0007 //extern int n_dotracking;
0008
0009 bool dotracking(clustEvent &vCluster, TH2F* h_dphi_nocut )
0010 {
0011   for (unsigned int i = 0; i < vCluster.vclus.size(); i++)
0012     {
0013       cluster c1 = vCluster.vclus[i];
0014       if (2 <= c1.layer) // inner
0015     continue;
0016
0017       for (unsigned int j = i + 1; j < vCluster.vclus.size(); j++)
0018         {
0019       cluster c2 = vCluster.vclus[j];
0020       if (c2.layer <= 1) // outer
0021         continue;
0022
0023       TVector3 beamspot = {vCluster.dca_mean[0], vCluster.dca_mean[1], 0};
0024       // TVector3 beamspot = {0, 0, 0};
0025
0026       TVector3 p1 = c1.pos - beamspot;
0027       TVector3 p2 = c2.pos - beamspot;
0028
0029       double p1_phi = atan2(p1.y(), p1.x());
0030       double p2_phi = atan2(p2.y(), p2.x());
0031       double d_phi = p2_phi - p1_phi;
0032       d_phi -= -2 * (M_PI * (int)(d_phi / (M_PI)));
0033
0034       double p1_r = sqrt(p1.x() * p1.x() + p1.y() * p1.y());
0035       double p2_r = sqrt(p2.x() * p2.x() + p2.y() * p2.y());
0036       double p1_theta = atan2(p1_r, p1.z());
0037       double p2_theta = atan2(p2_r, p2.z());
0038       double d_theta = p2_theta - p1_theta;
0039       if (n_dotracking == 1)
0040         h_dphi_nocut->Fill(p1_phi, d_phi);
0041
0042       if (fabs(d_phi) > 0.04)
0043         continue;
0044
0045       if( fabs(d_theta) > 0.1 ) // d_theta cut
0046         continue;
0047
0048       TVector3 u = p2 - p1;
0049       double unorm = sqrt(u.x() * u.x() + u.y() * u.y());
0050       // unit vector in 2D
0051       double ux = u.x() / unorm;
0052       double uy = u.y() / unorm;
0053       double uz = u.z() / unorm; // same normalization factor(2D) is multiplied
0054
0055       TVector3 p0 = beamspot - p1;
0056
0057       double dca_p0 = p0.x() * uy - p0.y() * ux; // cross product of p0 x u
0058       double len_p0 = p0.x() * ux + p0.y() * uy; // dot product of p0 . u
0059
0060       // beam center in X-Y plane
0061       double vx = len_p0 * ux + p1.x();
0062       double vy = len_p0 * uy + p1.y();
0063       double vz = len_p0 * uz + p1.z();
0064
0065       track *tr = new track();
0066       tr->phi_in = p1_phi;
0067       tr->phi_out = p2_phi;
0068       tr->theta_in = c1.gettheta_clus();
0069       tr->theta_out = c2.gettheta_clus();
0070       tr->d_phi = d_phi;
0071       tr->d_theta = d_theta;
0072       tr->charge = dca_p0 / fabs(dca_p0);
0073
0074       tr->dca[0] = vx;
0075       tr->dca[1] = vy;
0076       tr->dca[2] = vz;
0077
0078       tr->p1 = c1.pos;
0079       tr->p2 = c2.pos;
0080
0081       tr->zline.at(0) = c1.zline;
0082       tr->zline.at(1) = c2.zline;
0083
0084       tr->dca_2d = dca_p0;
0085
0086       tr->AddAssociatedClusterId( i );
0087       tr->AddAssociatedClusterId( j );
0088
0089       vCluster.vtrack.push_back(tr);
0090         } // end of for (unsigned int j = i + 1; j < vCluster.vclus.size(); j++)
0091     } // end of   for (unsigned int i = 0; i < vCluster.vclus.size(); i++)
0092
0093   return true;
0094 }