plugins/tracking/SiStripFineDelayTOF.cc

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 //#define TIF_COSMIC_SETUP
0009
0010 double SiStripFineDelayTOF::timeOfFlight(
0011     bool cosmics, bool field, double* trackParameters, double* hit, double* phit, bool onDisk) {
0012   // case 1: cosmics with no field.
0013   if (cosmics && !field) {
0014     return timeOfFlightCosmic(hit, phit);
0015   }
0016   // case 2: cosmics with field.
0017   else if (cosmics && field) {
0018     return timeOfFlightCosmicB(trackParameters, hit, phit, onDisk);
0019   }
0020   // case 3: beam with no field.
0021   else if (!cosmics && !field) {
0022     return timeOfFlightBeam(hit, phit);
0023   }
0024   // case 4: beam with field.
0025   else {
0026     return timeOfFlightBeamB(trackParameters, hit, phit, onDisk);
0027   }
0028 }
0029
0030 double SiStripFineDelayTOF::timeOfFlightCosmic(double* hit, double* phit) {
0031   // constants
0032   const double c = 30;  // cm/ns
0033 #ifndef TIF_COSMIC_SETUP
0034   const double Rmu = 385;  // cm
0035   const double zmu = 560;  // cm
0036   // estimate the time for crossing the barrel
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   // estimate the time for crossing endcaps
0044   double t_endcap = fabs(((phit[2] / fabs(phit[2]) * zmu) + hit[2]) / c * pr.Mag() / pr.Pz());
0045   // take the minimum
0046   return t_barrel < t_endcap ? t_barrel : t_endcap;
0047 #else
0048   const double y_trigger = 100;  //cm
0049   double p = sqrt(phit[0] * phit[0] + phit[1] * phit[1] + phit[2] * phit[2]);
0050   //  LogDebug("TOF") << "momentum:" << phit[0] << " " << phit[1] << " " << phit[2];
0051   //  LogDebug("TOF") << "p/py=" << p/phit[1] << "  Y0,Y,dY = " << y_trigger << " " << hit[1] << " " << (y_trigger-hit[1]);
0052   //  LogDebug("TOF") << "d, t=d/c : " << ((y_trigger-hit[1])*(p/phit[1])) << " " << ((y_trigger-hit[1])*(p/phit[1])/c);
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   // constants
0059   const double Rmu = 385;  // cm
0060   const double zmu = 560;  // cm
0061   const double c = 30;     // cm/ns
0062   // track parameters
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   // computes the value of the track parameter that correspond to the hit, relative to phi_0
0070   //double phi = kappa*tan(theta)*(hit[2]-z_0);
0071   double phi = getPhi(trackParameters, hit, onDisk) - phi_0;
0072   // computes the value of the track parameter that correspond to the muon system, relative to phi_0
0073   // phi_mu = phi_mu0 - phi_0
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   // estimate the time for crossing the barrel
0078   double t_barrel = fabs(invkappa / (c * sin(theta)) * (phi - phi_mu_b));
0079   // estimate the time for crossing endcaps
0080   double t_endcap = fabs(invkappa / (c * sin(theta)) * (phi - phi_mu_e));
0081   // take the minimum
0082   return t_barrel < t_endcap ? t_barrel : t_endcap;
0083 }
0084
0085 double SiStripFineDelayTOF::timeOfFlightBeam(double* hit, double*) {
0086   // constants
0087   const double c = 30;  // cm/ns
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   // constants
0094   const double c = 30;  // cm/ns
0095   // track parameters
0096   double& theta = trackParameters[1];
0097   // returns the time of flight from the origin
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)  //use z coordinate to find phi
0116     return trackParameters[2] + trackParameters[0] * tan(trackParameters[1]) * (hit[2] - trackParameters[4]);
0117   else  // use x,y coordinate to find phi
0118   {
0119     double phi = 0;
0120     if (trackParameters[2] > -2.35 && trackParameters[2] < -0.78) {
0121       // first use y
0122       phi = acos(-trackParameters[0] *
0123                  (hit[1] - (trackParameters[3] - 1 / trackParameters[0]) * (-cos(trackParameters[2]))));
0124       // use x to resolve the ambiguity
0125       if (fabs(x(trackParameters, phi) - hit[0]) > fabs(x(trackParameters, -phi) - hit[0]))
0126         phi *= -1.;
0127     } else {
0128       // first use x
0129       phi =
0130           asin(trackParameters[0] * (hit[0] - (trackParameters[3] - 1 / trackParameters[0]) * sin(trackParameters[2])));
0131       // use y to resolve the ambiguity
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   // computes the center of curvature
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   // the transverse IP
0151   trackParameters[3] = position.rho() - fabs(1 / trackParameters[0]);
0152   if (trackParameters[0] > 0)
0153     trackParameters[3] *= -1.;
0154   // phi_0
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   // z_0
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 }