Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DFDCSegment_factory.cc
Go to the documentation of this file.
1 //************************************************************************
2 // DFDCSegment_factory.cc - factory producing track segments from pseudopoints
3 //************************************************************************
4 
5 #include "DFDCSegment_factory.h"
6 #include "DANA/DApplication.h"
7 //#include "HDGEOMETRY/DLorentzMapCalibDB.h"
8 #include <math.h>
9 #include <cmath>
10 using namespace std;
11 
12 // some of these aren't being used anymore, so I've commented them out - sdobbs, 6/3/14
13 //#define HALF_CELL 0.5
14 //#define MAX_DEFLECTION 0.15
15 #define EPS 1e-8
16 //#define KILL_RADIUS 5.0
17 #define MATCH_RADIUS 2.0
18 #define ADJACENT_MATCH_DISTANCE 0.3
19 //#define SIGN_CHANGE_CHISQ_CUT 10.0
20 //#define USED_IN_SEGMENT 0x8
21 //#define CORRECTED 0x10
22 //#define MAX_ITER 10
23 //#define MIN_TANL 2.0
24 #define ONE_THIRD 0.33333333333333333
25 #define SQRT3 1.73205080756887719
26 
27 
28 
29 // Variance for position along wire using PHENIX angle dependence, transverse
30 // diffusion, and an intrinsic resolution of 127 microns.
31 // Deprecated!
32 //inline double fdc_y_variance(double alpha,double x){
33 // return 0.00060+0.0064*tan(alpha)*tan(alpha)+0.0004*fabs(x);
34 //}
35 
36 bool DFDCSegment_package_cmp(const DFDCPseudo* a, const DFDCPseudo* b) {
37  if(a->wire->layer == b->wire->layer) {
38  if (fabs(a->time-b->time)<EPS){
39  double tsum_1=a->t_u+a->t_v;
40  double tsum_2=b->t_u+b->t_v;
41  return (tsum_1<tsum_2);
42  }
43  return(a->time<b->time);
44  }
45  return a->wire->layer>b->wire->layer;
46 }
47 
48 
50  _log = new JStreamLog(std::cout, "FDCSEGMENT >>");
51  *_log << "File initialized." << endMsg;
52 }
53 
54 
55 ///
56 /// DFDCSegment_factory::~DFDCSegment_factory():
57 /// default destructor -- closes log file
58 ///
60  delete _log;
61 }
62 ///
63 /// DFDCSegment_factory::brun():
64 /// Initialization: read in deflection map, get magnetic field map
65 ///
66 jerror_t DFDCSegment_factory::brun(JEventLoop* eventLoop, int32_t runnumber) {
67  DApplication* dapp=dynamic_cast<DApplication*>(eventLoop->GetJApplication());
68  const DMagneticFieldMap *bfield = dapp->GetBfield(runnumber);
69  double Bz=bfield->GetBz(0.,0.,65.);
70  RotationSenseToCharge=(Bz>0.)?-1.:1.;
71 
72  // get the geometry
73  const DGeometry *geom = dapp->GetDGeometry(runnumber);
74 
75  geom->GetTargetZ(TARGET_Z);
76  /*
77  bfield = dapp->GetBfield();
78  lorentz_def=dapp->GetLorentzDeflections();
79 
80  *_log << "Table of Lorentz deflections initialized." << endMsg;
81  */
82  DEBUG_LEVEL=0;
83  gPARMS->SetDefaultParameter("FDC:DEBUG_LEVEL", DEBUG_LEVEL);
84 
85  BEAM_VARIANCE=1.0;
86  gPARMS->SetDefaultParameter("FDC:BEAM_VARIANCE",BEAM_VARIANCE);
87 
88 
89  return NOERROR;
90 }
91 
92 ///
93 /// DFDCSegment_factory::evnt():
94 /// Routine where pseudopoints are combined into track segments
95 ///
96 jerror_t DFDCSegment_factory::evnt(JEventLoop* eventLoop, uint64_t eventNo) {
97  myeventno=eventNo;
98 
99  vector<const DFDCPseudo*>pseudopoints;
100  eventLoop->Get(pseudopoints);
101 
102  // Skip segment finding if there aren't enough points to form a sensible
103  // segment
104  if (pseudopoints.size()>=3){
105  // Group pseudopoints by package
106  vector<const DFDCPseudo*>package[4];
107  for (vector<const DFDCPseudo*>::const_iterator i=pseudopoints.begin();
108  i!=pseudopoints.end();i++){
109  package[((*i)->wire->layer-1)/6].push_back(*i);
110  }
111 
112  // Find the segments in each package
113  for (int j=0;j<4;j++){
114  std::stable_sort(package[j].begin(),package[j].end(),DFDCSegment_package_cmp);
115  // We need at least 3 points to define a circle, so skip if we don't
116  // have enough points.
117  if (package[j].size()>2) FindSegments(package[j]);
118  }
119  } // pseudopoints>2
120 
121  return NOERROR;
122 }
123 
124 // Riemann Line fit: linear regression of s on z to determine the tangent of
125 // the dip angle and the z position of the closest approach to the beam line.
126 // Also returns predicted positions along the helical path.
127 //
128 jerror_t DFDCSegment_factory::RiemannLineFit(vector<const DFDCPseudo *>&points,
129  DMatrix &CR,vector<xyz_t>&XYZ){
130  unsigned int n=points.size()+1;
131  vector<int>bad(n); // Keep track of "bad" intersection points
132  bool got_bad_intersection=false;
133  // Fill matrix of intersection points
134  for (unsigned int m=0;m<n-1;m++){
135  double r2=points[m]->xy.Mod2();
136  double denom= N[0]*N[0]+N[1]*N[1];
137  double numer=dist_to_origin+r2*N[2];
138  double ratio=numer/denom;
139 
140  DVector2 xy_int0(-N[0]*ratio,-N[1]*ratio);
141  double temp=denom*r2-numer*numer;
142  if (temp<0){
143  got_bad_intersection=true;
144  bad[m]=1;
145  XYZ[m].xy=points[m]->xy;
146  continue;
147  }
148  temp=sqrt(temp)/denom;
149 
150  // Choose sign of square root based on proximity to actual measurements
151  DVector2 delta(N[1]*temp,-N[0]*temp);
152  DVector2 xy1=xy_int0+delta;
153  DVector2 xy2=xy_int0-delta;
154  double diff1=(xy1-points[m]->xy).Mod2();
155  double diff2=(xy2-points[m]->xy).Mod2();
156  if (diff1>diff2){
157  XYZ[m].xy=xy2;
158  }
159  else{
160  XYZ[m].xy=xy1;
161  }
162  }
163  // Fake target point
164  XYZ[n-1].xy.Set(0.,0.);
165 
166  // Linear regression to find z0, tanl
167  double sumv=0.,sumx=0.;
168  double sumy=0.,sumxx=0.,sumxy=0.;
169  double sperp=0.,sperp_old=0.,ratio,Delta;
170  double z=0,zlast=0;
171  double two_rc=2.*rc;
172  DVector2 oldxy=XYZ[0].xy;
173  unsigned int num_z=0;
174  for (unsigned int k=0;k<n;k++){
175  zlast=z;
176  z=XYZ[k].z;
177  if (fabs(z-zlast)<0.01) continue;
178  num_z++;
179 
180  DVector2 diffxy=XYZ[k].xy-oldxy;
181  double var=CR(k,k);
182  if (bad[k]) var=XYZ[k].covr;
183 
184  sperp_old=sperp;
185  ratio=diffxy.Mod()/(two_rc);
186  // Make sure the argument for the arcsin does not go out of range...
187  sperp=sperp_old+(ratio>1?two_rc*(M_PI_2):two_rc*asin(ratio));
188  //z=XYZ(k,2);
189 
190  // Assume errors in s dominated by errors in R
191  double inv_var=1./var;
192  sumv+=inv_var;
193  sumy+=sperp*inv_var;
194  sumx+=z*inv_var;
195  sumxx+=z*z*inv_var;
196  sumxy+=sperp*z*inv_var;
197 
198  // Save the current x and y coordinates
199  //oldx=XYZ(k,0);
200  //oldy=XYZ(k,1);
201  oldxy=XYZ[k].xy;
202  }
203 
204  Delta=sumv*sumxx-sumx*sumx;
205  double denom=sumv*sumxy-sumy*sumx;
206  if (fabs(Delta)>EPS && fabs(denom)>EPS){
207  // Track parameters z0 and tan(lambda)
208  tanl=-Delta/denom;
209  //z0=(sumxx*sumy-sumx*sumxy)/Delta*tanl;
210  // Error in tanl
211  var_tanl=sumv/Delta*(tanl*tanl*tanl*tanl);
212 
213  // Vertex position
214  sperp-=sperp_old;
215  zvertex=zlast-tanl*sperp;
216  }
217  else{
218  // assume that the particle came from the target
219  zvertex=TARGET_Z;
220  if (sperp<EPS){
221  // This is a rare failure mode where it seems that the arc length data
222  // cannot be reliably used. Provide some default value of tanl to
223  // avoid division by zero errors
224  tanl=57.3; // theta=1 degree
225  }
226  else tanl=(zlast-zvertex)/sperp;
227  var_tanl=100.; // guess for now
228  }
229  if (num_z<3){
230  //_DBG_ << "Too few z points..." << endl;
231  return VALUE_OUT_OF_RANGE;
232  }
233  if (got_bad_intersection) return VALUE_OUT_OF_RANGE;
234  return NOERROR;
235 }
236 
237 // Use predicted positions (propagating from plane 1 using a helical model) to
238 // update the R and RPhi covariance matrices.
239 //
241  double r1sq,vector<xyz_t>&XYZ,DMatrix &CRPhi,DMatrix &CR){
242  double delta_x=XYZ[ref_plane].xy.X()-xc;
243  double delta_y=XYZ[ref_plane].xy.Y()-yc;
244  double r1=sqrt(r1sq);
245  double denom=delta_x*delta_x+delta_y*delta_y;
246 
247  // Auxiliary matrices for correcting for non-normal track incidence to FDC
248  // The correction is
249  // CRPhi'= C*CRPhi*C+S*CR*S, where S(i,i)=R_i*kappa/2
250  // and C(i,i)=sqrt(1-S(i,i)^2)
251  DMatrix C(n,n);
252  DMatrix S(n,n);
253 
254  // Predicted positions
255  Phi1=atan2(delta_y,delta_x);
256  double z1=XYZ[ref_plane].z;
257  double y1=XYZ[ref_plane].xy.X();
258  double x1=XYZ[ref_plane].xy.Y();
259  double var_R1=CR(ref_plane,ref_plane);
260  for (unsigned int k=0;k<n;k++){
261  double sperp=rotation_sense*(XYZ[k].z-z1)/tanl;
262  double phi_s=Phi1+sperp/rc;
263  double sinp=sin(phi_s);
264  double cosp=cos(phi_s);
265  XYZ[k].xy.Set(xc+rc*cosp,yc+rc*sinp);
266 
267  // Error analysis. We ignore errors in N because there doesn't seem to
268  // be any obvious established way to estimate errors in eigenvalues for
269  // small eigenvectors. We assume that most of the error comes from
270  // the measurement for the reference plane radius and the dip angle anyway.
271  double Phi=XYZ[k].xy.Phi();
272  double sinPhi=sin(Phi);
273  double cosPhi=cos(Phi);
274  double dRPhi_dx=Phi*cosPhi-sinPhi;
275  double dRPhi_dy=Phi*sinPhi+cosPhi;
276 
277  double rc_sinphi_over_denom=rc*sinp/denom;
278  //double dx_dx1=rc*sinp*delta_y/denom;
279  //double dx_dy1=-rc*sinp*delta_x/denom;
280  double dx_dx1=rc_sinphi_over_denom*delta_y;
281  double dx_dy1=-rc_sinphi_over_denom*delta_x;
282  double dx_dtanl=sinp*sperp/tanl;
283 
284  double rc_cosphi_over_denom=rc*cosp/denom;
285  //double dy_dx1=-rc*cosp*delta_y/denom;
286  //double dy_dy1=rc*cosp*delta_x/denom;
287  double dy_dx1=-rc_cosphi_over_denom*delta_y;
288  double dy_dy1=rc_cosphi_over_denom*delta_x;
289  double dy_dtanl=-cosp*sperp/tanl;
290 
291  double dRPhi_dx1=dRPhi_dx*dx_dx1+dRPhi_dy*dy_dx1;
292  double dRPhi_dy1=dRPhi_dx*dx_dy1+dRPhi_dy*dy_dy1;
293  double dRPhi_dtanl=dRPhi_dx*dx_dtanl+dRPhi_dy*dy_dtanl;
294 
295  double dR_dx1=cosPhi*dx_dx1+sinPhi*dy_dx1;
296  double dR_dy1=cosPhi*dx_dy1+sinPhi*dy_dy1;
297  double dR_dtanl=cosPhi*dx_dtanl+sinPhi*dy_dtanl;
298 
299  double cdist=dist_to_origin+r1sq*N[2];
300  double n2=N[0]*N[0]+N[1]*N[1];
301 
302  double ydenom=y1*n2+N[1]*cdist;
303  double dy1_dr1=-r1*(2.*N[1]*N[2]*y1+2.*N[2]*cdist-N[0]*N[0])/ydenom;
304  double var_y1=dy1_dr1*dy1_dr1*var_R1;
305 
306  double xdenom=x1*n2+N[0]*cdist;
307  double dx1_dr1=-r1*(2.*N[0]*N[2]*x1+2.*N[2]*cdist-N[1]*N[1])/xdenom;
308  double var_x1=dx1_dr1*dx1_dr1*var_R1;
309 
310  CRPhi(k,k)=dRPhi_dx1*dRPhi_dx1*var_x1+dRPhi_dy1*dRPhi_dy1*var_y1
311  +dRPhi_dtanl*dRPhi_dtanl*var_tanl;
312  CR(k,k)=dR_dx1*dR_dx1*var_x1+dR_dy1*dR_dy1*var_y1
313  +dR_dtanl*dR_dtanl*var_tanl;
314 
315  // double stemp=sqrt(XYZ(k,0)*XYZ(k,0)+XYZ(k,1)*XYZ(k,1))/(4.*rc);
316  double stemp=XYZ[k].xy.Mod()/(4.*rc);
317  double ctemp=1.-stemp*stemp;
318  if (ctemp>0){
319  S(k,k)=stemp;
320  C(k,k)=sqrt(ctemp);
321  }
322  else{
323  S(k,k)=0.;
324  C(k,k)=1.;
325  }
326  }
327 
328  // Correction for non-normal incidence of track on FDC
329  CRPhi=C*CRPhi*C+S*CR*S;
330 
331  return NOERROR;
332 }
333 
334 // Riemann Circle fit: points on a circle in x,y project onto a plane cutting
335 // the circular paraboloid surface described by (x,y,x^2+y^2). Therefore the
336 // task of fitting points in (x,y) to a circle is transormed to the taks of
337 // fitting planes in (x,y, w=x^2+y^2) space
338 //
339 jerror_t DFDCSegment_factory::RiemannCircleFit(vector<const DFDCPseudo *>&points,
340  DMatrix &CRPhi){
341  unsigned int n=points.size()+1;
342  DMatrix X(n,3);
343  DMatrix Xavg(1,3);
344  DMatrix A(3,3);
345  // vector of ones
346  DMatrix OnesT(1,n);
347  double W_sum=0.;
348  DMatrix W(n,n);
349 
350  // The goal is to find the eigenvector corresponding to the smallest
351  // eigenvalue of the equation
352  // lambda=n^T (X^T W X - W_sum Xavg^T Xavg)n
353  // where n is the normal vector to the plane slicing the cylindrical
354  // paraboloid described by the parameterization (x,y,w=x^2+y^2),
355  // and W is the weight matrix, assumed for now to be diagonal.
356  // In the absence of multiple scattering, W_sum is the sum of all the
357  // diagonal elements in W.
358  // At this stage we ignore the multiple scattering.
359  unsigned int last_index=n-1;
360  for (unsigned int i=0;i<last_index;i++){
361  X(i,0)=points[i]->xy.X();
362  X(i,1)=points[i]->xy.Y();
363  X(i,2)=points[i]->xy.Mod2();
364  OnesT(0,i)=1.;
365  W(i,i)=1./CRPhi(i,i);
366  W_sum+=W(i,i);
367  }
368  OnesT(0,last_index)=1.;
369  W(last_index,last_index)=1./CRPhi(last_index,last_index);
370  W_sum+=W(last_index,last_index);
371  var_avg=1./W_sum;
372  Xavg=var_avg*(OnesT*(W*X));
373  // Store in private array for use in other routines
374  xavg[0]=Xavg(0,0);
375  xavg[1]=Xavg(0,1);
376  xavg[2]=Xavg(0,2);
377 
378  A=DMatrix(DMatrix::kTransposed,X)*(W*X)
379  -W_sum*(DMatrix(DMatrix::kTransposed,Xavg)*Xavg);
380  if(!A.IsValid())return UNRECOVERABLE_ERROR;
381 
382  // The characteristic equation is
383  // lambda^3+B2*lambda^2+lambda*B1+B0=0
384  //
385  double B2=-(A(0,0)+A(1,1)+A(2,2));
386  double B1=A(0,0)*A(1,1)-A(1,0)*A(0,1)+A(0,0)*A(2,2)-A(2,0)*A(0,2)
387  +A(1,1)*A(2,2)-A(2,1)*A(1,2);
388  double B0=-A.Determinant();
389  if(B0==0 || !isfinite(B0))return UNRECOVERABLE_ERROR;
390 
391  // The roots of the cubic equation are given by
392  // lambda1= -B2/3 + S+T
393  // lambda2= -B2/3 - (S+T)/2 + i sqrt(3)/2. (S-T)
394  // lambda3= -B2/3 - (S+T)/2 - i sqrt(3)/2. (S-T)
395  // where we define some temporary variables:
396  // S= (R+sqrt(Q^3+R^2))^(1/3)
397  // T= (R-sqrt(Q^3+R^2))^(1/3)
398  // Q=(3*B1-B2^2)/9
399  // R=(9*B2*B1-27*B0-2*B2^3)/54
400  // sum=S+T;
401  // diff=i*(S-T)
402  // We divide Q and R by a safety factor to prevent multiplying together
403  // enormous numbers that cause unreliable results.
404 
405  double Q=(3.*B1-B2*B2)/9.e4;
406  double R=(9.*B2*B1-27.*B0-2.*B2*B2*B2)/54.e6;
407  double Q1=Q*Q*Q+R*R;
408  if (Q1<0) Q1=sqrt(-Q1);
409  else{
410  return VALUE_OUT_OF_RANGE;
411  }
412 
413  // DeMoivre's theorem for fractional powers of complex numbers:
414  // (r*(cos(theta)+i sin(theta)))^(p/q)
415  // = r^(p/q)*(cos(p*theta/q)+i sin(p*theta/q))
416  //
417  //double temp=100.*pow(R*R+Q1*Q1,0.16666666666666666667);
418  double temp=100.*sqrt(cbrt(R*R+Q1*Q1));
419  double theta1=ONE_THIRD*atan2(Q1,R);
420  double sum_over_2=temp*cos(theta1);
421  double diff_over_2=-temp*sin(theta1);
422  // Third root
423  double lambda_min=-ONE_THIRD*B2-sum_over_2+SQRT3*diff_over_2;
424 
425  // Normal vector to plane
426  N[0]=1.;
427  double A11_minus_lambda_min=A(1,1)-lambda_min;
428  N[1]=(A(1,0)*A(0,2)-(A(0,0)-lambda_min)*A(1,2))
429  /(A(0,1)*A(2,1)-(A11_minus_lambda_min)*A(0,2));
430  N[2]=(A(2,0)*A11_minus_lambda_min-A(1,0)*A(2,1))
431  /(A(1,2)*A(2,1)-(A(2,2)-lambda_min)*A11_minus_lambda_min);
432 
433  // Normalize: n1^2+n2^2+n3^2=1
434  double denom=sqrt(N[0]*N[0]+N[1]*N[1]+N[2]*N[2]);
435  for (int i=0;i<3;i++){
436  N[i]/=denom;
437  }
438 
439  // Distance to origin
440  dist_to_origin=-(N[0]*Xavg(0,0)+N[1]*Xavg(0,1)+N[2]*Xavg(0,2));
441 
442  // Center and radius of the circle
443  double two_N2=2.*N[2];
444  xc=-N[0]/two_N2;
445  yc=-N[1]/two_N2;
446  rc=sqrt(1.-N[2]*N[2]-4.*dist_to_origin*N[2])/fabs(two_N2);
447 
448  return NOERROR;
449 }
450 
451 // Riemann Helical Fit based on transforming points on projection to x-y plane
452 // to a circular paraboloid surface combined with a linear fit of the arc
453 // length versus z. Uses RiemannCircleFit and RiemannLineFit.
454 //
455 jerror_t
456 DFDCSegment_factory::RiemannHelicalFit(vector<const DFDCPseudo*>&points){
457  double Phi;
458  unsigned int num_measured=points.size();
459  unsigned int last_index=num_measured;
460  unsigned int num_points=num_measured+1;
461  Ndof=num_points-3;
462  // list of points on track and corresponding covariance matrices
463  vector<xyz_t>XYZ(num_points);
464  DMatrix CR(num_points,num_points);
465  DMatrix CRPhi(num_points,num_points);
466 
467 
468  // Fill initial matrices for R and RPhi measurements
469  for (unsigned int m=0;m<num_measured;m++){
470  XYZ[m].z=points[m]->wire->origin.z();
471 
472  //Phi=atan2(points[m]->y,points[m]->x);
473  Phi=points[m]->xy.Phi();
474 
475  double u=points[m]->w;
476  double v=points[m]->s;
477  double temp1=u*Phi-v;
478  double temp2=v*Phi+u;
479  double var_u=0.08333;
480  double var_v=0.0075;
481  double one_over_R2=1./points[m]->xy.Mod2();
482  CRPhi(m,m)=one_over_R2*(var_v*temp1*temp1+var_u*temp2*temp2);
483  CR(m,m)=one_over_R2*(var_u*u*u+var_v*v*v);
484 
485  XYZ[m].covr=CR(m,m);
486  XYZ[m].covrphi=CRPhi(m,m);
487  }
488  XYZ[last_index].z=TARGET_Z;
489  CR(last_index,last_index)=BEAM_VARIANCE;
490  CRPhi(last_index,last_index)=BEAM_VARIANCE;
491 
492  // Reference track:
493  jerror_t error=NOERROR;
494  // First find the center and radius of the projected circle
495  error=RiemannCircleFit(points,CRPhi);
496  if (error!=NOERROR){
497  if (DEBUG_LEVEL>0) _DBG_ << "Circle fit failed..." << endl;
498  return error;
499  }
500 
501  // Get reference track estimates for z0 and tanl and intersection points
502  // (stored in XYZ)
503  error=RiemannLineFit(points,CR,XYZ);
504  if (error!=NOERROR){
505  // The circle fit is inconsistent with at least one of the measured points
506  // (i.e., the calcuation of the intersection point failed). Relax the
507  // emphasis on the "target" point and refit.
508  for (unsigned int i=0;i<last_index;i++){
509  CR(i,i)=XYZ[i].covr;
510  CRPhi(i,i)=XYZ[i].covrphi;
511  }
512  CRPhi(last_index,last_index)=1e6;
513  CR(last_index,last_index)=1e6;
514 
515  // First find the center and radius of the projected circle
516  error=RiemannCircleFit(points,CRPhi);
517  if (error!=NOERROR){
518  if (DEBUG_LEVEL>0) _DBG_ << "Circle fit failed..." << endl;
519  return error;
520  }
521 
522  // Get reference track estimates for z0 and tanl and intersection points
523  // (stored in XYZ)
524  error=RiemannLineFit(points,CR,XYZ);
525  if (error!=NOERROR){
526  if (DEBUG_LEVEL>0) _DBG_ << "Line fit failed..." << endl;
527  return error;
528  }
529  }
530 
531  // Guess particle sense of rotation (+/-1);
532  rotation_sense=GetRotationSense(points.size(),XYZ,CR,CRPhi,points);
533 
534  double r1sq=XYZ[ref_plane].xy.Mod2();
535  UpdatePositionsAndCovariance(num_points,r1sq,XYZ,CRPhi,CR);
536 
537  // Compute the chi-square value for this iteration
538  double chisq0=0.;
539  double rc0=rc,xc0=xc,yc0=yc,tanl0=tanl,zvertex0=zvertex,Phi1_0=Phi1;
540  double rotation_sense0=rotation_sense;
541  for (unsigned int m=0;m<points.size();m++){
542  double sperp=rotation_sense*(XYZ[m].z-XYZ[ref_plane].z)/tanl;
543  double phi_s=Phi1+sperp/rc;
544  DVector2 XY(xc+rc*cos(phi_s),yc+rc*sin(phi_s));
545  chisq0+=(XY-points[m]->xy).Mod2()/CR(m,m);
546  }
547 
548  // Preliminary circle fit
549  error=RiemannCircleFit(points,CRPhi);
550  if (error!=NOERROR){
551  if (DEBUG_LEVEL>0) _DBG_ << "Circle fit failed..." << endl;
552  return error;
553  }
554 
555  // Preliminary line fit
556  error=RiemannLineFit(points,CR,XYZ);
557  if (error!=NOERROR){
558  if (DEBUG_LEVEL>0) _DBG_ << "Line fit failed..." << endl;
559  return error;
560  }
561 
562  // Guess particle sense of rotation (+/-1);
563  rotation_sense=GetRotationSense(points.size(),XYZ,CR,CRPhi,points);
564 
565  r1sq=XYZ[ref_plane].xy.Mod2();
566  UpdatePositionsAndCovariance(num_points,r1sq,XYZ,CRPhi,CR);
567 
568  // Compute the chi-square value for this iteration
569  double chisq_=0.;
570  double rotation_sense_=rotation_sense;
571  double rc_=rc,xc_=xc,yc_=yc,tanl_=tanl,zvertex_=zvertex,Phi1_=Phi1;
572  for (unsigned int m=0;m<points.size();m++){
573  double sperp=rotation_sense*(XYZ[m].z-XYZ[ref_plane].z)/tanl;
574  double phi_s=Phi1+sperp/rc;
575  DVector2 XY(xc+rc*cos(phi_s),yc+rc*sin(phi_s));
576  chisq_+=(XY-points[m]->xy).Mod2()/CR(m,m);
577  }
578  // If we did not improve the fit, bail from this routine...
579  if (chisq_>chisq0){
580  chisq=chisq0;
581  rotation_sense=rotation_sense0;
582  rc=rc0,xc=xc0,yc=yc0,tanl=tanl0,zvertex=zvertex0,Phi1=Phi1_0;
583 
584  return NOERROR;
585  }
586 
587  // Final circle fit
588  error=RiemannCircleFit(points,CRPhi);
589  if (error!=NOERROR){
590  if (DEBUG_LEVEL>0) _DBG_ << "Circle fit failed..." << endl;
591  return error;
592  }
593 
594  // Final line fit
595  error=RiemannLineFit(points,CR,XYZ);
596  if (error!=NOERROR){
597  if (DEBUG_LEVEL>0) _DBG_ << "Line fit failed..." << endl;
598  return error;
599  }
600 
601  // Guess particle sense of rotation (+/-1)
602  rotation_sense=GetRotationSense(num_measured,XYZ,CR,CRPhi,points);
603 
604  // Final update to covariance matrices
605  r1sq=XYZ[ref_plane].xy.Mod2();
606  UpdatePositionsAndCovariance(num_points,r1sq,XYZ,CRPhi,CR);
607 
608  // Store residuals and path length for each measurement
609  chisq=0.;
610  for (unsigned int m=0;m<num_points;m++){
611  double sperp=rotation_sense*(XYZ[m].z-XYZ[ref_plane].z)/tanl;
612  double phi_s=Phi1+sperp/rc;
613  DVector2 XY(xc+rc*cos(phi_s),yc+rc*sin(phi_s));
614 
615  if (m<num_measured){
616  chisq+=(XY-points[m]->xy).Mod2()/CR(m,m);
617  }
618  else{
619  chisq+=XY.Mod2()/CR(m,m);
620  }
621  }
622 
623  // If the last iteration did not improve the results, restore the previously
624  // saved values for the parameters and chi-square.
625  if (chisq>chisq_){
626  chisq=chisq_;
627  rotation_sense=rotation_sense_;
628  rc=rc_,xc=xc_,yc=yc_,tanl=tanl_,zvertex=zvertex_,Phi1=Phi1_;
629  }
630 
631  return NOERROR;
632 }
633 
634 
635 // DFDCSegment_factory::FindSegments
636 // Associate nearest neighbors within a package with track segment candidates.
637 // Provide guess for (seed) track parameters
638 //
639 jerror_t DFDCSegment_factory::FindSegments(vector<const DFDCPseudo*>&points){
640  // Create a vector to store whether or not a hit has been used
641  vector<bool>used(points.size());
642  unsigned int total_num_used=0;
643 
644  // Put indices for the first point in each plane before the most downstream
645  // plane in the vector x_list.
646  double old_z=points[0]->wire->origin.z();
647  vector<unsigned int>x_list;
648  x_list.push_back(0);
649  for (unsigned int i=0;i<points.size();i++){
650  if (points[i]->wire->origin.z()!=old_z){
651  x_list.push_back(i);
652  }
653  old_z=points[i]->wire->origin.z();
654  }
655  x_list.push_back(points.size());
656 
657  unsigned int start=0;
658  // loop over the start indices, starting with the first plane
659  while (start<x_list.size()-1){
660  int num_used=0;
661  // For each iteration, count how many hits we have used already in segments
662  for (unsigned int i=0;i<points.size();i++){
663  if(used[i]==true) num_used++;
664  }
665  // Break out of loop if we don't have enough hits left to fit a circle
666  if (points.size()-num_used<3) break;
667 
668  // Now loop over the list of track segment start points
669  for (unsigned int i=x_list[start];i<x_list[start+1];i++){
670  if (used[i]==false){
671  used[i]=true;
672  total_num_used++;
673  chisq=1.e8;
674 
675  // Clear track parameters
676  tanl=D=z0=phi0=Phi1=xc=yc=rc=0.;
677  rotation_sense=1.;
678 
679  // Point in the current plane in the package
680  DVector2 XY=points[i]->xy;
681 
682  // Create list of nearest neighbors
683  vector<const DFDCPseudo*>neighbors;
684  neighbors.push_back(points[i]);
685  unsigned int match=0;
686  double delta,delta_min=1000.;
687  for (unsigned int k=0;k<x_list.size()-1;k++){
688  delta_min=1000.;
689  match=0;
690  if (x_list[k]>0){
691  for (unsigned int m=x_list[k];m<x_list[k+1];m++){
692  delta=(XY-points[m]->xy).Mod();
693  if (delta<delta_min && delta<MATCH_RADIUS){
694  delta_min=delta;
695  match=m;
696  }
697  }
698  }
699  if (match!=0
700  && used[match]==false
701  ){
702  XY=points[match]->xy;
703  used[match]=true;
704  total_num_used++;
705  neighbors.push_back(points[match]);
706  }
707  }
708  unsigned int num_neighbors=neighbors.size();
709 
710  // Skip to next segment seed if we don't have enough points to fit a
711  // circle
712  if (num_neighbors<3) continue;
713 
714  bool do_sort=false;
715  // Look for hits adjacent to the ones we have in our segment candidate
716  for (unsigned int k=0;k<points.size();k++){
717  if (!used[k]){
718  for (unsigned int j=0;j<num_neighbors;j++){
719  delta=fabs(points[k]->xy.Y()-neighbors[j]->xy.Y());
720 
721  if (delta<ADJACENT_MATCH_DISTANCE &&
722  abs(neighbors[j]->wire->wire-points[k]->wire->wire)<=1
723  && neighbors[j]->wire->origin.z()==points[k]->wire->origin.z()){
724  used[k]=true;
725  total_num_used++;
726  neighbors.push_back(points[k]);
727  do_sort=true;
728  }
729  }
730  }
731  }
732  // Sort if we added another point
733  if (do_sort)
734  std::stable_sort(neighbors.begin(),neighbors.end(),DFDCSegment_package_cmp);
735 
736  // Arc lengths in helical model are referenced relative to the plane
737  // ref_plane within a segment. For a 6 hit segment, ref_plane=2 is
738  // roughly in the center of the segment.
739  ref_plane=0;
740 
741  // Perform the Riemann Helical Fit on the track segment
742  jerror_t error=RiemannHelicalFit(neighbors);
743 
744  if (error==NOERROR){
745  // Since the cell size is 0.5 cm and the Lorentz deflection can be
746  // mm scale, a circle radius on the order of 1 cm is at the level of
747  // how well the points are defined at this stage. Use a simple circle
748  // fit assuming the circle goes through the origin.
749  if (rc<1.0) {
750  if (CircleFit(neighbors)==NOERROR){
751  if (LineFit(neighbors)==NOERROR){
752  chisq=ComputeCircleChiSq(neighbors);
753  }
754  }
755  }
756 
757  // Create a new segment
758  DFDCSegment *segment = new DFDCSegment;
759  segment->hits=neighbors;
760  segment->package=(neighbors[0]->wire->layer-1)/6;
761  FillSegmentData(segment);
762 
763  _data.push_back(segment);
764  }
765  else {
766  // Fit assuming particle came from (x,y)=(0,0)
767  if (CircleFit(neighbors)==NOERROR){
768  if (LineFit(neighbors)==NOERROR){
769  chisq=ComputeCircleChiSq(neighbors);
770  }
771  }
772 
773  // Create a new segment
774  DFDCSegment *segment = new DFDCSegment;
775  segment->hits=neighbors;
776  segment->package=(neighbors[0]->wire->layer-1)/6;
777  FillSegmentData(segment);
778 
779  _data.push_back(segment);
780  }
781 
782  }
783  }
784 
785  // Move on to the next plane to start looking for segments
786  start++;
787  } //while loop over x_list planes
788 
789  // Use fitted segment results to look for additional stray hits to add to
790  // segments.
791  if (total_num_used<used.size()){
792  for (unsigned int k=0;k<_data.size();k++){
793  if (total_num_used==used.size()) break;
794  DFDCSegment *segment=_data[k];
795  bool added_hits=false;
796  for (unsigned int i=0;i<used.size();i++){
797  if (total_num_used==used.size()) break;
798  if (used[i]==0){
799  unsigned int packNo=(points[i]->wire->layer-1)/6;
800  if (packNo==segment->package){
801  double z=points[i]->wire->origin.z();
802  double z0=segment->hits[0]->wire->origin.z();
803  double phi_s=segment->Phi1
804  +rotation_sense*(z-z0)/(segment->rc*segment->tanl);
805 
806  double dx=segment->xc+segment->rc*cos(phi_s)-points[i]->xy.X();
807  double dy=segment->yc+segment->rc*sin(phi_s)-points[i]->xy.Y();
808  double dr=sqrt(dx*dx+dy*dy);
809  if (dr<MATCH_RADIUS){
810  used[i]=1;
811  added_hits=true;
812  total_num_used++;
813  segment->hits.push_back(points[i]);
814  }
815  } // check z position of hit
816 
817  } // check if hit used already
818  }// loop over hits
819  if (added_hits){
820  stable_sort(segment->hits.begin(),segment->hits.end(),DFDCSegment_package_cmp);
821  if (RiemannHelicalFit(segment->hits)==NOERROR){
822  // Since the cell size is 0.5 cm and the Lorentz deflection can be
823  // mm scale, a circle radius on the order of 1 cm is at the level of
824  // how well the points are defined at this stage. Use a simple circle
825  // fit assuming the circle goes through the origin.
826  if (rc<1.0) {
827  if (CircleFit(segment->hits)==NOERROR){
828  if (LineFit(segment->hits)==NOERROR){
829  chisq=ComputeCircleChiSq(segment->hits);
830  }
831  }
832  }
833 
834  FillSegmentData(segment);
835  }
836  }
837  } // loop over segments
838  }
839 
840  return NOERROR;
841 }
842 
843 // Linear regression to find charge
844 double DFDCSegment_factory::GetRotationSense(unsigned int n,vector<xyz_t>&XYZ,
845  DMatrix &CR,
846  DMatrix &CRPhi,
847  vector<const DFDCPseudo *>&points){
848  double Phi1=atan2(points[0]->xy.Y()-yc,points[0]->xy.X()-xc);
849  double z0=points[0]->wire->origin.z();
850  double plus_sum=0.,minus_sum=0.;
851  for (unsigned int j=1;j<points.size();j++){
852  double dphi=(points[j]->wire->origin.z()-z0)/(rc*tanl);
853  double my_x=points[j]->xy.X();
854  double my_y=points[j]->xy.Y();
855 
856  double phiplus=Phi1+dphi;
857  double dxplus=xc+rc*cos(phiplus)-my_x;
858  double dyplus=yc+rc*sin(phiplus)-my_y;
859  double dxplus2=dxplus*dxplus;
860  double dyplus2=dyplus*dyplus;
861  double d2plus=dxplus2+dyplus2;
862  double varplus=(dxplus2*points[j]->covxx+dyplus2*points[j]->covyy
863  +2.*dxplus*dyplus*points[j]->covxy)/d2plus;
864 
865  plus_sum+=d2plus/varplus;
866 
867  double phiminus=Phi1-dphi;
868  double dxminus=xc+rc*cos(phiminus)-my_x;
869  double dyminus=yc+rc*sin(phiminus)-my_y;
870  double dxminus2=dxminus*dxminus;
871  double dyminus2=dyminus*dyminus;
872  double d2minus=dxminus2+dyminus2;
873  double varminus=(dxminus2*points[j]->covxx+dyminus2*points[j]->covyy
874  +2.*dxminus*dyminus*points[j]->covxy)/d2minus;
875 
876  minus_sum+=d2minus/varminus;
877  }
878 
879  // Look for smallest sum to determine q
880  if (minus_sum<plus_sum) return -1.;
881  return 1.;
882 }
883 
884 // Linear regression to find z0, tanl
885 jerror_t DFDCSegment_factory::LineFit(vector<const DFDCPseudo *>&points){
886  double sumv=0.,sumx=0.;
887  double sumy=0.,sumxx=0.,sumxy=0.;
888  double sperp=0.,sperp_old=0.,ratio,Delta;
889  double z=0,zlast=0;
890  double two_rc=2.*rc;
891  DVector2 oldxy=points[0]->xy;
892  for (unsigned int k=1;k<points.size();k++){
893  zlast=z;
894  z=points[k]->wire->origin.z();
895  if (fabs(z-zlast)<0.01) continue;
896 
897  DVector2 diffxy=points[k]->xy-oldxy;
898  double Phi=points[k]->xy.Phi();
899  double cosPhi=cos(Phi);
900  double sinPhi=sin(Phi);
901  double var=cosPhi*cosPhi*points[k]->covxx
902  +sinPhi*sinPhi*points[k]->covyy
903  +2.*sinPhi*cosPhi*points[k]->covxy;
904 
905  sperp_old=sperp;
906  ratio=diffxy.Mod()/(two_rc);
907  // Make sure the argument for the arcsin does not go out of range...
908  sperp=sperp_old+(ratio>1?two_rc*(M_PI_2):two_rc*asin(ratio));
909 
910  // Assume errors in s dominated by errors in R
911  double inv_var=1./var;
912  sumv+=inv_var;
913  sumy+=sperp*inv_var;
914  sumx+=z*inv_var;
915  sumxx+=z*z*inv_var;
916  sumxy+=sperp*z*inv_var;
917 
918  // Save the current x and y coordinates
919  oldxy=points[k]->xy;
920  }
921  // last point is at "vertex" --> give it a large error...
922  sperp+=oldxy.Mod()/two_rc;
923  double inv_var=0.01;
924  sumv+=inv_var;
925  sumy+=sperp*inv_var;
926  sumx+=TARGET_Z*inv_var;
927  sumxx+=TARGET_Z*TARGET_Z*inv_var;
928  sumxy+=sperp*TARGET_Z*inv_var;
929 
930  Delta=sumv*sumxx-sumx*sumx;
931  double denom=sumv*sumxy-sumy*sumx;
932  if (fabs(Delta)>EPS && fabs(denom)>EPS){
933  // Track parameters z0 and tan(lambda)
934  tanl=-Delta/denom;
935  // Vertex position
936  sperp-=sperp_old;
937  zvertex=zlast-tanl*sperp;
938 
939  return NOERROR;
940  }
941 
942  return VALUE_OUT_OF_RANGE;
943 }
944 
945 // Simple least squares circle fit forcing the circle to go through (0,0)
946 jerror_t DFDCSegment_factory::CircleFit(vector<const DFDCPseudo *>&points){
947  double alpha=0.0, beta=0.0, gamma=0.0, deltax=0.0, deltay=0.0;
948  // Loop over hits to calculate alpha, beta, gamma, and delta
949  for(unsigned int i=0;i<points.size();i++){
950  const DFDCPseudo *hit = points[i];
951  double x=hit->xy.X();
952  double y=hit->xy.Y();
953  double x_sq=x*x;
954  double y_sq=y*y;
955  double x_sq_plus_y_sq=x_sq+y_sq;
956  alpha += x_sq;
957  beta += y_sq;
958  gamma += x*y;
959  deltax += 0.5*x*x_sq_plus_y_sq;
960  deltay += 0.5*y*x_sq_plus_y_sq;
961  }
962 
963  // Calculate xc,yc - the center of the circle
964  double denom = alpha*beta-gamma*gamma;
965  if(fabs(denom)<1.0E-20)return UNRECOVERABLE_ERROR;
966  xc = (deltax*beta-deltay*gamma)/denom;
967  yc = (deltay*alpha-deltax*gamma)/denom;
968  rc = sqrt(xc*xc + yc*yc);
969  Ndof=points.size()-2;
970 
971  return NOERROR;
972 }
973 
974 // Crude chi^2 calculation for circle fit that was forced to go through (0,0)
975 double DFDCSegment_factory::ComputeCircleChiSq(vector<const DFDCPseudo *>&neighbors) {
976  chisq=0.;
977  Phi1=atan2(neighbors[0]->xy.Y()-yc,neighbors[0]->xy.X()-xc);
978  double z0=neighbors[0]->wire->origin.z();
979  for (unsigned int m=0;m<neighbors.size();m++){
980  double sperp=rotation_sense*(neighbors[m]->wire->origin.z()-z0)/tanl;
981  double phi_s=Phi1+sperp/rc;
982  DVector2 XY(xc+rc*cos(phi_s),yc+rc*sin(phi_s));
983  chisq+=(XY-neighbors[m]->xy).Mod2();
984  // Very crude estimate: we assumed errors of size 1. in circle fit...
985  }
986 
987  return chisq;
988 }
989 
990 
991 // Put fit results into the segment
993  // the particle's charge
994  double charge=RotationSenseToCharge*rotation_sense;
995 
996  // Estimate for azimuthal angle
997  phi0=atan2(-xc,yc);
998  if (rotation_sense<0) phi0+=M_PI;
999 
1000  // Look for distance of closest approach nearest to target
1001  D=-rotation_sense*rc-xc/sin(phi0);
1002 
1003  segment->q=charge; //charge
1004  segment->phi0=phi0; // Phi
1005  segment->D=D; // D=distance of closest approach to origin
1006  segment->tanl=tanl; // tan(lambda), lambda=dip angle
1007  segment->z_vertex=zvertex;// z-position at closest approach to origin
1008  segment->xc=xc;
1009  segment->yc=yc;
1010  segment->rc=rc;
1011  segment->Phi1=Phi1;
1012  segment->chisq=chisq;
1013  segment->Ndof=Ndof;
1014 }
1015 
1016 
1017 //----------------------------------------------------------------------------
1018 // The following routine is no longer in use:
1019 // Correct avalanche position along wire and incorporate drift data for
1020 // coordinate away from the wire using results of preliminary hit-based fit
1021 //
1022 //#define R_START 7.6
1023 //#define Z_TOF 617.4
1024 //#include "HDGEOMETRY/DLorentzMapCalibDB.h
1025 //#define SC_V_EFF 15.
1026 //#define SC_LIGHT_GUIDE 140.
1027 //#define SC_CENTER 38.75
1028 //#define TOF_BAR_LENGTH 252.0
1029 //#define TOF_V_EFF 15.
1030 //#define FDC_X_RESOLUTION 0.028
1031 //#define FDC_Y_RESOLUTION 0.02 //cm
1032 /*
1033 jerror_t DFDCSegment_factory::CorrectPoints(vector<DFDCPseudo*>points,
1034  DMatrix XYZ){
1035  // dip angle
1036  double lambda=atan(tanl);
1037  double alpha=M_PI/2.-lambda;
1038 
1039  if (alpha == 0. || rc==0.) return VALUE_OUT_OF_RANGE;
1040 
1041  // Get Bfield, needed to guess particle momentum
1042  double Bx,By,Bz,B;
1043  double x=XYZ(ref_plane,0);
1044  double y=XYZ(ref_plane,1);
1045  double z=XYZ(ref_plane,2);
1046 
1047  bfield->GetField(x,y,z,Bx,By,Bz);
1048  B=sqrt(Bx*Bx+By*By+Bz*Bz);
1049 
1050  // Momentum and beta
1051  double p=0.002998*B*rc/cos(lambda);
1052  double beta=p/sqrt(p*p+0.140*0.140);
1053 
1054  // Correct start time for propagation from (0,0)
1055  double my_start_time=0.;
1056  if (use_tof){
1057  //my_start_time=ref_time-(Z_TOF-TARGET_Z)/sin(lambda)/beta/29.98;
1058  // If we need to use the tof, the angle relative to the beam line is
1059  // small enough that sin(lambda) ~ 1.
1060  my_start_time=ref_time-(Z_TOF-TARGET_Z)/beta/29.98;
1061  //my_start_time=0;
1062  }
1063  else{
1064  double ratio=R_START/2./rc;
1065  if (ratio<=1.)
1066  my_start_time=ref_time
1067  -2.*rc*asin(R_START/2./rc)*(1./cos(lambda)/beta/29.98);
1068  else
1069  my_start_time=ref_time
1070  -rc*M_PI*(1./cos(lambda)/beta/29.98);
1071 
1072  }
1073  //my_start_time=0.;
1074 
1075  for (unsigned int m=0;m<points.size();m++){
1076  DFDCPseudo *point=points[m];
1077 
1078  double cosangle=point->wire->udir(1);
1079  double sinangle=point->wire->udir(0);
1080  x=XYZ(m,0);
1081  y=XYZ(m,1);
1082  z=point->wire->origin.z();
1083  double delta_x=0,delta_y=0;
1084  // Variances based on expected resolution
1085  double sigx2=FDC_X_RESOLUTION*FDC_X_RESOLUTION;
1086 
1087  // Find difference between point on helical path and wire
1088  double w=x*cosangle-y*sinangle-point->w;
1089  // .. and use it to determine which sign to use for the drift time data
1090  double sign=(w>0?1.:-1.);
1091 
1092  // Correct the drift time for the flight path and convert to distance units
1093  // assuming the particle is a pion
1094  delta_x=sign*(point->time-fdc_track[m].s/beta/29.98-my_start_time)*55E-4;
1095 
1096  // Variance for position along wire. Includes angle dependence from PHENIX
1097  // and transverse diffusion
1098  double sigy2=fdc_y_variance(alpha,delta_x);
1099 
1100  // Next find correction to y from table of deflections
1101  delta_y=lorentz_def->GetLorentzCorrection(x,y,z,alpha,delta_x);
1102 
1103  // Fill x and y elements with corrected values
1104  point->ds =-delta_y;
1105  point->dw =delta_x;
1106  point->x=(point->w+point->dw)*cosangle+(point->s+point->ds)*sinangle;
1107  point->y=-(point->w+point->dw)*sinangle+(point->s+point->ds)*cosangle;
1108  point->covxx=sigx2*cosangle*cosangle+sigy2*sinangle*sinangle;
1109  point->covyy=sigx2*sinangle*sinangle+sigy2*cosangle*cosangle;
1110  point->covxy=(sigy2-sigx2)*sinangle*cosangle;
1111  point->status|=CORRECTED;
1112  }
1113  return NOERROR;
1114 }
1115 */
DVector2 xy
rough x,y coordinates in lab coordinate system
Definition: DFDCPseudo.h:100
jerror_t CircleFit(vector< const DFDCPseudo * > &points)
DApplication * dapp
double chisq
Definition: DFDCSegment.h:55
double rc
Definition: DFDCSegment.h:59
TMatrixD DMatrix
Definition: DMatrix.h:14
double ComputeCircleChiSq(vector< const DFDCPseudo * > &neighbors)
double z_vertex
Definition: DFDCSegment.h:65
double GetRotationSense(unsigned int n, vector< xyz_t > &XYZ, DMatrix &CR, DMatrix &CRPhi, vector< const DFDCPseudo * > &points)
TVector2 DVector2
Definition: DVector2.h:9
jerror_t LineFit(vector< const DFDCPseudo * > &points)
jerror_t brun(JEventLoop *eventLoop, int32_t runnumber)
DFDCSegment_factory::brun():
DFDCSegment_factory()
DFDCSegment_factory::DFDCSegment_factory(): default constructor – initializes log file...
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
#define EPS
double phi0
Definition: DFDCSegment.h:65
double t_v
time of the two cathode clusters
Definition: DFDCPseudo.h:86
#define y
DMagneticFieldMap * GetBfield(unsigned int run_number=1)
const DFDCWire * wire
DFDCWire for this wire.
Definition: DFDCPseudo.h:92
bool DFDCSegment_package_cmp(const DFDCPseudo *a, const DFDCPseudo *b)
#define B0
Definition: src/md5.cpp:23
jerror_t RiemannLineFit(vector< const DFDCPseudo * > &points, DMatrix &CR, vector< xyz_t > &XYZ)
class DFDCPseudo: definition for a reconstructed point in the FDC
Definition: DFDCPseudo.h:74
jerror_t evnt(JEventLoop *eventLoop, uint64_t eventNo)
DFDCSegment_factory::evnt(): this is the place that finds track segments and converts pseudopoints in...
#define X(str)
Definition: hddm-c.cpp:83
double xc
Definition: DFDCSegment.h:59
int layer
1-24
Definition: DFDCWire.h:16
unsigned int package
Definition: DFDCSegment.h:73
#define ONE_THIRD
double t_u
Definition: DFDCPseudo.h:86
const double alpha
vector< const DFDCPseudo * > hits
Definition: DFDCSegment.h:76
~DFDCSegment_factory()
DFDCSegment_factory::~DFDCSegment_factory(): default destructor – closes log file.
jerror_t RiemannCircleFit(vector< const DFDCPseudo * > &points, DMatrix &CRPhi)
double Phi1
Definition: DFDCSegment.h:63
DGeometry * GetDGeometry(unsigned int run_number)
jerror_t RiemannHelicalFit(vector< const DFDCPseudo * > &points)
#define _DBG_
Definition: HDEVIO.h:12
jerror_t UpdatePositionsAndCovariance(unsigned int n, double r1sq, vector< xyz_t > &XYZ, DMatrix &CRPhi, DMatrix &CR)
#define MATCH_RADIUS
#define ADJACENT_MATCH_DISTANCE
double time
time corresponding to this pseudopoint.
Definition: DFDCPseudo.h:93
double sqrt(double)
virtual double GetBz(double x, double y, double z) const =0
double tanl
Definition: DFDCSegment.h:67
double D
Definition: DFDCSegment.h:71
double sin(double)
#define S(str)
Definition: hddm-c.cpp:84
jerror_t FindSegments(vector< const DFDCPseudo * > &points)
double q
Definition: DFDCSegment.h:69
class DFDCSegment: definition for a track segment in the FDC
Definition: DFDCSegment.h:31
TH2F * temp
bool GetTargetZ(double &z_target) const
z-location of center of target
Definition: DGeometry.cc:1933
void FillSegmentData(DFDCSegment *segment)
union @6::@8 u
#define SQRT3
double yc
Definition: DFDCSegment.h:59