11 #include <JANA/JCalibration.h>
23 #define MAX_TB_PASSES 20
24 #define MAX_WB_PASSES 20
27 #define CHISQ_DELTA 0.01
37 if (fabs(a->
z-b->
z)<
EPS){
38 if (fabs(a->
t-b->
t)<
EPS){
41 if (fabs(tsum_1-tsum_2)<
EPS){
44 return (tsum_1<tsum_2);
51 if (a==NULL || b==NULL){
52 cout <<
"Null pointer in CDC hit list??" << endl;
69 int ascnd=(xx[n-1]>=xx[0]);
72 if ( (x>=xx[jm])==ascnd)
78 else if (x==xx[n-1]) *j=n-2;
87 if (x==xx[0])
return 0;
88 else if (x==xx[n-1])
return n-2;
92 int ascnd=(xx[n-1]>=xx[0]);
95 if ( (x>=xx[jm])==ascnd)
115 double &d,
double &V,
double &tcorr){
148 double my_t=0.001*tcorr;
149 double sqrt_t=
sqrt(my_t);
150 double t3=my_t*my_t*my_t;
151 double delta_mag=fabs(delta);
152 double a=a1+a2*delta_mag;
153 double b=b1+b2*delta_mag;
154 double c=c1+c2*delta_mag+c3*delta*delta;
155 f_delta=a*sqrt_t+b*my_t+c*t3;
156 f_0=a1*sqrt_t+b1*my_t+c1*t3;
158 dd_dt=0.001*(0.5*a/sqrt_t+b+3.*c*my_t*my_t);
161 double my_t=0.001*tcorr;
162 double sqrt_t=
sqrt(my_t);
163 double delta_mag=fabs(delta);
173 double delta_sq=delta*delta;
174 double a=a1+a2*delta_mag+a3*delta_sq;
175 double b=b1+b2*delta_mag+b3*delta_sq;
176 f_delta=a*sqrt_t+b*my_t;
177 f_0=a1*sqrt_t+b1*my_t;
179 dd_dt=0.001*(0.5*a/sqrt_t+b);
186 V=sigma*sigma+
mVarT0*dd_dt*dd_dt;
192 unsigned int index=0;
195 double frac=(tcorr-cdc_drift_table[
index])/dt;
196 double d_0=0.01*(double(index)+frac);
198 if (fabs(delta) <
EPS2){
200 V=sigma*sigma+
mVarT0*dd_dt*dd_dt;
208 d=f_delta*(d_0/f_0*P+1.-P);
209 V=sigma*sigma+
mVarT0*dd_dt*dd_dt;
216 V=sigma*sigma+
mVarT0*0.0001/(dt*dt);
222 #define FDC_T0_OFFSET 17.6
244 if (time<0.)
return 0.;
247 double tsq=time*time;
255 double t_high_sq=t_high*t_high;
273 double endplate_rmin,endplate_rmax;
282 vector<double>cdc_center;
283 vector<double>cdc_upstream_endplate_pos;
284 vector<double>cdc_endplate_dim;
286 geom->
Get(
"//posXYZ[@volume='centralDC']/@X_Y_Z",cdc_center);
287 geom->
Get(
"//posXYZ[@volume='CDPU']/@X_Y_Z",cdc_upstream_endplate_pos);
288 geom->
Get(
"//tubs[@name='CDPU']/@Rio_Z",cdc_endplate_dim);
289 cdc_origin[2]+=cdc_center[2]+cdc_upstream_endplate_pos[2]
290 +0.5*cdc_endplate_dim[2];
295 for(uint ring=0; ring<
cdcwires.size(); ring++)
301 vector<double>tof_face;
302 geom->
Get(
"//section/composition/posXYZ[@volume='ForwardTOF']/@X_Y_Z",
304 vector<double>tof_plane;
305 geom->
Get(
"//composition[@name='ForwardTOF']/posXYZ[@volume='forwardTOF']/@X_Y_Z/plane[@value='0']", tof_plane);
306 dTOFz=tof_face[2]+tof_plane[2];
307 geom->
Get(
"//composition[@name='ForwardTOF']/posXYZ[@volume='forwardTOF']/@X_Y_Z/plane[@value='1']", tof_plane);
308 dTOFz+=tof_face[2]+tof_plane[2];
316 for (
int i=0;i<30;i++){
317 vector<DVector3>
temp;
318 for (
unsigned int j=0;j<
sc_pos[i].size()-1;j++){
322 temp.push_back(
DVector3(dx/dz,dy/dz,1.));
343 gPARMS->SetDefaultParameter(
"KALMAN:THETA_CUT",
THETA_CUT);
346 gPARMS->SetDefaultParameter(
"KALMAN:RING_TO_SKIP",
RING_TO_SKIP);
349 gPARMS->SetDefaultParameter(
"KALMAN:PLANE_TO_SKIP",
PLANE_TO_SKIP);
356 gPARMS->SetDefaultParameter(
"KALMAN:PHOTON_ENERGY_CUTOFF",
360 gPARMS->SetDefaultParameter(
"TRKFIT:USE_FDC_HITS",
USE_FDC_HITS);
362 gPARMS->SetDefaultParameter(
"TRKFIT:USE_CDC_HITS",
USE_CDC_HITS);
366 gPARMS->SetDefaultParameter(
"TRKFIT:ALIGNMENT",
ALIGNMENT);
377 gPARMS->SetDefaultParameter(
"KALMAN:DEBUG_HISTS",
DEBUG_HISTS);
380 gPARMS->SetDefaultParameter(
"KALMAN:DEBUG_LEVEL",
DEBUG_LEVEL);
386 gPARMS->SetDefaultParameter(
"KALMAN:ESTIMATE_T0_TB",
ESTIMATE_T0_TB);
389 gPARMS->SetDefaultParameter(
"GEOM:ENABLE_BOUNDARY_CHECK",
393 gPARMS->SetDefaultParameter(
"TRKFIT:USE_MULS_COVARIANCE",
400 gPARMS->SetDefaultParameter(
"TRKFIT:USE_FDC_DRIFT_TIMES",
409 "maximum distance in number of sigmas away from projection to accept cdc hit");
411 "maximum distance in number of sigmas away from projection to accept fdc hit");
415 gPARMS->SetDefaultParameter(
"KALMAN:ANNEAL_SCALE",
ANNEAL_SCALE,
416 "Scale factor for annealing");
418 "Annealing parameter");
421 gPARMS->SetDefaultParameter(
"KALMAN:FORWARD_ANNEAL_SCALE",
423 "Scale factor for annealing");
424 gPARMS->SetDefaultParameter(
"KALMAN:FORWARD_ANNEAL_POW_CONST",
426 "Annealing parameter");
444 gPARMS->SetDefaultParameter(
"KALMAN:COVARIANCE_SCALE_FACTOR_CENTRAL",
448 gPARMS->SetDefaultParameter(
"KALMAN:COVARIANCE_SCALE_FACTOR_FORWARD",
452 gPARMS->SetDefaultParameter(
"KALMAN:WRITE_ML_TRAINING_OUTPUT",
456 mlfile.open(
"mltraining.dat");
461 JCalibration *jcalib = dapp->GetJCalibration((loop->GetJEvent()).GetRunNumber());
462 vector< map<string, double> > tvals;
464 if (jcalib->Get(
"CDC/cdc_drift_table", tvals)==
false){
465 for(
unsigned int i=0; i<tvals.size(); i++){
466 map<string, double> &row = tvals[i];
471 jerr <<
" CDC time-to-distance table not available... bailing..." << endl;
475 int straw_number[28]={42,42,54,54,66,66,80,80,93,93,106,106,
476 123,123,135,135,146,146,158,158,170,170,
477 182,182,197,197,209,209};
480 int straw_count=0,ring_count=0;
481 if (jcalib->Get(
"CDC/sag_parameters", tvals)==
false){
482 vector<double>
temp,temp2;
483 for(
unsigned int i=0; i<tvals.size(); i++){
484 map<string, double> &row = tvals[i];
486 temp.push_back(row[
"offset"]);
487 temp2.push_back(row[
"phi"]);
490 if (straw_count==straw_number[ring_count]){
501 if (jcalib->Get(
"CDC/drift_parameters", tvals)==
false){
502 map<string, double> &row = tvals[0];
529 map<string, double> cdc_drift_parms;
530 jcalib->Get(
"CDC/cdc_drift_parms", cdc_drift_parms);
534 map<string, double> cdc_res_parms;
535 jcalib->Get(
"CDC/cdc_resolution_parms_v2", cdc_res_parms);
541 map<string,double>lorentz_parms;
542 jcalib->Get(
"FDC/lorentz_deflection_parms",lorentz_parms);
549 map<string,double>drift_res_parms;
550 jcalib->Get(
"FDC/drift_resolution_parms",drift_res_parms);
556 map<string,double>drift_func_parms;
557 jcalib->Get(
"FDC/drift_function_parms",drift_func_parms);
564 map<string,double>drift_func_ext;
565 if (jcalib->Get(
"FDC/drift_function_ext",drift_func_ext)==
false){
570 map<string, double> fdc_drift_parms;
571 jcalib->Get(
"FDC/fdc_drift_parms", fdc_drift_parms);
590 for (
unsigned int i=0;i<5;i++)
I5x5(i,i)=1.;
594 map<string, double> targetparms;
595 if (jcalib->Get(
"TARGET/target_parms",targetparms)==
false){
596 TARGET_Z = targetparms[
"TARGET_Z_POSITION"];
602 gPARMS->SetDefaultParameter(
"KALMAN:VERTEX_POSITION",
TARGET_Z);
606 map<string, double> beam_vals;
607 jcalib->Get(
"PHOTON_BEAM/beam_spot",beam_vals);
609 beam_dir.Set(beam_vals[
"dxdz"],beam_vals[
"dydz"]);
611 <<
" z=" << beam_vals[
"z"]
616 static bool config_printed =
false;
618 config_printed =
true;
620 ss <<
"vertex constraint: " ;
622 ss <<
"z = " <<
TARGET_Z <<
"cm" << endl;
624 ss <<
"<none>" << endl;
630 for (
auto i=0; i < 46; i++){
631 double min = -10.,
max=10.;
632 if(i%23<12) {min=-5;
max=5;}
633 if(i<23)
alignDerivHists[i]=
new TH1I(Form(
"CentralDeriv%i",i),Form(
"CentralDeriv%i",i),200, min,
max);
634 else alignDerivHists[i]=
new TH1I(Form(
"ForwardDeriv%i",i%23),Form(
"ForwardDeriv%i",i%23),200, min,
max);
636 brentCheckHists[0]=
new TH2I(
"ForwardBrentCheck",
"DOCA vs ds", 100, -5., 5., 100, 0.0, 1.5);
637 brentCheckHists[1]=
new TH2I(
"CentralBrentCheck",
"DOCA vs ds", 100, -5., 5., 100, 0.0, 1.5);
658 for (
unsigned int i=0;i<
my_cdchits.size();i++){
661 for (
unsigned int i=0;i<
my_fdchits.size();i++){
743 unsigned int num_good_cdchits=
my_cdchits.size();
744 unsigned int num_good_fdchits=
my_fdchits.size();
751 if (num_good_cdchits>0){
758 for (
unsigned int i=0;i<
my_cdchits.size()-1;i++){
773 if (num_good_fdchits>0){
777 for (
unsigned int i=0;i<
my_fdchits.size()-1;i++){
800 if (num_good_cdchits==0 && num_good_fdchits<4)
return kFitNotDone;
801 if(num_good_cdchits+num_good_fdchits < 6)
return kFitNotDone;
855 _DBG_ <<
"------Starting "
857 <<
" Fit with " <<
my_fdchits.size() <<
" FDC hits and "
858 <<
my_cdchits.size() <<
" CDC hits.-------" <<endl;
860 _DBG_ <<
" Using t0=" <<
mT0 <<
" from DET="
868 _DBG_ <<
"Fit failed with Error = " << error <<endl;
884 _DBG_ <<
"----- Pass: "
887 <<
" p=" << mom.Mag()
888 <<
" theta=" << 90.0-180./M_PI*atan(
tanl_)
889 <<
" vertex=(" <<
x_ <<
"," <<
y_ <<
"," <<
z_<<
")"
894 for (
unsigned int iPull = 0; iPull <
pulls.size(); iPull++){
896 _DBG_ <<
" ring: " <<
pulls[iPull].cdc_hit->wire->ring
897 <<
" straw: " <<
pulls[iPull].cdc_hit->wire->straw
898 <<
" Residual: " <<
pulls[iPull].resi
899 <<
" Err: " <<
pulls[iPull].err
900 <<
" tdrift: " <<
pulls[iPull].tdrift
901 <<
" doca: " <<
pulls[iPull].d
902 <<
" docaphi: " <<
pulls[iPull].docaphi
903 <<
" z: " <<
pulls[iPull].z
904 <<
" cos(theta_rel): " <<
pulls[iPull].cosThetaRel
905 <<
" tcorr: " <<
pulls[iPull].tcorr
917 for (
unsigned int i=0;i<5;i++){
918 for (
unsigned int j=0;j<5;j++){
919 errMatrix(i,j)=
fcov[i][j];
937 else if (
cov.size()!=0){
942 for (
unsigned int i=0;i<5;i++){
943 for (
unsigned int j=0;j<5;j++){
944 errMatrix(i,j)=
cov[i][j];
953 locTrackingCovarianceMatrix->ResizeTo(5, 5);
954 for(
unsigned int loc_i = 0; loc_i < 5; ++loc_i)
956 for(
unsigned int loc_j = 0; loc_j < 5; ++loc_j)
957 (*locTrackingCovarianceMatrix)(loc_i, loc_j) = errMatrix(loc_i, loc_j);
966 set<const DCDCWire *> expected_hit_straws;
967 set<int> expected_hit_fdc_planes;
973 for(; ring<
cdc_rmid.size(); ring++) {
983 -
cdcwires[ring][0]->origin).Mag());
985 for(uint straw=1; straw<
cdcwires[ring].size(); straw++) {
989 double ds = dz*tan(
cdcwires[ring][straw]->stereo);
992 - wire_position).Mag());
993 if( diff < best_dist_diff )
997 expected_hit_straws.insert(
cdcwires[ring][best_straw]);
1007 for(uint plane=0; plane<
fdc_z_wires.size(); plane++) {
1008 int package = plane/6;
1011 expected_hit_fdc_planes.insert(plane);
1040 unsigned int ndf =
GetNDF();
1042 if(chisq_ptr)*chisq_ptr =
chisq;
1043 if(dof_ptr)*dof_ptr = int(ndf);
1044 if(pulls_ptr)*pulls_ptr =
pulls;
1046 return chisq/double(ndf);
1064 hit->
t=fdchit->
time;
1077 hit->
dE=1e6*fdchit->
dE;
1093 double one_over_uz=1./cdchit->
wire->
udir.z();
1095 one_over_uz*cdchit->
wire->
udir.y());
1121 double kq_over_p=
qBr2p*q_over_p;
1125 double one_plus_tx2=1.+tx2;
1126 double dsdz=
sqrt(one_plus_tx2+ty2);
1127 double dtx_Bfac=ty*
Bz+txty*
Bx-one_plus_tx2*
By;
1128 double dty_Bfac=
Bx*(1.+ty2)-txty*By-tx*
Bz;
1129 double kq_over_p_dsdz=kq_over_p*dsdz;
1134 D(
state_tx)=kq_over_p_dsdz*dtx_Bfac;
1135 D(
state_ty)=kq_over_p_dsdz*dty_Bfac;
1139 double q_over_p_sq=q_over_p*q_over_p;
1153 bool stepped_to_boundary=
false;
1162 for (
unsigned int i=0;i<
my_cdchits.size();i++){
1168 while(z<endplate_z_downstream && z>
cdc_origin[2] &&
1174 !=NOERROR)
return UNRECOVERABLE_ERROR;
1190 for (
int j=0;j<forward_traj_length-i;j++){
1196 if (
forward_traj.size()<2)
return RESOURCE_UNAVAILABLE;
1200 cout <<
"--- Forward cdc trajectory ---" <<endl;
1205 double phi=atan2(ty,tx);
1206 double cosphi=cos(phi);
1207 double sinphi=
sin(phi);
1209 double tanl=1./
sqrt(tx*tx+ty*ty);
1210 double sinl=
sin(atan(tanl));
1211 double cosl=cos(atan(tanl));
1213 << setiosflags(ios::fixed)<<
"pos: " << setprecision(4)
1217 << p*cosphi*cosl<<
", " << p*sinphi*cosl <<
", "
1218 << p*sinl <<
" -> " << p
1219 <<
" s: " << setprecision(3)
1221 <<
" t: " << setprecision(3)
1242 double &z,
double &r2,
1244 bool &stepped_to_boundary){
1251 double s_to_boundary=1e6;
1266 double one_over_beta2=1.+
mass2*q_over_p_sq;
1267 if (one_over_beta2>
BIG) one_over_beta2=
BIG;
1277 &s_to_boundary)!=NOERROR){
1278 return UNRECOVERABLE_ERROR;
1288 return UNRECOVERABLE_ERROR;
1315 double max_step_size
1319 if (!stepped_to_boundary){
1320 stepped_to_boundary=
false;
1321 if (fabs(dEdx)>
EPS){
1324 if (ds>max_step_size) ds=max_step_size;
1325 if (s_to_boundary<ds){
1326 ds=s_to_boundary+
EPS3;
1327 stepped_to_boundary=
true;
1333 stepped_to_boundary=
false;
1336 double newz=z+ds*dz_ds;
1388 double &var_t_factor,
1390 bool &stepped_to_boundary){
1397 double s_to_boundary=1e6;
1401 double q_over_p_sq=q_over_p*q_over_p;
1402 double one_over_beta2=1.+
mass2*q_over_p_sq;
1403 if (one_over_beta2>
BIG) one_over_beta2=
BIG;
1428 return UNRECOVERABLE_ERROR;
1436 return UNRECOVERABLE_ERROR;
1462 if (stepped_to_boundary){
1464 stepped_to_boundary=
false;
1467 if (fabs(dEdx)>
EPS){
1471 if (s_to_boundary<step_size){
1472 step_size=s_to_boundary+
EPS3;
1473 stepped_to_boundary=
true;
1477 double r2=my_xy.Mod2();
1487 double dt=step_size*
sqrt(one_over_beta2);
1488 double one_minus_beta2=1.-1./one_over_beta2;
1490 var_ftime+=dt*dt*one_minus_beta2*one_minus_beta2*0.0004;
1491 var_t_factor=dt*dt*one_minus_beta2*one_minus_beta2;
1542 bool stepped_to_boundary=
false;
1546 double var_t_factor=0;
1556 double r2=xy.Mod2(),z=
z_;
1559 for (
unsigned int j=0;j<
my_cdchits.size();j++){
1568 if (
PropagateCentral(central_traj_length,i,my_xy,var_t_factor,Sc,stepped_to_boundary)
1569 !=NOERROR)
return UNRECOVERABLE_ERROR;
1578 for (
int j=0;j<central_traj_length-i;j++){
1592 if (
central_traj.size()<2)
return RESOURCE_UNAVAILABLE;
1596 cout <<
"---------" <<
central_traj.size() <<
" entries------" <<endl;
1606 << setiosflags(ios::fixed)<<
"pos: " << setprecision(4)
1610 << pt*cosphi <<
", " << pt*sinphi <<
", "
1611 << pt*tanl <<
" -> " << pt/cos(atan(tanl))
1612 <<
" s: " << setprecision(3)
1614 <<
" t: " << setprecision(3)
1630 double &z,
double zhit,
1632 bool &stepped_to_boundary,
1633 bool &stepped_to_endplate){
1641 double s_to_boundary=1e6;
1646 double r2=pos.Perp2();
1655 double one_over_beta2=1.+
mass2*q_over_p_sq;
1656 if (one_over_beta2>
BIG) one_over_beta2=
BIG;
1668 return UNRECOVERABLE_ERROR;
1678 return UNRECOVERABLE_ERROR;
1707 double max_step_size
1711 if (!stepped_to_boundary){
1712 stepped_to_boundary=
false;
1713 if (fabs(dEdx)>
EPS){
1716 if (ds>max_step_size) ds=max_step_size;
1717 if (s_to_boundary<ds){
1718 ds=s_to_boundary+
EPS3;
1719 stepped_to_boundary=
true;
1726 stepped_to_boundary=
false;
1730 double dz=stepped_to_endplate ?
endplate_dz : ds*dz_ds;
1736 stepped_to_endplate=
true;
1809 double zhit=0.,old_zhit=0.;
1810 bool stepped_to_boundary=
false;
1811 bool stepped_to_endplate=
false;
1818 || z>400. || z<
Z_MIN
1824 if (fabs(old_zhit-zhit)>
EPS){
1831 || z>400. || z<
Z_MIN
1837 stepped_to_boundary,stepped_to_endplate)
1839 return UNRECOVERABLE_ERROR;
1847 if (m<2)
return UNRECOVERABLE_ERROR;
1854 stepped_to_boundary,stepped_to_endplate)
1856 return UNRECOVERABLE_ERROR;
1858 stepped_to_boundary,stepped_to_endplate)
1860 return UNRECOVERABLE_ERROR;
1868 for (
int j=0;j<mylen-i;j++){
1878 if (zhit<z) my_id=m;
1884 if (z<zhit+
EPS2)
break;
1902 double ds=
Step(z,newz,0.,S);
1909 double one_over_beta2=1.+
mass2*q_over_p_sq;
1910 if (one_over_beta2>
BIG) one_over_beta2=
BIG;
1924 if (
forward_traj.size()<2)
return RESOURCE_UNAVAILABLE;
1929 unsigned int hit_id=my_id-1;
1973 cout <<
"--- Forward fdc trajectory ---" <<endl;
1977 double phi=atan2(ty,tx);
1978 double cosphi=cos(phi);
1979 double sinphi=
sin(phi);
1981 double tanl=1./
sqrt(tx*tx+ty*ty);
1982 double sinl=
sin(atan(tanl));
1983 double cosl=cos(atan(tanl));
1985 << setiosflags(ios::fixed)<<
"pos: " << setprecision(4)
1989 << p*cosphi*cosl<<
", " << p*sinphi*cosl <<
", "
1990 << p*sinl <<
" -> " << p
1991 <<
" s: " << setprecision(3)
1993 <<
" t: " << setprecision(3)
2013 double delta_z=newz-oldz;
2014 if (fabs(delta_z)<
EPS)
return 0.;
2019 double ds=
sqrt(1.+tx*tx+ty*ty)*delta_z;
2021 double delta_z_over_2=0.5*delta_z;
2022 double midz=oldz+delta_z_over_2;
2029 double Bx0=
Bx,By0=
By,Bz0=
Bz;
2044 S1=S+delta_z_over_2*D2;
2092 double delta_z=newz-oldz;
2093 if (fabs(delta_z)<
EPS)
return 0.;
2098 double ds=
sqrt(1.+tx*tx+ty*ty)*delta_z;
2100 double delta_z_over_2=0.5*delta_z;
2101 double midz=oldz+delta_z_over_2;
2103 double Bx0=
Bx,By0=
By,Bz0=
Bz;
2121 S1=S+delta_z_over_2*D2;
2176 double delta_z=newz-oldz;
2177 if (fabs(delta_z)<
EPS)
return NOERROR;
2198 double kq_over_p=
qBr2p*q_over_p;
2200 double twotx2=2.*tx2;
2202 double twoty2=2.*ty2;
2204 double one_plus_tx2=1.+tx2;
2205 double one_plus_ty2=1.+ty2;
2206 double one_plus_twotx2_plus_ty2=one_plus_ty2+twotx2;
2207 double one_plus_twoty2_plus_tx2=one_plus_tx2+twoty2;
2208 double dsdz=
sqrt(1.+tx2+ty2);
2209 double ds=dsdz*delta_z;
2210 double kds=
qBr2p*ds;
2211 double kqdz_over_p_over_dsdz=kq_over_p*delta_z/dsdz;
2212 double kq_over_p_ds=kq_over_p*ds;
2213 double dtx_Bdep=ty*
Bz+txty*
Bx-one_plus_tx2*
By;
2214 double dty_Bdep=
Bx*one_plus_ty2-txty*By-tx*
Bz;
2217 double Bztxty=Bz*txty;
2225 J(
state_tx,
state_tx)+=kqdz_over_p_over_dsdz*(Bxty*(one_plus_twotx2_plus_ty2)
2226 -Bytx*(3.*one_plus_tx2+twoty2)
2230 -Bytx*(one_plus_twoty2_plus_tx2)
2234 *((Bxtx+
Bz)*(one_plus_twoty2_plus_tx2)-Byty*one_plus_tx2);
2237 -Bxtx*one_plus_ty2);
2240 double one_over_p_sq=q_over_p*q_over_p;
2243 double temp=-(q_over_p*one_over_p_sq/dsdz)*E*dEdx*delta_z;
2261 double cosphi=cos(phi);
2262 double sinphi=
sin(phi);
2263 double lambda=atan(tanl);
2264 double sinl=
sin(lambda);
2265 double cosl=cos(lambda);
2268 double pt=fabs(1./q_over_pt);
2274 double kq_over_pt=
qBr2p*q_over_pt;
2277 double By_cosphi_minus_Bx_sinphi=
By*cosphi-
Bx*sinphi;
2279 =kq_over_pt*q_over_pt*sinl*By_cosphi_minus_Bx_sinphi;
2280 double one_over_cosl=1./cosl;
2282 double p=pt*one_over_cosl;
2290 D1(
state_tanl)=kq_over_pt*By_cosphi_minus_Bx_sinphi*one_over_cosl;
2294 dpos.Set(cosl*cosphi,cosl*sinphi);
2313 double cosphi=cos(phi);
2314 double sinphi=
sin(phi);
2315 double lambda=atan(tanl);
2316 double sinl=
sin(lambda);
2317 double cosl=cos(lambda);
2318 double cosl2=cosl*cosl;
2319 double cosl3=cosl*cosl2;
2320 double one_over_cosl=1./cosl;
2323 double pt=fabs(1./q_over_pt);
2324 double q=pt*q_over_pt;
2330 double kq_over_pt=
qBr2p*q_over_pt;
2338 double By_cosphi_minus_Bx_sinphi=By*cosphi-
Bx*sinphi;
2339 double By_sinphi_plus_Bx_cosphi=By*sinphi+
Bx*cosphi;
2340 D1(
state_q_over_pt)=kq_over_pt*q_over_pt*sinl*By_cosphi_minus_Bx_sinphi;
2341 D1(
state_phi)=kq_over_pt*(By_sinphi_plus_Bx_cosphi*sinl-Bz*cosl);
2342 D1(
state_tanl)=kq_over_pt*By_cosphi_minus_Bx_sinphi*one_over_cosl;
2346 dxy.Set(cosl*cosphi,cosl*sinphi);
2351 =
qBr2p*(By_sinphi_plus_Bx_cosphi*sinl-Bz*cosl);
2355 =kq_over_pt*(dBxdz*cosphi*sinl+dBydz*sinphi*sinl-
dBzdz*cosl);
2362 =-kq_over_pt*q_over_pt*sinl*By_sinphi_plus_Bx_cosphi;
2364 =2.*kq_over_pt*sinl*By_cosphi_minus_Bx_sinphi;
2366 =kq_over_pt*q_over_pt*cosl3*By_cosphi_minus_Bx_sinphi;
2368 double p=pt*one_over_cosl;
2370 double m2_over_p2=
mass2/p_sq;
2378 =kq_over_pt*q_over_pt*sinl*(dBydz*cosphi-dBxdz*sinphi);
2393 if (fabs(ds)<
EPS)
return NOERROR;
2399 double Bx0=
Bx,By0=
By,Bz0=
Bz;
2415 double dx=ds_2*dxy1.X();
2416 double dy=ds_2*dxy1.Y();
2417 Bx=Bx0+dBxdx*dx+dBxdy*dy+dBxdz*dz;
2418 By=By0+dBydx*dx+dBydy*dy+dBydz*dz;
2419 Bz=Bz0+dBzdx*dx+dBzdy*dy+
dBzdz*dz;
2429 Bx=Bx0+dBxdx*dx+dBxdy*dy+dBxdz*dz;
2430 By=By0+dBydx*dx+dBydy*dy+dBydz*dz;
2431 Bz=Bz0+dBzdx*dx+dBzdy*dy+
dBzdz*dz;
2441 Bx=Bx0+dBxdx*dx+dBxdy*dy+dBxdz*dz;
2442 By=By0+dBydx*dx+dBydy*dy+dBydz*dz;
2443 Bz=Bz0+dBzdx*dx+dBzdy*dy+
dBzdz*dz;
2454 dxy+=ds_over_3*dxy2;
2455 dxy+=ds_over_3*dxy3;
2456 dxy+=ds_over_6*dxy4;
2462 double B=
sqrt(Bx0*Bx0+By0*By0+Bz0*Bz0);
2465 double q=(qpt>0.)?1.:-1.;
2466 double qrc_old=qpt/(
qBr2p*B);
2469 double qrc_plus_D=
S(
state_D)+qrc_old;
2472 double dx_sinphi_minus_dy_cosphi=dx*sinphi-dy*cosphi;
2473 double rc=
sqrt(dxy.Mod2()
2474 +2.*qrc_plus_D*(dx_sinphi_minus_dy_cosphi)
2475 +qrc_plus_D*qrc_plus_D);
2476 double q_over_rc=q/rc;
2478 J(
state_D,
state_D)=q_over_rc*(dx_sinphi_minus_dy_cosphi+qrc_plus_D);
2515 if (fabs(ds)<
EPS)
return NOERROR;
2522 double Bx0=
Bx,By0=
By,Bz0=
Bz;
2533 double dx=ds_2*dxy1.X();
2534 double dy=ds_2*dxy1.Y();
2598 if (fabs(ds)<
EPS)
return NOERROR;
2604 double Bx0=
Bx,By0=
By,Bz0=
Bz;
2619 double dx=ds_2*dxy1.X();
2620 double dy=ds_2*dxy1.Y();
2621 Bx=Bx0+dBxdx*dx+dBxdy*dy+dBxdz*dz;
2622 By=By0+dBydx*dx+dBydy*dy+dBydz*dz;
2623 Bz=Bz0+dBzdx*dx+dBzdy*dy+
dBzdz*dz;
2633 Bx=Bx0+dBxdx*dx+dBxdy*dy+dBxdz*dz;
2634 By=By0+dBydx*dx+dBydy*dy+dBydz*dz;
2635 Bz=Bz0+dBzdx*dx+dBzdy*dy+
dBzdz*dz;
2645 Bx=Bx0+dBxdx*dx+dBxdy*dy+dBxdz*dz;
2646 By=By0+dBydx*dx+dBydy*dy+dBydz*dz;
2647 Bz=Bz0+dBzdx*dx+dBzdy*dy+
dBzdz*dz;
2687 double p=1./one_over_p;
2689 double denom=
sqrt(1.+tx*tx+ty*ty);
2690 double px=p*tx/denom;
2691 double py=p*ty/denom;
2695 double ds_over_p=ds*one_over_p;
2696 double factor=k_q*(0.25*ds_over_p);
2697 double Ax=factor*
Bx,Ay=factor*
By,Az=factor*
Bz;
2698 double Ax2=Ax*Ax,Ay2=Ay*Ay,Az2=Az*Az;
2699 double AxAy=Ax*Ay,AxAz=Ax*Az,AyAz=Ay*Az;
2700 double one_plus_Ax2=1.+Ax2;
2701 double scale=ds_over_p/(one_plus_Ax2+Ay2+Az2);
2704 double dx=scale*(px*one_plus_Ax2+py*(AxAy+Az)+pz*(AxAz-Ay));
2705 double dy=scale*(px*(AxAy-Az)+py*(1.+Ay2)+pz*(AyAz+Ax));
2706 double dz=scale*(px*(AxAz+Ay)+py*(AyAz-Ax)+pz*(1.+Az2));
2712 px+=k_q*(Bz*dy-By*dz);
2713 py+=k_q*(Bx*dz-Bz*dx);
2714 pz+=k_q*(By*dx-Bx*dy);
2717 if (fabs(dEdx)>
EPS){
2718 double one_over_p_sq=one_over_p*one_over_p;
2736 double ds_over_p=ds*one_over_p;
2737 double factor=k_q*(0.25*ds_over_p);
2738 double Ax=factor*
Bx,Ay=factor*
By,Az=factor*
Bz;
2739 double Ax2=Ax*Ax,Ay2=Ay*Ay,Az2=Az*Az;
2740 double AxAy=Ax*Ay,AxAz=Ax*Az,AyAz=Ay*Az;
2741 double one_plus_Ax2=1.+Ax2;
2742 double scale=ds_over_p/(one_plus_Ax2+Ay2+Az2);
2745 double dx=scale*(px*one_plus_Ax2+py*(AxAy+Az)+pz*(AxAz-Ay));
2746 double dy=scale*(px*(AxAy-Az)+py*(1.+Ay2)+pz*(AyAz+Ax));
2747 double dz=scale*(px*(AxAz+Ay)+py*(AyAz-Ax)+pz*(1.+Az2));
2748 xy.Set(xy.X()+dx,xy.Y()+dy);
2752 px+=k_q*(Bz*dy-By*dz);
2753 py+=k_q*(Bx*dz-Bz*dx);
2754 pz+=k_q*(By*dx-Bx*dy);
2755 pt=
sqrt(px*px+py*py);
2759 if (fabs(dEdx)>
EPS){
2760 double one_over_p_sq=one_over_p*one_over_p;
2777 if (fabs(ds)<
EPS)
return NOERROR;
2792 double sinl=
sin(lambda);
2793 double cosl=cos(lambda);
2794 double cosl2=cosl*cosl;
2795 double cosl3=cosl*cosl2;
2796 double one_over_cosl=1./cosl;
2797 double pt=fabs(1./q_over_pt);
2803 double kds=
qBr2p*ds;
2804 double kq_ds_over_pt=kds*q_over_pt;
2805 double By_cosphi_minus_Bx_sinphi=
By*cosphi-
Bx*sinphi;
2806 double By_sinphi_plus_Bx_cosphi=
By*sinphi+
Bx*cosphi;
2821 =-kq_ds_over_pt*q_over_pt*sinl*By_sinphi_plus_Bx_cosphi;
2823 +=2.*kq_ds_over_pt*sinl*By_cosphi_minus_Bx_sinphi;
2825 =kq_ds_over_pt*q_over_pt*cosl3*By_cosphi_minus_Bx_sinphi;
2827 double p=pt*one_over_cosl;
2829 double m2_over_p2=
mass2/p_sq;
2831 double dE_over_E=dEdx*ds/E;
2837 =kq_ds_over_pt*q_over_pt*sinl*(
dBydz*cosphi-
dBxdz*sinphi);
2844 double qrc_old=qpt/(
qBr2p*B);
2845 double qrc_plus_D=D+qrc_old;
2848 double dx_sinphi_minus_dy_cosphi=dx*sinphi-dy*cosphi;
2849 double rc=
sqrt(dxy.Mod2()
2850 +2.*qrc_plus_D*(dx_sinphi_minus_dy_cosphi)
2851 +qrc_plus_D*qrc_plus_D);
2854 J(
state_D,
state_D)=q_over_rc*(dx_sinphi_minus_dy_cosphi+qrc_plus_D);
2873 if (fabs(ds)<
EPS)
return NOERROR;
2902 double qrc_old=qpt/(
qBr2p*B);
2903 double qrc_plus_D=D+qrc_old;
2906 double dx_sinphi_minus_dy_cosphi=dx*sinphi-dy*cosphi;
2907 double rc=
sqrt(dxy.Mod2()
2908 +2.*qrc_plus_D*(dx_sinphi_minus_dy_cosphi)
2909 +qrc_plus_D*qrc_plus_D);
2910 double q_over_rc=q/rc;
2912 J(
state_D,
state_D)=q_over_rc*(dx_sinphi_minus_dy_cosphi+qrc_plus_D);
2922 double chi2c_factor,
2923 double chi2a_factor,
2931 double tanl2=tanl*tanl;
2932 double one_plus_tanl2=1.+tanl2;
2934 double my_ds=fabs(ds);
2935 double my_ds_2=0.5*my_ds;
2945 double sign_D=(Sc(
state_D)>0.0)?1.:-1.;
2950 double one_over_p_sq=one_over_pt*one_over_pt/one_plus_tanl2;
2951 double one_over_beta2=1.+
mass2*one_over_p_sq;
2952 double chi2c_p_sq=chi2c_factor*my_ds*one_over_beta2;
2953 double chi2a_p_sq=chi2a_factor*(1.+chi2a_corr*one_over_beta2);
2957 double one_plus_nu=1.+nu;
2960 *(one_plus_nu/nu*log(one_plus_nu)-1.);
2971 double chi2c_factor,
2972 double chi2a_factor,
2982 double my_ds=fabs(ds);
2983 double my_ds_2=0.5*my_ds;
2986 double one_plus_tx2=1.+tx2;
2987 double one_plus_ty2=1.+ty2;
2988 double tsquare=tx2+ty2;
2989 double one_plus_tsquare=1.+tsquare;
2998 = my_ds_2*
sqrt(one_plus_tsquare*one_plus_ty2);
3000 = my_ds_2*
sqrt(one_plus_tsquare*one_plus_tx2);
3002 double one_over_beta2=1.+one_over_p_sq*
mass2;
3003 double chi2c_p_sq=chi2c_factor*my_ds*one_over_beta2;
3004 double chi2a_p_sq=chi2a_factor*(1.+chi2a_corr*one_over_beta2);
3008 double one_plus_nu=1.+nu;
3010 *(one_plus_nu/nu*log(one_plus_nu)-1.);
3025 double rho_Z_over_A,
double LnI,
double Z){
3026 if (rho_Z_over_A<=0.)
return 0.;
3029 double p=fabs(1./q_over_p);
3030 double betagamma=p/
MASS;
3031 double betagamma2=betagamma*betagamma;
3032 double gamma2=1.+betagamma2;
3033 double beta2=betagamma2/gamma2;
3034 if (beta2<
EPS) beta2=
EPS;
3042 double two_Me_betagamma_sq=
two_m_e*betagamma2;
3046 dEdx=K_rho_Z_over_A/beta2*(-log(two_Me_betagamma_sq*Tmax)
3047 +2.*(LnI + beta2)+delta);
3051 double tau=
sqrt(gamma2)-1.;
3052 double tau_sq=tau*tau;
3053 double tau_plus_1=tau+1.;
3054 double tau_plus_2=tau+2.;
3055 double tau_plus_2_sq=tau_plus_2*tau_plus_2;
3058 f=1.-beta2+(0.125*tau_sq-(2.*tau+1.)*log(2.))/(tau_plus_1*tau_plus_1);
3061 f=2.*log(2.)-(beta2/12.)*(23.+14./tau_plus_2+10./tau_plus_2_sq
3062 +4./(tau_plus_2*tau_plus_2_sq));
3067 =-K_rho_Z_over_A/beta2*(log(0.5*tau_sq*tau_plus_2*
m_e_sq)-LnI+f-delta);
3074 double epsilon2=epsilon*epsilon;
3075 double f_Z=a2*(1./(1.+a2)+0.20206-0.0369*a2+0.0083*a4-0.002*a2*a4);
3079 double dEdx_rad=-K_rho_Z_over_A*tau_plus_1*(2.*ALPHA/M_PI)*(Z+1.)
3080 *((log(183.*pow(Z,-1./3.))-f_Z)
3081 *(1.-epsilon-(1./3.)*(epsilon2-epsilon*epsilon2))
3082 +1./18.*(1.-epsilon2));
3087 dEdx=dEdx_coll+dEdx_rad;
3100 double one_over_beta2,
3101 double K_rho_Z_over_A){
3102 if (K_rho_Z_over_A<=0.)
return 0.;
3104 double betagamma2=1./(one_over_beta2-1.);
3105 double gamma2=betagamma2*one_over_beta2;
3106 double two_Me_betagamma_sq=
two_m_e*betagamma2;
3107 double Tmax=two_Me_betagamma_sq/(1.+2.*
sqrt(gamma2)*m_ratio+
m_ratio_sq);
3108 double var=K_rho_Z_over_A*one_over_beta2*fabs(ds)*Tmax*(1.-0.5/one_over_beta2);
3114 if (
z_<
Z_MIN)
return VALUE_OUT_OF_RANGE;
3117 vector<const DCDCTrackHit*>forward_cdc_used_in_fit;
3134 double dpt_over_pt=0.1;
3153 double sig_lambda=0.02;
3155 =dpt_over_pt*dpt_over_pt+
tanl_*
tanl_*sig_lambda*sig_lambda;
3162 double p_mag=pvec.Mag();
3183 sperp=(my_z-z0)/
tanl_;
3186 double one_over_2k=1./twokappa;
3188 for (
unsigned int i=
my_cdchits.size()-1;i>1;i--){
3191 double tworc=2.*fabs(one_over_2k);
3192 double ratio=(
my_cdchits[i]->hit->wire->origin
3194 sperp=(ratio<1.)?tworc*asin(ratio):tworc*M_PI_2;
3195 if (sperp<25.) sperp=25.;
3200 double twoks=twokappa*sperp;
3201 double cosphi=cos(phi0);
3202 double sinphi=
sin(phi0);
3203 double sin2ks=
sin(twoks);
3204 double cos2ks=cos(twoks);
3205 double one_minus_cos2ks=1.-cos2ks;
3206 double myx=x0+one_over_2k*(cosphi*sin2ks-sinphi*one_minus_cos2ks);
3207 double myy=y0+one_over_2k*(sinphi*sin2ks+cosphi*one_minus_cos2ks);
3208 double mypx=px*cos2ks-py*sin2ks;
3209 double mypy=py*cos2ks+px*sin2ks;
3210 double myphi=atan2(mypy,mypx);
3221 if (!isfinite(x0) || !isfinite(y0) || !isfinite(q_over_p0)){
3223 return UNRECOVERABLE_ERROR;
3227 double tx0=
tx_=px/pz;
3228 double ty0=
ty_=py/pz;
3232 double fdc_prob=0.,fdc_chisq=-1.;
3233 unsigned int fdc_ndf=0;
3236 (isfinite(tx0) && isfinite(ty0))
3239 _DBG_ <<
"Using forward parameterization." <<endl;
3252 double temp=sig_lambda*one_plus_tsquare;
3269 return UNRECOVERABLE_ERROR;
3283 map<DetectorSystem_t,vector<Extrapolation_t> >saved_extrapolations;
3291 double forward_prob=0.;
3292 double chisq_forward=-1.;
3293 unsigned int ndof_forward=0;
3297 vector< vector <double> > fcov_save;
3298 vector<pull_t>pulls_save;
3299 pulls_save.assign(
pulls.begin(),
pulls.end());
3301 fcov_save.assign(
fcov.begin(),
fcov.end());
3310 _DBG_ <<
"Using forward parameterization." <<endl;
3325 double temp=sig_lambda*one_plus_tsquare;
3345 if (fdc_ndf==0 || forward_prob>fdc_prob){
3358 if (!fcov_save.empty()){
3359 fcov.assign(fcov_save.begin(),fcov_save.end());
3361 if (!saved_extrapolations.empty()){
3365 pulls.assign(pulls_save.begin(),pulls_save.end());
3380 fcov_save.assign(
fcov.begin(),
fcov.end());
3381 pulls_save.assign(
pulls.begin(),
pulls.end());
3396 _DBG_ <<
"Using central parameterization." <<endl;
3418 double one_plus_tanl2=1.+tanl2;
3420 *sig_lambda*sig_lambda;
3442 double central_prob=TMath::Prob(
chisq_,
ndf_);
3444 if (central_prob<forward_prob){
3453 fcov.assign(fcov_save.begin(),fcov_save.end());
3454 pulls.assign(pulls_save.begin(),pulls_save.end());
3456 if (!saved_extrapolations.empty()){
3460 cdchits_used_in_fit.assign(forward_cdc_used_in_fit.begin(),forward_cdc_used_in_fit.end());
3480 if (!saved_extrapolations.empty()){
3484 fcov.assign(fcov_save.begin(),fcov_save.end());
3485 pulls.assign(pulls_save.begin(),pulls_save.end());
3486 cdchits_used_in_fit.assign(forward_cdc_used_in_fit.begin(),forward_cdc_used_in_fit.end());
3493 else return UNRECOVERABLE_ERROR;
3496 if (
ndf_==0)
return UNRECOVERABLE_ERROR;
3502 #define CGOLD 0.3819660
3503 #define ZEPS 1.0e-10
3504 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
3505 #define SIGN(a,b) ((b)>=0.0?fabs(a):-fabs(a))
3510 const double z0wire,
3521 double a=(ax<cx?ax:cx);
3522 double b=(ax>cx?ax:cx);
3523 double x=bx,w=bx,v=bx;
3534 Step(pos,x,Sc,dedx);
3542 return VALUE_OUT_OF_RANGE;
3550 double fw=(pos-wirepos).Mod2();
3554 for (
unsigned int iter=1;iter<=
ITMAX;iter++){
3555 double xm=0.5*(a+b);
3557 double tol2=2.0*tol1;
3559 if (fabs(x-xm)<=(tol2-0.5*(b-a))){
3561 unsigned int iter2=0;
3574 return VALUE_OUT_OF_RANGE;
3578 Step(pos,u_old-u,Sc,dedx);
3587 unsigned int iter2=0;
3601 return VALUE_OUT_OF_RANGE;
3605 Step(pos,u_old-u,Sc,dedx);
3618 double x_minus_w=x-w;
3619 double x_minus_v=x-v;
3620 double r=x_minus_w*(fx-fv);
3621 double q=x_minus_v*(fx-fw);
3622 double p=x_minus_v*q-x_minus_w*r;
3628 if (fabs(p)>=fabs(0.5*q*etemp) || p<=q*(a-x) || p>=q*(b-x))
3630 d=
CGOLD*(e=(x>=xm?a-x:b-
x));
3635 if (u-a<tol2 || b-u <tol2)
3639 d=
CGOLD*(e=(x>=xm?a-x:b-
x));
3641 u=(fabs(d)>=tol1 ? x+d: x+
SIGN(tol1,d));
3650 return VALUE_OUT_OF_RANGE;
3654 Step(pos,u_old-u,Sc,dedx);
3658 wirepos+=(Sc(
state_z)-z0wire)*dir;
3659 double fu=(pos-wirepos).Mod2();
3664 if (u>=x) a=
x;
else b=
x;
3669 if (u<x) a=
u;
else b=
u;
3670 if (fu<=fw || w==x){
3676 else if (fu<=fv || v==x || v==w){
3691 const double z0wire,
3702 double a=(ax<cx?ax:cx);
3703 double b=(ax>cx?ax:cx);
3704 double x=bx,w=bx,v=bx;
3711 Step(z,z_new,dedx,S);
3719 return VALUE_OUT_OF_RANGE;
3722 double dz0wire=z-z0wire;
3723 DVector2 wirepos=origin+(dz0wire+
x)*dir;
3728 double fw=(pos-wirepos).Mod2();
3733 for (
unsigned int iter=1;iter<=
ITMAX;iter++){
3734 double xm=0.5*(a+b);
3736 double tol2=2.0*tol1;
3737 if (fabs(x-xm)<=(tol2-0.5*(b-a))){
3748 return VALUE_OUT_OF_RANGE;
3752 return VALUE_OUT_OF_RANGE;
3761 double x_minus_w=x-w;
3762 double x_minus_v=x-v;
3763 double r=x_minus_w*(fx-fv);
3764 double q=x_minus_v*(fx-fw);
3765 double p=x_minus_v*q-x_minus_w*r;
3771 if (fabs(p)>=fabs(0.5*q*etemp) || p<=q*(a-x) || p>=q*(b-x))
3773 d=
CGOLD*(e=(x>=xm?a-x:b-
x));
3778 if (
u-a<tol2 || b-
u <tol2)
3782 d=
CGOLD*(e=(x>=xm?a-x:b-
x));
3784 u=(fabs(d)>=tol1 ? x+d: x+
SIGN(tol1,d));
3796 return VALUE_OUT_OF_RANGE;
3799 Step(z_old,z_new,dedx,S);
3803 wirepos+=(dz0wire+
u)*dir;
3805 double fu=(pos-wirepos).Mod2();
3810 if (u>=x) a=
x;
else b=
x;
3815 if (u<x) a=
u;
else b=
u;
3816 if (fu<=fw || w==x){
3822 else if (fu<=fv || v==x || v==w){
3838 unsigned int &my_ndf
3895 double doca2,old_doca2=(xy-wirexy).Mod2();
3901 bool more_measurements=
true;
3907 unsigned int k_minus_1=k-1;
3948 _DBG_<<
"Went outside of tracking volume at z="<<Sc(
state_z)
3949 <<
" r="<<xy.Mod()<<endl;
3972 wirexy+=(Sc(
state_z)-z0w)*dir;
3975 doca2=(xy-wirexy).Mod2();
3981 _DBG_ <<
" Good Hit Ring " <<
my_cdchits[cdc_index]->hit->wire->ring <<
" Straw " <<
my_cdchits[cdc_index]->hit->wire->straw << endl;
3982 _DBG_ <<
" doca " <<
sqrt(doca2) << endl;
3999 double qrc_plus_D=D+qrc_old;
4004 double cosstereo=
my_cdchits[cdc_index]->cosstereo;
4031 double rc=
sqrt(dxy1.Mod2()
4032 +2.*qrc_plus_D*(dxy1.X()*sinphi-dxy1.Y()*cosphi)
4033 +qrc_plus_D*qrc_plus_D);
4038 wirexy+=(Sc(
state_z)-z0w)*dir;
4042 double doca=diff.Mod()+
EPS;
4043 double prediction=doca*cosstereo;
4046 double measurement=0.39,tdrift=0.,tcorr=0.,
dDdt0=0.;
4051 int ring_index=mywire->
ring-1;
4052 int straw_index=mywire->
straw-1;
4054 double phi_d=diff.Phi();
4056 =
max_sag[ring_index][straw_index]*(1.-dz*dz/5625.)
4058 double dphi=phi_d-mywire->
origin.Phi();
4059 while (dphi>M_PI) dphi-=2*M_PI;
4060 while (dphi<-M_PI) dphi+=2*M_PI;
4061 if (mywire->
origin.Y()<0) dphi*=-1.;
4074 dDdt0=(d_shifted-measurement)/dt;
4086 double cosstereo_over_doca=cosstereo/doca;
4087 H_T(
state_D)=(dy*cosphi-dx*sinphi)*cosstereo_over_doca;
4089 =-Sc(
state_D)*cosstereo_over_doca*(dx*cosphi+dy*sinphi);
4090 H_T(
state_z)=-cosstereo_over_doca*(dx*ux+dy*uy);
4100 double Vproj=H*Cc*H_T;
4102 double dm=measurement-prediction;
4106 _DBG_ << k <<
" "<<
central_traj.size()-1<<
" Negative variance??? " << V <<
" " << H*(Cc*H_T) <<endl;
4118 double chi2check=dm*dm*InvV;
4119 if (
DEBUG_LEVEL>9)
_DBG_ <<
" Prediction " << prediction <<
" Measurement " << measurement <<
" Chi2 " << chi2check << endl;
4120 if (chi2check<chi2cut)
4129 Ctest=Cc.
SubSym(K*(H*Cc));
4134 if (!Ctest.IsPosDef()){
4148 if (
DEBUG_LEVEL>9)
_DBG_ <<
" Marked Trajectory central_traj[k].h_id=cdc_index+1 (k cdc_index)" << k <<
" " << cdc_index << endl;
4151 double scale=(skip_ring)?1.:(1.-H*K);
4162 chisq+=scale*dm*dm/V;
4167 <<
"ring " <<
my_cdchits[cdc_index]->hit->wire->ring
4168 <<
" t " <<
my_cdchits[cdc_index]->hit->tdrift
4169 <<
" Dm-Dpred " << dm
4170 <<
" chi2 " << (1.-H*K)*dm*dm/V
4206 more_measurements=
false;
4210 wirexy=origin+(Sc(
state_z)-z0w)*dir;
4214 doca2=(xy-wirexy).Mod2();
4234 _DBG_ <<
"Unphysical momentum: P = " << p <<endl;
4248 unsigned int num_good=0;
4249 for (
unsigned int j=0;j<num_cdc;j++){
4262 double cdc_anneal_factor,
4266 unsigned int &numdof){
4314 unsigned int max_num_fdc_used_in_fit=num_fdc_hits;
4316 unsigned int cdc_index=0;
4317 if (num_cdc_hits>0) cdc_index=num_cdc_hits-1;
4318 bool more_cdc_measurements=(num_cdc_hits>0);
4319 double old_doca2=1e6;
4326 if (more_cdc_measurements){
4336 jout <<
"Entering DTrackFitterKalmanSIMD::KalmanForward ================================" << endl;
4337 jout <<
" We have the following starting parameters for our fit. S = State vector, C = Covariance matrix" << endl;
4340 jout << setprecision(6);
4341 jout <<
" There are " << num_cdc <<
" CDC Updates and " << num_fdc <<
" FDC Updates on this track" << endl;
4342 jout <<
" There are " << num_cdc_hits <<
" CDC Hits and " << num_fdc_hits <<
" FDC Hits on this track" << endl;
4343 jout <<
" With NUM_FDC_SIGMA_CUT = " << NUM_FDC_SIGMA_CUT <<
" and NUM_CDC_SIGMA_CUT = " << NUM_CDC_SIGMA_CUT << endl;
4344 jout <<
" fdc_anneal_factor = " << fdc_anneal_factor <<
" cdc_anneal_factor = " << cdc_anneal_factor << endl;
4345 jout <<
" yields fdc_chi2cut = " << fdc_chi2cut <<
" cdc_chi2cut = " << cdc_chi2cut << endl;
4347 jout <<
" S0_ which is the state vector of the reference trajectory at the break point step:" << endl;
4349 jout <<
" ===== Beginning pass over reference trajectory ======== " << endl;
4353 unsigned int k_minus_1=k-1;
4399 jout <<
" At reference trajectory index " << k <<
" at z=" << z << endl;
4400 jout <<
" The State vector from the reference trajectory" << endl;
4402 jout <<
" The Jacobian matrix " << endl;
4404 jout <<
" The Q matrix "<< endl;
4406 jout <<
" The updated State vector S=S0+J*(S-S0_)" << endl;
4408 jout <<
" The updated Covariance matrix C=J*(C*J_T)+Q;" << endl;
4413 if (num_fdc_hits>0){
4432 double vpred_uncorrected=x*sina+y*cosa;
4435 double upred=x*cosa-y*sina;
4438 double tu=tx*cosa-ty*sina;
4439 double alpha=atan(tu);
4440 double cosalpha=cos(alpha);
4441 double cosalpha2=cosalpha*cosalpha;
4442 double sinalpha=
sin(alpha);
4446 double doca=du*cosalpha;
4451 double nz_sinalpha_plus_nr_cosalpha=nz*sinalpha+nr*cosalpha;
4454 V(1,1)=
my_fdchits[id]->vvar*fdc_anneal_factor;
4457 double tv=tx*sina+ty*cosa;
4458 Mdiff(1)=v-(vpred_uncorrected+doca*(nz_sinalpha_plus_nr_cosalpha
4469 Mdiff(0)=drift-doca;
4478 double temp2=nz_sinalpha_plus_nr_cosalpha-tv*sinalpha;
4479 H_T(
state_x,1)=sina+cosa*cosalpha*temp2;
4480 H_T(
state_y,1)=cosa-sina*cosalpha*temp2;
4482 double cos2_minus_sin2=cosalpha2-sinalpha*sinalpha;
4483 double fac=nz*cos2_minus_sin2-2.*nr*cosalpha*sinalpha;
4484 double doca_cosalpha=doca*cosalpha;
4485 double temp=doca_cosalpha*fac;
4487 -doca_cosalpha*(tu*sina+tv*cosa*cos2_minus_sin2)
4490 -doca_cosalpha*(tu*cosa-tv*sina*cos2_minus_sin2)
4494 H_T(
state_y,0)=-sina*cosalpha;
4495 double one_plus_tu2=1.+tu*tu;
4496 double factor=du*tu/
sqrt(one_plus_tu2)/one_plus_tu2;
4509 vector<DMatrix5x2> Klist;
4510 vector<DMatrix2x1> Mlist;
4511 vector<DMatrix2x5> Hlist;
4512 vector<DMatrix5x2> HTlist;
4513 vector<DMatrix2x2> Vlist;
4514 vector<double>probs;
4515 vector<unsigned int>used_ids;
4523 double chi2_hit=Vtemp.
Chi2(Mdiff);
4524 double prob_hit=exp(-0.5*chi2_hit)
4531 probs.push_back(prob_hit);
4534 HTlist.push_back(H_T);
4535 Mlist.push_back(Mdiff);
4536 Klist.push_back(C*H_T*InvV);
4538 used_ids.push_back(
id);
4544 unsigned int my_id=
id-m;
4557 Mdiff(1)=v-(vpred_uncorrected+doca*(nz_sinalpha_plus_nr_cosalpha
4567 Mdiff(0)=drift-doca;
4575 doca_cosalpha=doca*cosalpha;
4576 temp=doca_cosalpha*fac;
4578 -doca_cosalpha*(tu*sina+tv*cosa*cos2_minus_sin2)
4581 -doca_cosalpha*(tu*cosa-tv*sina*cos2_minus_sin2)
4583 factor=du*tu/
sqrt(one_plus_tu2)/one_plus_tu2;
4599 double chi2_hit=Vtemp.
Chi2(Mdiff);
4600 double prob_hit=exp(-0.5*chi2_hit)
4604 if(chi2_hit<fdc_chi2cut){
4605 probs.push_back(prob_hit);
4606 Mlist.push_back(Mdiff);
4609 HTlist.push_back(H_T);
4610 Klist.push_back(C*H_T*InvV);
4612 used_ids.push_back(my_id);
4618 double prob_tot=1
e-100;
4619 for (
unsigned int m=0;m<probs.size();m++){
4626 if (skip_plane==
false){
4629 for (
unsigned int m=0;m<Klist.size();m++){
4630 double my_prob=probs[m]/prob_tot;
4631 S+=my_prob*(Klist[m]*Mlist[m]);
4632 sum-=my_prob*(Klist[m]*Hlist[m]);
4633 sum2+=(my_prob*my_prob)*(Klist[m]*Vlist[m]*
Transpose(Klist[m]));
4637 R=Mlist[m]-HK*Mlist[m];
4638 RC=Vlist[m]-HK*Vlist[m];
4639 chisq+=my_prob*RC.
Chi2(R);
4641 unsigned int my_id=used_ids[m];
4645 jout <<
" Adjusting state vector for FDC hit " << m <<
" with prob " << my_prob <<
" S:" << endl;
4655 for (
unsigned int m=0;m<Hlist.size();m++){
4656 unsigned int my_id=used_ids[m];
4670 if (skip_plane==
false){
4675 if (
DEBUG_LEVEL > 25) jout <<
" == There is a single FDC hit on this plane" << endl;
4682 double chi2_hit=Vtemp.
Chi2(Mdiff);
4683 if (chi2_hit<fdc_chi2cut){
4688 if (skip_plane==
false){
4697 jout <<
"S Update: " << endl; S.
Print();
4698 jout <<
"C Uodate: " << endl; C.Print();
4711 if (skip_plane==
false){
4727 printf(
"hit %d p %5.2f t %f dm %5.2f sig %f chi2 %5.2f z %5.2f\n",
4746 else if (more_cdc_measurements ){
4749 wirepos+=(z-z0w)*dir;
4754 double doca2=dx*dx+dy*dy;
4757 if (doca2>old_doca2 ){
4767 double tanl=1./
sqrt(tx*tx+ty*ty);
4768 double sinl=
sin(atan(tanl));
4776 double denom=my_ux*my_ux+my_uy*my_uy;
4788 bool do_brent=
false;
4796 double two_step=step1+step2;
4799 *two_step/sinl)<0.05
4803 dz=-((
S(
state_x)-origin.X()-ux*dzw)*my_ux
4804 +(
S(
state_y)-origin.Y()-uy*dzw)*my_uy)/denom;
4806 if (fabs(dz)>two_step || dz<0){
4827 if (
BrentForward(z,dedx,z0w,origin,dir,S,dz)!=NOERROR){
4835 num_steps=int(fabs(dz/my_dz));
4836 dz3=dz-num_steps*my_dz;
4839 for (
int m=0;m<num_steps;m++){
4851 Step(my_z,newz,dedx,S0);
4867 Step(my_z,newz,dedx,S0);
4881 Step(z,newz,dedx,S0);
4887 wirepos+=(newz-z0w)*dir;
4889 double xw=wirepos.X();
4890 double yw=wirepos.Y();
4893 if (do_brent==
false){
4906 double cosstereo=
my_cdchits[cdc_index]->cosstereo;
4907 double d=
sqrt(dx*dx+dy*dy)*cosstereo+
EPS;
4910 double cosstereo2_over_d=cosstereo*cosstereo/d;
4911 Hc_T(
state_x)=dx*cosstereo2_over_d;
4913 Hc_T(
state_y)=dy*cosstereo2_over_d;
4915 if (do_brent==
false){
4931 double dm=0.39,tdrift=0.,tcorr=0.,
dDdt0=0.;
4937 int ring_index=mywire->
ring-1;
4938 int straw_index=mywire->
straw-1;
4940 double phi_d=atan2(dy,dx);
4942 =
max_sag[ring_index][straw_index]*(1.-dz*dz/5625.)
4944 double dphi=phi_d-mywire->
origin.Phi();
4945 while (dphi>M_PI) dphi-=2*M_PI;
4946 while (dphi<-M_PI) dphi+=2*M_PI;
4947 if (mywire->
origin.Y()<0) dphi*=-1.;
4953 Vc*=cdc_anneal_factor;
4960 if (tdrift < 0.) d_shifted = dm;
4962 dDdt0=(d_shifted-dm)/dt;
4965 if (max_num_fdc_used_in_fit>4)
4980 double Vproj=Hc*C*Hc_T;
4981 double InvV1=1./(Vc+Vproj);
4984 _DBG_ <<
"Negative variance???" << endl;
4989 double chi2_hit=res*res*InvV1;
4990 if (chi2_hit<cdc_chi2cut){
4997 Ctest=C.
SubSym(Kc*(Hc*C));
5000 if (!Ctest.IsPosDef()){
5012 jout <<
" == Adding CDC Hit in Ring " <<
my_cdchits[cdc_index]->hit->wire->ring << endl;
5013 jout <<
" New S: " << endl; S.
Print();
5014 jout <<
" New C: " << endl; C.Print();
5022 double scale=(skip_ring)?1.:(1.-Hc*Kc);
5036 chisq+=scale*res*res/Vc;
5041 jout <<
"Ring " <<
my_cdchits[cdc_index]->hit->wire->ring
5042 <<
" Straw " <<
my_cdchits[cdc_index]->hit->wire->straw
5043 <<
" Pred " << d <<
" Meas " << dm
5044 <<
" Sigma meas " <<
sqrt(Vc)
5046 <<
" Chi2 " << (1.-Hc*Kc)*res*res/Vc << endl;
5064 Step(newz,z,dedx,S);
5067 Step(newz,z,dedx,S0);
5071 for (
int m=0;m<num_steps;m++){
5080 Step(my_z,z,dedx,S);
5083 Step(my_z,z,dedx,S0);
5095 Step(my_z,z,dedx,S);
5098 Step(my_z,z,dedx,S0);
5112 else more_cdc_measurements=
false;
5115 wirepos=origin+(z-z0w)*dir;
5133 double my_dz=
mStepSizeZ*(dz_to_target>0?1.:-1.);
5134 int num_steps=int(fabs(dz_to_target/my_dz));
5136 for (
int k=0;k<num_steps;k++){
5137 double newz=z_+my_dz;
5162 double d=
sqrt(dx*dx+dy*dy);
5165 double one_over_d=1./d;
5178 double InvV1=1./(Vc+Hc*(C*Hc_T));
5182 _DBG_ <<
"Negative variance???" << endl;
5190 double res=0.1466666667-d;
5197 chisq+=(1.-Hc*Kc)*res*res/Vc;
5207 unsigned int new_index=(3*num_fdc)/4;
5211 unsigned int new_index=(3*num_fdc)/4;
5230 cout <<
"Position after forward filter: " <<
x_ <<
", " <<
y_ <<
", " <<
z_ <<endl;
5240 unsigned int new_index=(3*num_fdc)/4;
5262 unsigned int new_index=3*(num_fdc)/4;
5271 unsigned int num_good=0;
5272 unsigned int num_hits=num_cdc+max_num_fdc_used_in_fit;
5273 for (
unsigned int j=0;j<num_cdc;j++){
5276 for (
unsigned int j=0;j<num_fdc;j++){
5282 unsigned int new_index=(3*num_fdc)/4;
5286 unsigned int new_index=(3*num_fdc)/4;
5309 unsigned int &numdof){
5347 if (cdc_index<num_cdc-1){
5348 num_cdc=cdc_index+1;
5356 DVector2 wirepos=origin+(z-z0w)*dir;
5357 bool more_measurements=
true;
5362 double doca2=0.,old_doca2=dx*dx+dy*dy;
5367 unsigned int k_minus_1=k-1;
5397 _DBG_<<
"Went outside of tracking volume at x=" <<
S(
state_x)
5398 <<
" y=" <<
S(
state_y) <<
" z="<<z<<endl;
5422 wirepos+=(z-z0w)*dir;
5438 double dz=0.,newz=z;
5453 double tanl=1./
sqrt(tx*tx+ty*ty);
5454 double sinl=
sin(atan(tanl));
5462 double denom=my_ux*my_ux+my_uy*my_uy;
5472 bool do_brent=
false;
5474 double two_step=step1+step2;
5478 *two_step/sinl)<0.05
5481 dz=-((
S(
state_x)-origin.X()-ux*dzw)*my_ux
5482 +(
S(
state_y)-origin.Y()-uy*dzw)*my_uy)/denom;
5484 if (fabs(dz)>two_step || dz<0.0){
5515 num_steps=int(fabs(dz/my_dz));
5516 dz3=dz-num_steps*my_dz;
5519 for (
int m=0;m<num_steps;m++){
5531 Step(my_z,newz,dedx,S0);
5547 Step(my_z,newz,dedx,S0);
5560 Step(z,newz,dedx,S0);
5566 wirepos+=(newz-z0w)*dir;
5568 double xw=wirepos.X();
5569 double yw=wirepos.Y();
5572 if (do_brent==
false){
5585 double cosstereo=
my_cdchits[cdc_index]->cosstereo;
5586 double d=
sqrt(dx*dx+dy*dy)*cosstereo+
EPS;
5591 double cosstereo2_over_d=cosstereo*cosstereo/d;
5592 H_T(
state_x)=dx*cosstereo2_over_d;
5594 H_T(
state_y)=dy*cosstereo2_over_d;
5596 if (do_brent==
false){
5612 double dm=0.39,tdrift=0.,tcorr=0.,
dDdt0=0.;
5617 int ring_index=mywire->
ring-1;
5618 int straw_index=mywire->
straw-1;
5620 double phi_d=atan2(dy,dx);
5622 =
max_sag[ring_index][straw_index]*(1.-dz*dz/5625.)
5624 double dphi=phi_d-mywire->
origin.Phi();
5625 while (dphi>M_PI) dphi-=2*M_PI;
5626 while (dphi<-M_PI) dphi+=2*M_PI;
5627 if (mywire->
origin.Y()<0) dphi*=-1.;
5639 if (tdrift < 0.) d_shifted = dm;
5641 dDdt0=(d_shifted-dm)/dt;
5652 double Vproj=H*C*H_T;
5653 double InvV=1./(V+Vproj);
5656 _DBG_ <<
"Negative variance???" << endl;
5662 double chi2check=res*res*InvV;
5663 if (chi2check<chi2cut/*&&dm>0.*/)
5674 if (!Ctest.IsPosDef()){
5690 double scale=(skip_ring)?1.:(1.-H*K);
5701 chisq+=scale*res*res/V;
5708 printf(
"Ring %d straw %d pred %f meas %f chi2 %f useBrent %i \n",
5711 d,dm,(1.-H*K)*res*res/V,do_brent);
5726 Step(newz,z,dedx,S);
5729 Step(newz,z,dedx,S0);
5733 for (
int m=0;m<num_steps;m++){
5742 Step(my_z,z,dedx,S);
5745 Step(my_z,z,dedx,S0);
5757 Step(my_z,z,dedx,S);
5760 Step(my_z,z,dedx,S0);
5768 if (cdc_index>0) cdc_index--;
5781 more_measurements=
false;
5786 wirepos+=(z-z0w)*dir;
5831 unsigned int num_good=0;
5832 for (
unsigned int j=0;j<num_cdc;j++){
5853 double z=
z_,newz=
z_;
5862 double rho_Z_over_A=0.,LnI=0.,K_rho_Z_over_A=0.,Z=0.;
5863 double chi2c_factor=0.,chi2a_factor=0.,chi2a_corr=0.;
5873 chi2c_factor,chi2a_factor,chi2a_corr,
5877 _DBG_ <<
"Material error in ExtrapolateToVertex at (x,y,z)=("
5878 << pos.X() <<
"," << pos.y()<<
","<< pos.z()<<
")"<< endl;
5879 return UNRECOVERABLE_ERROR;
5885 if (fabs(dEdx)>
EPS){
5907 Step(z,z-dz,dEdx,S1);
5915 return UNRECOVERABLE_ERROR;
5920 Step(z,z+dz,dEdx,S2);
5928 return UNRECOVERABLE_ERROR;
5933 if (r2plus>r2_old && r2minus>r2_old){
5942 return UNRECOVERABLE_ERROR;
5964 if (r2minus<r2_old && r2_old<r2plus){
5983 if (r2minus>r2_old && r2_old>r2plus){
6010 return UNRECOVERABLE_ERROR;
6018 double s_to_boundary=1.e6;
6022 chi2c_factor,chi2a_factor,chi2a_corr,
6025 _DBG_ <<
"Material error in ExtrapolateToVertex! " << endl;
6026 return UNRECOVERABLE_ERROR;
6031 chi2c_factor,chi2a_factor,chi2a_corr,
6034 _DBG_ <<
"Material error in ExtrapolateToVertex! " << endl;
6047 if (fabs(dEdx)>EPS){
6056 if (ds>s_to_boundary) ds=s_to_boundary;
6066 double one_over_beta2=1.+
mass2*q_over_p_sq;
6081 Step(z,newz,dEdx,S);
6088 double two_step=dz+dz_old;
6095 return UNRECOVERABLE_ERROR;
6135 double z=
z_,newz=
z_;
6146 double rho_Z_over_A=0.,LnI=0.,K_rho_Z_over_A=0.,Z=0.;
6147 double chi2c_factor=0.,chi2a_factor=0.,chi2a_corr=0.;
6159 return UNRECOVERABLE_ERROR;
6168 chi2c_factor,chi2a_factor,chi2a_corr,
6171 _DBG_ <<
"Material error in ExtrapolateToVertex! " << endl;
6182 if (fabs(dEdx)>EPS){
6187 double dz=-ds*dz_ds;
6193 Step(z,newz,dEdx,S);
6200 double two_step=dz+dz_old;
6205 return UNRECOVERABLE_ERROR;
6225 return VALUE_OUT_OF_RANGE;
6248 double r2=(xy-beam_pos).Mod2();
6260 Step(xy0,ds,S0,dedx);
6268 return UNRECOVERABLE_ERROR;
6271 r2=(xy0-beam_pos).Mod2();
6272 if (r2>r2_old) ds*=-1.;
6279 return UNRECOVERABLE_ERROR;
6295 return UNRECOVERABLE_ERROR;
6299 double rho_Z_over_A=0.,LnI=0.,K_rho_Z_over_A=0.,Z=0.;
6300 double chi2c_factor=0.,chi2a_factor=0.,chi2a_corr=0.;
6303 chi2c_factor,chi2a_factor,chi2a_corr,
6306 _DBG_ <<
"Material error in ExtrapolateToVertex! " << endl;
6313 dedx=-
GetdEdx(q_over_p,K_rho_Z_over_A,rho_Z_over_A,LnI,Z);
6316 double sign=(ds>0.0)?1.:-1.;
6317 if (fabs(dedx)>
EPS){
6327 double q_over_p_sq=q_over_p*q_over_p;
6328 double one_over_beta2=1.+
mass2*q_over_p*q_over_p;
6341 Cc=(sign*Q).AddSym(Cc);
6344 r2=(xy-beam_pos).Mod2();
6352 return UNRECOVERABLE_ERROR;
6392 double r2=(xy-beam_pos).Mod2();
6402 Step(xy1,ds,S1,dedx);
6404 r2=(xy1-beam_pos).Mod2();
6405 if (r2>r2_old) ds*=-1.;
6412 double rho_Z_over_A=0.,LnI=0.,K_rho_Z_over_A=0.,Z=0;
6413 double chi2c_factor=0.,chi2a_factor=0.,chi2a_corr=0.;
6416 chi2c_factor,chi2a_factor,chi2a_corr,
6419 _DBG_ <<
"Material error in ExtrapolateToVertex! " << endl;
6426 dedx=
GetdEdx(q_over_p,K_rho_Z_over_A,rho_Z_over_A,LnI,Z);
6429 double sign=(ds>0.0)?1.:-1.;
6430 if (fabs(dedx)>
EPS){
6437 Step(xy,ds,Sc,dedx);
6440 r2=(xy-beam_pos).Mod2();
6461 return VALUE_OUT_OF_RANGE;
6474 C7x7->ResizeTo(7, 7);
6479 double tanl2=tanl*tanl;
6480 double lambda=atan(tanl);
6481 double sinl=
sin(lambda);
6482 double sinl3=sinl*sinl*sinl;
6496 double frac=tanl2*tanl2;
6501 *C7x7=C.Similarity(J);
6512 C7x7->ResizeTo(7, 7);
6520 double cosphi=cos(
phi_);
6542 *C7x7=C.Similarity(J);
6568 unsigned int &numdof){
6569 if (
DEBUG_LEVEL>1)
_DBG_ <<
"Attempting to recover broken track ... " <<endl;
6572 double refit_chisq=-1.;
6573 unsigned int refit_ndf=0;
6582 vector<bool>old_cdc_used_status(num_hits);
6583 for (
unsigned int j=0;j<num_hits;j++){
6603 if (r2next>r2 && r2>r2_cdc)
break;
6608 if (
DEBUG_LEVEL>0)
_DBG_ <<
"Invalid reference trajectory in track recovery..." << endl;
6615 unsigned int refit_iter=0;
6638 refit_chisq,refit_ndf);
6661 for (
unsigned int k=0;k<num_hits;k++){
6687 unsigned int &numdof
6689 if (
DEBUG_LEVEL>1)
_DBG_ <<
"Attempting to recover broken track ... " <<endl;
6695 double refit_chisq=-1.;
6696 unsigned int refit_ndf=0;
6702 vector<bool>old_cdc_used_status(num_cdchits);
6703 vector<bool>old_fdc_used_status(num_fdchits);
6704 for (
unsigned int j=0;j<num_fdchits;j++){
6707 for (
unsigned int j=0;j<num_cdchits;j++){
6738 double r2=x1*x1+y1*y1;
6742 double r2next=x2*x2+y2*y2;
6743 double r2cdc=(origin+(
forward_traj[k].z-z0)*dir).Mod2();
6745 if (r2next>r2 && r2>r2cdc)
break;
6750 if (
DEBUG_LEVEL>0)
_DBG_ <<
"Invalid reference trajectory in track recovery..." << endl;
6759 unsigned int refit_iter=0;
6798 for (
unsigned int k=0;k<num_cdchits;k++){
6801 for (
unsigned int k=0;k<num_fdchits;k++){
6817 unsigned int &numdof){
6819 _DBG_ <<
"Attempting to recover broken track ... " <<endl;
6824 double refit_chisq=-1.;
6825 unsigned int refit_ndf=0;
6831 vector<int>old_cdc_used_status(num_cdchits);
6832 vector<int>old_fdc_used_status(num_fdchits);
6833 for (
unsigned int j=0;j<num_fdchits;j++){
6836 for (
unsigned int j=0;j<num_cdchits;j++){
6871 for (
unsigned int j=0;j<num_cdchits;j++){
6874 for (
unsigned int j=0;j<=break_id;j++){
6877 for (
unsigned int j=break_id+1;j<num_fdchits;j++){
6883 for (
unsigned int j=0;j<num_cdchits+break_id;j++){
6886 for (
unsigned int j=num_cdchits+break_id;j<num_cdchits;j++){
6889 for (
unsigned int j=0;j<num_fdchits;j++){
6902 refit_error=
KalmanForward(fdc_anneal,cdc_anneal,S1,C1,refit_chisq,refit_ndf);
6915 if (break_id>0) break_id--;
6929 for (
unsigned int k=0;k<num_cdchits;k++){
6932 for (
unsigned int k=0;k<num_fdchits;k++){
6945 unsigned int max_fdc_index=num_fdchits-1;
6949 vector<bool>last_fdc_used_in_fit(num_fdchits);
6950 vector<bool>last_cdc_used_in_fit(num_cdchits);
6951 vector<pull_t>forward_pulls;
6952 vector<pull_t>last_forward_pulls;
6971 double chisq=-1.,chisq_forward=-1.;
6972 unsigned int my_ndf=0;
6973 unsigned int last_ndf=1;
7003 if (ref_track_error!=NOERROR){
7009 for (
unsigned int j=0;j<num_cdchits;j++){
7015 bool did_second_refit=
false;
7016 error=
KalmanForward(fdc_anneal,cdc_anneal,S,C,chisq,my_ndf);
7018 _DBG_ <<
"Iter: " << iter+1 <<
" Chi2=" << chisq <<
" Ndf=" << my_ndf <<
" Error code: " << error << endl;
7030 unsigned int temp_ndf=my_ndf;
7031 double temp_chi2=
chisq;
7046 did_second_refit=
true;
7071 if (iter>
MIN_ITER && did_second_refit==
false){
7072 double new_reduced_chisq=chisq/my_ndf;
7073 double old_reduced_chisq=chisq_forward/last_ndf;
7074 double new_prob=TMath::Prob(chisq,my_ndf);
7075 double old_prob=TMath::Prob(chisq_forward,last_ndf);
7076 if (new_prob<old_prob
7077 || fabs(new_reduced_chisq-old_reduced_chisq)<
CHISQ_DELTA){
7082 chisq_forward=
chisq;
7095 forward_pulls.clear();
7098 pulls.assign(forward_pulls.begin(),forward_pulls.end());
7108 for (
unsigned int m=0;m<last_cdc_used_in_fit.size();m++){
7109 if (last_cdc_used_in_fit[m]){
7114 for (
unsigned int m=0;m<last_fdc_used_in_fit.size();m++){
7115 if (last_fdc_used_in_fit[m]){
7157 double cosl=cos(atan(
tanl_));
7162 double dx=
x_-beam_pos.X();
7163 double dy=
y_-beam_pos.Y();
7166 double cosphi=cos(
phi_);
7168 if ((dx>0.0 && sinphi>0.0) || (dy<0.0 && cosphi>0.0)
7169 || (dy>0.0 && cosphi<0.0) || (dx<0.0 && sinphi<0.0))
D_*=-1.;
7173 vector<double>dummy;
7174 for (
unsigned int i=0;i<5;i++){
7176 for(
unsigned int j=0;j<5;j++){
7177 dummy.push_back(Clast(i,j));
7179 fcov.push_back(dummy);
7188 unsigned int max_cdc_index=num_cdchits-1;
7204 vector<pull_t>cdc_pulls;
7205 vector<pull_t>last_cdc_pulls;
7206 vector<bool>last_cdc_used_in_fit;
7212 double chisq=-1.,chisq_forward=-1.;
7213 unsigned int my_ndf=0;
7214 unsigned int last_ndf=1;
7224 _DBG_ <<
"-------- iteration " << iter2+1 <<
" -----------" <<endl;
7251 if (ref_track_error!=NOERROR){
7257 for (
unsigned int j=0;j<num_cdchits;j++){
7263 bool did_second_refit=
false;
7272 unsigned int temp_ndf=my_ndf;
7273 double temp_chi2=
chisq;
7278 unsigned int new_index=(3*max_cdc_index)/4;
7288 unsigned int new_index=(max_cdc_index)/2;
7293 unsigned int new_index=(3*max_cdc_index)/4;
7300 anneal_factor=1000.;
7303 did_second_refit=
true;
7325 if (iter2==0)
return error;
7331 <<
" Prob: " << TMath::Prob(chisq,my_ndf)
7334 if (iter2>
MIN_ITER && did_second_refit==
false){
7335 double new_reduced_chisq=chisq/my_ndf;
7336 double old_reduced_chisq=chisq_forward/last_ndf;
7337 double new_prob=TMath::Prob(chisq,my_ndf);
7338 double old_prob=TMath::Prob(chisq_forward,last_ndf);
7339 if (new_prob<old_prob
7340 || fabs(new_reduced_chisq-old_reduced_chisq)<
CHISQ_DELTA)
break;
7343 chisq_forward=
chisq;
7359 pulls.assign(cdc_pulls.begin(),cdc_pulls.end());
7369 for (
unsigned int m=0;m<last_cdc_used_in_fit.size();m++){
7370 if (last_cdc_used_in_fit[m]){
7401 double cosl=cos(atan(
tanl_));
7406 double dx=
x_-beam_pos.X();
7407 double dy=
y_-beam_pos.Y();
7410 double cosphi=cos(
phi_);
7412 if ((dx>0.0 && sinphi>0.0) || (dy<0.0 && cosphi>0.0)
7413 || (dy>0.0 && cosphi<0.0) || (dx<0.0 && sinphi<0.0))
D_*=-1.;
7417 vector<double>dummy;
7420 for (
unsigned int i=0;i<5;i++){
7422 for(
unsigned int j=0;j<5;j++){
7423 dummy.push_back(Clast(i,j));
7425 fcov.push_back(dummy);
7451 unsigned int max_cdc_index=num_cdchits-1;
7456 vector<pull_t>last_cdc_pulls;
7457 vector<pull_t>cdc_pulls;
7458 vector<bool>last_cdc_used_in_fit(num_cdchits);
7463 double chisq_iter=-1.,
chisq=-1.;
7464 unsigned int my_ndf=0;
7466 unsigned int last_ndf=1;
7474 _DBG_ <<
"-------- iteration " << iter2+1 <<
" -----------" <<endl;
7499 if (ref_track_error!=NOERROR){
7505 for (
unsigned int j=0;j<num_cdchits;j++){
7511 bool did_second_refit=
false;
7521 unsigned int temp_ndf=my_ndf;
7522 double temp_chi2=
chisq;
7526 unsigned int new_index=(3*max_cdc_index)/4;
7536 unsigned int new_index=(max_cdc_index)/2;
7540 unsigned int new_index=(3*max_cdc_index)/4;
7548 anneal_factor=1000.;
7551 did_second_refit=
true;
7562 if (
DEBUG_LEVEL > 1)
_DBG_ <<
" Refit did not succeed, but restoring old values" << endl;
7574 if (iter2==0)
return error;
7579 <<
" Prob: " << TMath::Prob(
chisq,my_ndf)
7582 if (iter2>
MIN_ITER && did_second_refit==
false){
7583 double new_reduced_chisq=
chisq/my_ndf;
7584 double old_reduced_chisq=chisq_iter/last_ndf;
7585 double new_prob=TMath::Prob(
chisq,my_ndf);
7586 double old_prob=TMath::Prob(chisq_iter,last_ndf);
7587 if (new_prob<old_prob
7588 || fabs(new_reduced_chisq-old_reduced_chisq)<
CHISQ_DELTA)
break;
7608 pulls.assign(cdc_pulls.begin(),cdc_pulls.end());
7621 if (last_pos.Mod()>0.001){
7627 for (
unsigned int m=0;m<last_cdc_used_in_fit.size();m++){
7628 if (last_cdc_used_in_fit[m]){
7644 double dx=
x_-beam_pos.X();
7645 double dy=
y_-beam_pos.Y();
7648 double cosphi=cos(
phi_);
7650 if ((dx>0.0 && sinphi>0.0) || (dy <0.0 && cosphi>0.0)
7651 || (dy>0.0 && cosphi<0.0) || (dx<0.0 && sinphi<0.0))
D_*=-1.;
7659 if (!isfinite(
x_) || !isfinite(
y_) || !isfinite(
z_) || !isfinite(
phi_)
7662 _DBG_ <<
"At least one parameter is NaN or +-inf!!" <<endl;
7671 vector<double>dummy;
7672 for (
unsigned int i=0;i<5;i++){
7674 for(
unsigned int j=0;j<5;j++){
7675 dummy.push_back(Cclast(i,j));
7677 cov.push_back(dummy);
7693 if (
forward_traj.size()<2)
return RESOURCE_UNAVAILABLE;
7704 jout <<
"---- Smoothed residuals ----" <<endl;
7705 jout << setprecision(4);
7708 for (
unsigned int m=max-1;m>0;m--){
7712 for (
unsigned int k=0;k<5;k++){
7715 for (
unsigned int k=0;k<5;k++){
7716 mlfile <<
" " << Cs(k,k);
7717 for (
unsigned int j=k+1;j<5;j++){
7734 if (
DEBUG_LEVEL>1)
_DBG_ <<
"Invalid values for smoothed parameters..." << endl;
7735 return VALUE_OUT_OF_RANGE;
7741 _DBG_ <<
"Covariance Matrix not PosDef..." << endl;
7742 return VALUE_OUT_OF_RANGE;
7760 double vpred_uncorrected=x*sina+y*cosa;
7763 double upred=x*cosa-y*sina;
7766 double tu=tx*cosa-ty*sina;
7767 double one_plus_tu2=1.+tu*tu;
7768 double alpha=atan(tu);
7769 double cosalpha=cos(alpha);
7771 double sinalpha=
sin(alpha);
7775 double doca=du*cosalpha;
7780 double nz_sinalpha_plus_nr_cosalpha=nz*sinalpha+nr*cosalpha;
7783 double tv=tx*sina+ty*cosa;
7784 double resi=v-(vpred_uncorrected+doca*(nz_sinalpha_plus_nr_cosalpha
7793 double resi_a=drift-doca;
7800 double temp2=nz_sinalpha_plus_nr_cosalpha-tv*sinalpha;
7801 H_T(
state_x,1)=sina+cosa*cosalpha*temp2;
7802 H_T(
state_y,1)=cosa-sina*cosalpha*temp2;
7804 double cos2_minus_sin2=cosalpha*cosalpha-sinalpha*sinalpha;
7805 double fac=nz*cos2_minus_sin2-2.*nr*cosalpha*sinalpha;
7806 double doca_cosalpha=doca*cosalpha;
7807 double temp=doca_cosalpha*fac;
7809 -doca_cosalpha*(tu*sina+tv*cosa*cos2_minus_sin2)
7812 -doca_cosalpha*(tu*cosa-tv*sina*cos2_minus_sin2)
7816 H_T(
state_y,0)=-sina*cosalpha;
7817 double factor=du*tu/
sqrt(one_plus_tu2)/one_plus_tu2;
7834 vector<double> alignmentDerivatives;
7849 double cos2a = cos(2.*
my_fdchits[
id]->hit->wire->angle);
7851 double cos3a = cos(3.*
my_fdchits[
id]->hit->wire->angle);
7856 pow(1 + pow(tx*cosa - ty*sina,2),1.5);
7858 pow(1 + pow(tx*cosa - ty*sina,2),1.5));
7860 cosa*(-(tx*ty*x) + y + pow(tx,2)*y + (tx - ty)*(tx + ty)*u*sina))/
7861 pow(1 + pow(tx*cosa - ty*sina,2),1.5);
7865 double drift_shift = 0.0;
7867 if (drift_time < 0.) drift_shift = drift;
7879 alignmentDerivatives[
FDCTrackD::dDOCAW_dtx] = (cosa*(tx*cosa - ty*sina)*(u - x*cosa + y*sina))/pow(1 + pow(tx*cosa - ty*sina,2),1.5);
7883 pow(1 + pow(tx*cosa - ty*sina,2),1.5);
7889 (-nr + (-nz + ty*cosa + tx*sina)*(tx*cosa - ty*sina))/(1 + pow(tx*cosa - ty*sina,2));
7893 nz + (-nz + (nr*tx + ty)*cosa + (tx - nr*ty)*sina)/(1 + pow(tx*cosa - ty*sina,2));
7897 (-2*y*cosa*sina*(tx*cosa - ty*sina) - 2*x*pow(sina,2)*(tx*cosa - ty*sina) -
7898 (u - x*cosa + y*sina)*(-(nz*sina) + sina*(ty*cosa + tx*sina) -
7899 cosa*(tx*cosa - ty*sina)))/(1 + pow(tx*cosa - ty*sina,2)) +
7900 (2*sina*(tx*cosa - ty*sina)*(-((u - x*cosa + y*sina)*
7901 (nr + nz*(tx*cosa - ty*sina) - (ty*cosa + tx*sina)*(tx*cosa - ty*sina))) +
7902 y*cosa*(1 + pow(tx*cosa - ty*sina,2)) + x*sina*(1 + pow(tx*cosa - ty*sina,2))))/
7903 pow(1 + pow(tx*cosa - ty*sina,2),2);
7908 (u - x*cosa + y*sina)*(-(nz*cosa) + cosa*(ty*cosa + tx*sina) +
7909 sina*(tx*cosa - ty*sina)))/(1 + pow(tx*cosa - ty*sina,2)) +
7910 (2*cosa*(tx*cosa - ty*sina)*(-((u - x*cosa + y*sina)*
7911 (nr + nz*(tx*cosa - ty*sina) - (ty*cosa + tx*sina)*(tx*cosa - ty*sina))) +
7912 y*cosa*(1 + pow(tx*cosa - ty*sina,2)) + x*sina*(1 + pow(tx*cosa - ty*sina,2))))/
7913 pow(1 + pow(tx*cosa - ty*sina,2),2);
7917 (-((u - x*cosa + y*sina)*(nr + nz*(tx*cosa - ty*sina) -
7918 (ty*cosa + tx*sina)*(tx*cosa - ty*sina))) +
7919 y*cosa*(1 + pow(tx*cosa - ty*sina,2)) + x*sina*(1 + pow(tx*cosa - ty*sina,2))))/
7920 pow(1 + pow(tx*cosa - ty*sina,2),2) +
7921 (2*y*cosa*(ty*cosa + tx*sina)*(tx*cosa - ty*sina) +
7922 2*x*sina*(ty*cosa + tx*sina)*(tx*cosa - ty*sina) -
7923 (-(y*cosa) - x*sina)*(nr + nz*(tx*cosa - ty*sina) -
7924 (ty*cosa + tx*sina)*(tx*cosa - ty*sina)) -
7925 x*cosa*(1 + pow(tx*cosa - ty*sina,2)) + y*sina*(1 + pow(tx*cosa - ty*sina,2)) -
7926 (u - x*cosa + y*sina)*(nz*(ty*cosa + tx*sina) - pow(ty*cosa + tx*sina,2) -
7927 (tx*cosa - ty*sina)*(-(tx*cosa) + ty*sina)))/(1 + pow(tx*cosa - ty*sina,2));
7930 alignmentDerivatives[
FDCTrackD::dDOCAC_dx] = (cosa*(nr - tx*ty + nz*tx*cosa) + sina + ty*(ty - nz*cosa)*sina)/
7931 (1 + pow(tx*cosa - ty*sina,2));
7934 alignmentDerivatives[
FDCTrackD::dDOCAC_dy] = ((1 + pow(tx,2))*cosa - (nr + tx*ty + nz*tx*cosa)*sina + nz*ty*pow(sina,2))/
7935 (1 + pow(tx*cosa - ty*sina,2));
7938 alignmentDerivatives[
FDCTrackD::dDOCAC_dtx] = ((u - x*cosa + y*sina)*(4*nr*tx - 2*ty*(pow(tx,2) + pow(ty,2)) + nz*(-4 + 3*pow(tx,2) + pow(ty,2))*cosa +
7939 2*(2*nr*tx + ty*(2 - pow(tx,2) + pow(ty,2)))*cos2a + nz*(tx - ty)*(tx + ty)*cos3a - 2*nz*tx*ty*sina +
7940 4*(tx - nr*ty + tx*pow(ty,2))*sin2a - 2*nz*tx*ty*sin3a))/
7941 pow(2 + pow(tx,2) + pow(ty,2) + (tx - ty)*(tx + ty)*cos2a - 2*tx*ty*sin2a,2);
7944 alignmentDerivatives[
FDCTrackD::dDOCAC_dty] = -(((u - x*cosa + y*sina)*(-2*(pow(tx,3) + 2*nr*ty + tx*pow(ty,2)) - 2*nz*tx*ty*cosa -
7945 2*(2*tx + pow(tx,3) - 2*nr*ty - tx*pow(ty,2))*cos2a + 2*nz*tx*ty*cos3a +
7946 nz*(-4 + pow(tx,2) + 3*pow(ty,2))*sina + 4*(ty + tx*(nr + tx*ty))*sin2a + nz*(tx - ty)*(tx + ty)*sin3a))
7947 /pow(2 + pow(tx,2) + pow(ty,2) + (tx - ty)*(tx + ty)*cos2a - 2*tx*ty*sin2a,2));
7952 jout <<
"Layer " <<
my_fdchits[id]->hit->wire->layer
7953 <<
": t " << drift_time <<
" x "<< x <<
" y " << y
7954 <<
" coordinate along wire " << v <<
" resi_c " <<resi
7955 <<
" coordinate transverse to wire " << drift
7956 <<
" resi_a " << resi_a
7960 double scale=1./
sqrt(1.+tx*tx+ty*ty);
7961 double cosThetaRel=
my_fdchits[id]->hit->wire->udir.Dot(
DVector3(scale*tx,scale*ty,scale));
7972 forward_pulls.push_back(thisPull);
7989 if (!Cs.IsPosDef()){
7991 _DBG_ <<
"Covariance Matrix not PosDef..." << endl;
7992 return VALUE_OUT_OF_RANGE;
7995 if (
DEBUG_LEVEL>5)
_DBG_ <<
"Invalid values for smoothed parameters..." << endl;
7996 return VALUE_OUT_OF_RANGE;
8001 cdc_updates[
id],forward_pulls) != NOERROR)
return VALUE_OUT_OF_RANGE;
8027 if (
central_traj.size()<2)
return RESOURCE_UNAVAILABLE;
8038 _DBG_ <<
" S C JT at start of smoothing " << endl;
8042 for (
unsigned int m=max-1;m>0;m--){
8045 if (
DEBUG_LEVEL>1)
_DBG_ <<
" Encountered Hit ID " <<
id <<
" At trajectory position " << m <<
"/" << max << endl;
8047 if (
DEBUG_LEVEL>1)
_DBG_ <<
" SmoothCentral CDC Hit ID " <<
id <<
" used in fit " << endl;
8056 _DBG_ <<
"Invalid values for smoothed parameters..." << endl;
8057 return VALUE_OUT_OF_RANGE;
8059 if (!Cs.IsPosDef()){
8061 _DBG_ <<
"Covariance Matrix not PosDef... Ckk dC A" << endl;
8064 return VALUE_OUT_OF_RANGE;
8083 if(
BrentCentral(dEdx,xy,z0wire,origin,dir,myS,myds)!=NOERROR)
return VALUE_OUT_OF_RANGE;
8089 double d=cosstereo*diff.Mod()+
EPS;
8093 vector<double> alignmentDerivatives;
8095 double dscut_min=0., dscut_max=1.;
8097 double cosstereo_shifted;
8100 alignmentDerivatives.resize(12);
8103 double wposShift=0.025;
8104 double wdirShift=0.00005;
8107 double shiftFraction=0.01;
8109 double shift_phi=0.0001;
8110 double shift_tanl=shiftFraction*Ss(
state_tanl);
8111 double shift_D=0.01;
8112 double shift_z=0.01;
8116 DVector2 shift, origin_shifted, dir_shifted, wirepos_shifted, diff_shifted, xy_shifted;
8119 double d_dOriginX=0., d_dOriginY=0., d_dOriginZ=0.;
8120 double d_dDirX=0., d_dDirY=0., d_dDirZ=0.;
8121 double d_dS0=0., d_dS1=0., d_dS2=0., d_dS3=0., d_dS4=0.;
8127 shift.Set(wposShift, 0.);
8128 origin_shifted=origin+shift;
8133 if (
BrentCentral(dEdx,xy_shifted,z0_shifted,origin_shifted,
8134 dir_shifted,alignS,alignds)!=NOERROR)
return VALUE_OUT_OF_RANGE;
8135 if (alignds < dscut_min || alignds > dscut_max)
return VALUE_OUT_OF_RANGE;
8138 wirepos_shifted=origin_shifted+(alignS(
state_z)-z0_shifted)*dir_shifted;
8139 diff_shifted=xy_shifted-wirepos_shifted;
8140 d_dOriginX=cosstereo*diff_shifted.Mod()+
EPS;
8151 shift.Set(0.,wposShift);
8152 origin_shifted=origin+shift;
8157 if (
BrentCentral(dEdx,xy_shifted,z0_shifted,origin_shifted,
8158 dir_shifted,alignS,alignds)!=NOERROR)
return VALUE_OUT_OF_RANGE;
8159 if (alignds < dscut_min || alignds > dscut_max)
return VALUE_OUT_OF_RANGE;
8162 wirepos_shifted=origin_shifted+(alignS(
state_z)-z0_shifted)*dir_shifted;
8163 diff_shifted=xy_shifted-wirepos_shifted;
8164 d_dOriginY=cosstereo*diff_shifted.Mod()+
EPS;
8175 origin_shifted=origin;
8177 z0_shifted=z0wire+wposShift;
8180 if (
BrentCentral(dEdx,xy_shifted,z0_shifted,origin_shifted,
8181 dir_shifted,alignS,alignds)!=NOERROR)
return VALUE_OUT_OF_RANGE;
8182 if (alignds < dscut_min || alignds > dscut_max)
return VALUE_OUT_OF_RANGE;
8185 wirepos_shifted=origin_shifted+(alignS(
state_z)-z0_shifted)*dir_shifted;
8186 diff_shifted=xy_shifted-wirepos_shifted;
8187 d_dOriginZ=cosstereo*diff_shifted.Mod()+
EPS;
8198 shift.Set(wdirShift,0.);
8199 origin_shifted=origin;
8203 dir_shifted=dir+shift;
8204 cosstereo_shifted = cos((wireDir+
DVector3(wdirShift,0.,0.)).Angle(
DVector3(0.,0.,1.)));
8205 if (
BrentCentral(dEdx,xy_shifted,z0_shifted,origin_shifted,
8206 dir_shifted,alignS,alignds)!=NOERROR)
return VALUE_OUT_OF_RANGE;
8207 if (alignds < dscut_min || alignds > dscut_max)
return VALUE_OUT_OF_RANGE;
8210 wirepos_shifted=origin_shifted+(alignS(
state_z)-z0_shifted)*dir_shifted;
8211 diff_shifted=xy_shifted-wirepos_shifted;
8212 d_dDirX=cosstereo_shifted*diff_shifted.Mod()+
EPS;
8222 shift.Set(0.,wdirShift);
8223 origin_shifted=origin;
8227 dir_shifted=dir+shift;
8228 cosstereo_shifted = cos((wireDir+
DVector3(0.,wdirShift,0.)).Angle(
DVector3(0.,0.,1.)));
8229 if (
BrentCentral(dEdx,xy_shifted,z0_shifted,origin_shifted,
8230 dir_shifted,alignS,alignds)!=NOERROR)
return VALUE_OUT_OF_RANGE;
8231 if (alignds < dscut_min || alignds > dscut_max)
return VALUE_OUT_OF_RANGE;
8234 wirepos_shifted=origin_shifted+(alignS(
state_z)-z0_shifted)*dir_shifted;
8235 diff_shifted=xy_shifted-wirepos_shifted;
8236 d_dDirY=cosstereo_shifted*diff_shifted.Mod()+
EPS;
8246 origin_shifted=origin;
8247 dir_shifted.Set(wireDir.X()/(wireDir.Z()+wdirShift), wireDir.Y()/(wireDir.Z()+wdirShift));
8251 cosstereo_shifted = cos((wireDir+
DVector3(0.,0.,wdirShift)).Angle(
DVector3(0.,0.,1.)));
8252 if (
BrentCentral(dEdx,xy_shifted,z0_shifted,origin_shifted,
8253 dir_shifted,alignS,alignds)!=NOERROR)
return VALUE_OUT_OF_RANGE;
8254 if (alignds < dscut_min || alignds > dscut_max)
return VALUE_OUT_OF_RANGE;
8257 wirepos_shifted=origin_shifted+(alignS(
state_z)-z0_shifted)*dir_shifted;
8258 diff_shifted=xy_shifted-wirepos_shifted;
8259 d_dDirZ=cosstereo_shifted*diff_shifted.Mod()+
EPS;
8269 DMatrix5x1 trackShiftS0(shift_q_over_pt, 0., 0., 0., 0.);
8270 DMatrix5x1 trackShiftS1(0., shift_phi, 0., 0., 0.);
8271 DMatrix5x1 trackShiftS2(0., 0., shift_tanl, 0., 0.);
8272 DMatrix5x1 trackShiftS3(0., 0., 0., shift_D, 0.);
8273 DMatrix5x1 trackShiftS4(0., 0., 0., 0., shift_z);
8276 alignS=Ss+trackShiftS0;
8280 if(
BrentCentral(dEdx,xy_shifted,z0wire,origin,dir,alignS,alignds) != NOERROR)
return VALUE_OUT_OF_RANGE;
8281 if (alignds < dscut_min || alignds > dscut_max)
return VALUE_OUT_OF_RANGE;
8283 wirepos_shifted=origin+(alignS(
state_z)-z0wire)*dir;
8284 diff_shifted=xy_shifted-wirepos_shifted;
8285 d_dS0=cosstereo*diff_shifted.Mod()+
EPS;
8293 alignS=Ss+trackShiftS1;
8297 if(
BrentCentral(dEdx,xy_shifted,z0wire,origin,dir,alignS,alignds) != NOERROR)
return VALUE_OUT_OF_RANGE;
8298 if (alignds < dscut_min || alignds > dscut_max)
return VALUE_OUT_OF_RANGE;
8300 wirepos_shifted=origin+(alignS(
state_z)-z0wire)*dir;
8301 diff_shifted=xy_shifted-wirepos_shifted;
8302 d_dS1=cosstereo*diff_shifted.Mod()+
EPS;
8310 alignS=Ss+trackShiftS2;
8314 if(
BrentCentral(dEdx,xy_shifted,z0wire,origin,dir,alignS,alignds) != NOERROR)
return VALUE_OUT_OF_RANGE;
8315 if (alignds < dscut_min || alignds > dscut_max)
return VALUE_OUT_OF_RANGE;
8317 wirepos_shifted=origin+(alignS(
state_z)-z0wire)*dir;
8318 diff_shifted=xy_shifted-wirepos_shifted;
8319 d_dS2=cosstereo*diff_shifted.Mod()+
EPS;
8327 alignS=Ss+trackShiftS3;
8331 if(
BrentCentral(dEdx,xy_shifted,z0wire,origin,dir,alignS,alignds) != NOERROR)
return VALUE_OUT_OF_RANGE;
8332 if (alignds < dscut_min || alignds > dscut_max)
return VALUE_OUT_OF_RANGE;
8334 wirepos_shifted=origin+(alignS(
state_z)-z0wire)*dir;
8335 diff_shifted=xy_shifted-wirepos_shifted;
8336 d_dS3=cosstereo*diff_shifted.Mod()+
EPS;
8344 alignS=Ss+trackShiftS4;
8348 if(
BrentCentral(dEdx,xy_shifted,z0wire,origin,dir,alignS,alignds) != NOERROR)
return VALUE_OUT_OF_RANGE;
8349 if (alignds < dscut_min || alignds > dscut_max)
return VALUE_OUT_OF_RANGE;
8351 wirepos_shifted=origin+(alignS(
state_z)-z0wire)*dir;
8352 diff_shifted=xy_shifted-wirepos_shifted;
8353 d_dS4=cosstereo*diff_shifted.Mod()+
EPS;
8376 double cosstereo2_over_doca=cosstereo*cosstereo/d;
8377 H_T(
state_D)=(dy*cosphi-dx*sinphi)*cosstereo2_over_doca;
8379 =-myS(
state_D)*cosstereo2_over_doca*(dx*cosphi+dy*sinphi);
8380 H_T(
state_z)=-cosstereo2_over_doca*(dx*dir.X()+dy*dir.Y());
8389 double Vtrack=H*Cs*H_T;
8393 if (skip_ring) VRes = Vhit + Vtrack;
8394 else VRes = Vhit - Vtrack;
8396 if (
DEBUG_LEVEL>1 && (!isfinite(VRes) || VRes < 0.0) )
_DBG_ <<
" SmoothCentral Problem: VRes is " << VRes <<
" = " << Vhit <<
" - " << Vtrack << endl;
8399 double sinl=
sin(lambda);
8400 double cosl=cos(lambda);
8407 diff.Phi(),myS(
state_z),cosThetaRel,
8411 cdc_pulls.push_back(thisPull);
8442 if (
forward_traj.size()<2)
return RESOURCE_UNAVAILABLE;
8451 for (
unsigned int m=max-1;m>0;m--){
8456 _DBG_ <<
" Smoothing CDC index " << cdc_index <<
" ring " <<
my_cdchits[cdc_index]->hit->wire->ring
8457 <<
" straw " <<
my_cdchits[cdc_index]->hit->wire->straw << endl;
8464 _DBG_ <<
"Invalid values for smoothed parameters..." << endl;
8465 return VALUE_OUT_OF_RANGE;
8472 _DBG_ <<
"Covariance Matrix not Pos Def..." << endl;
8473 _DBG_ <<
" cdc_updates[cdc_index].C A C_ Cs " << endl;
8479 return VALUE_OUT_OF_RANGE;
8482 cdc_updates[cdc_index],cdc_pulls) != NOERROR)
return VALUE_OUT_OF_RANGE;
8511 vector<pull_t>&my_pulls){
8525 double z0wire=hit->
z0wire;
8526 if(
BrentForward(z,dEdx,z0wire,origin,dir,myS,mydz) != NOERROR)
return VALUE_OUT_OF_RANGE;
8528 double new_z=z+mydz;
8529 DVector2 wirepos=origin+(new_z-z0wire)*dir;
8534 double d=cosstereo*diff.Mod();
8538 vector<double> alignmentDerivatives;
8540 double dzcut_min=0., dzcut_max=1.;
8543 double cosstereo_shifted;
8545 alignmentDerivatives.resize(12);
8551 double wposShift=0.025;
8552 double wdirShift=0.00005;
8555 double shiftFraction=0.01;
8556 double shift_x=0.01;
8557 double shift_y=0.01;
8558 double shift_tx=shiftFraction*Ss(
state_tx);
8559 double shift_ty=shiftFraction*Ss(
state_ty);;
8563 double z0_shifted, new_z_shifted;
8564 DVector2 shift, origin_shifted, dir_shifted, wirepos_shifted, diff_shifted, xy_shifted;
8567 double d_dOriginX=0., d_dOriginY=0., d_dOriginZ=0.;
8568 double d_dDirX=0., d_dDirY=0., d_dDirZ=0.;
8569 double d_dS0=0., d_dS1=0., d_dS2=0., d_dS3=0., d_dS4=0.;
8575 shift.Set(wposShift, 0.);
8576 origin_shifted=origin+shift;
8579 if(
BrentForward(z,dEdx,z0_shifted,origin_shifted,dir_shifted,alignS,aligndz) != NOERROR)
return VALUE_OUT_OF_RANGE;
8580 if(aligndz < dzcut_min || aligndz > dzcut_max)
return VALUE_OUT_OF_RANGE;
8582 new_z_shifted=z+aligndz;
8583 wirepos_shifted=origin_shifted+(new_z_shifted-z0_shifted)*dir_shifted;
8585 diff_shifted=xy_shifted-wirepos_shifted;
8586 d_dOriginX=cosstereo*diff_shifted.Mod()+
EPS;
8597 shift.Set(0.,wposShift);
8598 origin_shifted=origin+shift;
8601 if(
BrentForward(z,dEdx,z0_shifted,origin_shifted,dir_shifted,alignS,aligndz) != NOERROR)
return VALUE_OUT_OF_RANGE;
8602 if(aligndz < dzcut_min || aligndz > dzcut_max)
return VALUE_OUT_OF_RANGE;
8604 new_z_shifted=z+aligndz;
8605 wirepos_shifted=origin_shifted+(new_z_shifted-z0_shifted)*dir_shifted;
8607 diff_shifted=xy_shifted-wirepos_shifted;
8608 d_dOriginY=cosstereo*diff_shifted.Mod()+
EPS;
8619 origin_shifted=origin;
8621 z0_shifted=z0wire+wposShift;
8622 if(
BrentForward(z,dEdx,z0_shifted,origin_shifted,dir_shifted,alignS,aligndz) != NOERROR)
return VALUE_OUT_OF_RANGE;
8623 if(aligndz < dzcut_min || aligndz > dzcut_max)
return VALUE_OUT_OF_RANGE;
8625 new_z_shifted=z+aligndz;
8626 wirepos_shifted=origin_shifted+(new_z_shifted-z0_shifted)*dir_shifted;
8628 diff_shifted=xy_shifted-wirepos_shifted;
8629 d_dOriginZ=cosstereo*diff_shifted.Mod()+
EPS;
8640 shift.Set(wdirShift,0.);
8641 origin_shifted=origin;
8643 dir_shifted=dir+shift;
8644 cosstereo_shifted = cos((wireDir+
DVector3(wdirShift,0.,0.)).Angle(
DVector3(0.,0.,1.)));
8645 if(
BrentForward(z,dEdx,z0_shifted,origin_shifted,dir_shifted,alignS,aligndz) != NOERROR)
return VALUE_OUT_OF_RANGE;
8646 if(aligndz < dzcut_min || aligndz > dzcut_max)
return VALUE_OUT_OF_RANGE;
8648 new_z_shifted=z+aligndz;
8649 wirepos_shifted=origin_shifted+(new_z_shifted-z0_shifted)*dir_shifted;
8651 diff_shifted=xy_shifted-wirepos_shifted;
8652 d_dDirX=cosstereo_shifted*diff_shifted.Mod()+
EPS;
8662 shift.Set(0.,wdirShift);
8663 origin_shifted=origin;
8665 dir_shifted=dir+shift;
8666 cosstereo_shifted = cos((wireDir+
DVector3(0.,wdirShift,0.)).Angle(
DVector3(0.,0.,1.)));
8667 if(
BrentForward(z,dEdx,z0_shifted,origin_shifted,dir_shifted,alignS,aligndz) != NOERROR)
return VALUE_OUT_OF_RANGE;
8668 if(aligndz < dzcut_min || aligndz > dzcut_max)
return VALUE_OUT_OF_RANGE;
8670 new_z_shifted=z+aligndz;
8671 wirepos_shifted=origin_shifted+(new_z_shifted-z0_shifted)*dir_shifted;
8673 diff_shifted=xy_shifted-wirepos_shifted;
8674 d_dDirY=cosstereo_shifted*diff_shifted.Mod()+
EPS;
8684 origin_shifted=origin;
8685 dir_shifted.Set(wireDir.X()/(wireDir.Z()+wdirShift), wireDir.Y()/(wireDir.Z()+wdirShift));
8687 cosstereo_shifted = cos((wireDir+
DVector3(0.,0.,wdirShift)).Angle(
DVector3(0.,0.,1.)));
8688 if(
BrentForward(z,dEdx,z0_shifted,origin_shifted,dir_shifted,alignS,aligndz) != NOERROR)
return VALUE_OUT_OF_RANGE;
8689 if(aligndz < dzcut_min || aligndz > dzcut_max)
return VALUE_OUT_OF_RANGE;
8691 new_z_shifted=z+aligndz;
8692 wirepos_shifted=origin_shifted+(new_z_shifted-z0_shifted)*dir_shifted;
8694 diff_shifted=xy_shifted-wirepos_shifted;
8695 d_dDirZ=cosstereo_shifted*diff_shifted.Mod()+
EPS;
8704 DMatrix5x1 trackShiftS0(shift_x, 0., 0., 0., 0.);
8705 DMatrix5x1 trackShiftS1(0., shift_y, 0., 0., 0.);
8706 DMatrix5x1 trackShiftS2(0., 0., shift_tx, 0., 0.);
8707 DMatrix5x1 trackShiftS3(0., 0., 0., shift_ty, 0.);
8708 DMatrix5x1 trackShiftS4(0., 0., 0., 0., shift_q_over_p);
8711 alignS=Ss+trackShiftS0;
8713 if(
BrentForward(z,dEdx,z0wire,origin,dir,alignS,aligndz) != NOERROR)
return VALUE_OUT_OF_RANGE;
8714 if(aligndz < dzcut_min || aligndz > dzcut_max)
return VALUE_OUT_OF_RANGE;
8716 new_z_shifted=z+aligndz;
8717 wirepos_shifted=origin+(new_z_shifted-z0wire)*dir;
8719 diff_shifted=xy_shifted-wirepos_shifted;
8720 d_dS0=cosstereo*diff_shifted.Mod()+
EPS;
8728 alignS=Ss+trackShiftS1;
8730 if(
BrentForward(z,dEdx,z0wire,origin,dir,alignS,aligndz) != NOERROR)
return VALUE_OUT_OF_RANGE;
8731 if(aligndz < dzcut_min || aligndz > dzcut_max)
return VALUE_OUT_OF_RANGE;
8733 new_z_shifted=z+aligndz;
8734 wirepos_shifted=origin+(new_z_shifted-z0wire)*dir;
8736 diff_shifted=xy_shifted-wirepos_shifted;
8737 d_dS1=cosstereo*diff_shifted.Mod()+
EPS;
8745 alignS=Ss+trackShiftS2;
8747 if(
BrentForward(z,dEdx,z0wire,origin,dir,alignS,aligndz) != NOERROR)
return VALUE_OUT_OF_RANGE;
8748 if(aligndz < dzcut_min || aligndz > dzcut_max)
return VALUE_OUT_OF_RANGE;
8750 new_z_shifted=z+aligndz;
8751 wirepos_shifted=origin+(new_z_shifted-z0wire)*dir;
8753 diff_shifted=xy_shifted-wirepos_shifted;
8754 d_dS2=cosstereo*diff_shifted.Mod()+
EPS;
8762 alignS=Ss+trackShiftS3;
8764 if(
BrentForward(z,dEdx,z0wire,origin,dir,alignS,aligndz) != NOERROR)
return VALUE_OUT_OF_RANGE;
8765 if(aligndz < dzcut_min || aligndz > dzcut_max)
return VALUE_OUT_OF_RANGE;
8767 new_z_shifted=z+aligndz;
8768 wirepos_shifted=origin+(new_z_shifted-z0wire)*dir;
8770 diff_shifted=xy_shifted-wirepos_shifted;
8771 d_dS3=cosstereo*diff_shifted.Mod()+
EPS;
8779 alignS=Ss+trackShiftS4;
8781 if(
BrentForward(z,dEdx,z0wire,origin,dir,alignS,aligndz) != NOERROR)
return VALUE_OUT_OF_RANGE;
8782 if(aligndz < dzcut_min || aligndz > dzcut_max)
return VALUE_OUT_OF_RANGE;
8784 new_z_shifted=z+aligndz;
8785 wirepos_shifted=origin+(new_z_shifted-z0wire)*dir;
8787 diff_shifted=xy_shifted-wirepos_shifted;
8788 d_dS4=cosstereo*diff_shifted.Mod()+
EPS;
8806 double cosstereo2_over_d=cosstereo*cosstereo/d;
8807 H_T(
state_x)=diff.X()*cosstereo2_over_d;
8808 H_T(
state_y)=diff.Y()*cosstereo2_over_d;
8819 if (skip_ring) V+=H*myC*H_T;
8822 if (
DEBUG_LEVEL>1 && (!isfinite(V) || V < 0.0) )
_DBG_ <<
" Problem: V is " << V << endl;
8826 double scale=1./
sqrt(1.+tx*tx+ty*ty);
8830 NULL,diff.Phi(),new_z,cosThetaRel,update.
tcorr);
8832 my_pulls.push_back(thisPull);
8841 double cosl=cos(atan(
tanl_));
8843 double tanl3=tanl2*
tanl_;
8844 double factor=1./
sqrt(1.+tsquare);
8848 double frac=1./(tsquare*tsquare);
8869 wirepos+=(z-z0w)*dir;
8872 double doca2 = dx*dx+dy*dy;
8875 return VALUE_OUT_OF_RANGE;
8879 unsigned int maxSteps=5;
8880 unsigned int stepCounter=0;
8883 double old_doca2=doca2;
8887 Step(ztemp,newz,dedx,S);
8890 wirepos+=(newz-z0w)*dir;
8897 while(doca2<old_doca2 && stepCounter<maxSteps){
8902 Step(ztemp,newz,dedx,S);
8907 wirepos+=(newz-z0w)*dir;
8918 if (
BrentsAlgorithm(newz,-mStepSizeZ,dedx,z0w,origin,dir,S,dz2)!=NOERROR){
8919 return VALUE_OUT_OF_RANGE;
8933 wirepos+=(ztemp-z0w)*dir;
8936 double old_doca2=doca2;
8942 while(doca2<old_doca2 && stepCounter<10*maxSteps){
8947 Step(ztemp,newz,dedx,S);
8952 wirepos+=(newz-z0w)*dir;
8964 return VALUE_OUT_OF_RANGE;
8977 wirexy+=(Sc(
state_z)-z0w)*dir;
8980 double doca2=(xy-wirexy).Mod2();
8981 double old_doca2=doca2;
8984 origin,dir,Sc,ds)!=NOERROR){
8985 return VALUE_OUT_OF_RANGE;
8988 unsigned int maxSteps=3;
8989 unsigned int stepCounter=0;
8997 wirexy+=(Sc(
state_z)-z0w)*dir;
8998 doca2=(xy-wirexy).Mod2();
8999 while(doca2<old_doca2 && stepCounter<maxSteps){
9003 return VALUE_OUT_OF_RANGE;
9011 wirexy+=(Sc(
state_z)-z0w)*dir;
9012 doca2=(xy-wirexy).Mod2();
9019 origin,dir,Sc,ds2)!=NOERROR){
9020 return VALUE_OUT_OF_RANGE;
9029 wirexy+=(Sc(
state_z)-z0w)*dir;
9033 doca2=(xy-wirexy).Mod2();
9035 while(doca2<old_doca2 && stepCounter<maxSteps){
9040 return VALUE_OUT_OF_RANGE;
9049 wirexy+=(Sc(
state_z)-z0w)*dir;
9050 doca2=(xy-wirexy).Mod2();
9057 origin,dir,Sc,ds2)!=NOERROR){
9058 return VALUE_OUT_OF_RANGE;
9068 if (
forward_traj.size()<2)
return RESOURCE_UNAVAILABLE;
9073 unsigned int index_beyond_start_counter=inner_index;
9075 bool intersected_start_counter=
false;
9079 double d_old=1000.,d=1000.,z=0.;
9080 unsigned int index=0;
9081 for (
unsigned int m=0;m<12;m++){
9082 unsigned int k=inner_index;
9088 if (dphi<0) dphi+=2.*M_PI;
9089 index=int(floor(dphi/(2.*M_PI/30.)));
9090 if (index>29) index=0;
9094 if (m==0) index_beyond_start_counter=k;
9101 if (z>
sc_pos[index][m+1].z()&&m<11){
9117 while (fabs(d)>0.001 && count<20){
9128 double one_over_beta2=1.+
mass2*q_over_p_sq;
9129 if (one_over_beta2>
BIG) one_over_beta2=
BIG;
9130 t-=ds*
sqrt(one_over_beta2);
9135 if (dphi<0) dphi+=2.*M_PI;
9136 index=int(floor(dphi/(2.*M_PI/30.)));
9142 if (fabs(d)<fabs(dmin)){
9153 double tanl=1./
sqrt(tsquare);
9154 double cosl=cos(atan(tanl));
9158 DVector3 momentum(pt*cos(phi),pt*
sin(phi),pt*tanl);
9164 intersected_start_counter=
true;
9170 double s_theta_ms_sum=0.;
9171 double theta2ms_sum=0.;
9172 if (intersected_start_counter){
9173 for (
unsigned int k=inner_index;k>index_beyond_start_counter;k--){
9175 s_theta_ms_sum+=
sqrt(ds_theta_ms_sq);
9177 theta2ms_sum+=ds_theta_ms_sq/(ds*ds);
9182 unsigned int fdc_plane=0;
9185 for (
int k=intersected_start_counter?index_beyond_start_counter:inner_index;k>=0;k--){
9191 double tanl=1./
sqrt(tsquare);
9192 double cosl=cos(atan(tanl));
9217 s_theta_ms_sum+=
sqrt(ds_theta_ms_sq);
9219 theta2ms_sum+=ds_theta_ms_sq/(ds*ds);
9224 DVector3 momentum(pt*cos(phi),pt*
sin(phi),pt*tanl);
9227 s_theta_ms_sum,theta2ms_sum));
9240 Step(z,fdc_z_wires[fdc_plane],0.,S);
9242 tanl=1./
sqrt(tsquare);
9243 cosl=cos(atan(tanl));
9248 DVector3 momentum(pt*cos(phi),pt*
sin(phi),pt*tanl);
9251 s_theta_ms_sum,theta2ms_sum));
9274 double rho_Z_over_A=0.,LnI=0.,K_rho_Z_over_A=0.,Z=0.;
9275 double chi2c_factor=0.,chi2a_factor=0.,chi2a_corr=0.;
9279 double newz=z,dz=0.;
9287 const double z_outer_max=650.;
9288 const double x_max=130.;
9289 const double y_max=130.;
9291 bool hit_dirc=
false;
9301 return VALUE_OUT_OF_RANGE;
9306 if (r2>89.*89. && z<400.)
return VALUE_OUT_OF_RANGE;
9307 if (r2>64.9*64.9 && r2<89.*89.){
9309 return VALUE_OUT_OF_RANGE;
9313 double tanl=1./
sqrt(tsquare);
9314 double cosl=cos(atan(tanl));
9318 DVector3 momentum(pt*cos(phi),pt*
sin(phi),pt*tanl);
9337 double s_to_boundary=0.;
9339 chi2c_factor,chi2a_factor,chi2a_corr,
9344 _DBG_ <<
"Material error in ExtrapolateForwardToOuterDetectors!"<< endl;
9345 _DBG_ <<
" Position (x,y,z)=("<<pos.x()<<
","<<pos.y()<<
","<<pos.z()<<
")"
9348 return VALUE_OUT_OF_RANGE;
9358 if (fabs(dEdx)>
EPS){
9362 if (s_to_boundary<ds) ds=s_to_boundary;
9364 if (ds<0.5 && z<406. && r2>65.*65.) ds=0.5;
9367 if (hit_tof==
false && newz>
dTOFz){
9371 if (hit_dirc==
false && newz>
dDIRCz){
9379 bool got_fdc_hit=
false;
9389 double one_over_beta2=1.+
mass2*q_over_p_sq;
9390 if (one_over_beta2>
BIG) one_over_beta2=
BIG;
9391 t+=ds*
sqrt(one_over_beta2);
9397 s_theta_ms_sum+=
sqrt(ds_theta_ms_sq);
9398 theta2ms_sum+=ds_theta_ms_sq/(ds*ds);
9401 Step(z,newz,dEdx,S);
9405 double tanl=1./
sqrt(tsquare);
9406 double cosl=cos(atan(tanl));
9410 DVector3 momentum(pt*cos(phi),pt*
sin(phi),pt*tanl);
9413 s_theta_ms_sum,theta2ms_sum));
9417 if (hit_dirc==
false && newz>
dDIRCz){
9421 double tanl=1./
sqrt(tsquare);
9422 double cosl=cos(atan(tanl));
9426 DVector3 momentum(pt*cos(phi),pt*
sin(phi),pt*tanl);
9430 if (hit_tof==
false && newz>
dTOFz){
9434 double tanl=1./
sqrt(tsquare);
9435 double cosl=cos(atan(tanl));
9439 DVector3 momentum(pt*cos(phi),pt*
sin(phi),pt*tanl);
9445 double tanl=1./
sqrt(tsquare);
9446 double cosl=cos(atan(tanl));
9450 DVector3 momentum(pt*cos(phi),pt*
sin(phi),pt*tanl);
9455 double zend=newz+45.;
9456 Step(newz,zend,dEdx,S);
9457 ds=(zend-newz)/dz_ds;
9458 t+=ds*
sqrt(one_over_beta2);
9461 tanl=1./
sqrt(tsquare);
9462 cosl=cos(atan(tanl));
9466 momentum.SetXYZ(pt*cos(phi),pt*
sin(phi),pt*tanl);
9468 t*TIME_UNIT_CONVERSION,s));
9480 if (
central_traj.size()<2)
return RESOURCE_UNAVAILABLE;
9485 unsigned int index_beyond_start_counter=inner_index;
9490 double d_old=1000.,d=1000.,z=0.;
9491 unsigned int index=0;
9492 for (
unsigned int m=0;m<12;m++){
9493 unsigned int k=inner_index;
9500 if (dphi<0) dphi+=2.*M_PI;
9501 index=int(floor(dphi/(2.*M_PI/30.)));
9502 if (index>29) index=0;
9509 if (m==0) index_beyond_start_counter=k;
9516 if (z>
sc_pos[index][m+1].z()&&m<11){
9533 while (fabs(d)>0.05 && count<20){
9544 double q_over_p_sq=q_over_p*q_over_p;
9545 double one_over_beta2=1.+
mass2*q_over_p_sq;
9546 if (one_over_beta2>
BIG) one_over_beta2=
BIG;
9547 t-=ds*
sqrt(one_over_beta2);
9553 if (dphi<0) dphi+=2.*M_PI;
9554 index=int(floor(dphi/(2.*M_PI/30.)));
9555 if (index>29) index=0;
9561 if (fabs(d)<fabs(dmin)){
9574 DVector3 momentum(pt*cos(phi),pt*
sin(phi),pt*tanl);
9586 double s_theta_ms_sum=0.,theta2ms_sum=0.;
9587 for (
unsigned int k=inner_index;k>index_beyond_start_counter;k--){
9589 s_theta_ms_sum+=
sqrt(ds_theta_ms_sq);
9591 theta2ms_sum+=ds_theta_ms_sq/(ds*ds);
9597 for (
int k=index_beyond_start_counter;k>=0;k--){
9619 s_theta_ms_sum+=
sqrt(ds_theta_ms_sq);
9621 theta2ms_sum+=ds_theta_ms_sq/(ds*ds);
9625 DVector3 momentum(pt*cos(phi),pt*
sin(phi),pt*tanl);
9627 s_theta_ms_sum,theta2ms_sum));
9639 double r2=xy.Mod2();
9663 return VALUE_OUT_OF_RANGE;
9667 double rho_Z_over_A=0.,LnI=0.,K_rho_Z_over_A=0.,Z=0.;
9668 double chi2c_factor=0.,chi2a_factor=0.,chi2a_corr=0.;
9670 double s_to_boundary=0.;
9673 chi2c_factor,chi2a_factor,chi2a_corr,
9676 _DBG_ <<
"Material error in ExtrapolateToVertex! " << endl;
9677 _DBG_ <<
" Position (x,y,z)=("<<pos3d.x()<<
","<<pos3d.y()<<
","
9686 dedx=
GetdEdx(q_over_p,K_rho_Z_over_A,rho_Z_over_A,LnI,Z);
9689 if (fabs(dedx)>
EPS){
9694 if (s_to_boundary<ds) ds=s_to_boundary;
9696 if (ds<0.5 &&
S(
state_z)<400. && pos3d.Perp()>65.) ds=0.5;
9700 double q_over_p_sq=q_over_p*q_over_p;
9701 double one_over_beta2=1.+
mass2*q_over_p_sq;
9702 if (one_over_beta2>
BIG) one_over_beta2=
BIG;
9703 t+=ds*
sqrt(one_over_beta2);
9710 theta2ms_sum+=ds_theta_ms_sq/(ds*ds);
9719 return VALUE_OUT_OF_RANGE;
9726 DVector3 momentum(pt*cos(phi),pt*
sin(phi),pt*tanl);
unsigned int break_point_fdc_index
jerror_t ExtrapolateToVertex(DVector2 &xy, DMatrix5x1 &Sc, DMatrix5x5 &Cc)
void setMomentum(const DVector3 &aMomentum)
double fdc_drift_variance(double t) const
jerror_t SetCDCForwardReferenceTrajectory(DMatrix5x1 &S)
double CalcDensityEffect(double p, double mass, double density, double Z_over_A, double I)
void setTime(double locTime)
jerror_t ExtrapolateCentralToOtherDetectors(void)
bool GetFDCZ(vector< double > &z_wires) const
z-locations for each of the FDC wire planes in cm
vector< vector< double > > fcov
jerror_t BrentsAlgorithm(double ds1, double ds2, double dedx, DVector2 &pos, const double z0wire, const DVector2 &origin, const DVector2 &dir, DMatrix5x1 &Sc, double &ds_out)
vector< double > fdc_z_wires
virtual jerror_t SmoothForward(vector< pull_t > &mypulls)
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.
double MINIMUM_HIT_FRACTION
double CDC_DRIFT_BSCALE_PAR1
vector< vector< double > > max_sag
jerror_t CalcDerivAndJacobian(double z, double dz, const DMatrix5x1 &S, double dEdx, DMatrix5x5 &J, DMatrix5x1 &D)
void ClearExtrapolations(void)
double fdc_drift_distance(double t, double Bz) const
TH2I * brentCheckHists[2]
double t_v
time of the two cathode clusters
virtual void GetFieldAndGradient(double x, double y, double z, double &Bx, double &By, double &Bz, double &dBxdx, double &dBxdy, double &dBxdz, double &dBydx, double &dBydy, double &dBydz, double &dBzdx, double &dBzdy, double &dBzdz) const =0
deque< DKalmanForwardTrajectory_t > forward_traj
jerror_t PropagateForward(int length, int &index, double &z, double zhit, DMatrix5x1 &S, bool &done, bool &stepped_to_boundary, bool &stepped_to_endplate)
bool Get(string xpath, string &sval) const
bool GetFDCRmin(vector< double > &rmin_packages) const
beam hole size for each FDC package in cm
bool RECOVER_BROKEN_TRACKS
const DVector3 & position(void) const
void setTrackingErrorMatrix(const shared_ptr< const TMatrixFSym > &aMatrix)
double FORWARD_ANNEAL_SCALE
const DFDCWire * wire
DFDCWire for this wire.
virtual jerror_t ExtrapolateForwardToOtherDetectors(void)
vector< vector< DVector3 > > sc_norm
vector< double > cdc_drift_table
void setT0(double at0, double at0_err, DetectorSystem_t at0_detector)
void setForwardParmFlag(bool aFlag)
vector< double > cdc_origin
kalman_error_t ForwardFit(const DMatrix5x1 &S, const DMatrix5x5 &C0)
unsigned int break_point_step_index
shared_ptr< DResourcePool< TMatrixFSym > > dResourcePool_TMatrixFSym
vector< DKalmanUpdate_t > fdc_updates
kalman_error_t CentralFit(const DVector2 &startpos, const DMatrix5x1 &Sc, const DMatrix5x5 &C0)
class DFDCPseudo: definition for a reconstructed point in the FDC
map< DetectorSystem_t, vector< Extrapolation_t > > extrapolations
vector< const DFDCPseudo * > fdchits_used_in_fit
vector< DKalmanSIMDFDCHit_t * > my_fdchits
unsigned int potential_cdc_hits_on_track
bool GetFCALZ(double &z_fcal) const
z-location of front face of CCAL in cm
static char index(char c)
#define FDC_CATHODE_VARIANCE
double PHOTON_ENERGY_CUTOFF
vector< vector< DCDCWire * > > cdcwires
vector< vector< double > > sag_phi_offset
double ds
local coordinate of pseudopoint in the direction along the wire and its uncertainty ...
bool GetFDCRmax(double &rmax_active_fdc) const
outer radius of FDC active area in cm
vector< double > cdc_rmid
jerror_t SetCDCReferenceTrajectory(const DVector2 &xy, DMatrix5x1 &Sc)
static bool DKalmanSIMDFDCHit_cmp(DKalmanSIMDFDCHit_t *a, DKalmanSIMDFDCHit_t *b)
bool GetCDCWires(vector< vector< DCDCWire * > > &cdcwires) const
DTrackFitterKalmanSIMD(JEventLoop *loop)
double COVARIANCE_SCALE_FACTOR_FORWARD
double Chi2(const DMatrix2x1 &R) const
double long_drift_func[3][3]
unsigned int break_point_cdc_index
double FORWARD_ANNEAL_POW_CONST
vector< double > fdc_rmin_packages
const DMagneticFieldMap * bfield
DMatrix2x5 Transpose(const DMatrix5x2 &M)
jerror_t BrentForward(double z, double dedx, const double z0w, const DVector2 &origin, const DVector2 &dir, DMatrix5x1 &S, double &dz)
void setTrackingStateVector(double a1, double a2, double a3, double a4, double a5)
fit_status_t FitTrack(void)
void setErrorMatrix(const shared_ptr< const TMatrixFSym > &aMatrix)
double FDC_DRIFT_BSCALE_PAR2
double DRIFT_RES_PARMS[3]
jerror_t AddCDCHit(const DCDCTrackHit *cdchit)
vector< const DCDCTrackHit * > cdchits
void GetPosition(DVector3 &pos)
void locate(const double *xx, int n, double x, int *j)
unsigned int MIN_HITS_FOR_REFIT
double endplate_z_downstream
jerror_t PropagateForwardCDC(int length, int &index, double &z, double &r2, DMatrix5x1 &S, bool &stepped_to_boundary)
vector< vector< double > > cov
DMatrix5x5 SubSym(const DMatrix5x5 &m2) const
kalman_error_t ForwardCDCFit(const DMatrix5x1 &S, const DMatrix5x5 &C0)
double mCDCInternalStepSize
double dE
energy deposition, from integral
unsigned int GetNDF(void)
#define TIME_UNIT_CONVERSION
double FactorForSenseOfRotation
kalman_error_t KalmanCentral(double anneal_factor, DMatrix5x1 &S, DMatrix5x5 &C, DVector2 &xy, double &chisq, unsigned int &myndf)
kalman_error_t RecoverBrokenTracks(double anneal_factor, DMatrix5x1 &S, DMatrix5x5 &C, const DMatrix5x5 &C0, double &chisq, unsigned int &numdof)
double mT0MinimumDriftTime
double s
< wire position computed from cathode data, assuming the avalanche occurs at the wire ...
TH1I * alignDerivHists[46]
Double_t sigma[NCHANNELS]
shared_ptr< TMatrixFSym > Get7x7ErrorMatrixForward(DMatrixDSym C)
double CDC_VAR_SCALE_FACTOR
vector< DKalmanSIMDCDCHit_t * > my_cdchits
vector< const DCDCTrackHit * > cdchits_used_in_fit
void FastStep(double &z, double ds, double dEdx, DMatrix5x1 &S)
double charge(void) const
double COVARIANCE_SCALE_FACTOR_CENTRAL
jerror_t SmoothForwardCDC(vector< pull_t > &mypulls)
void TransformCovariance(DMatrix5x5 &C)
DetectorSystem_t mT0Detector
double time
time corresponding to this pseudopoint.
void setPID(Particle_t locPID)
kalman_error_t KalmanForwardCDC(double anneal, DMatrix5x1 &S, DMatrix5x5 &C, double &chisq, unsigned int &numdof)
void cdc_hit(Int_t &, Int_t &, Int_t &, Int_t[], Int_t, Int_t, Int_t, Int_t, Int_t, Int_t)
virtual double GetBz(double x, double y, double z) const =0
vector< const DFDCPseudo * > fdchits
double short_drift_func[3][3]
double ChiSq(fit_type_t fit_type, DReferenceTrajectory *rt, double *chisq_ptr=NULL, int *dof_ptr=NULL, vector< pull_t > *pulls_ptr=NULL)
jerror_t StepJacobian(double oldz, double newz, const DMatrix5x1 &S, double dEdx, DMatrix5x5 &J)
unsigned int last_material_map
double short_drift_Bscale_par1
jerror_t KalmanLoop(void)
double FDC_DRIFT_BSCALE_PAR1
double long_drift_Bscale_par2
void AddTrackDerivatives(vector< double > d)
void ComputeCDCDrift(double dphi, double delta, double t, double B, double &d, double &V, double &tcorr)
const DVector3 & momentum(void) const
jerror_t SmoothCentral(vector< pull_t > &cdc_pulls)
double DRIFT_FUNC_PARMS[6]
double long_drift_Bscale_par1
DTrackingData input_params
jerror_t StepStateAndCovariance(DVector2 &xy, double ds, double dEdx, DMatrix5x1 &S, DMatrix5x5 &J, DMatrix5x5 &C)
double CDC_DRIFT_BSCALE_PAR2
vector< vector< DVector3 > > sc_pos
DetectorSystem_t t0_detector(void) const
jerror_t GetProcessNoise(double z, double ds, double chi2c_factor, double chi2a_factor, double chi2a_corr, const DMatrix5x1 &S, DMatrix5x5 &Q)
kalman_error_t RecoverBrokenForwardTracks(double fdc_anneal_factor, double cdc_anneal_factor, DMatrix5x1 &S, DMatrix5x5 &C, const DMatrix5x5 &C0, double &chisq, unsigned int &numdof)
double Step(double oldz, double newz, double dEdx, DMatrix5x1 &S)
jerror_t CalcDeriv(double z, const DMatrix5x1 &S, double dEdx, DMatrix5x1 &D)
bool WRITE_ML_TRAINING_OUTPUT
shared_ptr< TMatrixFSym > Get7x7ErrorMatrix(DMatrixDSym C)
jerror_t AddFDCHit(const DFDCPseudo *fdchit)
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
vector< bool > cdc_used_in_fit
bool GetCDCEndplate(double &z, double &dz, double &rmin, double &rmax) const
static bool DKalmanSIMDCDCHit_cmp(DKalmanSIMDCDCHit_t *a, DKalmanSIMDCDCHit_t *b)
jerror_t SetReferenceTrajectory(DMatrix5x1 &S)
bool GetDIRCZ(double &z_dirc) const
z-location of DIRC in cm
void setPosition(const DVector3 &aPosition)
bool ENABLE_BOUNDARY_CHECK
jerror_t GetProcessNoiseCentral(double ds, double chi2c_factor, double chi2a_factor, double chi2a_corr, const DMatrix5x1 &S, DMatrix5x5 &Q)
virtual void GetField(const DVector3 &pos, DVector3 &Bout) const =0
static Particle_t IDTrack(float locCharge, float locMass)
jerror_t FillPullsVectorEntry(const DMatrix5x1 &Ss, const DMatrix5x5 &Cs, const DKalmanForwardTrajectory_t &traj, const DKalmanSIMDCDCHit_t *hit, const DKalmanUpdate_t &update, vector< pull_t > &mypulls)
printf("string=%s", string)
void ResetKalmanSIMD(void)
bool GetTargetZ(double &z_target) const
z-location of center of target
vector< DKalmanUpdate_t > cdc_updates
DMatrix5x5 AddSym(const DMatrix5x5 &m2) const
jerror_t BrentCentral(double dedx, DVector2 &xy, const double z0w, const DVector2 &origin, const DVector2 &dir, DMatrix5x1 &Sc, double &ds)
unsigned int potential_fdc_hits_on_track
vector< vector< DVector3 > > sc_dir
void GetMomentum(DVector3 &mom)
double GetEnergyVariance(double ds, double beta2, double K_rho_Z_over_A)
double short_drift_Bscale_par2
jerror_t PropagateCentral(int length, int &index, DVector2 &my_xy, double &var_t_factor, DMatrix5x1 &Sc, bool &stepped_to_boundary)
double GetdEdx(double q_over_p, double K_rho_Z_over_A, double rho_Z_over_A, double rho_Z_over_A_LnI, double Z)
vector< bool > fdc_used_in_fit
bool GetStartCounterGeom(vector< vector< DVector3 > > &pos, vector< vector< DVector3 > > &norm) const
deque< DKalmanCentralTrajectory_t > central_traj
virtual kalman_error_t KalmanForward(double fdc_anneal, double cdc_anneal, DMatrix5x1 &S, DMatrix5x5 &C, double &chisq, unsigned int &numdof)
double FasterStep(double oldz, double newz, double dEdx, DMatrix5x1 &S)