CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SiStripFineDelayTOF.cc
Go to the documentation of this file.
4 #include <iostream>
5 #include "TVector3.h"
6 #include "TMath.h"
7 
8 //#define TIF_COSMIC_SETUP
9 
10 double SiStripFineDelayTOF::timeOfFlight(bool cosmics, bool field, double* trackParameters, double* hit, double* phit, bool onDisk)
11 {
12 
13  // case 1: cosmics with no field.
14  if(cosmics && !field) {
15  return timeOfFlightCosmic(hit,phit);
16  }
17  // case 2: cosmics with field.
18  else if(cosmics && field) {
19  return timeOfFlightCosmicB(trackParameters,hit,phit,onDisk);
20  }
21  // case 3: beam with no field.
22  else if(!cosmics && !field) {
23  return timeOfFlightBeam(hit,phit);
24  }
25  // case 4: beam with field.
26  else {
27  return timeOfFlightBeamB(trackParameters,hit,phit,onDisk);
28  }
29 }
30 
31 double SiStripFineDelayTOF::timeOfFlightCosmic(double* hit, double* phit)
32 {
33  // constants
34  const double c = 30; // cm/ns
35 #ifndef TIF_COSMIC_SETUP
36  const double Rmu = 385; // cm
37  const double zmu = 560; // cm
38  // estimate the time for crossing the barrel
39  TVector3 r0(hit[0],hit[1],hit[2]);
40  TVector3 pr(phit[0],phit[1],phit[2]);
41  TVector2 r0_xy = r0.XYvector();
42  TVector2 pr_xy = pr.Unit().XYvector();
43  double t_barrel = ((r0_xy*pr_xy)+sqrt( (r0_xy*pr_xy)*(r0_xy*pr_xy) - r0_xy.Mod2() + Rmu*Rmu ))/c*pr.Mag()/pr.Perp();
44  // estimate the time for crossing endcaps
45  double t_endcap = fabs(((phit[2]/fabs(phit[2])*zmu)+hit[2])/c*pr.Mag()/pr.Pz());
46  // take the minimum
47  return t_barrel<t_endcap ? t_barrel : t_endcap;
48 #else
49  const double y_trigger = 100; //cm
50  double p = sqrt(phit[0]*phit[0]+phit[1]*phit[1]+phit[2]*phit[2]);
51 // LogDebug("TOF") << "momentum:" << phit[0] << " " << phit[1] << " " << phit[2];
52 // LogDebug("TOF") << "p/py=" << p/phit[1] << " Y0,Y,dY = " << y_trigger << " " << hit[1] << " " << (y_trigger-hit[1]);
53 // LogDebug("TOF") << "d, t=d/c : " << ((y_trigger-hit[1])*(p/phit[1])) << " " << ((y_trigger-hit[1])*(p/phit[1])/c);
54  return fabs((y_trigger-hit[1])*(p/phit[1])/c);
55 #endif
56 }
57 
58 double SiStripFineDelayTOF::timeOfFlightCosmicB(double* trackParameters, double* hit, double* phit, bool onDisk)
59 {
60  // constants
61  const double Rmu = 385; // cm
62  const double zmu = 560; // cm
63  const double c = 30; // cm/ns
64  // track parameters
65  double& kappa = trackParameters[0];
66  double& theta = trackParameters[1];
67  double& phi_0 = trackParameters[2];
68  double& d_0 = trackParameters[3];
69  double& z_0 = trackParameters[4];
70  double invkappa = 1/kappa;
71  // computes the value of the track parameter that correspond to the hit, relative to phi_0
72  //double phi = kappa*tan(theta)*(hit[2]-z_0);
73  double phi = getPhi(trackParameters,hit,onDisk) - phi_0;
74  // computes the value of the track parameter that correspond to the muon system, relative to phi_0
75  // phi_mu = phi_mu0 - phi_0
76  double phi_mu_b = (kappa>0 ? -1 : 1 ) * acos((Rmu*Rmu-d_0*d_0-2*invkappa*invkappa+2*invkappa*d_0)/(2*invkappa*d_0-2*invkappa*invkappa));
77  double phi_mu_e = kappa*tan(theta)*((phit[2]<0 ? 1 : -1)*zmu-z_0);
78  // estimate the time for crossing the barrel
79  double t_barrel = fabs(invkappa/(c*sin(theta))*(phi-phi_mu_b));
80  // estimate the time for crossing endcaps
81  double t_endcap = fabs(invkappa/(c*sin(theta))*(phi-phi_mu_e));
82  // take the minimum
83  return t_barrel<t_endcap ? t_barrel : t_endcap;
84 }
85 
87 {
88  // constants
89  const double c = 30; // cm/ns
90  TVector3 r0(hit[0],hit[1],hit[2]);
91  return r0.Mag()/c;
92 }
93 
94 double SiStripFineDelayTOF::timeOfFlightBeamB(double* trackParameters, double* hit, double* phit, bool onDisk)
95 {
96  // constants
97  const double c = 30; // cm/ns
98  // track parameters
99  double& theta = trackParameters[1];
100  // returns the time of flight from the origin
101  return onDisk? fabs(hit[2]/(c*cos(theta))) : fabs(sqrt(hit[0]*hit[0]+hit[1]*hit[1])/(c*sin(theta)));
102 }
103 
104 double SiStripFineDelayTOF::x(double* trackParameters, double phi)
105 {
106  return trackParameters[3]*sin(trackParameters[2]) + (1/trackParameters[0])*(sin(phi)-sin(trackParameters[2]));
107 }
108 
109 double SiStripFineDelayTOF::y(double* trackParameters, double phi)
110 {
111  return -trackParameters[3]*cos(trackParameters[2]) - (1/trackParameters[0])*(cos(phi)-cos(trackParameters[2]));
112 }
113 
114 double SiStripFineDelayTOF::z(double* trackParameters, double phi)
115 {
116  return trackParameters[4] + (phi-trackParameters[2])/(trackParameters[0]*tan(trackParameters[1]));
117 }
118 
119 double SiStripFineDelayTOF::getPhi(double* trackParameters, double* hit, bool onDisk)
120 {
121  if(onDisk) //use z coordinate to find phi
122  return trackParameters[2] + trackParameters[0]*tan(trackParameters[1])*(hit[2]-trackParameters[4]);
123  else // use x,y coordinate to find phi
124  {
125  double phi =0;
126  if(trackParameters[2]>-2.35 && trackParameters[2]<-0.78) {
127  // first use y
128  phi = acos(-trackParameters[0]*(hit[1]-(trackParameters[3]-1/trackParameters[0])*(-cos(trackParameters[2]))));
129  // use x to resolve the ambiguity
130  if(fabs(x(trackParameters,phi)-hit[0])>fabs(x(trackParameters,-phi)-hit[0])) phi *= -1.;
131  } else {
132  // first use x
133  phi = asin(trackParameters[0]*(hit[0]-(trackParameters[3]-1/trackParameters[0])*sin(trackParameters[2])));
134  // use y to resolve the ambiguity
135  if((fabs(y(trackParameters,phi)-hit[1])>fabs(y(trackParameters,TMath::Pi()-phi)-hit[1])))
136  phi = phi>0 ? TMath::Pi()-phi: -TMath::Pi()-phi;
137  }
138  return phi;
139  }
140  return 0.;
141 }
142 
143 void SiStripFineDelayTOF::trackParameters(const reco::Track& tk,double* trackParameters)
144 {
146  math::XYZVector momentum = tk.outerMomentum();
147  LogDebug("trackParameters") << "outer position: " << position.x() << " " << position.y() << " " << position.z();
148  LogDebug("trackParameters") << "outer momentum: " << momentum.x() << " " << momentum.y() << " " << momentum.z();
149  math::XYZVector fieldDirection(0,0,1);
150  // computes the center of curvature
151  math::XYZVector radius = momentum.Cross(fieldDirection).Unit();
152  position -= radius / trackParameters[0];
153  LogDebug("trackParameters") << "center of curvature: " << position.x() << " " << position.y() << " " << position.z();
154  // the transverse IP
155  trackParameters[3] = position.rho() - fabs(1/trackParameters[0]);
156  if(trackParameters[0]>0) trackParameters[3] *= -1.;
157  // phi_0
158  double phi_out = trackParameters[2];
159  double phi_0 = position.phi() - TMath::Pi()/2;
160  if(trackParameters[0]<0) phi_0 -= TMath::Pi();
161  if(phi_0<-2*TMath::Pi()) phi_0 += 2*TMath::Pi();
162  if(phi_0>2*TMath::Pi()) phi_0 -= 2*TMath::Pi();
163  if(phi_0 > 0) phi_0 -= 2*TMath::Pi();
164  LogDebug("trackParameters") << "phi_0: " << phi_0;
165  trackParameters[2] = phi_0;
166  // z_0
167  trackParameters[4] = tk.outerPosition().z() - (phi_out-phi_0)/TMath::Tan(trackParameters[1])/trackParameters[0];
168  LogDebug("trackParameters") << "z_0: " << tk.outerPosition().z() << " " << trackParameters[4] << " " << tk.innerPosition().z();
169 }
170 
#define LogDebug(id)
const double Pi
static double timeOfFlightBeamB(double *trackParameters, double *hit, double *phit, bool onDisk)
static double x(double *trackParameters, double phi)
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Geom::Theta< T > theta() const
const math::XYZPoint & outerPosition() const
position of the outermost hit
Definition: Track.h:47
const math::XYZPoint & innerPosition() const
position of the innermost hit
Definition: Track.h:42
static double timeOfFlight(bool cosmics, bool field, double *trackParameters, double *hit, double *phit, bool onDisk)
T sqrt(T t)
Definition: SSEVec.h:48
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
static double timeOfFlightCosmicB(double *trackParameters, double *hit, double *phit, bool onDisk)
static void trackParameters(const reco::Track &tk, double *trackParameters)
static double getPhi(double *trackParameters, double *hit, bool onDisk)
static double timeOfFlightCosmic(double *hit, double *phit)
const math::XYZVector & outerMomentum() const
momentum vector at the outermost hit position
Definition: Track.h:49
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
static double y(double *trackParameters, double phi)
static double timeOfFlightBeam(double *hit, double *phit)
static double z(double *trackParameters, double phi)
static int position[264][3]
Definition: ReadPGInfo.cc:509
static const G4double kappa
Definition: DDAxes.h:10