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