StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Helper.cxx
1 /***************************************************************************
2  *
3  * $Id: Helper.cxx,v 1.5 2007/07/12 19:41:00 fisyak Exp $
4  *
5  * Author: Bum Choi, UT Austin, Apr 2002
6  *
7  ***************************************************************************
8  *
9  * Description: This class contains methods used by StHiMicroMaker to
10  * generate highpt uDST's from StEvent for highpt Analysis.
11  *
12  ***************************************************************************
13  *
14  * $Log: Helper.cxx,v $
15  * Revision 1.5 2007/07/12 19:41:00 fisyak
16  * Add includes for ROOT 5.16
17  *
18  * Revision 1.4 2003/09/02 17:58:36 perev
19  * gcc 3.2 updates + WarnOff
20  *
21  * Revision 1.3 2002/05/31 22:11:11 jklay
22  * Fixed a small bug
23  *
24  * Revision 1.2 2002/04/03 00:37:41 jklay
25  * Fixed some bugs, added new version of dcaz
26  *
27  * Revision 1.2 2002/04/03 00:23:27 jklay
28  * Fixed private member access bugs in analysis code
29  *
30  * Revision 1.1 2002/04/02 20:00:41 jklay
31  * Bums highpt uDST Maker
32  *
33  *
34  **************************************************************************/
35 #include "Helper.h"
36 #include <math.h>
37 #include "TMath.h"
38 
39 //#include "StPhysicalHelixD.hh"
40 #include "TRandom.h"
41 
42 // from ben norman.
43 // global X & Y coordinates of X unit vector in local coords, by sector
44 // (X => parallel to padrow in XY plane, right handed wrt Y)
45 static Double_t sectorX[] = {
46  0, 0,
47  0.866025403784439, -0.5,
48  0.5, -0.866025403784439,
49  0, -1,
50  -0.5, -0.866025403784439,
51  -0.866025403784439, -0.5,
52  -1, 0,
53  -0.866025403784439, 0.5,
54  -0.5, 0.866025403784438,
55  0, 1,
56  0.5, 0.866025403784439,
57  0.866025403784438, 0.5,
58  1, 0,
59  0.866025403784439, 0.5,
60  0.5, 0.866025403784438,
61  0, 1,
62  -0.5, 0.866025403784439,
63  -0.866025403784439, 0.5,
64  -1, 0,
65  -0.866025403784439, -0.5,
66  -0.5, -0.866025403784439,
67  0, -1,
68  0.5, -0.866025403784439,
69  0.866025403784438, -0.5,
70  1, 0
71 };
72 
73 // global X & Y coordinates of Y unit vector in local coords, by sector
74 // (Y => normal to padrow in XY plane, radially outward from Z axis)
75 static Double_t sectorY[] = {
76  0,0,
77  0.5, 0.866025403784439,
78  0.866025403784439, 0.5,
79  1, 0,
80  0.866025403784439, -0.5,
81  0.5, -0.866025403784439,
82  0, -1,
83  -0.5, -0.866025403784439,
84  -0.866025403784439, -0.5,
85  -1, 0,
86  -0.866025403784439, 0.5,
87  -0.5, 0.866025403784438,
88  0, 1,
89  -0.5, 0.866025403784439,
90  -0.866025403784439, 0.5,
91  -1, 0,
92  -0.866025403784439, -0.5,
93  -0.5, -0.866025403784439,
94  0, -1,
95  0.5, -0.866025403784439,
96  0.866025403784438, -0.5,
97  1, 0,
98  0.866025403784438, 0.5,
99  0.5, 0.866025403784439,
100  0, 1
101 };
102 
103 //--------------------------------------------------------------
104 /*
105  distance in 2d
106  */
107 
108 double distance(double a1, double b1, double a2, double b2){
109  return ::sqrt( (a1-b1)*(a1-b1) + (a2-b2)*(a2-b2) );
110 }
111 
112 /*
113  distance in 3d
114 */
115 double distance(const StThreeVectorF& point1,
116  const StThreeVectorF& point2)
117 {
118  return (point1-point2).mag();
119 }
120 
121 
122 double dca3d(const StPhysicalHelixD& helix, const StThreeVectorF& point)
123 {
124  return helix.distance(point);
125 
126 }
127 
128 double dca2d(const StPhysicalHelixD& helix, const StThreeVectorF& point,
129  const StThreeVectorF* origin)
130 {
131  if(fabs(helix.curvature())<=static_cast<double>(0)){
132  return lineDca2D(helix,point,*origin);
133  }
134  else{
135  return helixDca2D(helix,point);
136  }
137 }
138 
139 
140 //
141 // i like this one better. dip angle in s-z plane
142 // should be exact assuming a circle in x-y
143 // (redundant StTrack* for debugging)
144 double dcaz(const StPhysicalHelixD& helix, const StThreeVectorF& point,
145  const StTrack* track)
146 {
147  double z0 = helix.origin().z();
148  double phi = atan2(point.y()-helix.ycenter(),
149  point.x()-helix.xcenter());
150  int h = helix.h(); // -sign(q*B) (+ := ccw looking down from +z)
151  double dphi = h*(phi-helix.phase());
152 
153  // half circle assumption
154  dphi = (fabs(dphi) < M_PI ) ? dphi :
155  ((dphi<0) ? 2*M_PI + dphi : 2*M_PI - dphi);
156 
157  double arclength= (1./helix.curvature()) * dphi;
158 
159  double dcaZ = (point.z() - (z0 + arclength*tan(helix.dipAngle())));
160  /*
161  if(gRandom->Rndm(1)<0.1){
162  cout << "--------" << endl;
163  cout << "zpoint=" << point.z() << endl;
164  if(track)
165  cout << "pt=" << track->geometry()->momentum().perp()
166  << ",fit pts=" << track->fitTraits().numberOfFitPoints(kTpcId)
167  << endl;
168 
169  cout << ">>>zdca=" << dcaZ
170  << ", dphi=" << dphi*180./M_PI
171  << ", z0=" << z0 << ",arclength=" << arclength << endl
172  << ", dip=" << helix.dipAngle()*180./M_PI
173  << ", arc/tan(dip)="<< arclength/tan(helix.dipAngle()) << endl
174  << ">>>dca.z= " << point.z()-helix.at(helix.pathLength(point)).z()
175  << endl;
176  cout << "--------" << endl;
177  }
178  */
179  return dcaZ;
180 
181 }
182 /*
183 double dcaz(const StPhysicalHelixD& helix, const StThreeVectorF& point)
184 {
185  pairD path = helix.pathLength(point.perp());
186 
187  const StThreeVectorD& pos1 = helix.at(path.first);
188  const StThreeVectorD& pos2 = helix.at(path.second);
189  const StThreeVectorD dis1 = point - pos1;
190  const StThreeVectorD dis2 = point - pos2;
191 
192  double dcaZ = (dis1.mag() < dis2.mag()) ? dis1.z() : dis2.z();
193  if(isnan(dcaZ)) return 999;
194  return dcaZ;
195 }
196 */
197 
198 
199 double crossingAngle(const StPhysicalHelixD& helix,
200  const StTpcHit* hit, float bField)
201 {
202  if(fabs(helix.curvature())<=static_cast<double>(0)){
203  return lineCrossingAngle(helix,hit);
204  }
205  else{
206  return helixCrossingAngle(helix,hit,bField);
207  }
208 }
209 
210 double padrowDca(const StPhysicalHelixD& helix,
211  const StTpcHit* hit)
212 {
213  if(fabs(helix.curvature())<=static_cast<double>(0)){
214  return linePadrowDca(helix,hit);
215  }
216  else{
217  return helixPadrowDca(helix,hit);
218  }
219 }
220 
221 //-----------------------------------------------------
222 
223 
224 
225 
226 double
227 propagateToPadrow(const StPhysicalHelixD& helix,
228  const StTpcHit* hit)
229 {
230  // calculate intersection of helix (circle) with hit's padrow.
231  // There's a little trig to do, and 4 cases, but the net result is
232  // that if we define
233  //
234  // R = radius of helix's circular projection
235  // d = distance from center of circle to hit
236  // x = signed displacement from hit to track in sector's local x direction
237  // dHat = unit vector pointing from center to hit
238  // xHat = unit vector pointing along the padrow (x direction in sector's
239  // local coordinates)
240  // cos(theta) = dHat . xHat
241  //
242  // Then
243  //
244  // x = -d*cos(theta) +- ::sqrt(R^2 - d^2*(1-cos^2(theta)))
245  //
246  // '-' is used if cos(theta)<0, and '+' if cos(theta)>0. The overall sign
247  // is determined by whether or not the hit is in the circle, as above.
248 
249  double R = 1./helix.curvature();
250  double dX = hit->position().x() - helix.xcenter();
251  double dY = hit->position().y() - helix.ycenter();
252 
253  // before continuing, make sure that there *is* a solution, i.e., that the
254  // circle does intersect the padrow. We determine this by projecting the
255  // vector between the hit and the circle's center onto the padrow's radial
256  // normal. If this is >= R, we have no intersection and return the helix
257  // path length to the point nearest to the padrow.
258 
259  double d = TMath::Sqrt(dX*dX + dY*dY);
260  double dPerp = dX*sectorY[2*hit->sector()] +
261  dY*sectorY[2*hit->sector() + 1];
262 
263  double x;
264  if(dPerp >= R){
265  // this shouldn't happen for hi-pt, but...
266  //
267  // take cross product (Sin) to find displacement along padrow to point
268  // nearest circle
269  x = - dX*sectorY[2*hit->sector() + 1] + // padrow normal x d
270  dY*sectorY[2*hit->sector()];
271  }else{
272  // find analytic solution of intersection
273  double cosTheta = (dX*sectorX[2*hit->sector()] +
274  dY*sectorX[2*hit->sector() + 1])/d;
275 
276  x = -d*cosTheta + (cosTheta<0 ? -1. : 1.) *
277  TMath::Sqrt(R*R - d*d*(1 - cosTheta*cosTheta));
278  }
279 
280  // finally, propegate along local X direction x units from hit to track
281  // (reuse dX,dY)
282  dX = hit->position().x() + x*sectorX[2*hit->sector()];
283  dY = hit->position().y() + x*sectorX[2*hit->sector() + 1];
284 
285  /*
286  if(mDebug){
287  cout << " Position on helix & padrow in sector " << hit->sector()
288  << " closest to point '" << hit->position() << "' is (" << dX << ", "
289  << dY << "). ";
290  }
291  */
292  double s = helix.pathLength(dX, dY);
293 
294  /* non-analytic solution in helix, not used
295  // normal to plane
296  StThreeVectorD n(sectorY[2*hit->sector()], sectorY[2*hit->sector() + 1], 0);
297  // point on plane
298  double s = helix.pathLength(hit->position(), n);
299  */
300  return s;
301 
302 }
303 
304 double
305 helixPadrowDca(const StPhysicalHelixD& helix,
306  const StTpcHit* hit)
307 {
308  // let propegateToPadrow do the work of finding the intersection of
309  // the helix with the padrow
310  double s = propagateToPadrow(helix, hit);
311  double dX = hit->position().x() - helix.x(s);
312  double dY = hit->position().y() - helix.y(s);
313  double dca = TMath::Sqrt(dX*dX + dY*dY);
314 
315  // get our sign from the simple xyDca calculation
316  if (helixDca2D(helix, hit->position()) < 0) dca *= -1.;
317 
318  /*
319  if(mDebug){
320  cout << " Dca of helix '" << helix << "' to hit '"
321  << hit->position() << "' along padrow is " << dca << endl;
322  }
323  */
324  return dca;
325 }
326 
327 double
328 helixCrossingAngle(const StPhysicalHelixD& helix,
329  const StTpcHit* hit,
330  float bField){
331 
332  // first, we need the track momentum direction at the hit's padrow.
333  double s = propagateToPadrow(helix, hit);
334 
335  // to suppress warning messages
336  StPhysicalHelixD helixTemp(helix);
337 
338  StThreeVectorD p(helixTemp.momentumAt(s, bField));
339 
340  // calculate the cosine of the angle between the radially outward normal
341  // to the padrow (i.e. the local Y unit vector) and the track tangent.
342  double cosTheta = (p.x()*sectorY[2*hit->sector()] +
343  p.y()*sectorY[2*hit->sector() + 1])/p.perp();
344 
345  /*
346  if(mDebug){
347  cout << " Cos of rossing angle of helixTemp '" << helixTemp << "' in sector "
348  << hit->sector() << " is " << cosTheta << ", ";
349  }
350  */
351  // if the cosine is negative, take the inverse of the momentum vector
352  // to bring it between 0 and pi/2.
353  if (cosTheta < 0){
354  cosTheta *= -1.;
355  p.setX(-p.x());
356  p.setY(-p.y());
357  }
358 
359  //if(mDebug) cout << " changed to " << cosTheta << ", ";
360 
361  // use the cross product to determine whether this angle should be
362  // negative or positive.
363  // if cross product is along positive z, then p was clockwise from sectorY
364  double theta = TMath::ACos(cosTheta);
365  if (p.x()*sectorY[2*hit->sector() + 1] -
366  p.y()*sectorY[2*hit->sector()] > 0) theta *= -1.;
367 
368  // if(mDebug) cout << " finally theta=" << theta << endl;
369 
370  return theta;
371 } // crossingAngle
372 
373 double helixDca2D(const StPhysicalHelixD &helix, const StThreeVectorF& point)
374 {
375  return (distance(helix.xcenter(),point.x(),helix.ycenter(),point.y())
376  -(1./helix.curvature()));
377 }
378 
379 /*
380  return StThreeVectorF of point on the line closet to the input point
381  in xy plane. z coordinate is ignored.
382 
383 */
384 
386 lineAt2d(const StPhysicalHelixD& helix,
387  const StThreeVectorF& point)
388 {
389  // note: phase for straight lines is defined as phase=psi-pi/2.
390  // where psi is the 'direction' the line is pointing.
391  // the relationship b/psi and the slope is
392  // slope = tan(psi).
393 
394  const StThreeVectorD& origin = helix.origin();
395  double cosPhase = cos(helix.phase());
396  double sinPhase = sin(helix.phase());
397 
398  double px = point.x();
399  double py = point.y();
400  double dx = px-origin.x();
401  double dy = py-origin.y();
402 
403  double s = cosPhase*dy-sinPhase*dx;
404 
405  // point on line closest to point in xy plane.
406  double xDca = origin.x() - s*sinPhase;
407  double yDca = origin.y() + s*cosPhase;
408 
409  StThreeVectorF vec(xDca,yDca,0); // z coordinate ignored
410  return vec;
411 }
412 
413 
414 
415 
416 double
417 linePadrowDca(const StPhysicalHelixD& helix,
418  const StTpcHit* hit)
419 {
420  // note: phase for straight lines is defined as phase=psi-pi/2.
421  // where psi is the 'direction' the line is pointing.
422  // the relationship b/psi and the slope is
423  // slope = tan(psi).
424 
425  // point on line closest to hit in xy plane.
426 
427  StThreeVectorF at = lineAt2d(helix,hit->position());
428 
429  // vector perpendicular to the line in xy
430  StThreeVectorF perp = (hit->position() - at);
431 
432  // padplane dca
433  double cosTheta =
434  (perp.x()*sectorX[2*hit->sector()] +
435  perp.y()*sectorX[2*hit->sector()+1])/perp.perp();
436 
437  double padrowDca = perp.perp()/fabs(cosTheta);
438 
439  // get sign from lineDca2D
440  StThreeVectorF origin(0,0,0);
441  double sign = lineDca2D(helix,hit->position(),origin);
442 
443  return (sign>=0) ? padrowDca : -padrowDca;
444 }
445 
446 double
447 lineCrossingAngle(const StPhysicalHelixD& helix,
448  const StTpcHit* hit)
449 {
450  // psi
451  double psi = helix.phase() + TMath::Pi()/2.;
452  double psiX = cos(psi);
453  double psiY = sin(psi);
454 
455  double cosTheta =
456  psiX*sectorY[2*hit->sector()] + psiY*sectorY[2*hit->sector()+1];
457 
458  // should never be negative
459  if(cosTheta<0){
460  cout << "%%% line crossing angle negative? "
461  << cosTheta << endl;
462  return 99;
463  }
464 
465  // sign from the cross product (psi X sector Normal).
466  double sign = psiX*sectorY[2*hit->sector()+1]-psiY*sectorY[2*hit->sector()];
467 
468  return (sign>0) ? acos(cosTheta) : -acos(cosTheta);
469 }
470 
471 double
472 lineCrossingAngle(const StPhysicalHelixD& helix,
473  const int sector)
474 {
475  // psi
476  double psi = helix.phase() + TMath::Pi()/2.;
477  double psiX = cos(psi);
478  double psiY = sin(psi);
479 
480  double cosTheta =
481  psiX*sectorY[2*sector] + psiY*sectorY[2*sector+1];
482 
483  // should never be negative
484  if(cosTheta<0){
485  cout << "%%% line crossing angle negative? "
486  << cosTheta << endl;
487  return 99;
488  }
489 
490  // sign from the cross product (psi X sector Normal).
491  double sign = psiX*sectorY[2*sector+1]-psiY*sectorY[2*sector];
492 
493  return (sign>0) ? acos(cosTheta) : -acos(cosTheta);
494 }
495 
496 double
497 lineDca2D(const StPhysicalHelixD& helix,
498  const StThreeVectorF& hit,
499  const StThreeVectorF& origin)
500 {
501  // point on line closest to hit (in xy plane)
502  StThreeVectorF atHit = lineAt2d(helix,hit);
503 
504  // point on line closest to 'origin' in xy
505  StThreeVectorF atOrigin = lineAt2d(helix,origin);
506 
507  // vector from atOrigin to atHit
508  StThreeVectorF origin2atHit = (atHit-atOrigin);
509 
510  // vector fom atOrigin to hit
511  StThreeVectorF origin2hit = (hit-atOrigin);
512 
513  // vector perpendicular to the line to hit.
514  StThreeVectorF dca = (hit-atHit);
515 
516  // sign determined from origin2atHit X orgin2hit
517  double cross = origin2atHit.x()*origin2hit.y()
518  - origin2atHit.y()*origin2hit.x();
519 
520  return (cross>=0) ? dca.perp() : - dca.perp();
521 
522 }
int h() const
y-center of circle in xy-plane
Definition: StHelix.hh:174
double x(double s) const
coordinates of helix at point s
Definition: StHelix.hh:182
double distance(const StThreeVector< double > &p, bool scanPeriods=true) const
minimal distance between point and helix
Definition: StHelix.cc:240
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Definition: StHelix.cc:351
double ycenter() const
x-center of circle in xy-plane
Definition: StHelix.cc:215
const StThreeVector< double > & origin() const
-sign(q*B);
Definition: StHelix.hh:224
double xcenter() const
aziumth in xy-plane measured from ring center
Definition: StHelix.cc:207
double phase() const
1/R in xy-plane
Definition: StHelix.hh:180