File indexing completed on 2024-04-06 12:08:35
0001 #include "DQM/SiStripCommissioningSources/plugins/tracking/SiStripFineDelayTOF.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include <DataFormats/TrackReco/interface/Track.h>
0004 #include <iostream>
0005 #include "TVector3.h"
0006 #include "TMath.h"
0007
0008
0009
0010 double SiStripFineDelayTOF::timeOfFlight(
0011 bool cosmics, bool field, double* trackParameters, double* hit, double* phit, bool onDisk) {
0012
0013 if (cosmics && !field) {
0014 return timeOfFlightCosmic(hit, phit);
0015 }
0016
0017 else if (cosmics && field) {
0018 return timeOfFlightCosmicB(trackParameters, hit, phit, onDisk);
0019 }
0020
0021 else if (!cosmics && !field) {
0022 return timeOfFlightBeam(hit, phit);
0023 }
0024
0025 else {
0026 return timeOfFlightBeamB(trackParameters, hit, phit, onDisk);
0027 }
0028 }
0029
0030 double SiStripFineDelayTOF::timeOfFlightCosmic(double* hit, double* phit) {
0031
0032 const double c = 30;
0033 #ifndef TIF_COSMIC_SETUP
0034 const double Rmu = 385;
0035 const double zmu = 560;
0036
0037 TVector3 r0(hit[0], hit[1], hit[2]);
0038 TVector3 pr(phit[0], phit[1], phit[2]);
0039 TVector2 r0_xy = r0.XYvector();
0040 TVector2 pr_xy = pr.Unit().XYvector();
0041 double t_barrel =
0042 ((r0_xy * pr_xy) + sqrt((r0_xy * pr_xy) * (r0_xy * pr_xy) - r0_xy.Mod2() + Rmu * Rmu)) / c * pr.Mag() / pr.Perp();
0043
0044 double t_endcap = fabs(((phit[2] / fabs(phit[2]) * zmu) + hit[2]) / c * pr.Mag() / pr.Pz());
0045
0046 return t_barrel < t_endcap ? t_barrel : t_endcap;
0047 #else
0048 const double y_trigger = 100;
0049 double p = sqrt(phit[0] * phit[0] + phit[1] * phit[1] + phit[2] * phit[2]);
0050
0051
0052
0053 return fabs((y_trigger - hit[1]) * (p / phit[1]) / c);
0054 #endif
0055 }
0056
0057 double SiStripFineDelayTOF::timeOfFlightCosmicB(double* trackParameters, double* hit, double* phit, bool onDisk) {
0058
0059 const double Rmu = 385;
0060 const double zmu = 560;
0061 const double c = 30;
0062
0063 double& kappa = trackParameters[0];
0064 double& theta = trackParameters[1];
0065 double& phi_0 = trackParameters[2];
0066 double& d_0 = trackParameters[3];
0067 double& z_0 = trackParameters[4];
0068 double invkappa = 1 / kappa;
0069
0070
0071 double phi = getPhi(trackParameters, hit, onDisk) - phi_0;
0072
0073
0074 double phi_mu_b = (kappa > 0 ? -1 : 1) * acos((Rmu * Rmu - d_0 * d_0 - 2 * invkappa * invkappa + 2 * invkappa * d_0) /
0075 (2 * invkappa * d_0 - 2 * invkappa * invkappa));
0076 double phi_mu_e = kappa * tan(theta) * ((phit[2] < 0 ? 1 : -1) * zmu - z_0);
0077
0078 double t_barrel = fabs(invkappa / (c * sin(theta)) * (phi - phi_mu_b));
0079
0080 double t_endcap = fabs(invkappa / (c * sin(theta)) * (phi - phi_mu_e));
0081
0082 return t_barrel < t_endcap ? t_barrel : t_endcap;
0083 }
0084
0085 double SiStripFineDelayTOF::timeOfFlightBeam(double* hit, double*) {
0086
0087 const double c = 30;
0088 TVector3 r0(hit[0], hit[1], hit[2]);
0089 return r0.Mag() / c;
0090 }
0091
0092 double SiStripFineDelayTOF::timeOfFlightBeamB(double* trackParameters, double* hit, double* phit, bool onDisk) {
0093
0094 const double c = 30;
0095
0096 double& theta = trackParameters[1];
0097
0098 return onDisk ? fabs(hit[2] / (c * cos(theta))) : fabs(sqrt(hit[0] * hit[0] + hit[1] * hit[1]) / (c * sin(theta)));
0099 }
0100
0101 double SiStripFineDelayTOF::x(double* trackParameters, double phi) {
0102 return trackParameters[3] * sin(trackParameters[2]) + (1 / trackParameters[0]) * (sin(phi) - sin(trackParameters[2]));
0103 }
0104
0105 double SiStripFineDelayTOF::y(double* trackParameters, double phi) {
0106 return -trackParameters[3] * cos(trackParameters[2]) -
0107 (1 / trackParameters[0]) * (cos(phi) - cos(trackParameters[2]));
0108 }
0109
0110 double SiStripFineDelayTOF::z(double* trackParameters, double phi) {
0111 return trackParameters[4] + (phi - trackParameters[2]) / (trackParameters[0] * tan(trackParameters[1]));
0112 }
0113
0114 double SiStripFineDelayTOF::getPhi(double* trackParameters, double* hit, bool onDisk) {
0115 if (onDisk)
0116 return trackParameters[2] + trackParameters[0] * tan(trackParameters[1]) * (hit[2] - trackParameters[4]);
0117 else
0118 {
0119 double phi = 0;
0120 if (trackParameters[2] > -2.35 && trackParameters[2] < -0.78) {
0121
0122 phi = acos(-trackParameters[0] *
0123 (hit[1] - (trackParameters[3] - 1 / trackParameters[0]) * (-cos(trackParameters[2]))));
0124
0125 if (fabs(x(trackParameters, phi) - hit[0]) > fabs(x(trackParameters, -phi) - hit[0]))
0126 phi *= -1.;
0127 } else {
0128
0129 phi =
0130 asin(trackParameters[0] * (hit[0] - (trackParameters[3] - 1 / trackParameters[0]) * sin(trackParameters[2])));
0131
0132 if ((fabs(y(trackParameters, phi) - hit[1]) > fabs(y(trackParameters, TMath::Pi() - phi) - hit[1])))
0133 phi = phi > 0 ? TMath::Pi() - phi : -TMath::Pi() - phi;
0134 }
0135 return phi;
0136 }
0137 return 0.;
0138 }
0139
0140 void SiStripFineDelayTOF::trackParameters(const reco::Track& tk, double* trackParameters) {
0141 math::XYZPoint position = tk.outerPosition();
0142 const math::XYZVector& momentum = tk.outerMomentum();
0143 LogDebug("trackParameters") << "outer position: " << position.x() << " " << position.y() << " " << position.z();
0144 LogDebug("trackParameters") << "outer momentum: " << momentum.x() << " " << momentum.y() << " " << momentum.z();
0145 math::XYZVector fieldDirection(0, 0, 1);
0146
0147 math::XYZVector radius = momentum.Cross(fieldDirection).Unit();
0148 position -= radius / trackParameters[0];
0149 LogDebug("trackParameters") << "center of curvature: " << position.x() << " " << position.y() << " " << position.z();
0150
0151 trackParameters[3] = position.rho() - fabs(1 / trackParameters[0]);
0152 if (trackParameters[0] > 0)
0153 trackParameters[3] *= -1.;
0154
0155 double phi_out = trackParameters[2];
0156 double phi_0 = position.phi() - TMath::Pi() / 2;
0157 if (trackParameters[0] < 0)
0158 phi_0 -= TMath::Pi();
0159 if (phi_0 < -2 * TMath::Pi())
0160 phi_0 += 2 * TMath::Pi();
0161 if (phi_0 > 2 * TMath::Pi())
0162 phi_0 -= 2 * TMath::Pi();
0163 if (phi_0 > 0)
0164 phi_0 -= 2 * TMath::Pi();
0165 LogDebug("trackParameters") << "phi_0: " << phi_0;
0166 trackParameters[2] = phi_0;
0167
0168 trackParameters[4] = tk.outerPosition().z() - (phi_out - phi_0) / TMath::Tan(trackParameters[1]) / trackParameters[0];
0169 LogDebug("trackParameters") << "z_0: " << tk.outerPosition().z() << " " << trackParameters[4] << " "
0170 << tk.innerPosition().z();
0171 }