| File: | libraries/TRACKING/DQuickFit.cc |
| Location: | line 621, column 8 |
| Description: | Value stored to 'dzdphi' during its initialization is never read |
| 1 | // $Id: DQuickFit.cc 4274 2008-10-01 11:51:20Z staylor $ |
| 2 | // Author: David Lawrence June 25, 2004 |
| 3 | // |
| 4 | // |
| 5 | // |
| 6 | |
| 7 | #include <iostream> |
| 8 | #include <algorithm> |
| 9 | using namespace std; |
| 10 | |
| 11 | #include <math.h> |
| 12 | |
| 13 | #include "DQuickFit.h" |
| 14 | #include "DRiemannFit.h" |
| 15 | #define qBr2p0.003 0.003 // conversion for converting q*B*r to GeV/c |
| 16 | #define Z_VERTEX65.0 65.0 |
| 17 | |
| 18 | // The following is for sorting hits by z |
| 19 | class DQFHitLessThanZ{ |
| 20 | public: |
| 21 | bool operator()(DQFHit_t* const &a, DQFHit_t* const &b) const { |
| 22 | return a->z < b->z; |
| 23 | } |
| 24 | }; |
| 25 | |
| 26 | bool DQFHitLessThanZ_C(DQFHit_t* const &a, DQFHit_t* const &b) { |
| 27 | return a->z < b->z; |
| 28 | } |
| 29 | |
| 30 | //----------------- |
| 31 | // DQuickFit (Constructor) |
| 32 | //----------------- |
| 33 | DQuickFit::DQuickFit(void) |
| 34 | { |
| 35 | x0 = y0 = 0; |
| 36 | chisq = 0; |
| 37 | chisq_source = NOFIT; |
| 38 | bfield = NULL__null; |
| 39 | |
| 40 | c_origin=0.; |
| 41 | normal.SetXYZ(0.,0.,0.); |
| 42 | } |
| 43 | |
| 44 | //----------------- |
| 45 | // DQuickFit (Constructor) |
| 46 | //----------------- |
| 47 | DQuickFit::DQuickFit(const DQuickFit &fit) |
| 48 | { |
| 49 | Copy(fit); |
| 50 | } |
| 51 | //----------------- |
| 52 | // Copy |
| 53 | //----------------- |
| 54 | void DQuickFit::Copy(const DQuickFit &fit) |
| 55 | { |
| 56 | x0 = fit.x0; |
| 57 | y0 = fit.y0; |
| 58 | q = fit.q; |
| 59 | p = fit.p; |
| 60 | p_trans = fit.p_trans; |
| 61 | phi = fit.phi; |
| 62 | theta = fit.theta; |
| 63 | z_vertex = fit.z_vertex; |
| 64 | chisq = fit.chisq; |
| 65 | dzdphi = fit.dzdphi; |
| 66 | chisq_source = fit.chisq_source; |
| 67 | bfield = fit.GetMagneticFieldMap(); |
| 68 | Bz_avg = fit.GetBzAvg(); |
| 69 | z_mean = fit.GetZMean(); |
| 70 | phi_mean = fit.GetPhiMean(); |
| 71 | |
| 72 | const vector<DQFHit_t*> myhits = fit.GetHits(); |
| 73 | for(unsigned int i=0; i<myhits.size(); i++){ |
| 74 | DQFHit_t *a = new DQFHit_t; |
| 75 | *a = *myhits[i]; |
| 76 | hits.push_back(a); |
| 77 | } |
| 78 | } |
| 79 | |
| 80 | //----------------- |
| 81 | // operator= (Assignment operator) |
| 82 | //----------------- |
| 83 | DQuickFit& DQuickFit::operator=(const DQuickFit& fit) |
| 84 | { |
| 85 | if(this == &fit)return *this; |
| 86 | Copy(fit); |
| 87 | |
| 88 | return *this; |
| 89 | } |
| 90 | |
| 91 | //----------------- |
| 92 | // DQuickFit (Destructor) |
| 93 | //----------------- |
| 94 | DQuickFit::~DQuickFit() |
| 95 | { |
| 96 | Clear(); |
| 97 | } |
| 98 | |
| 99 | //----------------- |
| 100 | // AddHit |
| 101 | //----------------- |
| 102 | jerror_t DQuickFit::AddHit(float r, float phi, float z) |
| 103 | { |
| 104 | /// Add a hit to the list of hits using cylindrical coordinates |
| 105 | |
| 106 | return AddHitXYZ(r*cos(phi), r*sin(phi), z); |
| 107 | } |
| 108 | |
| 109 | //----------------- |
| 110 | // AddHitXYZ |
| 111 | //----------------- |
| 112 | jerror_t DQuickFit::AddHitXYZ(float x, float y, float z) |
| 113 | { |
| 114 | /// Add a hit to the list of hits useing Cartesian coordinates |
| 115 | DQFHit_t *hit = new DQFHit_t; |
| 116 | hit->x = x; |
| 117 | hit->y = y; |
| 118 | hit->z = z; |
| 119 | hit->chisq = 0.0; |
| 120 | hits.push_back(hit); |
| 121 | |
| 122 | return NOERROR; |
| 123 | } |
| 124 | |
| 125 | //----------------- |
| 126 | // PruneHit |
| 127 | //----------------- |
| 128 | jerror_t DQuickFit::PruneHit(int idx) |
| 129 | { |
| 130 | /// Remove the hit specified by idx from the list |
| 131 | /// of hits. The value of idx can be anywhere from |
| 132 | /// 0 to GetNhits()-1. |
| 133 | if(idx<0 || idx>=(int)hits.size())return VALUE_OUT_OF_RANGE; |
| 134 | |
| 135 | delete hits[idx]; |
| 136 | hits.erase(hits.begin() + idx); |
| 137 | |
| 138 | return NOERROR; |
| 139 | } |
| 140 | |
| 141 | //----------------- |
| 142 | // Clear |
| 143 | //----------------- |
| 144 | jerror_t DQuickFit::Clear(void) |
| 145 | { |
| 146 | /// Remove all hits |
| 147 | for(unsigned int i=0; i<hits.size(); i++)delete hits[i]; |
| 148 | hits.clear(); |
| 149 | |
| 150 | return NOERROR; |
| 151 | } |
| 152 | |
| 153 | //----------------- |
| 154 | // PrintChiSqVector |
| 155 | //----------------- |
| 156 | jerror_t DQuickFit::PrintChiSqVector(void) const |
| 157 | { |
| 158 | /// Dump the latest chi-squared vector to the screen. |
| 159 | /// This prints the individual hits' chi-squared |
| 160 | /// contributions in the order in which the hits were |
| 161 | /// added. See PruneHits() for more detail. |
| 162 | |
| 163 | cout<<"Chisq vector from DQuickFit: (source="; |
| 164 | switch(chisq_source){ |
| 165 | case NOFIT: cout<<"NOFIT"; break; |
| 166 | case CIRCLE: cout<<"CIRCLE"; break; |
| 167 | case TRACK: cout<<"TRACK"; break; |
| 168 | }; |
| 169 | cout<<")"<<endl; |
| 170 | cout<<"----------------------------"<<endl; |
| 171 | |
| 172 | for(unsigned int i=0;i<hits.size();i++){ |
| 173 | cout<<i<<" "<<hits[i]->chisq<<endl; |
| 174 | } |
| 175 | cout<<"Total: "<<chisq<<endl<<endl; |
| 176 | |
| 177 | return NOERROR; |
| 178 | } |
| 179 | |
| 180 | //----------------- |
| 181 | // FitCircle |
| 182 | //----------------- |
| 183 | jerror_t DQuickFit::FitCircle(void) |
| 184 | { |
| 185 | /// Fit the current set of hits to a circle |
| 186 | /// |
| 187 | /// This takes the hits which have been added thus far via one |
| 188 | /// of the AddHit() methods and fits them to a circle. |
| 189 | /// The alogorithm used here calculates the parameters directly using |
| 190 | /// a technique very much like linear regression. The key assumptions |
| 191 | /// are: |
| 192 | /// 1. The magnetic field is uniform and along z so that the projection |
| 193 | /// of the track onto the X/Y plane will fall on a circle |
| 194 | /// (this also implies no multiple-scattering or energy loss) |
| 195 | /// 2. The vertex is at the target center (i.e. 0,0 in the coordinate |
| 196 | /// system of the points passed to us. |
| 197 | /// |
| 198 | /// IMPORTANT: The value of phi which results from this assumes |
| 199 | /// the particle was positively charged. If the particle was |
| 200 | /// negatively charged, then phi will be 180 degrees off. To |
| 201 | /// correct this, one needs z-coordinate information to determine |
| 202 | /// the sign of the charge. |
| 203 | /// |
| 204 | /// ALSO IMPORTANT: This assumes a charge of +1 for the particle. If |
| 205 | /// the particle actually has a charge of +2, then the resulting |
| 206 | /// value of p_trans will be half of what it should be. |
| 207 | |
| 208 | float alpha=0.0, beta=0.0, gamma=0.0, deltax=0.0, deltay=0.0; |
| 209 | chisq_source = NOFIT; // in case we reutrn early |
| 210 | |
| 211 | // Loop over hits to calculate alpha, beta, gamma, and delta |
| 212 | // if a magnetic field map was given, use it to find average Z B-field |
| 213 | DQFHit_t *a = NULL__null; |
| 214 | for(unsigned int i=0;i<hits.size();i++){ |
| 215 | a = hits[i]; |
| 216 | float x=a->x; |
| 217 | float y=a->y; |
| 218 | alpha += x*x; |
| 219 | beta += y*y; |
| 220 | gamma += x*y; |
| 221 | deltax += x*(x*x+y*y)/2.0; |
| 222 | deltay += y*(x*x+y*y)/2.0; |
| 223 | } |
| 224 | |
| 225 | // Calculate x0,y0 - the center of the circle |
| 226 | double denom = alpha*beta-gamma*gamma; |
| 227 | if(fabs(denom)<1.0E-20)return UNRECOVERABLE_ERROR; |
| 228 | x0 = (deltax*beta-deltay*gamma)/denom; |
| 229 | //y0 = (deltay-gamma*x0)/beta; |
| 230 | y0 = (deltay*alpha-deltax*gamma)/denom; |
| 231 | |
| 232 | // Calculate the phi angle traversed by the track from the |
| 233 | // origin to the last hit. NOTE: this can be off by a multiple |
| 234 | // of 2pi! |
| 235 | double delta_phi=0.0; |
| 236 | if(a){ // a should be pointer to last hit from above loop |
| 237 | delta_phi = atan2(a->y-y0, a->x-x0); |
| 238 | if(delta_phi<0.0)delta_phi += 2.0*M_PI3.14159265358979323846; |
| 239 | } |
| 240 | |
| 241 | // Momentum depends on magnetic field. If bfield has been |
| 242 | // set, we should use it to determine an average value of Bz |
| 243 | // for this track. Otherwise, assume -2T. |
| 244 | // Also assume a singly charged track (i.e. q=+/-1) |
| 245 | // The sign of the charge will be determined below. |
| 246 | Bz_avg=-2.0; |
| 247 | q = +1.0; |
| 248 | r0 = sqrt(x0*x0 + y0*y0); |
| 249 | p_trans = q*Bz_avg*r0*qBr2p0.003; // qBr2p converts to GeV/c |
| 250 | phi = atan2(y0,x0) - M_PI_21.57079632679489661923; |
| 251 | if(p_trans<0.0){ |
| 252 | p_trans = -p_trans; |
| 253 | } |
| 254 | if(phi<0)phi+=2.0*M_PI3.14159265358979323846; |
| 255 | if(phi>=2.0*M_PI3.14159265358979323846)phi-=2.0*M_PI3.14159265358979323846; |
| 256 | |
| 257 | // Calculate the chisq |
| 258 | ChisqCircle(); |
| 259 | chisq_source = CIRCLE; |
| 260 | |
| 261 | return NOERROR; |
| 262 | } |
| 263 | |
| 264 | //----------------- |
| 265 | // ChisqCircle |
| 266 | //----------------- |
| 267 | double DQuickFit::ChisqCircle(void) |
| 268 | { |
| 269 | /// Calculate the chisq for the hits and current circle |
| 270 | /// parameters. |
| 271 | /// NOTE: This does not return the chi2/dof, just the |
| 272 | /// raw chi2 with errs set to 1.0 |
| 273 | chisq = 0.0; |
| 274 | for(unsigned int i=0;i<hits.size();i++){ |
| 275 | DQFHit_t *a = hits[i]; |
| 276 | float x = a->x - x0; |
| 277 | float y = a->y - y0; |
| 278 | float c = sqrt(x*x + y*y) - r0; |
| 279 | c *= c; |
| 280 | a->chisq = c; |
| 281 | chisq+=c; |
| 282 | } |
| 283 | |
| 284 | // Do NOT divide by DOF |
| 285 | return chisq; |
| 286 | } |
| 287 | |
| 288 | //----------------- |
| 289 | // FitCircleRiemann |
| 290 | //----------------- |
| 291 | jerror_t DQuickFit::FitCircleRiemann(double BeamRMS) |
| 292 | { |
| 293 | /// This is a temporary solution for doing a Riemann circle fit. |
| 294 | /// It uses Simon's DRiemannFit class. This is not very efficient |
| 295 | /// since it creates a new DRiemann fit object, copies the hits |
| 296 | /// into it, and then copies the reults back to this DQuickFit |
| 297 | /// object every time this is called. At some point, the DRiemannFit |
| 298 | /// class will be merged into DQuickFit. |
| 299 | |
| 300 | chisq_source = NOFIT; // in case we reutrn early |
| 301 | |
| 302 | DRiemannFit rfit; |
| 303 | for(unsigned int i=0; i<hits.size(); i++){ |
| 304 | DQFHit_t *hit = hits[i]; |
| 305 | rfit.AddHitXYZ(hit->x, hit->y, hit->z); |
| 306 | } |
| 307 | |
| 308 | // Fake point at origin |
| 309 | double beam_var=BeamRMS*BeamRMS; |
| 310 | rfit.AddHit(0.,0.,Z_VERTEX65.0,beam_var,beam_var,0.); |
| 311 | |
| 312 | jerror_t err = rfit.FitCircle(BeamRMS); |
| 313 | if(err!=NOERROR)return err; |
| 314 | |
| 315 | x0 = rfit.xc; |
| 316 | y0 = rfit.yc; |
| 317 | phi = rfit.phi; |
| 318 | |
| 319 | //r0 = sqrt(x0*x0+y0*y0); |
| 320 | r0=rfit.rc; |
| 321 | |
| 322 | rfit.GetPlaneParameters(c_origin,normal); |
| 323 | |
| 324 | // Momentum depends on magnetic field. If bfield has been |
| 325 | // set, we should use it to determine an average value of Bz |
| 326 | // for this track. Otherwise, assume -2T. |
| 327 | // Also assume a singly charged track (i.e. q=+/-1) |
| 328 | // The sign of the charge will be determined below. |
| 329 | Bz_avg=-2.0; |
| 330 | q = +1.0; |
| 331 | float r0 = sqrt(x0*x0 + y0*y0); |
| 332 | p_trans = q*Bz_avg*r0*qBr2p0.003; // qBr2p converts to GeV/c |
| 333 | phi = atan2(y0,x0) - M_PI_21.57079632679489661923; |
| 334 | if(p_trans<0.0){ |
| 335 | p_trans = -p_trans; |
| 336 | } |
| 337 | if(phi<0)phi+=2.0*M_PI3.14159265358979323846; |
| 338 | if(phi>=2.0*M_PI3.14159265358979323846)phi-=2.0*M_PI3.14159265358979323846; |
| 339 | |
| 340 | // Calculate the chisq |
| 341 | ChisqCircle(); |
| 342 | chisq_source = CIRCLE; |
| 343 | |
| 344 | return NOERROR; |
| 345 | } |
| 346 | |
| 347 | |
| 348 | //----------------- |
| 349 | // FitCircleStraightTrack |
| 350 | //----------------- |
| 351 | jerror_t DQuickFit::FitCircleStraightTrack(void) |
| 352 | { |
| 353 | /// This is a last resort! |
| 354 | /// The circle fit can fail for high momentum tracks that are nearly |
| 355 | /// straight tracks. In these cases (when pt~2GeV or more) there |
| 356 | /// is not enough position resolution with wire positions only |
| 357 | /// to measure the curvature of the track. |
| 358 | /// We can though, fit the X/Y points with a straight line in order |
| 359 | /// to get phi and project out the stereo angles. |
| 360 | /// |
| 361 | /// For the radius of the circle (i.e. p_trans) we loop over momenta |
| 362 | /// from 1 to 9GeV and just use the one with the smallest chisq. |
| 363 | |
| 364 | // Average the x's and y's of the individual hits. We should be |
| 365 | // able to do this via linear regression, but I can't get it to |
| 366 | // work right now and I'm under the gun to get ready for a review |
| 367 | // so I take the easy way. Note that we don't average phi's since |
| 368 | // that will cause problems at the 0/2pi boundary. |
| 369 | double X=0.0, Y=0.0; |
| 370 | DQFHit_t *a = NULL__null; |
| 371 | for(unsigned int i=0;i<hits.size();i++){ |
| 372 | a = hits[i]; |
| 373 | double r = sqrt(pow((double)a->x,2.0) + pow((double)a->y, 2.0)); |
| 374 | // weight by r to give outer points more influence. Note that |
| 375 | // we really are really weighting by r^2 since x and y already |
| 376 | // have a magnitude component. |
| 377 | X += a->x*r; |
| 378 | Y += a->y*r; |
| 379 | } |
| 380 | phi = atan2(Y,X); |
| 381 | if(phi<0)phi+=2.0*M_PI3.14159265358979323846; |
| 382 | if(phi>=2.0*M_PI3.14159265358979323846)phi-=2.0*M_PI3.14159265358979323846; |
| 383 | |
| 384 | // Search the chi2 space for values for p_trans, x0, ... |
| 385 | SearchPtrans(9.0, 0.5); |
| 386 | |
| 387 | #if 0 |
| 388 | // We do a simple linear regression here to find phi. This is |
| 389 | // simplified by the intercept being zero (i.e. the track |
| 390 | // passes through the beamline). |
| 391 | float Sxx=0.0, Syy=0.0, Sxy=0.0; |
| 392 | chisq_source = NOFIT; // in case we reutrn early |
| 393 | |
| 394 | // Loop over hits to calculate Sxx, Syy, and Sxy |
| 395 | DQFHit_t *a = NULL__null; |
| 396 | for(unsigned int i=0;i<hits.size();i++){ |
| 397 | a = hits[i]; |
| 398 | float x=a->x; |
| 399 | float y=a->y; |
| 400 | Sxx += x*x; |
| 401 | Syy += y*y; |
| 402 | Sxy += x*y; |
| 403 | } |
| 404 | double A = 2.0*Sxy; |
| 405 | double B = Sxx - Syy; |
| 406 | phi = B>A ? atan2(A,B)/2.0 : 1.0/atan2(B,A)/2.0; |
| 407 | if(phi<0)phi+=2.0*M_PI3.14159265358979323846; |
| 408 | if(phi>=2.0*M_PI3.14159265358979323846)phi-=2.0*M_PI3.14159265358979323846; |
| 409 | #endif |
| 410 | |
| 411 | return NOERROR; |
| 412 | } |
| 413 | |
| 414 | //----------------- |
| 415 | // SearchPtrans |
| 416 | //----------------- |
| 417 | void DQuickFit::SearchPtrans(double ptrans_max, double ptrans_step) |
| 418 | { |
| 419 | /// Search the chi2 space as a function of the transverse momentum |
| 420 | /// for the minimal chi2. The values corresponding to the minmal |
| 421 | /// chi2 are left in chisq, x0, y0, r0, q, and p_trans. |
| 422 | /// |
| 423 | /// This does NOT optimize the value of p_trans. It simply |
| 424 | /// does a straight forward chisq calculation on a grid |
| 425 | /// and keeps the best one. It is intended for use in finding |
| 426 | /// a reasonable value for p_trans for straight tracks that |
| 427 | /// are contained to less than p_trans_max which is presumably |
| 428 | /// chosen based on a known θ. |
| 429 | |
| 430 | // Loop over several values for p_trans and calculate the |
| 431 | // chisq for each. |
| 432 | double min_chisq=1.0E6; |
| 433 | double min_x0=0.0, min_y0=0.0, min_r0=0.0; |
| 434 | for(double my_p_trans=ptrans_step; my_p_trans<=ptrans_max; my_p_trans+=ptrans_step){ |
| 435 | |
| 436 | r0 = my_p_trans/(0.003*2.0); |
| 437 | double alpha, my_chisq; |
| 438 | |
| 439 | // q = +1 |
| 440 | alpha = phi + M_PI_21.57079632679489661923; |
| 441 | x0 = r0*cos(alpha); |
| 442 | y0 = r0*sin(alpha); |
| 443 | my_chisq = ChisqCircle(); |
| 444 | if(my_chisq<min_chisq){ |
| 445 | min_chisq=my_chisq; |
| 446 | min_x0 = x0; |
| 447 | min_y0 = y0; |
| 448 | min_r0 = r0; |
| 449 | q = +1.0; |
| 450 | p_trans = my_p_trans; |
| 451 | } |
| 452 | |
| 453 | // q = -1 |
| 454 | alpha = phi - M_PI_21.57079632679489661923; |
| 455 | x0 = r0*cos(alpha); |
| 456 | y0 = r0*sin(alpha); |
| 457 | my_chisq = ChisqCircle(); |
| 458 | if(my_chisq<min_chisq){ |
| 459 | min_chisq=my_chisq; |
| 460 | min_x0 = x0; |
| 461 | min_y0 = y0; |
| 462 | min_r0 = r0; |
| 463 | q = -1.0; |
| 464 | p_trans = my_p_trans; |
| 465 | } |
| 466 | } |
| 467 | |
| 468 | // Copy params from minimum chisq |
| 469 | x0 = min_x0; |
| 470 | y0 = min_y0; |
| 471 | r0 = min_r0; |
| 472 | } |
| 473 | |
| 474 | //----------------- |
| 475 | // QuickPtrans |
| 476 | //----------------- |
| 477 | void DQuickFit::QuickPtrans(void) |
| 478 | { |
| 479 | /// Quickly calculate a value of p_trans by looking |
| 480 | /// for the hit furthest out and the hit closest |
| 481 | /// to half that distance. Those 2 hits along with the |
| 482 | /// origin are used to define a circle from which |
| 483 | /// p_trans is calculated. |
| 484 | |
| 485 | // Find hit with largest R |
| 486 | double R2max = 0.0; |
| 487 | DQFHit_t *hit_max = NULL__null; |
| 488 | for(unsigned int i=0; i<hits.size(); i++){ |
| 489 | // use cross product to decide if hits is to left or right |
| 490 | double x = hits[i]->x; |
| 491 | double y = hits[i]->y; |
| 492 | double R2 = x*x + y*y; |
| 493 | if(R2>R2max){ |
| 494 | R2max = R2; |
| 495 | hit_max = hits[i]; |
| 496 | } |
| 497 | } |
| 498 | |
| 499 | // Bullet proof |
| 500 | if(!hit_max)return; |
| 501 | |
| 502 | // Find hit closest to half-way out |
| 503 | double Rmid = 0.0; |
| 504 | double Rmax = sqrt(R2max); |
| 505 | DQFHit_t *hit_mid = NULL__null; |
| 506 | for(unsigned int i=0; i<hits.size(); i++){ |
| 507 | // use cross product to decide if hits is to left or right |
| 508 | double x = hits[i]->x; |
| 509 | double y = hits[i]->y; |
| 510 | double R = sqrt(x*x + y*y); |
| 511 | if(fabs(R-Rmax/2.0) < fabs(Rmid-Rmax/2.0)){ |
| 512 | Rmid = R; |
| 513 | hit_mid = hits[i]; |
| 514 | } |
| 515 | } |
| 516 | |
| 517 | // Bullet proof |
| 518 | if(!hit_mid)return; |
| 519 | |
| 520 | // Calculate p_trans from 3 points |
| 521 | double x1 = hit_mid->x; |
| 522 | double y1 = hit_mid->y; |
| 523 | double x2 = hit_max->x; |
| 524 | double y2 = hit_max->y; |
| 525 | double r2 = sqrt(x2*x2 + y2*y2); |
| 526 | double cos_phi = x2/r2; |
| 527 | double sin_phi = y2/r2; |
| 528 | double u1 = x1*cos_phi + y1*sin_phi; |
| 529 | double v1 = -x1*sin_phi + y1*cos_phi; |
| 530 | double u2 = x2*cos_phi + y2*sin_phi; |
| 531 | double u0 = u2/2.0; |
| 532 | double v0 = (u1*u1 + v1*v1 - u1*u2)/(2.0*v1); |
| 533 | |
| 534 | x0 = u0*cos_phi - v0*sin_phi; |
| 535 | y0 = u0*sin_phi + v0*cos_phi; |
| 536 | r0 = sqrt(x0*x0 + y0*y0); |
| 537 | |
| 538 | double B=-2.0; |
| 539 | p_trans = fabs(0.003*B*r0); |
| 540 | q = v1>0.0 ? -1.0:+1.0; |
| 541 | } |
| 542 | |
| 543 | //----------------- |
| 544 | // GuessChargeFromCircleFit |
| 545 | //----------------- |
| 546 | jerror_t DQuickFit::GuessChargeFromCircleFit(void) |
| 547 | { |
| 548 | /// Adjust the sign of the charge (magnitude will stay 1) based on |
| 549 | /// whether the hits tend to be to the right or to the left of the |
| 550 | /// line going from the origin through the center of the circle. |
| 551 | /// If the sign is flipped, the phi angle will also be shifted by |
| 552 | /// +/- pi since the majority of hits are assumed to come from |
| 553 | /// outgoing tracks. |
| 554 | /// |
| 555 | /// This is just a guess since tracks can bend all the way |
| 556 | /// around and have hits that look exactly like tracks from an |
| 557 | /// outgoing particle of opposite charge. The final charge should |
| 558 | /// come from the sign of the dphi/dz slope. |
| 559 | |
| 560 | // Simply count the number of hits whose phi angle relative to the |
| 561 | // phi of the center of the circle are greater than pi. |
| 562 | unsigned int N=0; |
| 563 | for(unsigned int i=0; i<hits.size(); i++){ |
| 564 | // use cross product to decide if hits is to left or right |
| 565 | double x = hits[i]->x; |
| 566 | double y = hits[i]->y; |
| 567 | if((x*y0 - y*x0) < 0.0)N++; |
| 568 | } |
| 569 | |
| 570 | // Check if more hits are negative and make sign negative if so. |
| 571 | if(N>hits.size()/2.0){ |
| 572 | q = -1.0; |
| 573 | phi += M_PI3.14159265358979323846; |
| 574 | if(phi>2.0*M_PI3.14159265358979323846)phi-=2.0*M_PI3.14159265358979323846; |
| 575 | } |
| 576 | |
| 577 | return NOERROR; |
| 578 | } |
| 579 | |
| 580 | //----------------- |
| 581 | // FitTrack |
| 582 | //----------------- |
| 583 | jerror_t DQuickFit::FitTrack(void) |
| 584 | { |
| 585 | /// Find theta, sign of electric charge, total momentum and |
| 586 | /// vertex z position. |
| 587 | |
| 588 | // Points must be in order of increasing Z |
| 589 | sort(hits.begin(), hits.end(), DQFHitLessThanZ_C); |
| 590 | |
| 591 | // Fit to circle to get circle's center |
| 592 | FitCircle(); |
| 593 | |
| 594 | // The thing that is really needed is dphi/dz (where phi is the angle |
| 595 | // of the point as measured from the center of the circle, not the beam |
| 596 | // line). The relation between phi and z is linear so we use linear |
| 597 | // regression to find the slope (dphi/dz). The one complication is |
| 598 | // that phi is periodic so the value obtained via the x and y of a |
| 599 | // specific point can be off by an integral number of 2pis. Handle this |
| 600 | // by assuming the first point is measured before a full rotation |
| 601 | // has occurred. Subsequent points should always be within 2pi of |
| 602 | // the previous point so we just need to calculate the relative phi |
| 603 | // between succesive points and keep a running sum. We do this in |
| 604 | // the first loop were we find the mean z and phi values. The regression |
| 605 | // formulae are calculated in the second loop. |
| 606 | |
| 607 | // Calculate phi about circle center for each hit |
| 608 | Fill_phi_circle(hits, x0, y0); |
| 609 | |
| 610 | // Linear regression for z/phi relation |
| 611 | float Sxx=0.0, Syy=0.0, Sxy=0.0; |
| 612 | for(unsigned int i=0;i<hits.size();i++){ |
| 613 | DQFHit_t *a = hits[i]; |
| 614 | float deltaZ = a->z - z_mean; |
| 615 | float deltaPhi = a->phi_circle - phi_mean; |
| 616 | Syy += deltaZ*deltaZ; |
| 617 | Sxx += deltaPhi*deltaPhi; |
| 618 | Sxy += deltaZ*deltaPhi; |
| 619 | //cout<<__FILE__<<":"<<__LINE__<<" deltaZ="<<deltaZ<<" deltaPhi="<<deltaPhi<<" Sxy(i)="<<deltaZ*deltaPhi<<endl; |
| 620 | } |
| 621 | float dzdphi = Syy/Sxy; |
Value stored to 'dzdphi' during its initialization is never read | |
| 622 | dzdphi = Syy/Sxy; |
| 623 | z_vertex = z_mean - phi_mean*dzdphi; |
| 624 | //cout<<__FILE__<<":"<<__LINE__<<" z_mean="<<z_mean<<" phi_mean="<<phi_mean<<" dphidz="<<dphidz<<" Sxy="<<Sxy<<" Syy="<<Syy<<" z_vertex="<<z_vertex<<endl; |
| 625 | |
| 626 | // Fill in the rest of the parameters |
| 627 | return FillTrackParams(); |
| 628 | } |
| 629 | |
| 630 | //----------------- |
| 631 | // FitTrack_FixedZvertex |
| 632 | //----------------- |
| 633 | jerror_t DQuickFit::FitTrack_FixedZvertex(float z_vertex) |
| 634 | { |
| 635 | /// Fit the points, but hold the z_vertex fixed at the specified value. |
| 636 | /// |
| 637 | /// This just calls FitCircle and FitLine_FixedZvertex to find all parameters |
| 638 | /// of the track |
| 639 | |
| 640 | // Fit to circle to get circle's center |
| 641 | FitCircle(); |
| 642 | |
| 643 | // Fit to line in phi-z plane and return error |
| 644 | return FitLine_FixedZvertex(z_vertex); |
| 645 | } |
| 646 | |
| 647 | //----------------- |
| 648 | // FitLine_FixedZvertex |
| 649 | //----------------- |
| 650 | jerror_t DQuickFit::FitLine_FixedZvertex(float z_vertex) |
| 651 | { |
| 652 | /// Fit the points, but hold the z_vertex fixed at the specified value. |
| 653 | /// |
| 654 | /// This assumes FitCircle has already been called and the values |
| 655 | /// in x0 and y0 are valid. |
| 656 | /// |
| 657 | /// This just fits the phi-z angle by minimizing the chi-squared |
| 658 | /// using a linear regression technique. As it turns out, the |
| 659 | /// chi-squared weights points by their distances squared which |
| 660 | /// leads to a quadratic equation for cot(theta) (where theta is |
| 661 | /// the angle in the phi-z plane). To pick the right solution of |
| 662 | /// the quadratic equation, we pick the one closest to the linearly |
| 663 | /// weighted angle. One could argue that we should just use the |
| 664 | /// linear weighting here rather than the square weighting. The |
| 665 | /// choice depends though on how likely the "out-lier" points |
| 666 | /// are to really be from this track. If they are likely, to |
| 667 | /// be, then we would want to leverage their "longer" lever arms |
| 668 | /// with the squared weighting. If they are more likely to be bad |
| 669 | /// points, then we would want to minimize their influence with |
| 670 | /// a linear (or maybe even root) weighting. It is expected than |
| 671 | /// for our use, the points are all likely valid so we use the |
| 672 | /// square weighting. |
| 673 | |
| 674 | // Points must be in order of increasing Z |
| 675 | sort(hits.begin(), hits.end(), DQFHitLessThanZ_C); |
| 676 | |
| 677 | // Fit is being done for a fixed Z-vertex |
| 678 | this->z_vertex = z_vertex; |
| 679 | |
| 680 | // Calculate phi about circle center for each hit |
| 681 | Fill_phi_circle(hits, x0, y0); |
| 682 | |
| 683 | // Do linear regression on phi-Z |
| 684 | float Sx=0, Sy=0; |
| 685 | float Sxx=0, Syy=0, Sxy = 0; |
| 686 | float r0 = sqrt(x0*x0 + y0*y0); |
| 687 | for(unsigned int i=0; i<hits.size(); i++){ |
| 688 | DQFHit_t *a = hits[i]; |
| 689 | |
| 690 | float dz = a->z - z_vertex; |
| 691 | float dphi = a->phi_circle; |
| 692 | Sx += dz; |
| 693 | Sy += r0*dphi; |
| 694 | Syy += r0*dphi*r0*dphi; |
| 695 | Sxx += dz*dz; |
| 696 | Sxy += r0*dphi*dz; |
| 697 | } |
| 698 | |
| 699 | float k = (Syy-Sxx)/(2.0*Sxy); |
| 700 | float s = sqrt(k*k + 1.0); |
| 701 | float cot_theta1 = -k+s; |
| 702 | float cot_theta2 = -k-s; |
| 703 | float cot_theta_lin = Sx/Sy; |
| 704 | float cot_theta; |
| 705 | if(fabs(cot_theta_lin-cot_theta1) < fabs(cot_theta_lin-cot_theta2)){ |
| 706 | cot_theta = cot_theta1; |
| 707 | }else{ |
| 708 | cot_theta = cot_theta2; |
| 709 | } |
| 710 | |
| 711 | dzdphi = r0*cot_theta; |
| 712 | |
| 713 | // Fill in the rest of the paramters |
| 714 | return FillTrackParams(); |
| 715 | } |
| 716 | |
| 717 | //------------------------------------------------------------------ |
| 718 | // Fill_phi_circle |
| 719 | //------------------------------------------------------------------ |
| 720 | jerror_t DQuickFit::Fill_phi_circle(vector<DQFHit_t*> hits, float x0, float y0) |
| 721 | { |
| 722 | float x_last = -x0; |
| 723 | float y_last = -y0; |
| 724 | float phi_last = 0.0; |
| 725 | z_mean = phi_mean = 0.0; |
| 726 | for(unsigned int i=0; i<hits.size(); i++){ |
| 727 | DQFHit_t *a = hits[i]; |
| 728 | |
| 729 | float dx = a->x - x0; |
| 730 | float dy = a->y - y0; |
| 731 | float dphi = atan2f(dx*y_last - dy*x_last, dx*x_last + dy*y_last)atan2((double)dx*y_last - dy*x_last,(double)dx*x_last + dy*y_last ); |
| 732 | float my_phi = phi_last +dphi; |
| 733 | a->phi_circle = my_phi; |
| 734 | |
| 735 | z_mean += a->z; |
| 736 | phi_mean += my_phi; |
| 737 | |
| 738 | x_last = dx; |
| 739 | y_last = dy; |
| 740 | phi_last = my_phi; |
| 741 | } |
| 742 | z_mean /= (float)hits.size(); |
| 743 | phi_mean /= (float)hits.size(); |
| 744 | |
| 745 | return NOERROR; |
| 746 | } |
| 747 | |
| 748 | //------------------------------------------------------------------ |
| 749 | // FillTrackParams |
| 750 | //------------------------------------------------------------------ |
| 751 | jerror_t DQuickFit::FillTrackParams(void) |
| 752 | { |
| 753 | /// Fill in and tweak some parameters like q, phi, theta using |
| 754 | /// other values already set in the class. This is used by |
| 755 | /// both FitTrack() and FitTrack_FixedZvertex(). |
| 756 | |
| 757 | float r0 = sqrt(x0*x0 + y0*y0); |
| 758 | theta = atan(r0/fabs(dzdphi)); |
| 759 | |
| 760 | // The sign of the electric charge will be opposite that |
| 761 | // of dphi/dz. Also, the value of phi will be PI out of phase |
| 762 | if(dzdphi<0.0){ |
| 763 | phi += M_PI3.14159265358979323846; |
| 764 | if(phi<0)phi+=2.0*M_PI3.14159265358979323846; |
| 765 | if(phi>=2.0*M_PI3.14159265358979323846)phi-=2.0*M_PI3.14159265358979323846; |
| 766 | }else{ |
| 767 | q = -q; |
| 768 | } |
| 769 | |
| 770 | // Re-calculate chi-sq (needs to be done) |
| 771 | chisq_source = TRACK; |
| 772 | |
| 773 | // Up to now, the fit has assumed a forward going particle. In |
| 774 | // other words, if the particle is going backwards, the helix does |
| 775 | // actually still go through the points, but only if extended |
| 776 | // backward in time. We use the average z of the hits compared |
| 777 | // to the reconstructed vertex to determine if it was back-scattered. |
| 778 | if(z_mean<z_vertex){ |
| 779 | // back scattered particle |
| 780 | theta = M_PI3.14159265358979323846 - theta; |
| 781 | phi += M_PI3.14159265358979323846; |
| 782 | if(phi<0)phi+=2.0*M_PI3.14159265358979323846; |
| 783 | if(phi>=2.0*M_PI3.14159265358979323846)phi-=2.0*M_PI3.14159265358979323846; |
| 784 | q = -q; |
| 785 | } |
| 786 | |
| 787 | // There is a problem with theta for data generated with a uniform |
| 788 | // field. The following factor is emprical until I can figure out |
| 789 | // what the source of the descrepancy is. |
| 790 | //theta += 0.1*pow((double)sin(theta),2.0); |
| 791 | |
| 792 | p = fabs(p_trans/sin(theta)); |
| 793 | |
| 794 | return NOERROR; |
| 795 | } |
| 796 | |
| 797 | //------------------------------------------------------------------ |
| 798 | |
| 799 | //------------------------------------------------------------------ |
| 800 | jerror_t DQuickFit::Print(void) const |
| 801 | { |
| 802 | cout<<"-- DQuickFit Params ---------------"<<endl; |
| 803 | cout<<" x0 = "<<x0<<endl; |
| 804 | cout<<" y0 = "<<y0<<endl; |
| 805 | cout<<" q = "<<q<<endl; |
| 806 | cout<<" p = "<<p<<endl; |
| 807 | cout<<" p_trans = "<<p_trans<<endl; |
| 808 | cout<<" phi = "<<phi<<endl; |
| 809 | cout<<" theta = "<<theta<<endl; |
| 810 | cout<<" z_vertex = "<<z_vertex<<endl; |
| 811 | cout<<" chisq = "<<chisq<<endl; |
| 812 | cout<<"chisq_source = "; |
| 813 | switch(chisq_source){ |
| 814 | case NOFIT: cout<<"NOFIT"; break; |
| 815 | case CIRCLE: cout<<"CIRCLE"; break; |
| 816 | case TRACK: cout<<"TRACK"; break; |
| 817 | } |
| 818 | cout<<endl; |
| 819 | cout<<" Bz(avg) = "<<Bz_avg<<endl; |
| 820 | |
| 821 | return NOERROR; |
| 822 | } |
| 823 | |
| 824 | |
| 825 | //------------------------------------------------------------------ |
| 826 | // Dump |
| 827 | //------------------------------------------------------------------ |
| 828 | jerror_t DQuickFit::Dump(void) const |
| 829 | { |
| 830 | Print(); |
| 831 | |
| 832 | for(unsigned int i=0;i<hits.size();i++){ |
| 833 | DQFHit_t *v = hits[i]; |
| 834 | cout<<" x="<<v->x<<" y="<<v->y<<" z="<<v->z; |
| 835 | cout<<" phi_circle="<<v->phi_circle<<" chisq="<<v->chisq<<endl; |
| 836 | } |
| 837 | |
| 838 | return NOERROR; |
| 839 | } |
| 840 | |
| 841 | //------------------------------------------------------------------- |
| 842 | //------------------------------------------------------------------- |
| 843 | //--------------------------- UNUSED -------------------------------- |
| 844 | //------------------------------------------------------------------- |
| 845 | //------------------------------------------------------------------- |
| 846 | |
| 847 | |
| 848 | #if 0 |
| 849 | |
| 850 | //----------------- |
| 851 | // AddHits |
| 852 | //----------------- |
| 853 | jerror_t DQuickFit::AddHits(int N, DVector3 *v) |
| 854 | { |
| 855 | /// Append a list of hits to the current list of hits using |
| 856 | /// DVector3 objects. The DVector3 objects are copied internally |
| 857 | /// so it is safe to delete the objects after calling AddHits() |
| 858 | /// For 2D hits, the value of z will be ignored. |
| 859 | |
| 860 | for(int i=0; i<N; i++, v++){ |
| 861 | DVector3 *vec = new DVector3(*v); |
| 862 | if(vec) |
| 863 | hits.push_back(vec); |
| 864 | else |
| 865 | cerr<<__FILE__"DQuickFit.cc"<<":"<<__LINE__865<<" NULL vector in DQuickFit::AddHits. Hit dropped."<<endl; |
| 866 | } |
| 867 | |
| 868 | return NOERROR; |
| 869 | } |
| 870 | |
| 871 | //----------------- |
| 872 | // PruneHits |
| 873 | //----------------- |
| 874 | jerror_t DQuickFit::PruneHits(float chisq_limit) |
| 875 | { |
| 876 | /// Remove hits whose individual contribution to the chi-squared |
| 877 | /// value exceeds <i>chisq_limit</i>. The value of the individual |
| 878 | /// contribution is calculated like: |
| 879 | /// |
| 880 | /// \f$ r_0 = \sqrt{x_0^2 + y_0^2} \f$ <br> |
| 881 | /// \f$ r_i = \sqrt{(x_i-x_0)^2 + (y_i-y_0)^2} \f$ <br> |
| 882 | /// \f$ \chi_i^2 = (r_i-r_0)^2 \f$ <br> |
| 883 | |
| 884 | if(hits.size() != chisqv.size()){ |
| 885 | cerr<<__FILE__"DQuickFit.cc"<<":"<<__LINE__885<<" hits and chisqv do not have the same number of rows!"<<endl; |
| 886 | cerr<<"Call FitCircle() or FitTrack() method first!"<<endl; |
| 887 | |
| 888 | return NOERROR; |
| 889 | } |
| 890 | |
| 891 | // Loop over these backwards to make it easier to |
| 892 | // handle the deletes |
| 893 | for(int i=chisqv.size()-1; i>=0; i--){ |
| 894 | if(chisqv[i] > chisq_limit)PruneHit(i); |
| 895 | } |
| 896 | |
| 897 | return NOERROR; |
| 898 | } |
| 899 | |
| 900 | //----------------- |
| 901 | // PruneOutliers |
| 902 | //----------------- |
| 903 | jerror_t DQuickFit::PruneOutlier(void) |
| 904 | { |
| 905 | /// Remove the point which is furthest from the geometric |
| 906 | /// center of all the points in the X/Y plane. |
| 907 | /// |
| 908 | /// This just calculates the average x and y values of |
| 909 | /// registered hits. It finds the distance of every hit |
| 910 | /// with respect to the geometric mean and removes the |
| 911 | /// hit whose furthest from the mean. |
| 912 | |
| 913 | float X=0, Y=0; |
| 914 | for(unsigned int i=0;i<hits.size();i++){ |
| 915 | DVector3 *v = hits[i]; |
| 916 | X += v->x(); |
| 917 | Y += v->y(); |
| 918 | } |
| 919 | X /= (float)hits.size(); |
| 920 | Y /= (float)hits.size(); |
| 921 | |
| 922 | float max =0.0; |
| 923 | int idx = -1; |
| 924 | for(unsigned int i=0;i<hits.size();i++){ |
| 925 | DVector3 *v = hits[i]; |
| 926 | float x = v->x()-X; |
| 927 | float y = v->y()-Y; |
| 928 | float dist_sq = x*x + y*y; // we don't need to take sqrt just to find max |
| 929 | if(dist_sq>max){ |
| 930 | max = dist_sq; |
| 931 | idx = i; |
| 932 | } |
| 933 | } |
| 934 | if(idx>=0)PruneHit(idx); |
| 935 | |
| 936 | return NOERROR; |
| 937 | } |
| 938 | |
| 939 | //----------------- |
| 940 | // PruneOutliers |
| 941 | //----------------- |
| 942 | jerror_t DQuickFit::PruneOutliers(int n) |
| 943 | { |
| 944 | /// Remove the n points which are furthest from the geometric |
| 945 | /// center of all the points in the X/Y plane. |
| 946 | /// |
| 947 | /// This just calls PruneOutlier() n times. Since the mean |
| 948 | /// can change each time, it should be recalculated after |
| 949 | /// every hit is removed. |
| 950 | |
| 951 | for(int i=0;i<n;i++)PruneOutlier(); |
| 952 | |
| 953 | return NOERROR; |
| 954 | } |
| 955 | |
| 956 | //----------------- |
| 957 | // firstguess_curtis |
| 958 | //----------------- |
| 959 | jerror_t DCDC::firstguess_curtis(s_Cdc_trackhit_t *hits, int Npoints |
| 960 | , float &theta ,float &phi, float &p, float &q) |
| 961 | { |
| 962 | if(Npoints<3)return NOERROR; |
| 963 | // pick out 3 points to calculate the circle with |
| 964 | s_Cdc_trackhit_t *hit1, *hit2, *hit3; |
| 965 | hit1 = hits; |
| 966 | hit2 = &hits[Npoints/2]; |
| 967 | hit3 = &hits[Npoints-1]; |
| 968 | |
| 969 | float x1 = hit1->x, y1=hit1->y; |
| 970 | float x2 = hit2->x, y2=hit2->y; |
| 971 | float x3 = hit3->x, y3=hit3->y; |
| 972 | |
| 973 | float b = (x2*x2+y2*y2-x1*x1-y1*y1)/(2.0*(x2-x1)); |
| 974 | b -= (x3*x3+y3*y3-x1*x1-y1*y1)/(2.0*(x3-x1)); |
| 975 | b /= ((y1-y2)/(x1-x2)) - ((y1-y3)/(x1-x3)); |
| 976 | float a = (x2*x2-y2*y2-x1*x1-y1*y1)/(2.0*(x2-x1)); |
| 977 | a -= b*(y2-y1)/(x2-x1); |
| 978 | |
| 979 | // Above is the method from Curtis's crystal barrel note 93, pg 72 |
| 980 | // Below here is just a copy of the code from David's firstguess |
| 981 | // routine above (after the x0,y0 calculation) |
| 982 | float x0=a, y0=b; |
| 983 | |
| 984 | // Momentum depends on magnetic field. Assume 2T for now. |
| 985 | // Also assume a singly charged track (i.e. q=+/-1) |
| 986 | // The sign of the charge will be deterined below. |
| 987 | float B=-2.0*0.593; // The 0.61 is empirical |
| 988 | q = +1.0; |
| 989 | float r0 = sqrt(x0*x0 + y0*y0); |
| 990 | float hbarc = 197.326; |
| 991 | float p_trans = q*B*r0/hbarc; // are these the right units? |
| 992 | |
| 993 | // Assuming points are ordered in increasing z, the sign of the |
| 994 | // cross-product between successive points will be the opposite |
| 995 | // sign of the charge. Since it's possible to have a few "bad" |
| 996 | // points, we don't want to rely on any one to determine this. |
| 997 | // The method we use is to sum cross-products of the first and |
| 998 | // middle points, the next-to-first and next-to-middle, etc. |
| 999 | float xprod_sum = 0.0; |
| 1000 | int n_2 = Npoints/2; |
| 1001 | s_Cdc_trackhit_t *v = hits; |
| 1002 | s_Cdc_trackhit_t *v2=&hits[n_2]; |
| 1003 | for(int i=0;i<n_2;i++, v++, v2++){ |
| 1004 | xprod_sum += v->x*v2->y - v2->x*v->y; |
| 1005 | } |
| 1006 | if(xprod_sum>0.0)q = -q; |
| 1007 | |
| 1008 | // Phi is pi/2 out of phase with x0,y0. The sign of the phase difference |
| 1009 | // depends on the charge |
| 1010 | phi = atan2(y0,x0); |
| 1011 | phi += q>0.0 ? -M_PI_21.57079632679489661923:M_PI_21.57079632679489661923; |
| 1012 | |
| 1013 | // Theta is determined by extrapolating the helix back to the target. |
| 1014 | // To do this, we need dphi/dz and a phi,z point. The easiest way to |
| 1015 | // get these is by a simple average (I guess). |
| 1016 | v = hits; |
| 1017 | v2=&v[1]; |
| 1018 | float dzdphi =0.0; |
| 1019 | int Ndzdphipoints = 0; |
| 1020 | for(int i=0;i<Npoints-1;i++, v++, v2++){ |
| 1021 | float myphi1 = atan2(v->y-y0, v->x-x0); |
| 1022 | float myphi2 = atan2(v2->y-y0, v2->x-x0); |
| 1023 | float mydzdphi = (v2->z-v->z)/(myphi2-myphi1); |
| 1024 | if(finite(mydzdphi)){ |
| 1025 | dzdphi+=mydzdphi; |
| 1026 | Ndzdphipoints++; |
| 1027 | } |
| 1028 | } |
| 1029 | if(Ndzdphipoints){ |
| 1030 | dzdphi/=(float)Ndzdphipoints; |
| 1031 | } |
| 1032 | |
| 1033 | theta = atan(r0/fabs(dzdphi)); |
| 1034 | p = -p_trans/sin(theta); |
| 1035 | |
| 1036 | return NOERROR; |
| 1037 | } |
| 1038 | |
| 1039 | #endif |
| 1040 |