| File: | libraries/TRACKING/DHelicalFit.cc |
| Location: | line 311, column 20 |
| Description: | Value stored to 'delta_phi' is never read |
| 1 | // Set of routines for fitting tracks in both the cdc and fdc assuming a |
| 2 | // uniform magnetic field. The track therefore is assumed to follow a helical |
| 3 | // path. |
| 4 | |
| 5 | #include <iostream> |
| 6 | #include <algorithm> |
| 7 | using namespace std; |
| 8 | |
| 9 | #include "DANA/DApplication.h" |
| 10 | #include <TDecompLU.h> |
| 11 | #include <math.h> |
| 12 | |
| 13 | #include "DHelicalFit.h" |
| 14 | #define qBr2p0.003 0.003 // conversion for converting q*B*r to GeV/c |
| 15 | |
| 16 | #define ONE_THIRD0.33333333333333333 0.33333333333333333 |
| 17 | #define SQRT31.73205080756887719 1.73205080756887719 |
| 18 | #define EPS1e-8 1e-8 |
| 19 | #ifndef M_TWO_PI6.28318530717958647692 |
| 20 | #define M_TWO_PI6.28318530717958647692 6.28318530717958647692 |
| 21 | #endif |
| 22 | |
| 23 | // The following is for sorting hits by z |
| 24 | class DHFHitLessThanZ{ |
| 25 | public: |
| 26 | bool operator()(DHFHit_t* const &a, DHFHit_t* const &b) const { |
| 27 | return a->z < b->z; |
| 28 | } |
| 29 | }; |
| 30 | |
| 31 | inline bool DHFHitLessThanZ_C(DHFHit_t* const &a, DHFHit_t* const &b) { |
| 32 | return a->z < b->z; |
| 33 | } |
| 34 | inline bool RiemannFit_hit_cmp(DHFHit_t *a,DHFHit_t *b){ |
| 35 | return (a->z>b->z); |
| 36 | } |
| 37 | |
| 38 | |
| 39 | |
| 40 | //----------------- |
| 41 | // DHelicalFit (Constructor) |
| 42 | //----------------- |
| 43 | DHelicalFit::DHelicalFit(void) |
| 44 | { |
| 45 | x0 = y0 = r0 = tanl = z_vertex = p_trans = phi= theta = q = p = dzdphi =0; |
| 46 | chisq = 0; |
| 47 | ndof=0; |
| 48 | chisq_source = NOFIT; |
| 49 | bfield = NULL__null; |
| 50 | |
| 51 | // For Riemann Fit |
| 52 | CovR_=NULL__null; |
| 53 | CovRPhi_=NULL__null; |
| 54 | c_origin=0.; |
| 55 | normal.SetXYZ(0.,0.,0.); |
| 56 | } |
| 57 | |
| 58 | //----------------- |
| 59 | // DHelicalFit (Constructor) |
| 60 | //----------------- |
| 61 | DHelicalFit::DHelicalFit(const DHelicalFit &fit) |
| 62 | { |
| 63 | Copy(fit); |
| 64 | } |
| 65 | //----------------- |
| 66 | // Copy |
| 67 | //----------------- |
| 68 | void DHelicalFit::Copy(const DHelicalFit &fit) |
| 69 | { |
| 70 | x0 = fit.x0; |
| 71 | y0 = fit.y0; |
| 72 | r0 = fit.r0; |
| 73 | q = fit.q; |
| 74 | p = fit.p; |
| 75 | tanl=fit.tanl; |
| 76 | p_trans = fit.p_trans; |
| 77 | phi = fit.phi; |
| 78 | theta = fit.theta; |
| 79 | z_vertex = fit.z_vertex; |
| 80 | chisq = fit.chisq; |
| 81 | ndof=fit.ndof; |
| 82 | dzdphi = fit.dzdphi; |
| 83 | chisq_source = fit.chisq_source; |
| 84 | bfield = fit.GetMagneticFieldMap(); |
| 85 | Bz_avg = fit.GetBzAvg(); |
| 86 | z_mean = fit.GetZMean(); |
| 87 | phi_mean = fit.GetPhiMean(); |
| 88 | |
| 89 | const vector<DHFHit_t*> myhits = fit.GetHits(); |
| 90 | for(unsigned int i=0; i<myhits.size(); i++){ |
| 91 | DHFHit_t *a = new DHFHit_t; |
| 92 | *a = *myhits[i]; |
| 93 | hits.push_back(a); |
| 94 | } |
| 95 | if (fit.CovR_!=NULL__null){ |
| 96 | CovR_ = new DMatrix(hits.size(),hits.size()); |
| 97 | for (unsigned int i=0;i<hits.size();i++) |
| 98 | for (unsigned int j=0;j<hits.size();j++) |
| 99 | CovR_->operator()(i,j)=fit.CovR_->operator()(i,j); |
| 100 | } |
| 101 | else CovR_=NULL__null; |
| 102 | if (fit.CovRPhi_!=NULL__null){ |
| 103 | CovRPhi_ = new DMatrix(hits.size(),hits.size()); |
| 104 | for (unsigned int i=0;i<hits.size();i++) |
| 105 | for (unsigned int j=0;j<hits.size();j++) |
| 106 | CovRPhi_->operator()(i,j)=fit.CovRPhi_->operator()(i,j); |
| 107 | } |
| 108 | else CovRPhi_=NULL__null; |
| 109 | } |
| 110 | |
| 111 | //----------------- |
| 112 | // operator= (Assignment operator) |
| 113 | //----------------- |
| 114 | DHelicalFit& DHelicalFit::operator=(const DHelicalFit& fit) |
| 115 | { |
| 116 | if(this == &fit)return *this; |
| 117 | Copy(fit); |
| 118 | |
| 119 | return *this; |
| 120 | } |
| 121 | |
| 122 | //----------------- |
| 123 | // DHelicalFit (Destructor) |
| 124 | //----------------- |
| 125 | DHelicalFit::~DHelicalFit() |
| 126 | { |
| 127 | Clear(); |
| 128 | } |
| 129 | //----------------- |
| 130 | // AddHit -- add an FDC hit to the list |
| 131 | //----------------- |
| 132 | jerror_t DHelicalFit::AddHit(const DFDCPseudo *fdchit){ |
| 133 | DHFHit_t *hit = new DHFHit_t; |
| 134 | hit->x = fdchit->xy.X(); |
| 135 | hit->y = fdchit->xy.Y(); |
| 136 | hit->z = fdchit->wire->origin.z(); |
| 137 | hit->covx=fdchit->covxx; |
| 138 | hit->covy=fdchit->covyy; |
| 139 | hit->covxy=fdchit->covxy; |
| 140 | hit->chisq = 0.0; |
| 141 | hits.push_back(hit); |
| 142 | |
| 143 | return NOERROR; |
| 144 | |
| 145 | } |
| 146 | |
| 147 | |
| 148 | //----------------- |
| 149 | // AddHit |
| 150 | //----------------- |
| 151 | jerror_t DHelicalFit::AddHit(float r, float phi, float z) |
| 152 | { |
| 153 | /// Add a hit to the list of hits using cylindrical coordinates |
| 154 | |
| 155 | return AddHitXYZ(r*cos(phi), r*sin(phi), z); |
| 156 | } |
| 157 | |
| 158 | //----------------- |
| 159 | // AddHitXYZ |
| 160 | //----------------- |
| 161 | jerror_t DHelicalFit::AddHitXYZ(float x, float y, float z) |
| 162 | { |
| 163 | /// Add a hit to the list of hits using Cartesian coordinates |
| 164 | DHFHit_t *hit = new DHFHit_t; |
| 165 | hit->x = x; |
| 166 | hit->y = y; |
| 167 | hit->z = z; |
| 168 | hit->covx=1.; |
| 169 | hit->covy=1.; |
| 170 | hit->covxy=0.; |
| 171 | hit->chisq = 0.0; |
| 172 | hits.push_back(hit); |
| 173 | |
| 174 | return NOERROR; |
| 175 | } |
| 176 | |
| 177 | // Add a hit to the list of hits using Cartesian coordinates |
| 178 | jerror_t DHelicalFit::AddHitXYZ(float x,float y, float z,float covx, |
| 179 | float covy, float covxy){ |
| 180 | DHFHit_t *hit = new DHFHit_t; |
| 181 | hit->x = x; |
| 182 | hit->y = y; |
| 183 | hit->z = z; |
| 184 | hit->covx=covx; |
| 185 | hit->covy=covy; |
| 186 | hit->covxy=covxy; |
| 187 | hits.push_back(hit); |
| 188 | |
| 189 | return NOERROR; |
| 190 | } |
| 191 | |
| 192 | |
| 193 | //----------------- |
| 194 | // PruneHit |
| 195 | //----------------- |
| 196 | jerror_t DHelicalFit::PruneHit(int idx) |
| 197 | { |
| 198 | /// Remove the hit specified by idx from the list |
| 199 | /// of hits. The value of idx can be anywhere from |
| 200 | /// 0 to GetNhits()-1. |
| 201 | if(idx<0 || idx>=(int)hits.size())return VALUE_OUT_OF_RANGE; |
| 202 | |
| 203 | delete hits[idx]; |
| 204 | hits.erase(hits.begin() + idx); |
| 205 | |
| 206 | return NOERROR; |
| 207 | } |
| 208 | |
| 209 | //----------------- |
| 210 | // Clear |
| 211 | //----------------- |
| 212 | jerror_t DHelicalFit::Clear(void) |
| 213 | { |
| 214 | /// Remove covariance matrices |
| 215 | if (CovR_!=NULL__null) delete CovR_; |
| 216 | if (CovRPhi_!=NULL__null) delete CovRPhi_; |
| 217 | |
| 218 | /// Remove all hits |
| 219 | for(unsigned int i=0; i<hits.size(); i++)delete hits[i]; |
| 220 | hits.clear(); |
| 221 | |
| 222 | return NOERROR; |
| 223 | } |
| 224 | |
| 225 | //----------------- |
| 226 | // PrintChiSqVector |
| 227 | //----------------- |
| 228 | jerror_t DHelicalFit::PrintChiSqVector(void) const |
| 229 | { |
| 230 | /// Dump the latest chi-squared vector to the screen. |
| 231 | /// This prints the individual hits' chi-squared |
| 232 | /// contributions in the order in which the hits were |
| 233 | /// added. See PruneHits() for more detail. |
| 234 | |
| 235 | cout<<"Chisq vector from DHelicalFit: (source="; |
| 236 | switch(chisq_source){ |
| 237 | case NOFIT: cout<<"NOFIT"; break; |
| 238 | case CIRCLE: cout<<"CIRCLE"; break; |
| 239 | case TRACK: cout<<"TRACK"; break; |
| 240 | }; |
| 241 | cout<<")"<<endl; |
| 242 | cout<<"----------------------------"<<endl; |
| 243 | |
| 244 | for(unsigned int i=0;i<hits.size();i++){ |
| 245 | cout<<i<<" "<<hits[i]->chisq<<endl; |
| 246 | } |
| 247 | cout<<"Total: "<<chisq<<endl<<endl; |
| 248 | |
| 249 | return NOERROR; |
| 250 | } |
| 251 | |
| 252 | //----------------- |
| 253 | // FitCircle |
| 254 | //----------------- |
| 255 | jerror_t DHelicalFit::FitCircle(void) |
| 256 | { |
| 257 | /// Fit the current set of hits to a circle |
| 258 | /// |
| 259 | /// This takes the hits which have been added thus far via one |
| 260 | /// of the AddHit() methods and fits them to a circle. |
| 261 | /// The alogorithm used here calculates the parameters directly using |
| 262 | /// a technique very much like linear regression. The key assumptions |
| 263 | /// are: |
| 264 | /// 1. The magnetic field is uniform and along z so that the projection |
| 265 | /// of the track onto the X/Y plane will fall on a circle |
| 266 | /// (this also implies no multiple-scattering or energy loss) |
| 267 | /// 2. The vertex is at the target center (i.e. 0,0 in the coordinate |
| 268 | /// system of the points passed to us. |
| 269 | /// |
| 270 | /// IMPORTANT: The value of phi which results from this assumes |
| 271 | /// the particle was positively charged. If the particle was |
| 272 | /// negatively charged, then phi will be 180 degrees off. To |
| 273 | /// correct this, one needs z-coordinate information to determine |
| 274 | /// the sign of the charge. |
| 275 | /// |
| 276 | /// ALSO IMPORTANT: This assumes a charge of +1 for the particle. If |
| 277 | /// the particle actually has a charge of +2, then the resulting |
| 278 | /// value of p_trans will be half of what it should be. |
| 279 | |
| 280 | float alpha=0.0, beta=0.0, gamma=0.0, deltax=0.0, deltay=0.0; |
| 281 | chisq_source = NOFIT; // in case we return early |
| 282 | ndof=hits.size()-2; |
| 283 | |
| 284 | // Loop over hits to calculate alpha, beta, gamma, and delta |
| 285 | // if a magnetic field map was given, use it to find average Z B-field |
| 286 | DHFHit_t *a = NULL__null; |
| 287 | for(unsigned int i=0;i<hits.size();i++){ |
| 288 | a = hits[i]; |
| 289 | float x=a->x; |
| 290 | float y=a->y; |
| 291 | alpha += x*x; |
| 292 | beta += y*y; |
| 293 | gamma += x*y; |
| 294 | deltax += x*(x*x+y*y)/2.0; |
| 295 | deltay += y*(x*x+y*y)/2.0; |
| 296 | } |
| 297 | |
| 298 | // Calculate x0,y0 - the center of the circle |
| 299 | double denom = alpha*beta-gamma*gamma; |
| 300 | if(fabs(denom)<1.0E-20)return UNRECOVERABLE_ERROR; |
| 301 | x0 = (deltax*beta-deltay*gamma)/denom; |
| 302 | //y0 = (deltay-gamma*x0)/beta; |
| 303 | y0 = (deltay*alpha-deltax*gamma)/denom; |
| 304 | |
| 305 | // Calculate the phi angle traversed by the track from the |
| 306 | // origin to the last hit. NOTE: this can be off by a multiple |
| 307 | // of 2pi! |
| 308 | double delta_phi=0.0; |
| 309 | if(a){ // a should be pointer to last hit from above loop |
| 310 | delta_phi = atan2(a->y-y0, a->x-x0); |
| 311 | if(delta_phi<0.0)delta_phi += M_TWO_PI6.28318530717958647692; |
Value stored to 'delta_phi' is never read | |
| 312 | } |
| 313 | |
| 314 | // Momentum depends on magnetic field. If bfield has been |
| 315 | // set, we should use it to determine an average value of Bz |
| 316 | // for this track. Otherwise, assume -2T. |
| 317 | // Also assume a singly charged track (i.e. q=+/-1) |
| 318 | // The sign of the charge will be determined below. |
| 319 | Bz_avg=-2.0; |
| 320 | q = +1.0; |
| 321 | r0 = sqrt(x0*x0 + y0*y0); |
| 322 | p_trans = q*Bz_avg*r0*qBr2p0.003; // qBr2p converts to GeV/c |
| 323 | phi = atan2(y0,x0) - M_PI_21.57079632679489661923; |
| 324 | if(p_trans<0.0){ |
| 325 | p_trans = -p_trans; |
| 326 | } |
| 327 | if(phi<0)phi+=M_TWO_PI6.28318530717958647692; |
| 328 | if(phi>=M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
| 329 | |
| 330 | // Calculate the chisq |
| 331 | ChisqCircle(); |
| 332 | chisq_source = CIRCLE; |
| 333 | |
| 334 | return NOERROR; |
| 335 | } |
| 336 | |
| 337 | //----------------- |
| 338 | // ChisqCircle |
| 339 | //----------------- |
| 340 | double DHelicalFit::ChisqCircle(void) |
| 341 | { |
| 342 | /// Calculate the chisq for the hits and current circle |
| 343 | /// parameters. |
| 344 | /// NOTE: This does not return the chi2/dof, just the |
| 345 | /// raw chi2 with errs set to 1.0 |
| 346 | chisq = 0.0; |
| 347 | for(unsigned int i=0;i<hits.size();i++){ |
| 348 | DHFHit_t *a = hits[i]; |
| 349 | float x = a->x - x0; |
| 350 | float y = a->y - y0; |
| 351 | float c = sqrt(x*x + y*y) - r0; |
| 352 | c *= c; |
| 353 | a->chisq = c; |
| 354 | chisq+=c; |
| 355 | } |
| 356 | |
| 357 | // Do NOT divide by DOF |
| 358 | return chisq; |
| 359 | } |
| 360 | |
| 361 | //----------------- |
| 362 | // FitCircleRiemann |
| 363 | //----------------- |
| 364 | jerror_t DHelicalFit::FitCircleRiemann(float z_vertex,float BeamRMS) |
| 365 | { |
| 366 | chisq_source = NOFIT; // in case we return early |
| 367 | |
| 368 | // Fake point at origin |
| 369 | float beam_var=BeamRMS*BeamRMS; |
| 370 | AddHitXYZ(0.,0.,z_vertex,beam_var,beam_var,0.); |
| 371 | |
| 372 | jerror_t err = FitCircleRiemann(); |
| 373 | if(err!=NOERROR)return err; |
| 374 | |
| 375 | // Number of degrees of freedom |
| 376 | ndof=hits.size()-3; |
| 377 | |
| 378 | // Momentum depends on magnetic field. If bfield has been |
| 379 | // set, we should use it to determine an average value of Bz |
| 380 | // for this track. Otherwise, assume -2T. |
| 381 | // Also assume a singly charged track (i.e. q=+/-1) |
| 382 | // The sign of the charge will be determined below. |
| 383 | Bz_avg=-2.0; |
| 384 | q = +1.0; |
| 385 | p_trans = q*Bz_avg*r0*qBr2p0.003; // qBr2p converts to GeV/c |
| 386 | if(p_trans<0.0){ |
| 387 | p_trans = -p_trans; |
| 388 | } |
| 389 | phi=atan2(-x0,y0); |
| 390 | if(phi<0)phi+=M_TWO_PI6.28318530717958647692; |
| 391 | if(phi>=M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
| 392 | |
| 393 | // Normal vector for plane intersecting Riemann surface |
| 394 | normal.SetXYZ(N[0],N[1],N[2]); |
| 395 | |
| 396 | return NOERROR; |
| 397 | } |
| 398 | |
| 399 | //----------------- |
| 400 | // FitCircleRiemann |
| 401 | //----------------- |
| 402 | jerror_t DHelicalFit::FitCircleRiemann(void){ |
| 403 | /// Riemann Circle fit: points on a circle in x,y project onto a plane cutting |
| 404 | /// the circular paraboloid surface described by (x,y,x^2+y^2). Therefore the |
| 405 | /// task of fitting points in (x,y) to a circle is transormed to the task of |
| 406 | /// fitting planes in (x,y, w=x^2+y^2) space |
| 407 | /// |
| 408 | DMatrix X(hits.size(),3); |
| 409 | DMatrix Xavg(1,3); |
| 410 | DMatrix A(3,3); |
| 411 | // vector of ones |
| 412 | DMatrix OnesT(1,hits.size()); |
| 413 | double W_sum=0.; |
| 414 | DMatrix W(hits.size(),hits.size()); |
| 415 | |
| 416 | // Make sure hit list is ordered in z |
| 417 | std::sort(hits.begin(),hits.end(),RiemannFit_hit_cmp); |
| 418 | |
| 419 | // Covariance matrix |
| 420 | DMatrix CRPhi(hits.size(),hits.size()); |
| 421 | if (CovRPhi_==NULL__null){ |
| 422 | CovRPhi_ = new DMatrix(hits.size(),hits.size()); |
| 423 | for (unsigned int i=0;i<hits.size();i++){ |
| 424 | double Phi=atan2(hits[i]->y,hits[i]->x); |
| 425 | double cosPhi=cos(Phi); |
| 426 | double sinPhi=sin(Phi); |
| 427 | double Phi_cosPhi_minus_sinPhi=Phi*cosPhi-sinPhi; |
| 428 | double Phi_sinPhi_plus_cosPhi=Phi*sinPhi+cosPhi; |
| 429 | CovRPhi_->operator()(i,i) |
| 430 | =Phi_cosPhi_minus_sinPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covx |
| 431 | +Phi_sinPhi_plus_cosPhi*Phi_sinPhi_plus_cosPhi*hits[i]->covy |
| 432 | +2.*Phi_sinPhi_plus_cosPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covxy; |
| 433 | } |
| 434 | } |
| 435 | for (unsigned int i=0;i<hits.size();i++) |
| 436 | for (unsigned int j=0;j<hits.size();j++) |
| 437 | CRPhi(i,j)=CovRPhi_->operator()(i, j); |
| 438 | |
| 439 | // The goal is to find the eigenvector corresponding to the smallest |
| 440 | // eigenvalue of the equation |
| 441 | // lambda=n^T (X^T W X - W_sum Xavg^T Xavg)n |
| 442 | // where n is the normal vector to the plane slicing the cylindrical |
| 443 | // paraboloid described by the parameterization (x,y,w=x^2+y^2), |
| 444 | // and W is the weight matrix, assumed for now to be diagonal. |
| 445 | // In the absence of multiple scattering, W_sum is the sum of all the |
| 446 | // diagonal elements in W. |
| 447 | |
| 448 | for (unsigned int i=0;i<hits.size();i++){ |
| 449 | X(i,0)=hits[i]->x; |
| 450 | X(i,1)=hits[i]->y; |
| 451 | X(i,2)=hits[i]->x*hits[i]->x+hits[i]->y*hits[i]->y; |
| 452 | OnesT(0,i)=1.; |
| 453 | W(i,i)=1./CRPhi(i,i); |
| 454 | W_sum+=W(i,i); |
| 455 | } |
| 456 | Xavg=(1./W_sum)*(OnesT*(W*X)); |
| 457 | |
| 458 | A=DMatrix(DMatrix::kTransposed,X)*(W*X) |
| 459 | -W_sum*(DMatrix(DMatrix::kTransposed,Xavg)*Xavg); |
| 460 | if(!A.IsValid())return UNRECOVERABLE_ERROR; |
| 461 | |
| 462 | // The characteristic equation is |
| 463 | // lambda^3+B2*lambda^2+lambda*B1+B0=0 |
| 464 | // |
| 465 | double B2=-(A(0,0)+A(1,1)+A(2,2)); |
| 466 | 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) |
| 467 | +A(1,1)*A(2,2)-A(2,1)*A(1,2); |
| 468 | double B0=-A.Determinant(); |
| 469 | if(B0==0 || !finite(B0))return UNRECOVERABLE_ERROR; |
| 470 | |
| 471 | // The roots of the cubic equation are given by |
| 472 | // lambda1= -B2/3 + S+T |
| 473 | // lambda2= -B2/3 - (S+T)/2 + i sqrt(3)/2. (S-T) |
| 474 | // lambda3= -B2/3 - (S+T)/2 - i sqrt(3)/2. (S-T) |
| 475 | // where we define some temporary variables: |
| 476 | // S= (R+sqrt(Q^3+R^2))^(1/3) |
| 477 | // T= (R-sqrt(Q^3+R^2))^(1/3) |
| 478 | // Q=(3*B1-B2^2)/9 |
| 479 | // R=(9*B2*B1-27*B0-2*B2^3)/54 |
| 480 | // sum=S+T; |
| 481 | // diff=i*(S-T) |
| 482 | // We divide Q and R by a safety factor to prevent multiplying together |
| 483 | // enormous numbers that cause unreliable results. |
| 484 | |
| 485 | long double Q=(3.*B1-B2*B2)/9.e4; |
| 486 | long double R=(9.*B2*B1-27.*B0-2.*B2*B2*B2)/54.e6; |
| 487 | long double Q1=Q*Q*Q+R*R; |
| 488 | if (Q1<0) Q1=sqrt(-Q1); |
| 489 | else{ |
| 490 | return VALUE_OUT_OF_RANGE; |
| 491 | } |
| 492 | |
| 493 | // DeMoivre's theorem for fractional powers of complex numbers: |
| 494 | // (r*(cos(theta)+i sin(theta)))^(p/q) |
| 495 | // = r^(p/q)*(cos(p*theta/q)+i sin(p*theta/q)) |
| 496 | // |
| 497 | //double temp=100.*pow(R*R+Q1*Q1,0.16666666666666666667); |
| 498 | long double temp=100.*sqrt(cbrt(R*R+Q1*Q1)); |
| 499 | long double theta1=ONE_THIRD0.33333333333333333*atan2(Q1,R); |
| 500 | long double sum_over_2=temp*cos(theta1); |
| 501 | long double diff_over_2=-temp*sin(theta1); |
| 502 | // Third root |
| 503 | long double lambda_min=-ONE_THIRD0.33333333333333333*B2-sum_over_2+SQRT31.73205080756887719*diff_over_2; |
| 504 | |
| 505 | // Calculate the (normal) eigenvector corresponding to the eigenvalue lambda |
| 506 | N[0]=1.; |
| 507 | N[1]=(A(1,0)*A(0,2)-(A(0,0)-lambda_min)*A(1,2)) |
| 508 | /(A(0,1)*A(2,1)-(A(1,1)-lambda_min)*A(0,2)); |
| 509 | N[2]=(A(2,0)*(A(1,1)-lambda_min)-A(1,0)*A(2,1)) |
| 510 | /(A(1,2)*A(2,1)-(A(2,2)-lambda_min)*(A(1,1)-lambda_min)); |
| 511 | |
| 512 | // Normalize: n1^2+n2^2+n3^2=1 |
| 513 | double denom=sqrt(N[0]*N[0]+N[1]*N[1]+N[2]*N[2]); |
| 514 | for (int i=0;i<3;i++){ |
| 515 | N[i]/=denom; |
| 516 | } |
| 517 | |
| 518 | // Distance to origin |
| 519 | c_origin=-(N[0]*Xavg(0,0)+N[1]*Xavg(0,1)+N[2]*Xavg(0,2)); |
| 520 | |
| 521 | // Center and radius of the circle |
| 522 | double one_over_2Nz=1./(2.*N[2]); |
| 523 | //x0=-N[0]/2./N[2]; |
| 524 | //y0=-N[1]/2./N[2]; |
| 525 | //r0=sqrt(1.-N[2]*N[2]-4.*c_origin*N[2])/2./fabs(N[2]); |
| 526 | x0=-N[0]*one_over_2Nz; |
| 527 | y0=-N[1]*one_over_2Nz; |
| 528 | r0=sqrt(1.-N[2]*N[2]-4.*c_origin*N[2])*fabs(one_over_2Nz); |
| 529 | |
| 530 | // Phi value at "vertex" |
| 531 | phi=atan2(-x0,y0); |
| 532 | if(phi<0)phi+=M_TWO_PI6.28318530717958647692; |
| 533 | if(phi>=M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
| 534 | |
| 535 | // Calculate the chisq |
| 536 | ChisqCircle(); |
| 537 | chisq_source = CIRCLE; |
| 538 | |
| 539 | return NOERROR; |
| 540 | } |
| 541 | |
| 542 | |
| 543 | //----------------- |
| 544 | // FitCircleRiemannCorrected |
| 545 | //----------------- |
| 546 | // Riemann Circle fit with correction for non-normal track incidence |
| 547 | // |
| 548 | jerror_t DHelicalFit::FitCircleRiemannCorrected(float rc){ |
| 549 | // Covariance matrices |
| 550 | DMatrix CRPhi(hits.size(),hits.size()); |
| 551 | DMatrix CR(hits.size(),hits.size()); |
| 552 | // Auxiliary matrices for correcting for non-normal track incidence to FDC |
| 553 | // The correction is |
| 554 | // CRPhi'= C*CRPhi*C+S*CR*S, where S(i,i)=R_i*kappa/2 |
| 555 | // and C(i,i)=sqrt(1-S(i,i)^2) |
| 556 | DMatrix C(hits.size(),hits.size()); |
| 557 | DMatrix S(hits.size(),hits.size()); |
| 558 | for (unsigned int i=0;i<hits.size();i++){ |
| 559 | S(i,i)=0.; |
| 560 | C(i,i)=1.; |
| 561 | if (rc>0){ |
| 562 | double rtemp=sqrt(hits[i]->x*hits[i]->x+hits[i]->y*hits[i]->y); |
| 563 | double stemp=rtemp/4./rc; |
| 564 | double ctemp=1.-stemp*stemp; |
| 565 | if (ctemp>0){ |
| 566 | S(i,i)=stemp; |
| 567 | C(i,i)=sqrt(ctemp); |
| 568 | } |
| 569 | } |
| 570 | } |
| 571 | |
| 572 | // Covariance matrices |
| 573 | if (CovRPhi_==NULL__null){ |
| 574 | CovRPhi_ = new DMatrix(hits.size(),hits.size()); |
| 575 | for (unsigned int i=0;i<hits.size();i++){ |
| 576 | double Phi=atan2(hits[i]->y,hits[i]->x); |
| 577 | double cosPhi=cos(Phi); |
| 578 | double sinPhi=sin(Phi); |
| 579 | double Phi_sinPhi_plus_cosPhi=Phi*sinPhi+cosPhi; |
| 580 | double Phi_cosPhi_minus_sinPhi=Phi*cosPhi-sinPhi; |
| 581 | CovRPhi_->operator()(i,i) |
| 582 | =Phi_cosPhi_minus_sinPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covx |
| 583 | +Phi_sinPhi_plus_cosPhi*Phi_sinPhi_plus_cosPhi*hits[i]->covy |
| 584 | +2.*Phi_sinPhi_plus_cosPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covxy; |
| 585 | } |
| 586 | } |
| 587 | if (CovR_==NULL__null){ |
| 588 | CovR_= new DMatrix(hits.size(),hits.size()); |
| 589 | for (unsigned int m=0;m<hits.size();m++){ |
| 590 | double Phi=atan2(hits[m]->y,hits[m]->x); |
| 591 | double cosPhi=cos(Phi); |
| 592 | double sinPhi=sin(Phi); |
| 593 | CovR_->operator()(m,m)=cosPhi*cosPhi*hits[m]->covx |
| 594 | +sinPhi*sinPhi*hits[m]->covy |
| 595 | +2.*sinPhi*cosPhi*hits[m]->covxy; |
| 596 | } |
| 597 | } |
| 598 | for (unsigned int i=0;i<hits.size();i++){ |
| 599 | for (unsigned int j=0;j<hits.size();j++){ |
| 600 | CR(i,j)=CovR_->operator()(i, j); |
| 601 | CRPhi(i,j)=CovRPhi_->operator()(i, j); |
| 602 | } |
| 603 | } |
| 604 | |
| 605 | // Correction for non-normal incidence of track on FDC |
| 606 | CRPhi=C*CRPhi*C+S*CR*S; |
| 607 | for (unsigned int i=0;i<hits.size();i++) |
| 608 | for (unsigned int j=0;j<hits.size();j++) |
| 609 | CovRPhi_->operator()(i,j)=CRPhi(i,j); |
| 610 | return FitCircleRiemann(); |
| 611 | } |
| 612 | |
| 613 | |
| 614 | //----------------- |
| 615 | // GetChargeRiemann |
| 616 | //----------------- |
| 617 | // Charge-finding routine with corrected CRPhi (see above) |
| 618 | // |
| 619 | jerror_t DHelicalFit::GetChargeRiemann(float rc_input){ |
| 620 | // Covariance matrices |
| 621 | DMatrix CRPhi(hits.size(),hits.size()); |
| 622 | DMatrix CR(hits.size(),hits.size()); |
| 623 | // Auxiliary matrices for correcting for non-normal track incidence to FDC |
| 624 | // The correction is |
| 625 | // CRPhi'= C*CRPhi*C+S*CR*S, where S(i,i)=R_i*kappa/2 |
| 626 | // and C(i,i)=sqrt(1-S(i,i)^2) |
| 627 | DMatrix C(hits.size(),hits.size()); |
| 628 | DMatrix S(hits.size(),hits.size()); |
| 629 | for (unsigned int i=0;i<hits.size();i++){ |
| 630 | double rtemp=sqrt(hits[i]->x*hits[i]->x+hits[i]->y*hits[i]->y); |
| 631 | double stemp=rtemp/(4.*rc_input); |
| 632 | double ctemp=1.-stemp*stemp; |
| 633 | if (ctemp>0){ |
| 634 | S(i,i)=stemp; |
| 635 | C(i,i)=sqrt(ctemp); |
| 636 | } |
| 637 | else{ |
| 638 | S(i,i)=0.; |
| 639 | C(i,i)=1; |
| 640 | } |
| 641 | } |
| 642 | |
| 643 | // Covariance matrices |
| 644 | if (CovRPhi_==NULL__null){ |
| 645 | CovRPhi_ = new DMatrix(hits.size(),hits.size()); |
| 646 | for (unsigned int i=0;i<hits.size();i++){ |
| 647 | double Phi=atan2(hits[i]->y,hits[i]->x); |
| 648 | double cosPhi=cos(Phi); |
| 649 | double sinPhi=sin(Phi); |
| 650 | double Phi_sinPhi_plus_cosPhi=Phi*sinPhi+cosPhi; |
| 651 | double Phi_cosPhi_minus_sinPhi=Phi*cosPhi-sinPhi; |
| 652 | CovRPhi_->operator()(i,i) |
| 653 | =Phi_cosPhi_minus_sinPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covx |
| 654 | +Phi_sinPhi_plus_cosPhi*Phi_sinPhi_plus_cosPhi*hits[i]->covy |
| 655 | +2.*Phi_sinPhi_plus_cosPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covxy; |
| 656 | } |
| 657 | } |
| 658 | if (CovR_==NULL__null){ |
| 659 | CovR_= new DMatrix(hits.size(),hits.size()); |
| 660 | for (unsigned int m=0;m<hits.size();m++){ |
| 661 | double Phi=atan2(hits[m]->y,hits[m]->x); |
| 662 | double cosPhi=cos(Phi); |
| 663 | double sinPhi=sin(Phi); |
| 664 | CovR_->operator()(m,m)=cosPhi*cosPhi*hits[m]->covx |
| 665 | +sinPhi*sinPhi*hits[m]->covy |
| 666 | +2.*sinPhi*cosPhi*hits[m]->covxy; |
| 667 | } |
| 668 | } |
| 669 | for (unsigned int i=0;i<hits.size();i++){ |
| 670 | for (unsigned int j=0;j<hits.size();j++){ |
| 671 | CR(i,j)=CovR_->operator()(i, j); |
| 672 | CRPhi(i,j)=CovRPhi_->operator()(i, j); |
| 673 | } |
| 674 | } |
| 675 | // Correction for non-normal incidence of track on FDC |
| 676 | CRPhi=C*CRPhi*C+S*CR*S; |
| 677 | for (unsigned int i=0;i<hits.size();i++) |
| 678 | for (unsigned int j=0;j<hits.size();j++) |
| 679 | CovRPhi_->operator()(i,j)=CRPhi(i,j); |
| 680 | return GetChargeRiemann(); |
| 681 | } |
| 682 | |
| 683 | //----------------- |
| 684 | // GetChargeRiemann |
| 685 | //----------------- |
| 686 | // Linear regression to find charge |
| 687 | // |
| 688 | jerror_t DHelicalFit::GetChargeRiemann(){ |
| 689 | // Covariance matrices |
| 690 | DMatrix CRPhi(hits.size(),hits.size()); |
| 691 | DMatrix CR(hits.size(),hits.size()); |
| 692 | if (CovRPhi_==NULL__null){ |
| 693 | CovRPhi_ = new DMatrix(hits.size(),hits.size()); |
| 694 | for (unsigned int i=0;i<hits.size();i++){ |
| 695 | double Phi=atan2(hits[i]->y,hits[i]->x); |
| 696 | double cosPhi=cos(Phi); |
| 697 | double sinPhi=sin(Phi); |
| 698 | double Phi_cosPhi_minus_sinPhi=Phi*cosPhi-sinPhi; |
| 699 | double Phi_sinPhi_plus_cosPhi=Phi*sinPhi+cosPhi; |
| 700 | CovRPhi_->operator()(i,i) |
| 701 | =Phi_cosPhi_minus_sinPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covx |
| 702 | +Phi_sinPhi_plus_cosPhi*Phi_sinPhi_plus_cosPhi*hits[i]->covy |
| 703 | +2.*Phi_sinPhi_plus_cosPhi*Phi_cosPhi_minus_sinPhi*hits[i]->covxy; |
| 704 | } |
| 705 | } |
| 706 | if (CovR_==NULL__null){ |
| 707 | CovR_= new DMatrix(hits.size(),hits.size()); |
| 708 | for (unsigned int m=0;m<hits.size();m++){ |
| 709 | double Phi=atan2(hits[m]->y,hits[m]->x); |
| 710 | double cosPhi=cos(Phi); |
| 711 | double sinPhi=sin(Phi); |
| 712 | CovR_->operator()(m,m)=cosPhi*cosPhi*hits[m]->covx |
| 713 | +sinPhi*sinPhi*hits[m]->covy |
| 714 | +2.*sinPhi*cosPhi*hits[m]->covxy; |
| 715 | } |
| 716 | } |
| 717 | for (unsigned int i=0;i<hits.size();i++){ |
| 718 | for (unsigned int j=0;j<hits.size();j++){ |
| 719 | CR(i,j)=CovR_->operator()(i, j); |
| 720 | CRPhi(i,j)=CovRPhi_->operator()(i, j); |
| 721 | } |
| 722 | } |
| 723 | |
| 724 | double phi_old=atan2(hits[0]->y,hits[0]->x); |
| 725 | double sumv=0,sumy=0,sumx=0,sumxx=0,sumxy=0; |
| 726 | for (unsigned int k=0;k<hits.size();k++){ |
| 727 | DHFHit_t *hit=hits[k]; |
| 728 | double phi_z=atan2(hit->y,hit->x); |
| 729 | |
| 730 | // Check for problem regions near +pi and -pi |
| 731 | if (fabs(phi_z-phi_old)>M_PI3.14159265358979323846){ |
| 732 | if (phi_old<0) phi_z-=2.*M_PI3.14159265358979323846; |
| 733 | else phi_z+=2.*M_PI3.14159265358979323846; |
| 734 | } |
| 735 | double inv_var=(hit->x*hit->x+hit->y*hit->y) |
| 736 | /(CRPhi(k,k)+phi_z*phi_z*CR(k,k)); |
| 737 | sumv+=inv_var; |
| 738 | sumy+=phi_z*inv_var; |
| 739 | sumx+=hit->z*inv_var; |
| 740 | sumxx+=hit->z*hit->z*inv_var; |
| 741 | sumxy+=phi_z*hit->z*inv_var; |
| 742 | phi_old=phi_z; |
| 743 | } |
| 744 | double slope=(sumv*sumxy-sumy*sumx)/(sumv*sumxx-sumx*sumx); |
| 745 | |
| 746 | // Guess particle charge (+/-1); |
| 747 | q=+1.; |
| 748 | if (slope<0.) q= -1.; |
| 749 | |
| 750 | return NOERROR; |
| 751 | } |
| 752 | |
| 753 | //----------------- |
| 754 | // FitLineRiemann |
| 755 | //----------------- |
| 756 | jerror_t DHelicalFit::FitLineRiemann(){ |
| 757 | /// Riemann Line fit: linear regression of s on z to determine the tangent of |
| 758 | /// the dip angle and the z position of the closest approach to the beam line. |
| 759 | /// Computes intersection points along the helical path. |
| 760 | /// |
| 761 | // Get covariance matrix |
| 762 | DMatrix CR(hits.size(),hits.size()); |
| 763 | if (CovR_==NULL__null){ |
| 764 | CovR_= new DMatrix(hits.size(),hits.size()); |
| 765 | for (unsigned int m=0;m<hits.size();m++){ |
| 766 | double Phi=atan2(hits[m]->y,hits[m]->x); |
| 767 | double cosPhi=cos(Phi); |
| 768 | double sinPhi=sin(Phi); |
| 769 | CovR_->operator()(m,m)=cosPhi*cosPhi*hits[m]->covx |
| 770 | +sinPhi*sinPhi*hits[m]->covy |
| 771 | +2.*sinPhi*cosPhi*hits[m]->covxy; |
| 772 | } |
| 773 | } |
| 774 | for (unsigned int i=0;i<hits.size();i++) |
| 775 | for (unsigned int j=0;j<hits.size();j++) |
| 776 | CR(i,j)=CovR_->operator()(i, j); |
| 777 | |
| 778 | // Fill vector of intersection points |
| 779 | double denom= N[0]*N[0]+N[1]*N[1]; |
| 780 | vector<int>bad(hits.size()); |
| 781 | int numbad=0; |
| 782 | // projection vector |
| 783 | vector<DVector2>projections; |
| 784 | for (unsigned int m=0;m<hits.size();m++){ |
| 785 | double r2=hits[m]->x*hits[m]->x+hits[m]->y*hits[m]->y; |
| 786 | double numer=c_origin+r2*N[2]; |
| 787 | |
| 788 | if (r2==0){ |
| 789 | projections.push_back(DVector2(0.,0.)); |
| 790 | } |
| 791 | else{ |
| 792 | double ratio=numer/denom; |
| 793 | double x_int0=-N[0]*ratio; |
| 794 | double y_int0=-N[1]*ratio; |
| 795 | double temp=denom*r2-numer*numer; |
| 796 | if (temp<0){ // Skip point if the intersection gives nonsense |
| 797 | bad[m]=1; |
| 798 | numbad++; |
| 799 | projections.push_back(DVector2(x_int0,y_int0)); |
| 800 | continue; |
| 801 | } |
| 802 | temp=sqrt(temp)/denom; |
| 803 | |
| 804 | // Choose sign of square root based on proximity to actual measurements |
| 805 | double deltax=N[1]*temp; |
| 806 | double deltay=-N[0]*temp; |
| 807 | double x1=x_int0+deltax; |
| 808 | double y1=y_int0+deltay; |
| 809 | double x2=x_int0-deltax; |
| 810 | double y2=y_int0-deltay; |
| 811 | double diffx1=x1-hits[m]->x; |
| 812 | double diffy1=y1-hits[m]->y; |
| 813 | double diffx2=x2-hits[m]->x; |
| 814 | double diffy2=y2-hits[m]->y; |
| 815 | if (diffx1*diffx1+diffy1*diffy1 > diffx2*diffx2+diffy2*diffy2){ |
| 816 | //temphit->x=x2; |
| 817 | //temphit->y=y2; |
| 818 | projections.push_back(DVector2(x2,y2)); |
| 819 | } |
| 820 | else{ |
| 821 | //temphit->x=x1; |
| 822 | //temphit->y=y1; |
| 823 | projections.push_back(DVector2(x1,y1)); |
| 824 | } |
| 825 | } |
| 826 | //projections.push_back(temphit); |
| 827 | } |
| 828 | |
| 829 | // All arc lengths are measured relative to some reference plane with a hit. |
| 830 | // Don't use a "bad" hit for the reference... |
| 831 | unsigned int start=0; |
| 832 | for (unsigned int i=0;i<bad.size();i++){ |
| 833 | if (!bad[i]){ |
| 834 | start=i; |
| 835 | break; |
| 836 | } |
| 837 | } |
| 838 | |
| 839 | // Linear regression to find z0, tanl |
| 840 | unsigned int n=projections.size(); |
| 841 | double sumv=0.,sumx=0.,sumy=0.,sumxx=0.,sumxy=0.; |
| 842 | double sperp=0.,sperp_old=0., ratio=0, Delta; |
| 843 | double z_last=0.,z=0.; |
| 844 | DVector2 old_proj=projections[start]; |
| 845 | double two_r0=2.*r0; |
| 846 | for (unsigned int k=start;k<n;k++){ |
| 847 | if (!bad[k]){ |
| 848 | sperp_old=sperp; |
| 849 | z_last=z; |
| 850 | ratio=(projections[k]-old_proj).Mod()/(two_r0); |
| 851 | // Make sure the argument for the arcsin does not go out of range... |
| 852 | sperp=sperp_old+(ratio>1? two_r0*M_PI_21.57079632679489661923 : two_r0*asin(ratio)); |
| 853 | z=hits[k]->z; |
| 854 | |
| 855 | // Assume errors in s dominated by errors in R |
| 856 | double weight=1./CR(k,k); |
| 857 | sumv+=weight; |
| 858 | sumy+=sperp*weight; |
| 859 | sumx+=z*weight; |
| 860 | sumxx+=z*z*weight; |
| 861 | sumxy+=sperp*z*weight; |
| 862 | |
| 863 | // Store the current x and y projection values |
| 864 | old_proj=projections[k]; |
| 865 | } |
| 866 | } |
| 867 | Delta=sumv*sumxx-sumx*sumx; |
| 868 | double tanl_denom=sumv*sumxy-sumy*sumx; |
| 869 | if (fabs(Delta)<EPS1e-8 || fabs(tanl_denom)<EPS1e-8) return VALUE_OUT_OF_RANGE; |
| 870 | |
| 871 | // Track parameters tan(lambda) and z-vertex |
| 872 | tanl=-Delta/tanl_denom; |
| 873 | //z_vertex=(sumxx*sumy-sumx*sumxy)/Delta; |
| 874 | sperp-=sperp_old; |
| 875 | z_vertex=z_last-sperp*tanl; |
| 876 | |
| 877 | /* |
| 878 | if (z_vertex<Z_MIN){ |
| 879 | z_vertex=Z_MIN; |
| 880 | tanl=(z_last-Z_MIN)/sperp; |
| 881 | } |
| 882 | else if (z_vertex>Z_MAX){ |
| 883 | z_vertex=Z_MAX; |
| 884 | tanl=(z_last-Z_MAX)/sperp; |
| 885 | } |
| 886 | */ |
| 887 | theta=M_PI_21.57079632679489661923-atan(tanl); |
| 888 | |
| 889 | return NOERROR; |
| 890 | } |
| 891 | |
| 892 | //----------------- |
| 893 | // FitCircleStraightTrack |
| 894 | //----------------- |
| 895 | jerror_t DHelicalFit::FitCircleStraightTrack(void) |
| 896 | { |
| 897 | /// This is a last resort! |
| 898 | /// The circle fit can fail for high momentum tracks that are nearly |
| 899 | /// straight tracks. In these cases (when pt~2GeV or more) there |
| 900 | /// is not enough position resolution with wire positions only |
| 901 | /// to measure the curvature of the track. |
| 902 | /// We can though, fit the X/Y points with a straight line in order |
| 903 | /// to get phi and project out the stereo angles. |
| 904 | /// |
| 905 | /// For the radius of the circle (i.e. p_trans) we loop over momenta |
| 906 | /// from 1 to 9GeV and just use the one with the smallest chisq. |
| 907 | |
| 908 | // Average the x's and y's of the individual hits. We should be |
| 909 | // able to do this via linear regression, but I can't get it to |
| 910 | // work right now and I'm under the gun to get ready for a review |
| 911 | // so I take the easy way. Note that we don't average phi's since |
| 912 | // that will cause problems at the 0/2pi boundary. |
| 913 | double X=0.0, Y=0.0; |
| 914 | DHFHit_t *a = NULL__null; |
| 915 | for(unsigned int i=0;i<hits.size();i++){ |
| 916 | a = hits[i]; |
| 917 | double r = sqrt(pow((double)a->x,2.0) + pow((double)a->y, 2.0)); |
| 918 | // weight by r to give outer points more influence. Note that |
| 919 | // we really are really weighting by r^2 since x and y already |
| 920 | // have a magnitude component. |
| 921 | X += a->x*r; |
| 922 | Y += a->y*r; |
| 923 | } |
| 924 | phi = atan2(Y,X); |
| 925 | if(phi<0)phi+=M_TWO_PI6.28318530717958647692; |
| 926 | if(phi>=M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
| 927 | |
| 928 | // Search the chi2 space for values for p_trans, x0, ... |
| 929 | SearchPtrans(9.0, 0.5); |
| 930 | |
| 931 | #if 0 |
| 932 | // We do a simple linear regression here to find phi. This is |
| 933 | // simplified by the intercept being zero (i.e. the track |
| 934 | // passes through the beamline). |
| 935 | float Sxx=0.0, Syy=0.0, Sxy=0.0; |
| 936 | chisq_source = NOFIT; // in case we return early |
| 937 | |
| 938 | // Loop over hits to calculate Sxx, Syy, and Sxy |
| 939 | DHFHit_t *a = NULL__null; |
| 940 | for(unsigned int i=0;i<hits.size();i++){ |
| 941 | a = hits[i]; |
| 942 | float x=a->x; |
| 943 | float y=a->y; |
| 944 | Sxx += x*x; |
| 945 | Syy += y*y; |
| 946 | Sxy += x*y; |
| 947 | } |
| 948 | double A = 2.0*Sxy; |
| 949 | double B = Sxx - Syy; |
| 950 | phi = B>A ? atan2(A,B)/2.0 : 1.0/atan2(B,A)/2.0; |
| 951 | if(phi<0)phi+=M_TWO_PI6.28318530717958647692; |
| 952 | if(phi>=M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
| 953 | #endif |
| 954 | |
| 955 | return NOERROR; |
| 956 | } |
| 957 | |
| 958 | //----------------- |
| 959 | // SearchPtrans |
| 960 | //----------------- |
| 961 | void DHelicalFit::SearchPtrans(double ptrans_max, double ptrans_step) |
| 962 | { |
| 963 | /// Search the chi2 space as a function of the transverse momentum |
| 964 | /// for the minimal chi2. The values corresponding to the minmal |
| 965 | /// chi2 are left in chisq, x0, y0, r0, q, and p_trans. |
| 966 | /// |
| 967 | /// This does NOT optimize the value of p_trans. It simply |
| 968 | /// does a straight forward chisq calculation on a grid |
| 969 | /// and keeps the best one. It is intended for use in finding |
| 970 | /// a reasonable value for p_trans for straight tracks that |
| 971 | /// are contained to less than p_trans_max which is presumably |
| 972 | /// chosen based on a known θ. |
| 973 | |
| 974 | // Loop over several values for p_trans and calculate the |
| 975 | // chisq for each. |
| 976 | double min_chisq=1.0E6; |
| 977 | double min_x0=0.0, min_y0=0.0, min_r0=0.0; |
| 978 | for(double my_p_trans=ptrans_step; my_p_trans<=ptrans_max; my_p_trans+=ptrans_step){ |
| 979 | |
| 980 | r0 = my_p_trans/(0.003*2.0); |
| 981 | double alpha, my_chisq; |
| 982 | |
| 983 | // q = +1 |
| 984 | alpha = phi + M_PI_21.57079632679489661923; |
| 985 | x0 = r0*cos(alpha); |
| 986 | y0 = r0*sin(alpha); |
| 987 | my_chisq = ChisqCircle(); |
| 988 | if(my_chisq<min_chisq){ |
| 989 | min_chisq=my_chisq; |
| 990 | min_x0 = x0; |
| 991 | min_y0 = y0; |
| 992 | min_r0 = r0; |
| 993 | q = +1.0; |
| 994 | p_trans = my_p_trans; |
| 995 | } |
| 996 | |
| 997 | // q = -1 |
| 998 | alpha = phi - M_PI_21.57079632679489661923; |
| 999 | x0 = r0*cos(alpha); |
| 1000 | y0 = r0*sin(alpha); |
| 1001 | my_chisq = ChisqCircle(); |
| 1002 | if(my_chisq<min_chisq){ |
| 1003 | min_chisq=my_chisq; |
| 1004 | min_x0 = x0; |
| 1005 | min_y0 = y0; |
| 1006 | min_r0 = r0; |
| 1007 | q = -1.0; |
| 1008 | p_trans = my_p_trans; |
| 1009 | } |
| 1010 | } |
| 1011 | |
| 1012 | // Copy params from minimum chisq |
| 1013 | x0 = min_x0; |
| 1014 | y0 = min_y0; |
| 1015 | r0 = min_r0; |
| 1016 | } |
| 1017 | |
| 1018 | //----------------- |
| 1019 | // QuickPtrans |
| 1020 | //----------------- |
| 1021 | void DHelicalFit::QuickPtrans(void) |
| 1022 | { |
| 1023 | /// Quickly calculate a value of p_trans by looking |
| 1024 | /// for the hit furthest out and the hit closest |
| 1025 | /// to half that distance. Those 2 hits along with the |
| 1026 | /// origin are used to define a circle from which |
| 1027 | /// p_trans is calculated. |
| 1028 | |
| 1029 | // Find hit with largest R |
| 1030 | double R2max = 0.0; |
| 1031 | DHFHit_t *hit_max = NULL__null; |
| 1032 | for(unsigned int i=0; i<hits.size(); i++){ |
| 1033 | // use cross product to decide if hits is to left or right |
| 1034 | double x = hits[i]->x; |
| 1035 | double y = hits[i]->y; |
| 1036 | double R2 = x*x + y*y; |
| 1037 | if(R2>R2max){ |
| 1038 | R2max = R2; |
| 1039 | hit_max = hits[i]; |
| 1040 | } |
| 1041 | } |
| 1042 | |
| 1043 | // Bullet proof |
| 1044 | if(!hit_max)return; |
| 1045 | |
| 1046 | // Find hit closest to half-way out |
| 1047 | double Rmid = 0.0; |
| 1048 | double Rmax = sqrt(R2max); |
| 1049 | DHFHit_t *hit_mid = NULL__null; |
| 1050 | for(unsigned int i=0; i<hits.size(); i++){ |
| 1051 | // use cross product to decide if hits is to left or right |
| 1052 | double x = hits[i]->x; |
| 1053 | double y = hits[i]->y; |
| 1054 | double R = sqrt(x*x + y*y); |
| 1055 | if(fabs(R-Rmax/2.0) < fabs(Rmid-Rmax/2.0)){ |
| 1056 | Rmid = R; |
| 1057 | hit_mid = hits[i]; |
| 1058 | } |
| 1059 | } |
| 1060 | |
| 1061 | // Bullet proof |
| 1062 | if(!hit_mid)return; |
| 1063 | |
| 1064 | // Calculate p_trans from 3 points |
| 1065 | double x1 = hit_mid->x; |
| 1066 | double y1 = hit_mid->y; |
| 1067 | double x2 = hit_max->x; |
| 1068 | double y2 = hit_max->y; |
| 1069 | double r2 = sqrt(x2*x2 + y2*y2); |
| 1070 | double cos_phi = x2/r2; |
| 1071 | double sin_phi = y2/r2; |
| 1072 | double u1 = x1*cos_phi + y1*sin_phi; |
| 1073 | double v1 = -x1*sin_phi + y1*cos_phi; |
| 1074 | double u2 = x2*cos_phi + y2*sin_phi; |
| 1075 | double u0 = u2/2.0; |
| 1076 | double v0 = (u1*u1 + v1*v1 - u1*u2)/(2.0*v1); |
| 1077 | |
| 1078 | x0 = u0*cos_phi - v0*sin_phi; |
| 1079 | y0 = u0*sin_phi + v0*cos_phi; |
| 1080 | r0 = sqrt(x0*x0 + y0*y0); |
| 1081 | |
| 1082 | double B=-2.0; |
| 1083 | p_trans = fabs(0.003*B*r0); |
| 1084 | q = v1>0.0 ? -1.0:+1.0; |
| 1085 | } |
| 1086 | |
| 1087 | //----------------- |
| 1088 | // GuessChargeFromCircleFit |
| 1089 | //----------------- |
| 1090 | jerror_t DHelicalFit::GuessChargeFromCircleFit(void) |
| 1091 | { |
| 1092 | /// Adjust the sign of the charge (magnitude will stay 1) based on |
| 1093 | /// whether the hits tend to be to the right or to the left of the |
| 1094 | /// line going from the origin through the center of the circle. |
| 1095 | /// If the sign is flipped, the phi angle will also be shifted by |
| 1096 | /// +/- pi since the majority of hits are assumed to come from |
| 1097 | /// outgoing tracks. |
| 1098 | /// |
| 1099 | /// This is just a guess since tracks can bend all the way |
| 1100 | /// around and have hits that look exactly like tracks from an |
| 1101 | /// outgoing particle of opposite charge. The final charge should |
| 1102 | /// come from the sign of the dphi/dz slope. |
| 1103 | |
| 1104 | // Simply count the number of hits whose phi angle relative to the |
| 1105 | // phi of the center of the circle are greater than pi. |
| 1106 | unsigned int N=0; |
| 1107 | for(unsigned int i=0; i<hits.size(); i++){ |
| 1108 | // use cross product to decide if hits is to left or right |
| 1109 | double x = hits[i]->x; |
| 1110 | double y = hits[i]->y; |
| 1111 | if((x*y0 - y*x0) < 0.0)N++; |
| 1112 | } |
| 1113 | |
| 1114 | // Check if more hits are negative and make sign negative if so. |
| 1115 | if(N>hits.size()/2.0){ |
| 1116 | q = -1.0; |
| 1117 | phi += M_PI3.14159265358979323846; |
| 1118 | if(phi>M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
| 1119 | } |
| 1120 | |
| 1121 | return NOERROR; |
| 1122 | } |
| 1123 | |
| 1124 | //----------------- |
| 1125 | // FitTrack |
| 1126 | //----------------- |
| 1127 | jerror_t DHelicalFit::FitTrack(void) |
| 1128 | { |
| 1129 | /// Find theta, sign of electric charge, total momentum and |
| 1130 | /// vertex z position. |
| 1131 | |
| 1132 | // Points must be in order of increasing Z |
| 1133 | sort(hits.begin(), hits.end(), DHFHitLessThanZ_C); |
| 1134 | |
| 1135 | // Fit to circle to get circle's center |
| 1136 | FitCircle(); |
| 1137 | |
| 1138 | // The thing that is really needed is dphi/dz (where phi is the angle |
| 1139 | // of the point as measured from the center of the circle, not the beam |
| 1140 | // line). The relation between phi and z is linear so we use linear |
| 1141 | // regression to find the slope (dphi/dz). The one complication is |
| 1142 | // that phi is periodic so the value obtained via the x and y of a |
| 1143 | // specific point can be off by an integral number of 2pis. Handle this |
| 1144 | // by assuming the first point is measured before a full rotation |
| 1145 | // has occurred. Subsequent points should always be within 2pi of |
| 1146 | // the previous point so we just need to calculate the relative phi |
| 1147 | // between succesive points and keep a running sum. We do this in |
| 1148 | // the first loop were we find the mean z and phi values. The regression |
| 1149 | // formulae are calculated in the second loop. |
| 1150 | |
| 1151 | // Calculate phi about circle center for each hit |
| 1152 | Fill_phi_circle(hits, x0, y0); |
| 1153 | |
| 1154 | // Linear regression for z/phi relation |
| 1155 | float Sxx=0.0, Syy=0.0, Sxy=0.0; |
| 1156 | for(unsigned int i=0;i<hits.size();i++){ |
| 1157 | DHFHit_t *a = hits[i]; |
| 1158 | float deltaZ = a->z - z_mean; |
| 1159 | float deltaPhi = a->phi_circle - phi_mean; |
| 1160 | Syy += deltaZ*deltaZ; |
| 1161 | Sxx += deltaPhi*deltaPhi; |
| 1162 | Sxy += deltaZ*deltaPhi; |
| 1163 | //cout<<__FILE__<<":"<<__LINE__<<" deltaZ="<<deltaZ<<" deltaPhi="<<deltaPhi<<" Sxy(i)="<<deltaZ*deltaPhi<<endl; |
| 1164 | } |
| 1165 | float dzdphi = Syy/Sxy; |
| 1166 | dzdphi = Syy/Sxy; |
| 1167 | z_vertex = z_mean - phi_mean*dzdphi; |
| 1168 | //cout<<__FILE__<<":"<<__LINE__<<" z_mean="<<z_mean<<" phi_mean="<<phi_mean<<" dphidz="<<dphidz<<" Sxy="<<Sxy<<" Syy="<<Syy<<" z_vertex="<<z_vertex<<endl; |
| 1169 | |
| 1170 | // Fill in the rest of the parameters |
| 1171 | return FillTrackParams(); |
| 1172 | } |
| 1173 | |
| 1174 | //----------------- |
| 1175 | // FitTrackRiemann |
| 1176 | //----------------- |
| 1177 | jerror_t DHelicalFit::FitTrackRiemann(float rc_input){ |
| 1178 | jerror_t error=NOERROR; |
| 1179 | |
| 1180 | if (CovR_!=NULL__null) delete CovR_; |
| 1181 | if (CovRPhi_!=NULL__null) delete CovRPhi_; |
| 1182 | CovR_=NULL__null; |
| 1183 | CovRPhi_=NULL__null; |
| 1184 | |
| 1185 | error=FitCircleRiemannCorrected(rc_input); |
| 1186 | if (error!=NOERROR) return error; |
| 1187 | error=FitLineRiemann(); |
| 1188 | GetChargeRiemann(); |
| 1189 | |
| 1190 | // Shift phi by pi if the charge is negative |
| 1191 | if (q<0) phi+=M_PI3.14159265358979323846; |
| 1192 | |
| 1193 | return error; |
| 1194 | } |
| 1195 | |
| 1196 | |
| 1197 | jerror_t DHelicalFit::FitCircleAndLineRiemann(float rc_input){ |
| 1198 | jerror_t error=NOERROR; |
| 1199 | |
| 1200 | if (CovR_!=NULL__null) delete CovR_; |
| 1201 | if (CovRPhi_!=NULL__null) delete CovRPhi_; |
| 1202 | CovR_=NULL__null; |
| 1203 | CovRPhi_=NULL__null; |
| 1204 | |
| 1205 | error=FitCircleRiemannCorrected(rc_input); |
| 1206 | if (error!=NOERROR) return error; |
| 1207 | error=FitLineRiemann(); |
| 1208 | |
| 1209 | return error; |
| 1210 | } |
| 1211 | |
| 1212 | |
| 1213 | |
| 1214 | |
| 1215 | //----------------- |
| 1216 | // FitTrack_FixedZvertex |
| 1217 | //----------------- |
| 1218 | jerror_t DHelicalFit::FitTrack_FixedZvertex(float z_vertex) |
| 1219 | { |
| 1220 | /// Fit the points, but hold the z_vertex fixed at the specified value. |
| 1221 | /// |
| 1222 | /// This just calls FitCircle and FitLine_FixedZvertex to find all parameters |
| 1223 | /// of the track |
| 1224 | |
| 1225 | // Fit to circle to get circle's center |
| 1226 | FitCircle(); |
| 1227 | |
| 1228 | // Fit to line in phi-z plane and return error |
| 1229 | return FitLine_FixedZvertex(z_vertex); |
| 1230 | } |
| 1231 | |
| 1232 | //----------------- |
| 1233 | // FitLine_FixedZvertex |
| 1234 | //----------------- |
| 1235 | jerror_t DHelicalFit::FitLine_FixedZvertex(float z_vertex) |
| 1236 | { |
| 1237 | /// Fit the points, but hold the z_vertex fixed at the specified value. |
| 1238 | /// |
| 1239 | /// This assumes FitCircle has already been called and the values |
| 1240 | /// in x0 and y0 are valid. |
| 1241 | /// |
| 1242 | /// This just fits the phi-z angle by minimizing the chi-squared |
| 1243 | /// using a linear regression technique. As it turns out, the |
| 1244 | /// chi-squared weights points by their distances squared which |
| 1245 | /// leads to a quadratic equation for cot(theta) (where theta is |
| 1246 | /// the angle in the phi-z plane). To pick the right solution of |
| 1247 | /// the quadratic equation, we pick the one closest to the linearly |
| 1248 | /// weighted angle. One could argue that we should just use the |
| 1249 | /// linear weighting here rather than the square weighting. The |
| 1250 | /// choice depends though on how likely the "out-lier" points |
| 1251 | /// are to really be from this track. If they are likely, to |
| 1252 | /// be, then we would want to leverage their "longer" lever arms |
| 1253 | /// with the squared weighting. If they are more likely to be bad |
| 1254 | /// points, then we would want to minimize their influence with |
| 1255 | /// a linear (or maybe even root) weighting. It is expected than |
| 1256 | /// for our use, the points are all likely valid so we use the |
| 1257 | /// square weighting. |
| 1258 | |
| 1259 | // Points must be in order of increasing Z |
| 1260 | sort(hits.begin(), hits.end(), DHFHitLessThanZ_C); |
| 1261 | |
| 1262 | // Fit is being done for a fixed Z-vertex |
| 1263 | this->z_vertex = z_vertex; |
| 1264 | |
| 1265 | // Calculate phi about circle center for each hit |
| 1266 | Fill_phi_circle(hits, x0, y0); |
| 1267 | |
| 1268 | // Do linear regression on phi-Z |
| 1269 | float Sx=0, Sy=0; |
| 1270 | float Sxx=0, Syy=0, Sxy = 0; |
| 1271 | float r0 = sqrt(x0*x0 + y0*y0); |
| 1272 | for(unsigned int i=0; i<hits.size(); i++){ |
| 1273 | DHFHit_t *a = hits[i]; |
| 1274 | |
| 1275 | float dz = a->z - z_vertex; |
| 1276 | float dphi = a->phi_circle; |
| 1277 | Sx += dz; |
| 1278 | Sy += r0*dphi; |
| 1279 | Syy += r0*dphi*r0*dphi; |
| 1280 | Sxx += dz*dz; |
| 1281 | Sxy += r0*dphi*dz; |
| 1282 | } |
| 1283 | |
| 1284 | float k = (Syy-Sxx)/(2.0*Sxy); |
| 1285 | float s = sqrt(k*k + 1.0); |
| 1286 | float cot_theta1 = -k+s; |
| 1287 | float cot_theta2 = -k-s; |
| 1288 | float cot_theta_lin = Sx/Sy; |
| 1289 | float cot_theta; |
| 1290 | if(fabs(cot_theta_lin-cot_theta1) < fabs(cot_theta_lin-cot_theta2)){ |
| 1291 | cot_theta = cot_theta1; |
| 1292 | }else{ |
| 1293 | cot_theta = cot_theta2; |
| 1294 | } |
| 1295 | |
| 1296 | dzdphi = r0*cot_theta; |
| 1297 | |
| 1298 | // Fill in the rest of the paramters |
| 1299 | return FillTrackParams(); |
| 1300 | } |
| 1301 | |
| 1302 | //------------------------------------------------------------------ |
| 1303 | // Fill_phi_circle |
| 1304 | //------------------------------------------------------------------ |
| 1305 | jerror_t DHelicalFit::Fill_phi_circle(vector<DHFHit_t*> hits, float x0, float y0) |
| 1306 | { |
| 1307 | float x_last = -x0; |
| 1308 | float y_last = -y0; |
| 1309 | float phi_last = 0.0; |
| 1310 | z_mean = phi_mean = 0.0; |
| 1311 | for(unsigned int i=0; i<hits.size(); i++){ |
| 1312 | DHFHit_t *a = hits[i]; |
| 1313 | |
| 1314 | float dx = a->x - x0; |
| 1315 | float dy = a->y - y0; |
| 1316 | 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 ); |
| 1317 | float my_phi = phi_last +dphi; |
| 1318 | a->phi_circle = my_phi; |
| 1319 | |
| 1320 | z_mean += a->z; |
| 1321 | phi_mean += my_phi; |
| 1322 | |
| 1323 | x_last = dx; |
| 1324 | y_last = dy; |
| 1325 | phi_last = my_phi; |
| 1326 | } |
| 1327 | z_mean /= (float)hits.size(); |
| 1328 | phi_mean /= (float)hits.size(); |
| 1329 | |
| 1330 | return NOERROR; |
| 1331 | } |
| 1332 | |
| 1333 | //------------------------------------------------------------------ |
| 1334 | // FillTrackParams |
| 1335 | //------------------------------------------------------------------ |
| 1336 | jerror_t DHelicalFit::FillTrackParams(void) |
| 1337 | { |
| 1338 | /// Fill in and tweak some parameters like q, phi, theta using |
| 1339 | /// other values already set in the class. This is used by |
| 1340 | /// both FitTrack() and FitTrack_FixedZvertex(). |
| 1341 | |
| 1342 | float r0 = sqrt(x0*x0 + y0*y0); |
| 1343 | theta = atan(r0/fabs(dzdphi)); |
| 1344 | |
| 1345 | // The sign of the electric charge will be opposite that |
| 1346 | // of dphi/dz. Also, the value of phi will be PI out of phase |
| 1347 | if(dzdphi<0.0){ |
| 1348 | phi += M_PI3.14159265358979323846; |
| 1349 | if(phi<0)phi+=M_TWO_PI6.28318530717958647692; |
| 1350 | if(phi>=M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
| 1351 | }else{ |
| 1352 | q = -q; |
| 1353 | } |
| 1354 | |
| 1355 | // Re-calculate chi-sq (needs to be done) |
| 1356 | chisq_source = TRACK; |
| 1357 | |
| 1358 | // Up to now, the fit has assumed a forward going particle. In |
| 1359 | // other words, if the particle is going backwards, the helix does |
| 1360 | // actually still go through the points, but only if extended |
| 1361 | // backward in time. We use the average z of the hits compared |
| 1362 | // to the reconstructed vertex to determine if it was back-scattered. |
| 1363 | if(z_mean<z_vertex){ |
| 1364 | // back scattered particle |
| 1365 | theta = M_PI3.14159265358979323846 - theta; |
| 1366 | phi += M_PI3.14159265358979323846; |
| 1367 | if(phi<0)phi+=M_TWO_PI6.28318530717958647692; |
| 1368 | if(phi>=M_TWO_PI6.28318530717958647692)phi-=M_TWO_PI6.28318530717958647692; |
| 1369 | q = -q; |
| 1370 | } |
| 1371 | |
| 1372 | // There is a problem with theta for data generated with a uniform |
| 1373 | // field. The following factor is emprical until I can figure out |
| 1374 | // what the source of the descrepancy is. |
| 1375 | //theta += 0.1*pow((double)sin(theta),2.0); |
| 1376 | |
| 1377 | p = fabs(p_trans/sin(theta)); |
| 1378 | |
| 1379 | return NOERROR; |
| 1380 | } |
| 1381 | |
| 1382 | //------------------------------------------------------------------ |
| 1383 | |
| 1384 | //------------------------------------------------------------------ |
| 1385 | jerror_t DHelicalFit::Print(void) const |
| 1386 | { |
| 1387 | cout<<"-- DHelicalFit Params ---------------"<<endl; |
| 1388 | cout<<" x0 = "<<x0<<endl; |
| 1389 | cout<<" y0 = "<<y0<<endl; |
| 1390 | cout<<" q = "<<q<<endl; |
| 1391 | cout<<" p = "<<p<<endl; |
| 1392 | cout<<" p_trans = "<<p_trans<<endl; |
| 1393 | cout<<" phi = "<<phi<<endl; |
| 1394 | cout<<" theta = "<<theta<<endl; |
| 1395 | cout<<" z_vertex = "<<z_vertex<<endl; |
| 1396 | cout<<" chisq = "<<chisq<<endl; |
| 1397 | cout<<"chisq_source = "; |
| 1398 | switch(chisq_source){ |
| 1399 | case NOFIT: cout<<"NOFIT"; break; |
| 1400 | case CIRCLE: cout<<"CIRCLE"; break; |
| 1401 | case TRACK: cout<<"TRACK"; break; |
| 1402 | } |
| 1403 | cout<<endl; |
| 1404 | cout<<" Bz(avg) = "<<Bz_avg<<endl; |
| 1405 | |
| 1406 | return NOERROR; |
| 1407 | } |
| 1408 | |
| 1409 | |
| 1410 | //------------------------------------------------------------------ |
| 1411 | // Dump |
| 1412 | //------------------------------------------------------------------ |
| 1413 | jerror_t DHelicalFit::Dump(void) const |
| 1414 | { |
| 1415 | Print(); |
| 1416 | |
| 1417 | for(unsigned int i=0;i<hits.size();i++){ |
| 1418 | DHFHit_t *v = hits[i]; |
| 1419 | cout<<" x="<<v->x<<" y="<<v->y<<" z="<<v->z; |
| 1420 | cout<<" phi_circle="<<v->phi_circle<<" chisq="<<v->chisq<<endl; |
| 1421 | } |
| 1422 | |
| 1423 | return NOERROR; |
| 1424 | } |