Bug Summary

File:libraries/FDC/DFDCSegment_factory.cc
Location:line 555, column 25
Description:Dereference of null pointer

Annotated Source Code

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
10#define HALF_CELL0.5 0.5
11#define MAX_DEFLECTION0.15 0.15
12#define EPS1e-8 1e-8
13#define KILL_RADIUS5.0 5.0
14#define MATCH_RADIUS5.0 5.0
15#define ADJACENT_MATCH_RADIUS2.0 2.0
16#define SIGN_CHANGE_CHISQ_CUT10.0 10.0
17#define BEAM_VARIANCE0.01 0.01 // cm^2
18#define USED_IN_SEGMENT0x8 0x8
19#define CORRECTED0x10 0x10
20#define MAX_ITER10 10
21#define MIN_TANL2.0 2.0
22#define ONE_THIRD0.33333333333333333 0.33333333333333333
23#define SQRT31.73205080756887719 1.73205080756887719
24
25// Variance for position along wire using PHENIX angle dependence, transverse
26// diffusion, and an intrinsic resolution of 127 microns.
27inline double fdc_y_variance(double alpha,double x){
28 return 0.00060+0.0064*tan(alpha)*tan(alpha)+0.0004*fabs(x);
29}
30
31bool DFDCSegment_package_cmp(const DFDCPseudo* a, const DFDCPseudo* b) {
32 return a->wire->layer>b->wire->layer;
33}
34
35DFDCSegment_factory::DFDCSegment_factory() {
36 _log = new JStreamLog(std::cout, "FDCSEGMENT >>");
37 *_log << "File initialized." << endMsg;
38}
39
40
41///
42/// DFDCSegment_factory::~DFDCSegment_factory():
43/// default destructor -- closes log file
44///
45DFDCSegment_factory::~DFDCSegment_factory() {
46 delete _log;
47}
48///
49/// DFDCSegment_factory::brun():
50/// Initialization: read in deflection map, get magnetic field map
51///
52jerror_t DFDCSegment_factory::brun(JEventLoop* eventLoop, int runnumber) {
53 DApplication* dapp=dynamic_cast<DApplication*>(eventLoop->GetJApplication());
54 // get the geometry
55 const DGeometry *geom = dapp->GetDGeometry(runnumber);
56
57 geom->GetTargetZ(TARGET_Z);
58 /*
59 bfield = dapp->GetBfield();
60 lorentz_def=dapp->GetLorentzDeflections();
61
62 *_log << "Table of Lorentz deflections initialized." << endMsg;
63 */
64 return NOERROR;
65}
66
67///
68/// DFDCSegment_factory::evnt():
69/// Routine where pseudopoints are combined into track segments
70///
71jerror_t DFDCSegment_factory::evnt(JEventLoop* eventLoop, int eventNo) {
72 myeventno=eventNo;
73
74 vector<const DFDCPseudo*>pseudopoints;
75 eventLoop->Get(pseudopoints);
76
77 // Skip segment finding if there aren't enough points to form a sensible
78 // segment
79 if (pseudopoints.size()>=3){
80 // Group pseudopoints by package
81 vector<const DFDCPseudo*>package[4];
82 for (vector<const DFDCPseudo*>::const_iterator i=pseudopoints.begin();
83 i!=pseudopoints.end();i++){
84 package[((*i)->wire->layer-1)/6].push_back(*i);
85 }
86
87 // Find the segments in each package
88 for (int j=0;j<4;j++){
89 std::sort(package[j].begin(),package[j].end(),DFDCSegment_package_cmp);
90 // We need at least 3 points to define a circle, so skip if we don't
91 // have enough points.
92 if (package[j].size()>2) FindSegments(package[j]);
93 }
94 } // pseudopoints>2
95
96 return NOERROR;
97}
98
99// Riemann Line fit: linear regression of s on z to determine the tangent of
100// the dip angle and the z position of the closest approach to the beam line.
101// Also returns predicted positions along the helical path.
102//
103jerror_t DFDCSegment_factory::RiemannLineFit(vector<const DFDCPseudo *>points,
104 DMatrix &CR,vector<xyz_t>&XYZ){
105 unsigned int n=points.size()+1;
106 vector<int>bad(n); // Keep track of "bad" intersection points
107 // Fill matrix of intersection points
108 for (unsigned int m=0;m<n-1;m++){
109 double r2=points[m]->xy.Mod2();
110 double denom= N[0]*N[0]+N[1]*N[1];
111 double numer=dist_to_origin+r2*N[2];
112 double ratio=numer/denom;
113
114 DVector2 xy_int0(-N[0]*ratio,-N[1]*ratio);
115 double temp=denom*r2-numer*numer;
116 if (temp<0){
117 bad[m]=1;
118 XYZ[m].xy=xy_int0;
119 continue;
120 }
121 temp=sqrt(temp)/denom;
122
123 // Choose sign of square root based on proximity to actual measurements
124 DVector2 delta(N[1]*temp,-N[0]*temp);
125 DVector2 xy1=xy_int0+delta;
126 DVector2 xy2=xy_int0-delta;
127 double diff1=(xy1-points[m]->xy).Mod2();
128 double diff2=(xy2-points[m]->xy).Mod2();
129 if (diff1>diff2){
130 XYZ[m].xy=xy2;
131 }
132 else{
133 XYZ[m].xy=xy1;
134 }
135 }
136 // Fake target point
137 XYZ[n-1].xy.Set(0.,0.);
138
139 // All arc lengths are measured relative to some reference plane with a hit.
140 // Don't use a "bad" hit for the reference...
141 unsigned int start=0;
142 for (unsigned int i=0;i<bad.size();i++){
143 if (!bad[i]){
144 start=i;
145 break;
146 }
147 }
148 if ((start!=0 && ref_plane==0) || (start!=2&&ref_plane==2)) ref_plane=start;
149
150 // Linear regression to find z0, tanl
151 double sumv=0.,sumx=0.;
152 double sumy=0.,sumxx=0.,sumxy=0.;
153 double sperp=0.,sperp_old=0.,ratio,Delta;
154 double z=0,zlast=0;
155 double two_rc=2.*rc;
156 DVector2 oldxy=XYZ[start].xy;
157 for (unsigned int k=start;k<n;k++){
158 zlast=z;
159 sperp_old=sperp;
160 if (!bad[k]){
161 DVector2 diffxy=XYZ[k].xy-oldxy;
162 ratio=diffxy.Mod()/(two_rc);
163 // Make sure the argument for the arcsin does not go out of range...
164 sperp=sperp_old+(ratio>1?two_rc*(M_PI_21.57079632679489661923):two_rc*asin(ratio));
165 //z=XYZ(k,2);
166 z=XYZ[k].z;
167
168 // Assume errors in s dominated by errors in R
169 double inv_var=1./CR(k,k);
170 sumv+=inv_var;
171 sumy+=sperp*inv_var;
172 sumx+=z*inv_var;
173 sumxx+=z*z*inv_var;
174 sumxy+=sperp*z*inv_var;
175
176 // Save the current x and y coordinates
177 //oldx=XYZ(k,0);
178 //oldy=XYZ(k,1);
179 oldxy=XYZ[k].xy;
180 }
181 }
182
183 Delta=sumv*sumxx-sumx*sumx;
184 double denom=sumv*sumxy-sumy*sumx;
185 if (fabs(Delta)>EPS1e-8 && fabs(denom)>EPS1e-8){
186 // Track parameters z0 and tan(lambda)
187 tanl=-Delta/denom;
188 //z0=(sumxx*sumy-sumx*sumxy)/Delta*tanl;
189 // Error in tanl
190 var_tanl=sumv/Delta*(tanl*tanl*tanl*tanl);
191
192 // Vertex position
193 sperp-=sperp_old;
194 zvertex=zlast-tanl*sperp;
195 }
196 else{
197 // assume that the particle came from the target
198 zvertex=TARGET_Z;
199 if (sperp<EPS1e-8){
200 // This is a rare failure mode where it seems that the arc length data
201 // cannot be reliably used. Provide some default value of tanl to
202 // avoid division by zero errors
203 tanl=57.3; // theta=1 degree
204 }
205 else tanl=(zlast-zvertex)/sperp;
206 var_tanl=100.; // guess for now
207
208
209 }
210
211 return NOERROR;
212}
213
214// Use predicted positions (propagating from plane 1 using a helical model) to
215// update the R and RPhi covariance matrices.
216//
217jerror_t DFDCSegment_factory::UpdatePositionsAndCovariance(unsigned int n,
218 double r1sq,vector<xyz_t>&XYZ,DMatrix &CRPhi,DMatrix &CR){
219 double delta_x=XYZ[ref_plane].xy.X()-xc;
220 double delta_y=XYZ[ref_plane].xy.Y()-yc;
221 double r1=sqrt(r1sq);
222 double denom=delta_x*delta_x+delta_y*delta_y;
223
224 // Auxiliary matrices for correcting for non-normal track incidence to FDC
225 // The correction is
226 // CRPhi'= C*CRPhi*C+S*CR*S, where S(i,i)=R_i*kappa/2
227 // and C(i,i)=sqrt(1-S(i,i)^2)
228 DMatrix C(n,n);
229 DMatrix S(n,n);
230
231 // Predicted positions
232 Phi1=atan2(delta_y,delta_x);
233 double z1=XYZ[ref_plane].z;
234 double y1=XYZ[ref_plane].xy.X();
235 double x1=XYZ[ref_plane].xy.Y();
236 double var_R1=CR(ref_plane,ref_plane);
237 for (unsigned int k=0;k<n;k++){
238 double sperp=charge*(XYZ[k].z-z1)/tanl;
239 double phi_s=Phi1+sperp/rc;
240 double sinp=sin(phi_s);
241 double cosp=cos(phi_s);
242 XYZ[k].xy.Set(xc+rc*cosp,yc+rc*sinp);
243
244 // Error analysis. We ignore errors in N because there doesn't seem to
245 // be any obvious established way to estimate errors in eigenvalues for
246 // small eigenvectors. We assume that most of the error comes from
247 // the measurement for the reference plane radius and the dip angle anyway.
248 double Phi=XYZ[k].xy.Phi();
249 double sinPhi=sin(Phi);
250 double cosPhi=cos(Phi);
251 double dRPhi_dx=Phi*cosPhi-sinPhi;
252 double dRPhi_dy=Phi*sinPhi+cosPhi;
253
254 double rc_sinphi_over_denom=rc*sinp/denom;
255 //double dx_dx1=rc*sinp*delta_y/denom;
256 //double dx_dy1=-rc*sinp*delta_x/denom;
257 double dx_dx1=rc_sinphi_over_denom*delta_y;
258 double dx_dy1=-rc_sinphi_over_denom*delta_x;
259 double dx_dtanl=sinp*sperp/tanl;
260
261 double rc_cosphi_over_denom=rc*cosp/denom;
262 //double dy_dx1=-rc*cosp*delta_y/denom;
263 //double dy_dy1=rc*cosp*delta_x/denom;
264 double dy_dx1=-rc_cosphi_over_denom*delta_y;
265 double dy_dy1=rc_cosphi_over_denom*delta_x;
266 double dy_dtanl=-cosp*sperp/tanl;
267
268 double dRPhi_dx1=dRPhi_dx*dx_dx1+dRPhi_dy*dy_dx1;
269 double dRPhi_dy1=dRPhi_dx*dx_dy1+dRPhi_dy*dy_dy1;
270 double dRPhi_dtanl=dRPhi_dx*dx_dtanl+dRPhi_dy*dy_dtanl;
271
272 double dR_dx1=cosPhi*dx_dx1+sinPhi*dy_dx1;
273 double dR_dy1=cosPhi*dx_dy1+sinPhi*dy_dy1;
274 double dR_dtanl=cosPhi*dx_dtanl+sinPhi*dy_dtanl;
275
276 double cdist=dist_to_origin+r1sq*N[2];
277 double n2=N[0]*N[0]+N[1]*N[1];
278
279 double ydenom=y1*n2+N[1]*cdist;
280 double dy1_dr1=-r1*(2.*N[1]*N[2]*y1+2.*N[2]*cdist-N[0]*N[0])/ydenom;
281 double var_y1=dy1_dr1*dy1_dr1*var_R1;
282
283 double xdenom=x1*n2+N[0]*cdist;
284 double dx1_dr1=-r1*(2.*N[0]*N[2]*x1+2.*N[2]*cdist-N[1]*N[1])/xdenom;
285 double var_x1=dx1_dr1*dx1_dr1*var_R1;
286
287 CRPhi(k,k)=dRPhi_dx1*dRPhi_dx1*var_x1+dRPhi_dy1*dRPhi_dy1*var_y1
288 +dRPhi_dtanl*dRPhi_dtanl*var_tanl;
289 CR(k,k)=dR_dx1*dR_dx1*var_x1+dR_dy1*dR_dy1*var_y1
290 +dR_dtanl*dR_dtanl*var_tanl;
291
292 // double stemp=sqrt(XYZ(k,0)*XYZ(k,0)+XYZ(k,1)*XYZ(k,1))/(4.*rc);
293 double stemp=XYZ[k].xy.Mod()/(4.*rc);
294 double ctemp=1.-stemp*stemp;
295 if (ctemp>0){
296 S(k,k)=stemp;
297 C(k,k)=sqrt(ctemp);
298 }
299 else{
300 S(k,k)=0.;
301 C(k,k)=1.;
302 }
303 }
304
305 // Correction for non-normal incidence of track on FDC
306 CRPhi=C*CRPhi*C+S*CR*S;
307
308 return NOERROR;
309}
310
311// Riemann Circle fit: points on a circle in x,y project onto a plane cutting
312// the circular paraboloid surface described by (x,y,x^2+y^2). Therefore the
313// task of fitting points in (x,y) to a circle is transormed to the taks of
314// fitting planes in (x,y, w=x^2+y^2) space
315//
316jerror_t DFDCSegment_factory::RiemannCircleFit(vector<const DFDCPseudo *>points,
317 DMatrix &CRPhi){
318 unsigned int n=points.size()+1;
319 DMatrix X(n,3);
320 DMatrix Xavg(1,3);
321 DMatrix A(3,3);
322 // vector of ones
323 DMatrix OnesT(1,n);
324 double W_sum=0.;
325 DMatrix W(n,n);
326
327 // The goal is to find the eigenvector corresponding to the smallest
328 // eigenvalue of the equation
329 // lambda=n^T (X^T W X - W_sum Xavg^T Xavg)n
330 // where n is the normal vector to the plane slicing the cylindrical
331 // paraboloid described by the parameterization (x,y,w=x^2+y^2),
332 // and W is the weight matrix, assumed for now to be diagonal.
333 // In the absence of multiple scattering, W_sum is the sum of all the
334 // diagonal elements in W.
335 // At this stage we ignore the multiple scattering.
336 for (unsigned int i=0;i<n-1;i++){
337 X(i,0)=points[i]->xy.X();
338 X(i,1)=points[i]->xy.Y();
339 X(i,2)=points[i]->xy.Mod2();
340 OnesT(0,i)=1.;
341 W(i,i)=1./CRPhi(i,i);
342 W_sum+=W(i,i);
343 }
344 unsigned int last_index=n-1;
345 OnesT(0,last_index)=1.;
346 W(last_index,last_index)=1./CRPhi(last_index,last_index);
347 W_sum+=W(last_index,last_index);
348 var_avg=1./W_sum;
349 Xavg=var_avg*(OnesT*(W*X));
350 // Store in private array for use in other routines
351 xavg[0]=Xavg(0,0);
352 xavg[1]=Xavg(0,1);
353 xavg[2]=Xavg(0,2);
354
355 A=DMatrix(DMatrix::kTransposed,X)*(W*X)
356 -W_sum*(DMatrix(DMatrix::kTransposed,Xavg)*Xavg);
357 if(!A.IsValid())return UNRECOVERABLE_ERROR;
358
359 // The characteristic equation is
360 // lambda^3+B2*lambda^2+lambda*B1+B0=0
361 //
362 double B2=-(A(0,0)+A(1,1)+A(2,2));
363 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)
364 +A(1,1)*A(2,2)-A(2,1)*A(1,2);
365 double B0=-A.Determinant();
366 if(B0==0 || !finite(B0))return UNRECOVERABLE_ERROR;
367
368 // The roots of the cubic equation are given by
369 // lambda1= -B2/3 + S+T
370 // lambda2= -B2/3 - (S+T)/2 + i sqrt(3)/2. (S-T)
371 // lambda3= -B2/3 - (S+T)/2 - i sqrt(3)/2. (S-T)
372 // where we define some temporary variables:
373 // S= (R+sqrt(Q^3+R^2))^(1/3)
374 // T= (R-sqrt(Q^3+R^2))^(1/3)
375 // Q=(3*B1-B2^2)/9
376 // R=(9*B2*B1-27*B0-2*B2^3)/54
377 // sum=S+T;
378 // diff=i*(S-T)
379 // We divide Q and R by a safety factor to prevent multiplying together
380 // enormous numbers that cause unreliable results.
381
382 double Q=(3.*B1-B2*B2)/9.e4;
383 double R=(9.*B2*B1-27.*B0-2.*B2*B2*B2)/54.e6;
384 double Q1=Q*Q*Q+R*R;
385 if (Q1<0) Q1=sqrt(-Q1);
386 else{
387 return VALUE_OUT_OF_RANGE;
388 }
389
390 // DeMoivre's theorem for fractional powers of complex numbers:
391 // (r*(cos(theta)+i sin(theta)))^(p/q)
392 // = r^(p/q)*(cos(p*theta/q)+i sin(p*theta/q))
393 //
394 //double temp=100.*pow(R*R+Q1*Q1,0.16666666666666666667);
395 double temp=100.*sqrt(cbrt(R*R+Q1*Q1));
396 double theta1=ONE_THIRD0.33333333333333333*atan2(Q1,R);
397 double sum_over_2=temp*cos(theta1);
398 double diff_over_2=-temp*sin(theta1);
399 // Third root
400 double lambda_min=-ONE_THIRD0.33333333333333333*B2-sum_over_2+SQRT31.73205080756887719*diff_over_2;
401
402 // Normal vector to plane
403 N[0]=1.;
404 N[1]=(A(1,0)*A(0,2)-(A(0,0)-lambda_min)*A(1,2))
405 /(A(0,1)*A(2,1)-(A(1,1)-lambda_min)*A(0,2));
406 N[2]=(A(2,0)*(A(1,1)-lambda_min)-A(1,0)*A(2,1))
407 /(A(1,2)*A(2,1)-(A(2,2)-lambda_min)*(A(1,1)-lambda_min));
408
409 // Normalize: n1^2+n2^2+n3^2=1
410 double denom=sqrt(N[0]*N[0]+N[1]*N[1]+N[2]*N[2]);
411 for (int i=0;i<3;i++){
412 N[i]/=denom;
413 }
414
415 // Distance to origin
416 dist_to_origin=-(N[0]*Xavg(0,0)+N[1]*Xavg(0,1)+N[2]*Xavg(0,2));
417
418 // Center and radius of the circle
419 double two_N2=2.*N[2];
420 xc=-N[0]/two_N2;
421 yc=-N[1]/two_N2;
422 rc=sqrt(1.-N[2]*N[2]-4.*dist_to_origin*N[2])/fabs(two_N2);
423
424 return NOERROR;
425}
426
427// Riemann Helical Fit based on transforming points on projection to x-y plane
428// to a circular paraboloid surface combined with a linear fit of the arc
429// length versus z. Uses RiemannCircleFit and RiemannLineFit.
430//
431jerror_t DFDCSegment_factory::RiemannHelicalFit(vector<const DFDCPseudo*>points,
432 DMatrix &CR,vector<xyz_t>&XYZ){
433 double Phi;
434 unsigned int num_points=points.size()+1;
435 DMatrix CRPhi(num_points,num_points);
436
437 // Fill initial matrices for R and RPhi measurements
438 XYZ[num_points-1].z=TARGET_Z;
439 for (unsigned int m=0;m<points.size();m++){
440 XYZ[m].z=points[m]->wire->origin.z();
441
442 //Phi=atan2(points[m]->y,points[m]->x);
443 Phi=points[m]->xy.Phi();
444 double cosPhi=cos(Phi);
445 double sinPhi=sin(Phi);
446 double Phi_cosPhi=Phi*cosPhi;
447 double Phi_sinPhi=Phi*sinPhi;
448 double Phi_cosPhi_minus_sinPhi=Phi_cosPhi-sinPhi;
449 double Phi_sinPhi_plus_cosPhi=Phi_sinPhi+cosPhi;
450 CRPhi(m,m)
451 =Phi_cosPhi_minus_sinPhi*Phi_cosPhi_minus_sinPhi*points[m]->covxx
452 +Phi_sinPhi_plus_cosPhi*Phi_sinPhi_plus_cosPhi*points[m]->covyy
453 +2.*Phi_sinPhi_plus_cosPhi*Phi_cosPhi_minus_sinPhi*points[m]->covxy;
454
455 CR(m,m)=cosPhi*cosPhi*points[m]->covxx
456 +sinPhi*sinPhi*points[m]->covyy
457 +2.*sinPhi*cosPhi*points[m]->covxy;
458 }
459 CR(points.size(),points.size())=BEAM_VARIANCE0.01;
460 CRPhi(points.size(),points.size())=BEAM_VARIANCE0.01;
461
462 // Reference track:
463 jerror_t error=NOERROR;
464 // First find the center and radius of the projected circle
465 error=RiemannCircleFit(points,CRPhi);
466 if (error!=NOERROR) return error;
467
468 // Get reference track estimates for z0 and tanl and intersection points
469 // (stored in XYZ)
470 error=RiemannLineFit(points,CR,XYZ);
471 if (error!=NOERROR) return error;
472
473 // Guess particle charge (+/-1);
474 charge=GetCharge(points.size(),XYZ,CR,CRPhi);
475
476 double r1sq=XYZ[ref_plane].xy.Mod2();
477 UpdatePositionsAndCovariance(num_points,r1sq,XYZ,CRPhi,CR);
478
479 // Preliminary circle fit
480 error=RiemannCircleFit(points,CRPhi);
481 if (error!=NOERROR) return error;
482
483 // Preliminary line fit
484 error=RiemannLineFit(points,CR,XYZ);
485 if (error!=NOERROR) return error;
486
487 // Guess particle charge (+/-1);
488 charge=GetCharge(points.size(),XYZ,CR,CRPhi);
489
490 r1sq=XYZ[ref_plane].xy.Mod2();
491 UpdatePositionsAndCovariance(num_points,r1sq,XYZ,CRPhi,CR);
492
493 // Final circle fit
494 error=RiemannCircleFit(points,CRPhi);
495 if (error!=NOERROR) return error;
496
497 // Final line fit
498 error=RiemannLineFit(points,CR,XYZ);
499 if (error!=NOERROR) return error;
500
501 // Guess particle charge (+/-1)
502 charge=GetCharge(points.size(),XYZ,CR,CRPhi);
503
504 // Final update to covariance matrices
505 ref_plane=0;
506 r1sq=XYZ[ref_plane].xy.Mod2();
507 UpdatePositionsAndCovariance(num_points,r1sq,XYZ,CRPhi,CR);
508
509 // Store residuals and path length for each measurement
510 chisq=0.;
511 for (unsigned int m=0;m<points.size();m++){
512 double sperp=charge*(XYZ[m].z-XYZ[ref_plane].z)/tanl;
513 double phi_s=Phi1+sperp/rc;
514 DVector2 XY(xc+rc*cos(phi_s),yc+rc*sin(phi_s));
515 chisq+=(XY-points[m]->xy).Mod2()/CR(m,m);
516 }
517 return NOERROR;
518}
519
520
521// DFDCSegment_factory::FindSegments
522// Associate nearest neighbors within a package with track segment candidates.
523// Provide guess for (seed) track parameters
524//
525jerror_t DFDCSegment_factory::FindSegments(vector<const DFDCPseudo*>points){
526 // Clear all the "used" flags;
527 used.clear();
528
529 // Put indices for the first point in each plane before the most downstream
530 // plane in the vector x_list.
531 double old_z=points[0]->wire->origin.z();
532 vector<unsigned int>x_list;
533 x_list.push_back(0);
534 for (unsigned int i=0;i<points.size();i++){
1
Loop condition is false. Execution continues on line 541
→
535 used.push_back(false);
536 if (points[i]->wire->origin.z()!=old_z){
537 x_list.push_back(i);
538 }
539 old_z=points[i]->wire->origin.z();
540 }
541 x_list.push_back(points.size());
542
543 unsigned int start=0;
544 // loop over the start indices, starting with the first plane
545 while (start<x_list.size()-1){
2
←
Loop condition is true. Entering loop body
→
546 int num_used=0;
547 // For each iteration, count how many hits we have used already in segments
548 for (unsigned int i=0;i<points.size();i++){
3
←
Loop condition is false. Execution continues on line 552
→
549 if(used[i]==true) num_used++;
550 }
551 // Break out of loop if we don't have enough hits left to fit a circle
552 if (points.size()-num_used<3) break;
4
←
Taking false branch
→
553
554 // Now loop over the list of track segment start points
555 for (unsigned int i=x_list[start];i<x_list[start+1];i++){
5
←
Dereference of null pointer
556 if (used[i]==false){
557 used[i]=true;
558 chisq=1.e8;
559
560 // Clear track parameters
561 tanl=D=z0=phi0=Phi1=xc=yc=rc=0.;
562 charge=1.;
563
564 // Point in the current plane in the package
565 DVector2 XY=points[i]->xy;
566
567 // Create list of nearest neighbors
568 vector<const DFDCPseudo*>neighbors;
569 neighbors.push_back(points[i]);
570 unsigned int match=0;
571 double delta,delta_min=1000.;
572 for (unsigned int k=0;k<x_list.size()-1;k++){
573 delta_min=1000.;
574 match=0;
575 for (unsigned int m=x_list[k];m<x_list[k+1];m++){
576 delta=(XY-points[m]->xy).Mod();
577 if (delta<delta_min && delta<MATCH_RADIUS5.0){
578 delta_min=delta;
579 match=m;
580 }
581 }
582 if (match!=0
583 && used[match]==false
584 ){
585 XY=points[match]->xy;
586 used[match]=true;
587 neighbors.push_back(points[match]);
588 }
589 }
590 unsigned int num_neighbors=neighbors.size();
591
592 // Skip to next segment seed if we don't have enough points to fit a
593 // circle
594 if (num_neighbors<3) continue;
595
596 bool do_sort=false;
597 // Look for hits adjacent to the ones we have in our segment candidate
598 for (unsigned int k=0;k<points.size();k++){
599 if (!used[k]){
600 for (unsigned int j=0;j<num_neighbors;j++){
601 delta=(points[k]->xy-neighbors[j]->xy).Mod();
602
603 if (delta<ADJACENT_MATCH_RADIUS2.0 &&
604 abs(neighbors[j]->wire->wire-points[k]->wire->wire)<=1
605 && neighbors[j]->wire->origin.z()==points[k]->wire->origin.z()){
606 used[k]=true;
607 neighbors.push_back(points[k]);
608 do_sort=true;
609 }
610 }
611 }
612 }
613 // Sort if we added another point
614 if (do_sort)
615 std::sort(neighbors.begin(),neighbors.end(),DFDCSegment_package_cmp);
616
617 // list of points on track and the corresponding covariances
618 unsigned int num_points_on_segment=neighbors.size()+1;
619 vector<xyz_t>XYZ(num_points_on_segment);
620 DMatrix CR(num_points_on_segment,num_points_on_segment);
621
622 // Arc lengths in helical model are referenced relative to the plane
623 // ref_plane within a segment. For a 6 hit segment, ref_plane=2 is
624 // roughly in the center of the segment.
625 ref_plane=2;
626
627 // Perform the Riemann Helical Fit on the track segment
628 jerror_t error=RiemannHelicalFit(neighbors,CR,XYZ); /// initial hit-based fit
629
630 if (error==NOERROR){
631 // Estimate for azimuthal angle
632 phi0=atan2(-xc,yc);
633 if (charge<0) phi0+=M_PI3.14159265358979323846;
634
635 // Look for distance of closest approach nearest to target
636 D=-charge*rc-xc/sin(phi0);
637
638 // Creat a new segment
639 DFDCSegment *segment = new DFDCSegment;
640
641 // Initialize seed track parameters
642 segment->q=charge; //charge
643 segment->phi0=phi0; // Phi
644 segment->D=D; // D=distance of closest approach to origin
645 segment->tanl=tanl; // tan(lambda), lambda=dip angle
646 segment->z_vertex=zvertex;// z-position at closest approach to origin
647
648 segment->hits=neighbors;
649 segment->xc=xc;
650 segment->yc=yc;
651 segment->rc=rc;
652 segment->Phi1=Phi1;
653 segment->chisq=chisq;
654
655 //printf("tanl %f \n",tanl);
656 //printf("xc %f yc %f rc %f\n",xc,yc,rc);
657
658 _data.push_back(segment);
659 }
660 }
661 }
662 // Look for a new plane to start looking for a segment
663 while (start<x_list.size()-1){
664 if (used[x_list[start]]==false) break;
665 start++;
666 }
667 } //while loop over x_list planes
668
669 return NOERROR;
670}
671
672// Linear regression to find charge
673double DFDCSegment_factory::GetCharge(unsigned int n,vector<xyz_t>&XYZ,
674 DMatrix &CR,
675 DMatrix &CRPhi){
676 double inv_var=0.;
677 double sumv=0.;
678 double sumy=0.;
679 double sumx=0.;
680 double sumxx=0.,sumxy=0,Delta;
681 double slope,r2;
682 double phi_old=XYZ[0].xy.Phi();
683 for (unsigned int k=0;k<n;k++){
684 double tempz=XYZ[k].z;
685 double phi_z=XYZ[k].xy.Phi();
686 // Check for problem regions near +pi and -pi
687 if (fabs(phi_z-phi_old)>M_PI3.14159265358979323846){
688 if (phi_old<0) phi_z-=2.*M_PI3.14159265358979323846;
689 else phi_z+=2.*M_PI3.14159265358979323846;
690 }
691 r2=XYZ[k].xy.Mod2();
692 inv_var=r2/(CRPhi(k,k)+phi_z*phi_z*CR(k,k));
693 sumv+=inv_var;
694 sumy+=phi_z*inv_var;
695 sumx+=tempz*inv_var;
696 sumxx+=tempz*tempz*inv_var;
697 sumxy+=phi_z*tempz*inv_var;
698 phi_old=phi_z;
699 }
700 Delta=sumv*sumxx-sumx*sumx;
701 slope=(sumv*sumxy-sumy*sumx)/Delta;
702
703 // Guess particle charge (+/-1);
704 if (slope<0.) return -1.;
705 return 1.;
706}
707
708//----------------------------------------------------------------------------
709// The following routine is no longer in use:
710// Correct avalanche position along wire and incorporate drift data for
711// coordinate away from the wire using results of preliminary hit-based fit
712//
713//#define R_START 7.6
714//#define Z_TOF 617.4
715//#include "HDGEOMETRY/DLorentzMapCalibDB.h
716//#define SC_V_EFF 15.
717//#define SC_LIGHT_GUIDE 140.
718//#define SC_CENTER 38.75
719//#define TOF_BAR_LENGTH 252.0
720//#define TOF_V_EFF 15.
721//#define FDC_X_RESOLUTION 0.028
722//#define FDC_Y_RESOLUTION 0.02 //cm
723/*
724jerror_t DFDCSegment_factory::CorrectPoints(vector<DFDCPseudo*>points,
725 DMatrix XYZ){
726 // dip angle
727 double lambda=atan(tanl);
728 double alpha=M_PI/2.-lambda;
729
730 if (alpha == 0. || rc==0.) return VALUE_OUT_OF_RANGE;
731
732 // Get Bfield, needed to guess particle momentum
733 double Bx,By,Bz,B;
734 double x=XYZ(ref_plane,0);
735 double y=XYZ(ref_plane,1);
736 double z=XYZ(ref_plane,2);
737
738 bfield->GetField(x,y,z,Bx,By,Bz);
739 B=sqrt(Bx*Bx+By*By+Bz*Bz);
740
741 // Momentum and beta
742 double p=0.002998*B*rc/cos(lambda);
743 double beta=p/sqrt(p*p+0.140*0.140);
744
745 // Correct start time for propagation from (0,0)
746 double my_start_time=0.;
747 if (use_tof){
748 //my_start_time=ref_time-(Z_TOF-TARGET_Z)/sin(lambda)/beta/29.98;
749 // If we need to use the tof, the angle relative to the beam line is
750 // small enough that sin(lambda) ~ 1.
751 my_start_time=ref_time-(Z_TOF-TARGET_Z)/beta/29.98;
752 //my_start_time=0;
753 }
754 else{
755 double ratio=R_START/2./rc;
756 if (ratio<=1.)
757 my_start_time=ref_time
758 -2.*rc*asin(R_START/2./rc)*(1./cos(lambda)/beta/29.98);
759 else
760 my_start_time=ref_time
761 -rc*M_PI*(1./cos(lambda)/beta/29.98);
762
763 }
764 //my_start_time=0.;
765
766 for (unsigned int m=0;m<points.size();m++){
767 DFDCPseudo *point=points[m];
768
769 double cosangle=point->wire->udir(1);
770 double sinangle=point->wire->udir(0);
771 x=XYZ(m,0);
772 y=XYZ(m,1);
773 z=point->wire->origin.z();
774 double delta_x=0,delta_y=0;
775 // Variances based on expected resolution
776 double sigx2=FDC_X_RESOLUTION*FDC_X_RESOLUTION;
777
778 // Find difference between point on helical path and wire
779 double w=x*cosangle-y*sinangle-point->w;
780 // .. and use it to determine which sign to use for the drift time data
781 double sign=(w>0?1.:-1.);
782
783 // Correct the drift time for the flight path and convert to distance units
784 // assuming the particle is a pion
785 delta_x=sign*(point->time-fdc_track[m].s/beta/29.98-my_start_time)*55E-4;
786
787 // Variance for position along wire. Includes angle dependence from PHENIX
788 // and transverse diffusion
789 double sigy2=fdc_y_variance(alpha,delta_x);
790
791 // Next find correction to y from table of deflections
792 delta_y=lorentz_def->GetLorentzCorrection(x,y,z,alpha,delta_x);
793
794 // Fill x and y elements with corrected values
795 point->ds =-delta_y;
796 point->dw =delta_x;
797 point->x=(point->w+point->dw)*cosangle+(point->s+point->ds)*sinangle;
798 point->y=-(point->w+point->dw)*sinangle+(point->s+point->ds)*cosangle;
799 point->covxx=sigx2*cosangle*cosangle+sigy2*sinangle*sinangle;
800 point->covyy=sigx2*sinangle*sinangle+sigy2*cosangle*cosangle;
801 point->covxy=(sigy2-sigx2)*sinangle*cosangle;
802 point->status|=CORRECTED;
803 }
804 return NOERROR;
805}
806*/