11 #define qBr2p 0.003 // conversion for converting q*B*r to GeV/c
46 double covy,
double covxy){
62 jerror_t error=NOERROR;
82 N(1,0)=
N(0,0)*(A(1,0)*A(0,2)-(A(0,0)-lambda)*A(1,2))
83 /(A(0,1)*A(2,1)-(A(1,1)-lambda)*A(0,2));
84 N(2,0)=
N(0,0)*(A(2,0)*(A(1,1)-lambda)-A(1,0)*A(2,1))
85 /(A(1,2)*A(2,1)-(A(2,2)-lambda)*(A(1,1)-lambda));
88 for (
int i=0;i<3;i++){
91 for (
int i=0;i<3;i++){
109 for (
unsigned int i=0;i<
hits.size();i++){
111 double stemp=rtemp/4./
rc;
112 double ctemp=1.-stemp*stemp;
126 for (
unsigned int i=0;i<
hits.size();i++){
129 =(Phi*cos(Phi)-
sin(Phi))*(Phi*cos(Phi)-
sin(Phi))*
hits[i]->covx
130 +(Phi*
sin(Phi)+cos(Phi))*(Phi*
sin(Phi)+cos(Phi))*
hits[i]->covy
131 +2.*(Phi*
sin(Phi)+cos(Phi))*(Phi*cos(Phi)-
sin(Phi))*
hits[i]->covxy;
136 for (
unsigned int m=0;m<
hits.size();m++){
138 CovR_->operator()(m,m)=cos(Phi)*cos(Phi)*
hits[m]->covx
140 +2.*
sin(Phi)*cos(Phi)*
hits[m]->covxy;
143 for (
unsigned int i=0;i<
hits.size();i++){
144 for (
unsigned int j=0;j<
hits.size();j++){
145 CR(i,j)=
CovR_->operator()(i, j);
146 CRPhi(i,j)=
CovRPhi_->operator()(i, j);
151 CRPhi=C*CRPhi*C+
S*CR*
S;
152 for (
unsigned int i=0;i<
hits.size();i++)
153 for (
unsigned int j=0;j<
hits.size();j++)
154 CovRPhi_->operator()(i,j)=CRPhi(i,j);
166 if (
hits.size()==0)
return RESOURCE_UNAVAILABLE;
170 double B0,B1,B2,Q,Q1,R,
sum,diff;
171 double theta,lambda_min=0.;
186 for (
unsigned int i=0;i<
hits.size();i++){
189 =(Phi*cos(Phi)-
sin(Phi))*(Phi*cos(Phi)-
sin(Phi))*
hits[i]->covx
190 +(Phi*
sin(Phi)+cos(Phi))*(Phi*
sin(Phi)+cos(Phi))*
hits[i]->covy
191 +2.*(Phi*
sin(Phi)+cos(Phi))*(Phi*cos(Phi)-
sin(Phi))*
hits[i]->covxy;
194 for (
unsigned int i=0;i<
hits.size();i++)
195 for (
unsigned int j=0;j<
hits.size();j++)
196 CRPhi(i,j)=
CovRPhi_->operator()(i, j);
207 for (
unsigned int i=0;i<
hits.size();i++){
211 Ones(i,0)=OnesT(0,i)=1.;
216 if (lu.Decompose()==
false){
217 return UNRECOVERABLE_ERROR;
219 W=
DMatrix(DMatrix::kInverted,CRPhi);
220 W_sum=OnesT*(W*Ones);
221 Xavg=(1./W_sum(0,0))*(OnesT*(W*
X));
223 A=
DMatrix(DMatrix::kTransposed,X)*(W*
X)
224 -W_sum(0,0)*(
DMatrix(DMatrix::kTransposed,Xavg)*Xavg);
225 if(!A.IsValid())
return UNRECOVERABLE_ERROR;
230 B2=-(A(0,0)+A(1,1)+A(2,2));
231 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)+A(1,1)*A(2,2)
234 if(B0==0 || !isfinite(B0))
return UNRECOVERABLE_ERROR;
250 Q=(3.*B1-B2*B2)/9.e4;
251 R=(9.*B2*B1-27.*B0-2.*B2*B2*B2)/54.e6;
253 if (Q1<0) Q1=
sqrt(-Q1);
255 return VALUE_OUT_OF_RANGE;
262 double temp=100.*pow(R*R+Q1*Q1,0.16666666666666666667);
263 theta=atan2(Q1,R)/3.;
264 sum=2.*temp*cos(theta);
265 diff=-2.*temp*
sin(theta);
267 lambda_min=-B2/3.-sum/2.+
sqrt(3.)/2.*diff;
277 dist_to_origin=-(N1(0,0)*Xavg(0,0)+N1(1,0)*Xavg(0,1)+N1(2,0)*Xavg(0,2));
280 xc=-N1(0,0)/2./N1(2,0);
281 yc=-N1(1,0)/2./N1(2,0);
298 for (
unsigned int i=0;i<
hits.size();i++){
300 double stemp=rtemp/4./rc_input;
301 double ctemp=1.-stemp*stemp;
315 for (
unsigned int i=0;i<
hits.size();i++){
318 =(Phi*cos(Phi)-
sin(Phi))*(Phi*cos(Phi)-
sin(Phi))*
hits[i]->covx
319 +(Phi*
sin(Phi)+cos(Phi))*(Phi*
sin(Phi)+cos(Phi))*
hits[i]->covy
320 +2.*(Phi*
sin(Phi)+cos(Phi))*(Phi*cos(Phi)-
sin(Phi))*
hits[i]->covxy;
325 for (
unsigned int m=0;m<
hits.size();m++){
327 CovR_->operator()(m,m)=cos(Phi)*cos(Phi)*
hits[m]->covx
329 +2.*
sin(Phi)*cos(Phi)*
hits[m]->covxy;
332 for (
unsigned int i=0;i<
hits.size();i++){
333 for (
unsigned int j=0;j<
hits.size();j++){
334 CR(i,j)=
CovR_->operator()(i, j);
335 CRPhi(i,j)=
CovRPhi_->operator()(i, j);
339 CRPhi=C*CRPhi*C+
S*CR*
S;
340 for (
unsigned int i=0;i<
hits.size();i++)
341 for (
unsigned int j=0;j<
hits.size();j++)
342 CovRPhi_->operator()(i,j)=CRPhi(i,j);
354 for (
unsigned int i=0;i<
hits.size();i++){
357 =(Phi*cos(Phi)-
sin(Phi))*(Phi*cos(Phi)-
sin(Phi))*
hits[i]->covx
358 +(Phi*
sin(Phi)+cos(Phi))*(Phi*
sin(Phi)+cos(Phi))*
hits[i]->covy
359 +2.*(Phi*
sin(Phi)+cos(Phi))*(Phi*cos(Phi)-
sin(Phi))*
hits[i]->covxy;
364 for (
unsigned int m=0;m<
hits.size();m++){
366 CovR_->operator()(m,m)=cos(Phi)*cos(Phi)*
hits[m]->covx
368 +2.*
sin(Phi)*cos(Phi)*
hits[m]->covxy;
371 for (
unsigned int i=0;i<
hits.size();i++){
372 for (
unsigned int j=0;j<
hits.size();j++){
373 CR(i,j)=
CovR_->operator()(i, j);
374 CRPhi(i,j)=
CovRPhi_->operator()(i, j);
379 double sumv=0,sumy=0,sumx=0,sumxx=0,sumxy=0;
380 for (
unsigned int k=0;k<
hits.size();k++){
382 double phi_z=atan2(hit->
y,hit->
x);
385 if (fabs(phi_z-phi_old)>M_PI){
386 if (phi_old<0) phi_z-=2.*M_PI;
389 double inv_var=(hit->
x*hit->
x+hit->
y*hit->
y)
390 /(CRPhi(k,k)+phi_z*phi_z*CR(k,k));
393 sumx+=hit->
z*inv_var;
394 sumxx+=hit->
z*hit->
z*inv_var;
395 sumxy+=phi_z*hit->
z*inv_var;
398 double slope=(sumv*sumxy-sumy*sumx)/(sumv*sumxx-sumx*sumx);
401 if (slope<0.)
return -1.;
416 for (
unsigned int m=0;m<
hits.size();m++){
418 CovR_->operator()(m,m)=cos(Phi)*cos(Phi)*
hits[m]->covx
420 +2.*
sin(Phi)*cos(Phi)*
hits[m]->covxy;
423 for (
unsigned int i=0;i<
hits.size();i++)
424 for (
unsigned int j=0;j<
hits.size();j++)
425 CR(i,j)=
CovR_->operator()(i, j);
428 double x_int0,
temp,y_int0;
429 double denom=
N[0]*
N[0]+
N[1]*
N[1];
431 vector<int>bad(
hits.size());
435 for (
unsigned int m=0;m<
hits.size();m++){
439 temphit->
z=
hits[m]->z;
447 x_int0=-N[0]*numer/denom;
448 y_int0=-N[1]*numer/denom;
449 temp=denom*r2-numer*numer;
458 temp=
sqrt(temp)/denom;
461 double diffx1=x_int0+N[1]*temp-
hits[m]->x;
462 double diffy1=y_int0-N[0]*temp-
hits[m]->y;
463 double diffx2=x_int0-N[1]*temp-
hits[m]->x;
464 double diffy2=y_int0+N[0]*temp-
hits[m]->y;
465 if (diffx1*diffx1+diffy1*diffy1 > diffx2*diffx2+diffy2*diffy2){
466 temphit->
x=x_int0-N[1]*
temp;
467 temphit->
y=y_int0+N[0]*
temp;
470 temphit->
x=x_int0+N[1]*
temp;
471 temphit->
y=y_int0-N[0]*
temp;
479 unsigned int start=0;
480 for (
unsigned int i=0;i<bad.size();i++){
489 double sumv=0.,sumx=0.,sumy=0.,sumxx=0.,sumxy=0.;
490 double sperp=0., chord=0,ratio=0, Delta;
491 for (
unsigned int k=start;k<n;k++){
496 chord=
sqrt(diffx*diffx+diffy*diffy);
500 sperp=2.*
rc*(M_PI/2.);
502 sperp=2.*
rc*asin(ratio);
516 sperp=2.*
rc*(M_PI/2.);
518 sperp=2.*
rc*asin(ratio);
519 Delta=sumv*sumxx-sumx*sumx;
521 tanl=-Delta/(sumv*sumxy-sumy*sumx);
525 double zvertex_temp=
projections[start]->z-(2.*
rc*M_PI-sperp)*tanl;
jerror_t CalcNormal(DMatrix A, double lambda, DMatrix &N)
double z
point in lab coordinates
jerror_t DoFit(double rc)
vector< DRiemannHit_t * > projections
jerror_t AddHit(double r, double phi, double z)
Add a hit to the list of hits using cylindrical coordinates.
vector< DRiemannHit_t * > hits
bool DRiemannFit_hit_cmp(DRiemannHit_t *a, DRiemannHit_t *b)
jerror_t AddHitXYZ(double x, double y, double z)
Add a hit to the list of hits using Cartesian coordinates.
double covxy
error info for x and y coordinates