1 /*
2 * HelixParameterCalculator.java
3 *
4 * Created on July 7th, 2008, 11:09 AM
5 *
6 *
7 */
8 package org.lcsim.fit.helicaltrack;
9
10 import hep.physics.vec.BasicHep3Vector;
11 import hep.physics.vec.Hep3Vector;
12 import org.lcsim.constants.Constants;
13 import org.lcsim.event.MCParticle;
14 import org.lcsim.event.EventHeader;
15
16 /**
17 * Class used for calculating MC particle track paramters
18 * @author Pelham Keahey
19 *
20 */
21 public class HelixParamCalculator {
22
23 private double R, p, pt, theta, arclength;
24 //Some varibles are usesd in other calculations, thus they are global
25 /*
26 xc, yc --coordinates
27 mcdca -- MC distance of closest approach
28 mcphi0 --azimuthal angle
29 tanL -- Slope SZ plane ds/dz
30 x0,y0 are the position of the particle at the dca
31 */
32 private double mcdca,mcphi0,tanL,z0;
33 private double x0,y0;
34 /**
35 * Constructor that is fed a magnetic field and MCPARTICLE
36 * @param mcpc
37 * @param cBField
38 */
39 public HelixParamCalculator(MCParticle mcp, double BField)
40 {
41
42 //Calculate theta, the of the helix projected into an SZ plane, from the z axis
43 double px = mcp.getPX();
44 double py = mcp.getPY();
45 double pz = mcp.getPZ();
46 pt = Math.sqrt(px*px + py*py);
47 p = Math.sqrt(pt*pt + pz*pz);
48 double cth = pz / p;
49 theta = Math.acos(cth);
50
51 //Calculate Radius of the Helix
52 R = mcp.getCharge() * pt / (Constants.fieldConversion * BField);
53
54 //Slope in the Dz/Ds sense, tanL Calculation
55 tanL = pz / pt;
56
57 // Azimuthal direction at origin
58 double mcphi = Math.atan2(py, px);
59
60 //Distance of closest approach Calculation
61 double xc = mcp.getOriginX() + R * Math.sin(mcphi);
62 double yc = mcp.getOriginY() - R * Math.cos(mcphi);
63
64 double Rc = Math.sqrt(xc*xc + yc*yc);
65
66 if(mcp.getCharge()>0)
67 {
68 mcdca = R - Rc;
69 }
70 else
71 {
72 mcdca = R + Rc;
73 }
74
75
76 //azimuthal calculation of the momentum at the DCA, phi0, Calculation
77 mcphi0 = Math.atan2(xc/(R-mcdca), -yc/(R-mcdca));
78 if(mcphi0<0)
79 {
80 mcphi0 += 2*Math.PI;
81 }
82 //z0 Calculation, z position of the particle at dca
83 x0 = -mcdca*Math.sin(mcphi0);
84 y0 = mcdca*Math.cos(mcphi0);
85 arclength = (((mcp.getOriginX()-x0)*Math.cos(mcphi0))+((mcp.getOriginY()-y0)*Math.sin(mcphi0)));
86 z0 = mcp.getOriginZ() - arclength * tanL;
87
88 }
89
90
91 /*
92 * added 9/3/2012 by mgraham
93 * calculate the helix parameters given 3-momentum and a point (+ charge and field)
94 */
95 public HelixParamCalculator(Hep3Vector momentum, Hep3Vector origin, int charge,double BField)
96 {
97
98 //Calculate theta, the of the helix projected into an SZ plane, from the z axis
99 double px = momentum.x();
100 double py = momentum.y();
101 double pz = momentum.z();
102 pt = Math.sqrt(px*px + py*py);
103 p = Math.sqrt(pt*pt + pz*pz);
104 double cth = pz / p;
105 theta = Math.acos(cth);
106 // System.out.println("pt = "+pt+"; costh = "+cth);
107
108 //Calculate Radius of the Helix
109 R = charge* pt / (Constants.fieldConversion * BField);
110
111
112 //Slope in the Dz/Ds sense, tanL Calculation
113 tanL = pz / pt;
114
115 // Azimuthal direction at origin
116 double mcphi = Math.atan2(py, px);
117
118 //Distance of closest approach Calculation
119 double xc = origin.x() + R * Math.sin(mcphi);
120 double yc = origin.y() - R * Math.cos(mcphi);
121 // System.out.println("xc = "+xc+"; yc = "+yc);
122 double Rc = Math.sqrt(xc*xc + yc*yc);
123
124 if(charge>0)
125 {
126 mcdca = R - Rc;
127 }
128 else
129 {
130 mcdca = R + Rc;
131 }
132
133
134 //azimuthal calculation of the momentum at the DCA, phi0, Calculation
135 mcphi0 = Math.atan2(xc/(R-mcdca), -yc/(R-mcdca));
136 if(mcphi0<0)
137 {
138 mcphi0 += 2*Math.PI;
139 }
140 //z0 Calculation, z position of the particle at dca
141 x0 = -mcdca*Math.sin(mcphi0);
142 y0 = mcdca*Math.cos(mcphi0);
143 arclength = (((origin.x()-x0)*Math.cos(mcphi0))+((origin.y()-y0)*Math.sin(mcphi0)));
144 z0 =origin.z() - arclength * tanL;
145
146 }
147 /**
148 * Calculates the B-Field from event
149 * @param mcpc
150 * @param eventc
151 */
152 public HelixParamCalculator(MCParticle mcpc,EventHeader eventc)
153 {
154 this(mcpc,eventc.getDetector().getFieldMap().getField(new BasicHep3Vector(0.,0.,0.)).z());
155 }
156 /**
157 * Return the radius of the Helix track
158 * @return double R
159 */
160 public double getRadius()
161 {
162 return R;
163 }
164 /**
165 * Return the theta angle for the projection of the helix in the SZ plane
166 * from the z axis
167 * @return double theta
168 */
169 public double getTheta()
170 {
171 return theta;
172 }
173 /**
174 * Return the particle's momentum
175 * @return double mcp momentum
176 */
177 public double getMCMomentum()
178 {
179 return p;
180 }
181 /**
182 * Return the curvature (omega)
183 * @return double omega
184 */
185 public double getMCOmega()
186 {
187 return 1. / R;
188 }
189 /**
190 * Return the transvers momentum of the MC particle, Pt
191 * @return double Pt
192 */
193 public double getMCTransverseMomentum()
194 {
195 return pt;
196 }
197 /**
198 * Return the slope of the helix in the SZ plane, tan(lambda)
199 * @return double tanL
200 */
201 public double getSlopeSZPlane()
202 {
203 return tanL;
204 }
205 /**
206 * Return the distance of closest approach
207 * @return double mcdca
208 */
209 public double getDCA()
210 {
211 return mcdca;
212 }
213 /**
214 * Return the azimuthal angle of the momentum when at the position of closest approach
215 * @return double mcphi0
216 */
217 public double getPhi0()
218 {
219 return mcphi0;
220 }
221 /**
222 * Return the z position at the distance of closest approach
223 * @return double z0 position
224 */
225 public double getZ0()
226 {
227 return z0;
228 }
229 /**
230 * Return the arclength of the helix from the ORIGIN TO THE DCA
231 * @return double arclength
232 */
233 public double getArcLength()
234 {
235 return arclength;
236 }
237 /**
238 * Return the x position of the particle when at the dca
239 * @return double arclength
240 */
241 public double getX0()
242 {
243 return x0;
244 }
245 /**
246 * Return the y position of the particle at the dca
247 * @return double arclength
248 */
249 public double getY0()
250 {
251 return y0;
252 }
253
254 }