11 #include <TDecompLU.h>
16 #define MOLIERE_FRACTION 0.99
17 #define ONE_THIRD 0.33333333333333333
18 #define ONE_SIXTH 0.16666666666666667
19 #define TWO_THIRDS 0.66666666666666667
27 #define DIFFUSION_COEFF 1.1e-6 // cm^2/s --> 200 microns at 1 cm
28 #define DRIFT_SPEED .0055
29 #define FDC_CATHODE_VARIANCE 0.02*0.02
43 =(108.55 + 7.62391*x + 556.176*exp(-(1.12566)*pow(x,1.29645)))*1
e-4;
45 return sigma_d*sigma_d;
78 for (
unsigned int i=0;i<
cdchits.size();i++){
94 for (
unsigned int i=0;i<
fdchits.size();i++){
114 CRPhi.ResizeTo(ncirclehits,ncirclehits);
124 for (
unsigned int i=0;i<
cdchits.size();i++){
162 Cz.ResizeTo(nstereo,nstereo);
168 CR.ResizeTo(nhits,nhits);
181 double cosl=cos(atan(
tanl));
182 p=0.003*fabs(
B)*
rc/cosl;
206 for (
unsigned int i=0;i<nhits;i++){
215 double drw=dXY.Mod();
224 double zplus=
GetStereoZ(dXYtest.X(),dXYtest.Y(),hit);
225 double zminus=
GetStereoZ(-dXYtest.X(),-dXYtest.Y(),hit);
231 else if (zplus<zpred){
236 else if (zpred<=17.){
241 else if (zplus>zpred){
247 if (fabs(zplus-zpred)<fabs(zminus-zpred)){
256 hit->
covx=var*dir.X()*dir.X();
257 hit->
covy=var*dir.Y()*dir.Y();
258 hit->
covz*=var/0.2133;
281 double cosphi=cos(
phi0);
296 double pt=0.003*fabs(
B)*
rc;
311 unsigned int ndf = this->
Ndof;
313 if(chisq_ptr)*chisq_ptr =
chisq;
314 if(dof_ptr)*dof_ptr = int(ndf);
317 printf(
"reduced chi2 %f\n",chisq/
double(ndf));
318 return chisq/double(ndf);
329 double y0=
D*cos(
phi0);
338 sperp+=2.*
rc*((ratio>1.)?M_PI_2:asin(ratio));
345 double cosPhi=cos(Phi);
346 double sinPhi=
sin(Phi);
361 sperp+=2.*
rc*((ratio>1.)?M_PI_2:asin(ratio));
368 double cosPhi=cos(Phi);
369 double sinPhi=
sin(Phi);
384 if (!isfinite(pos.Mag()) || !isfinite(mom.Mag())){
385 _DBG_ <<
"Invalid seed data." <<endl;
386 return UNRECOVERABLE_ERROR;
418 for (
unsigned int i=0;i<
fdchits.size();i++){
445 double drw=dXY.Mod();
448 double sign=(drw<
rc)?1.:-1.;
449 double ratio=(XY-XYold).Mod()/(2.*
rc);
450 sperp+=2.*
rc*((ratio>1.)?M_PI_2:asin(ratio));
452 double d_drift=0.0055*(hit->
cdc->
tdrift-tflight);
453 hit->
XY=XY+sign*d_drift*
dir;
457 hit->
covx=var*dir.X()*dir.X();
458 hit->
covy=var*dir.Y()*dir.Y();
471 double denom=ux*ux+uy*uy;
478 double my_z=zwire0+(dxwire0*ux+dywire0*uy)/denom;
479 double temp=uy*dxwire0-ux*dywire0;
480 double rootz2=denom*
rc*
rc-temp*
temp;
482 double rootz=
sqrt(rootz2)/denom;
483 double z1=my_z+rootz;
484 double z2=my_z-rootz;
507 double denom=ux*ux+uy*uy;
513 hit->
z=zwire0+(dxwire0*ux+dywire0*uy)/denom;
514 double temp=uy*dxwire0-ux*dywire0;
515 double rootz2=denom*
rc*
rc-temp*
temp;
516 double dx0dz=ux/denom;
517 double dy0dz=uy/denom;
519 double rootz=
sqrt(rootz2)/denom;
520 double z1=hit->
z+rootz;
521 double z2=hit->
z-rootz;
523 double temp1=temp/(rootz*denom*denom);
540 if (hit->
z>167.3) hit->
z=167.3;
541 if (hit->
z<17.0) hit->
z=17.0;
542 double dzwire=hit->
z-zwire0;
547 double drw=dXY.Mod();
550 hit->
covx=0.2133*dir.X()*dir.X();
551 hit->
covy=0.2133*dir.Y()*dir.Y();
562 double twokappa=
q/
rc;
563 double twoks=twokappa*sperp;
564 double sin2ks=
sin(twoks);
565 double cos2ks=cos(twoks);
566 double cosphi=cos(
phi0);
568 double one_minus_cos2ks=1.-cos2ks;
569 double one_over_2k=1./(twokappa);
571 DVector2(-
D*sinphi+(cosphi*sin2ks-sinphi*one_minus_cos2ks)*one_over_2k,
572 D*cosphi+(sinphi*sin2ks+cosphi*one_minus_cos2ks)*one_over_2k);
591 double w=XYp.X()*cosa-XYp.Y()*sina-hit->
fdc->
w;
593 double sign=(w>0?1.:-1.);
602 double u=hit->
fdc->
w+delta_x;
603 double v=hit->
fdc->
s-delta_y;
604 hit->
XY.Set(u*cosa+v*sina,-u*sina+v*cosa);
607 double cosa2=cosa*cosa;
608 double sina2=sina*sina;
611 hit->
covx=sigx2*cosa2+sigy2*sina2;
612 hit->
covy=sigx2*sina2+sigy2*cosa2;
613 hit->
covxy=(sigy2-sigx2)*sina*cosa;
629 for (
unsigned int i=0;i<nhits;i++){
631 double cosPhi=cos(Phi);
632 double sinPhi=
sin(Phi);
633 double Phi_sinPhi_plus_cosPhi=Phi*sinPhi+cosPhi;
634 double Phi_cosPhi_minus_sinPhi=Phi*cosPhi-sinPhi;
636 +Phi_sinPhi_plus_cosPhi*Phi_sinPhi_plus_cosPhi*
my_circle_hits[i]->covy
637 +2.*Phi_sinPhi_plus_cosPhi*Phi_cosPhi_minus_sinPhi*
my_circle_hits[i]->covxy;
648 double lambda=atan(
tanl);
649 double cosl=cos(lambda);
651 double cosl2=cosl*cosl;
652 for (
unsigned int m=0;m<nhits;m++){
654 for (
unsigned int k=m;k<nhits;k++){
656 unsigned int imax=(k>m)?m:k;
657 for (
unsigned int i=0;i<imax;i++){
660 if (std::isnan(sigma2_ms)){
664 CRPhi_ms(m,k)+=sigma2_ms*(Rk-Ri)*(Rm-Ri)/cosl2;
666 CRPhi_ms(k,m)=CRPhi_ms(m,k);
679 double ctemp=1.-stemp*stemp;
705 for (
unsigned int i=0;i<nhits;i++){
707 double cosPhi=cos(Phi);
708 double sinPhi=
sin(Phi);
719 double lambda=atan(
tanl);
720 double sinl=
sin(lambda);
722 double sinl4=sinl*sinl*sinl*sinl;
723 for (
unsigned int m=0;m<nhits;m++){
725 for (
unsigned int k=m;k<nhits;k++){
727 unsigned int imax=(k>m)?m:k;
728 for (
unsigned int i=0;i<imax;i++){
731 if (std::isnan(sigma2_ms)){
734 CR_ms(m,k)+=sigma2_ms*(zk-zi)*(zm-zi)/sinl4;
736 CR_ms(k,m)=CR_ms(m,k);
754 for (
unsigned int i=0;i<n;i++){
760 double lambda=atan(
tanl);
761 double cosl=cos(lambda);
763 double cosl4=cosl*cosl*cosl*cosl;
764 for (
unsigned int m=0;m<n;m++){
765 for (
unsigned int k=m;k<n;k++){
766 unsigned int imax=(k>m)?m:k;
767 for (
unsigned int i=0;i<imax;i++){
770 if (std::isnan(sigma2_ms)){
773 Cz_ms(m,k)+=sigma2_ms*(
s[k]-
s[i])*(
s[m]-
s[i])/cosl4;
775 Cz_ms(k,m)=Cz_ms(m,k);
787 double rho_Z_over_A,K_rho_Z_over_A,LnI,Z;
788 double chi2c_factor,chi2a_factor,chi2a_corr;
790 unsigned int dummy=0;
792 chi2c_factor,chi2a_factor,chi2a_corr,
799 double one_over_beta2=1.+
mass2/
p2;
801 double chi2c=chi2c_factor*my_ds*one_over_beta2/
p2;
802 double chi2a=chi2a_factor*(1.+chi2a_corr*one_over_beta2)/p2;
803 double nu=0.5*chi2c/(chi2a*(1.-
F));
804 return (2.*chi2c*1
e-6/(1.+F*F)*((1.+nu)/nu*log(1.+nu)-1.));
821 double B0,B1,B2,Q,Q1,R,
sum,diff;
822 double angle,lambda_min=0.;
824 DMatrix Ones(nhits,1),OnesT(1,nhits);
837 for (
unsigned int i=0;i<nhits;i++){
841 Ones(i,0)=OnesT(0,i)=1.;
846 if (lu.Decompose()==
false){
847 return UNRECOVERABLE_ERROR;
850 W_sum=OnesT*(W*Ones);
851 Xavg=(1./W_sum(0,0))*(OnesT*(W*X));
853 A=
DMatrix(DMatrix::kTransposed,X)*(W*
X)
854 -W_sum(0,0)*(
DMatrix(DMatrix::kTransposed,Xavg)*Xavg);
856 return UNRECOVERABLE_ERROR;
862 B2=-(A(0,0)+A(1,1)+A(2,2));
863 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)
866 if(B0==0 || !isfinite(B0)){
867 return UNRECOVERABLE_ERROR;
884 Q=(3.*B1-B2*B2)/9.e4;
885 R=(9.*B2*B1-27.*B0-2.*B2*B2*B2)/54.e6;
888 if (fabs(Q1)<
EPS) Q1=0.;
894 double temp=100.*pow(R*R+Q1*Q1,0.16666666666666666667);
895 angle=atan2(Q1,R)/3.;
896 sum=2.*temp*cos(angle);
897 diff=-2.*temp*
sin(angle);
899 lambda_min=-B2/3.-sum/2.+
sqrt(3.)/2.*diff;
909 (A(1,0)*A(0,2)-(A(0,0)-lambda_min)*A(1,2))
910 /(A(0,1)*A(2,1)-(A(1,1)-lambda_min)*A(0,2)),
911 (A(2,0)*(A(1,1)-lambda_min)-A(1,0)*A(2,1))
912 /(A(1,2)*A(2,1)-(A(2,2)-lambda_min)*(A(1,1)-lambda_min)));
918 c_origin=-(
N.X()*Xavg(0,0)+
N.Y()*Xavg(0,1)+
N.Z()*Xavg(0,2));
921 double one_over_2Nz=1./(2.*
N.Z());
922 xc=-
N.X()*one_over_2Nz;
923 yc=-
N.Y()*one_over_2Nz;
946 double x_int0,
temp,y_int0;
947 double denom=
N.Perp();
955 double ratio=numer/denom;
959 temp=denom*r2-numer*numer;
968 temp=
sqrt(temp)/denom;
971 double deltax=
N.y()*
temp;
972 double deltay=-
N.x()*
temp;
973 DVector2 XY1(x_int0+deltax,y_int0+deltay);
974 DVector2 XY2(x_int0-deltax,y_int0-deltay);
989 my_s=my_s+(chord_ratio>1.?2.*
rc*M_PI_2:2.*
rc*asin(chord_ratio));
1006 double sumv=0.,sumx=0.,sumy=0.,sumxx=0.,sumxy=0.;
1012 for (
unsigned int k=0;k<n;k++){
1015 double weight=1./
Cz(k,k);
1019 sumxx+=
s[k]*
s[k]*weight;
1020 sumxy+=
s[k]*z*weight;
1023 Delta=(sumv*sumxx-sumx*sumx);
1027 tanl=(sumv*sumxy-sumx*sumy)/Delta;
1029 z_vertex=(sumxx*sumy-sumx*sumxy)/Delta;
1033 if (z_vertex<Z_MIN || z_vertex>
Z_MAX){
1034 double weight=1000.;
1037 Delta= sumv*sumxx-sumx*sumx;
1038 tanl=(sumv*sumxy-sumx*sumy)/Delta;
1040 z_vertex=(sumxx*sumy-sumx*sumxy)/Delta;
1044 for (
unsigned int k=0;k<n;k++){
1048 double weight=1./
CR(k,k);
1053 sumxy+=
s[k]*z*weight;
1056 Delta= sumv*sumxx-sumx*sumx;
1057 double Delta1=sumv*sumxy-sumy*sumx;
1062 z_vertex=-(sumxx*sumy-sumx*sumxy)/Delta1;
1066 if (z_vertex<Z_MIN || z_vertex>
Z_MAX){
1067 double weight=1000.;
1071 Delta= sumv*sumxx-sumx*sumx;
1072 Delta1=sumv*sumxy-sumy*sumx;
1076 z_vertex=-(sumxx*sumy-sumx*sumxy)/Delta1;
#define FDC_CATHODE_VARIANCE
void setMomentum(const DVector3 &aMomentum)
double GetProcessNoise(const DVector2 &XY, const double z)
double z
point in lab coordinates
jerror_t GetStereoPosition(double &sperp, DVector2 &XYold, DRiemannHit_t *hit)
The DTrackFitter class is a base class for different charged track fitting algorithms. It does not actually fit the track itself, but provides the interface and some common support features most algorthims will need to implement.
jerror_t ComputeIntersections()
const DVector3 & position(void) const
const DFDCWire * wire
DFDCWire for this wire.
vector< DRiemannHit_t * > my_line_hits
fit_status_t FitTrack(void)
class DFDCPseudo: definition for a reconstructed point in the FDC
jerror_t SetSeed(double my_q, const DVector3 &pos, const DVector3 &mom, double mass)
const DMagneticFieldMap * bfield
vector< const DCDCTrackHit * > cdchits
double cdc_variance(double x)
double fdc_y_variance(double my_tanl, double x)
DTrackFitterRiemann(JEventLoop *loop)
jerror_t GetAxialPosition(double &sperp, const DVector2 &XYold, DRiemannHit_t *hit)
double s
< wire position computed from cathode data, assuming the avalanche occurs at the wire ...
double charge(void) const
double covxy
error info for x and y coordinates
double time
time corresponding to this pseudopoint.
void setPID(Particle_t locPID)
vector< DRiemannHit_t * > my_circle_hits
virtual double GetBz(double x, double y, double z) const =0
vector< const DFDCPseudo * > fdchits
const DVector3 & momentum(void) const
DTrackingData input_params
vector< DVector2 > projections
jerror_t GetFDCPosition(DRiemannHit_t *hit)
jerror_t FindMatKalman(const DVector3 &pos, const DVector3 &mom, double &KrhoZ_overA, double &rhoZ_overA, double &LnI, double &Z, double &chi2c_factor, double &chi2a_factor, double &chi2a_factor2, unsigned int &last_index, double *s_to_boundary=NULL) const
double GetStereoZ(double dx, double dy, DRiemannHit_t *hit)
DVector2 GetHelicalPosition(double sperp)
void setPosition(const DVector3 &aPosition)
static Particle_t IDTrack(float locCharge, float locMass)
printf("string=%s", string)